# Topological defect engineering enables size and shape control in self-assembly

Lara Koehler, Markus Eder, Vincent Ouazan-Reboul, Christoph Karfusehr, Andrey Zelenskiy, Pierre Ronceray, Friedrich Simmel, Martin Lenz

► To cite this version:

Lara Koehler, Markus Eder, Vincent Ouazan-Reboul, Christoph Karfusehr, Andrey Zelenskiy, et al. 2025. hal-05337745

HAL Id: hal-05337745

https://hal.science/hal-05337745v1

Preprint submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Topological Defect Engineering enables Size and Shape Control in Self-Assembly

Lara Koehler$^{1,2*}$, Markus Eder$^{3*}$, Vincent Ouazan-Reboul$^1$, Christoph Karfusehr$^{3,4}$, Andrey Zelenskiy$^1$, Pierre Ronceray$^5$, Friedrich C. Simmel$^3$, Martin Lenz$^{1,6*}$

The self-assembly of complex structures from engineered subunits is a major goal of nanotechnology, but controlling their size becomes increasingly difficult in larger assemblies. Existing strategies present significant challenges, among which the use of multiple subunit types or the precise control of their shape and mechanics. Here we introduce an alternative approach based on identical subunits whose interactions promote crystals, but also favor crystalline defects. We theoretically show that topological restrictions on the scope of these defects in large assemblies imply that the assembly size is controlled by the magnitude of the defect-inducing interaction. Using DNA origami, we experimentally demonstrate

<footnote>$^1$Université Paris-Saclay, CNRS, LPTMS, F-91405, Orsay, France</footnote>

<footnote>$^2$Max Planck Institute for the Physics of Complex Systems, Dresden, Germany</footnote>

<footnote>$^3$Department of Bioscience, School of Natural Sciences, Technical University of Munich, Garching, Germany</footnote>

<footnote>$^4$Max Planck School Matter to Life, Jahnstraße 29, Heidelberg, D-69120, Germany</footnote>

<footnote>$^5$Aix Marseille Université, CNRS, CINAM, Turing Center for Living Systems, 13288 Marseille, France.</footnote>

<footnote>$^6$PMMH, CNRS, ESPCI Paris, PSL University, Sorbonne Université, Université Paris-Cité, F-75005, Paris, France.</footnote>

<footnote>*Corresponding author. Email: martin.lenz@universite-paris-saclay.fr</footnote>

<footnote>†These authors contributed equally to this work.</footnote>
both size and shape control in two-dimensional disk- and fiber-like assemblies. Our
basic concept of defect engineering could be generalized well beyond these simple
examples, and thus provide a broadly applicable scheme to control self-assembly.

Equilibrium self-assembly is a powerful way to build nano- and microscale structures out of interacting subunits (1, 2). Many biological and technological functions, however, require these structures to have a well-controlled size (3, 4), requiring specific strategies for size-controlled assembly (5). In a first approach known as programmable assembly, a set of several distinct subunits bind with one another according to a predefined pattern of interactions (6–9). The outer surface of a complete assembly is designed not to interact with any remaining subunits, leading to size control. A second strategy, self-closing assembly, uses identical subunits that bind at an angle. A collection of many such subunits forms a ring, a cylinder or a sphere whose closure terminates the assembly process (10–14). A last pathway relies on deformable subunits whose shapes change upon binding to one another. In properly designed systems, these deformations increase as more and more subunits are added to the assembly, which eventually inhibits further growth (15–18).

Several experimental platforms enable these strategies, including DNA origami (19–21), and DNA tiles at the sub-micrometer scale (22,23), as well as larger colloidal systems interacting through attached DNA strands (24–26) or through shape recognition mediated by depletion interactions (27–30). However, the formation of large, size-controlled assemblies remains challenging and requires pushing these techniques to their limits. Using large numbers of distinct subunits in programmable self-assembly can lead to high costs and low yields (31,32). In self-closing assembly, both flexibility and small uncertainties in the subunit shapes lead to the formation of diverse, potentially off-target structures (13, 33). Finally, deformable subunits are difficult to manufacture, and their elastic properties must be tightly controlled to avoid the formation of infinite assemblies.

Here we propose an alternative route to size-controlled assembly that circumvents these limi-
tations. Our approach relies on a single type of lattice-forming subunits and does not require fine
control over their shape or elasticity. Instead, the regulation of the assembly size and shape relies
on the competition between two binding energies, which in DNA-based systems can be controlled
to a high precision (34, 35). As illustrated in Fig. 1A, the first energy is associated with subunits
organizing into a crystalline lattice. The second describes defects in the crystal, specifically grain
<figure><img src="image_2.png" /><figcaption>**Figure 1: Crystals with energetically favorable defects may not grow to arbitrarily large sizes.**</figcaption></figure>

(A) Our proposed subunits form crystals with binding free energy $\tilde{\epsilon}_c < 0$, but also favor defective grain boundaries with $\tilde{\epsilon}_d < \tilde{\epsilon}_c$. Singular points at the end of these defects generically incur a cost $\tilde{\epsilon}_p > 0$ related to $\tilde{\epsilon}_d$ and $\tilde{\epsilon}_c$. Denoting the typical size scale of the assembly by $r$, the associated free energies respectively scale as $\tilde{\epsilon}_c r^2$, $\tilde{\epsilon}_d r^1$ and $\tilde{\epsilon}_p r^0$. (B) Because these quantities scale differently with $r$, small assemblies are dominated by energetically costly singular points and large ones by crystalline interactions. The favorable defect interactions may dominate in assemblies with intermediate sizes, making them the most stable. Expressing $r$ in units of the subunit size, this qualitative result is recovered by minimizing the free energy per subunit $\approx -|\tilde{\epsilon}_c| - |\tilde{\epsilon}_d|/r + |\tilde{\epsilon}_p|/r^2$, yielding a preferred radius $r^* \approx |\tilde{\epsilon}_p|/|\tilde{\epsilon}_d|$.

boundaries between domains with different orientations. Our strategy requires that unlike in usual crystals, the defect energy be more favorable than that of the crystal. In this setting, subunits self-assemble into structures with a set number of defects imposed by topological constraints. The crystalline interaction favors the formation of large assemblies, while forming a large number of small structures maximizes the defect interaction. As illustrated in Fig. 1B, the balance of these effects results in the emergence of an optimal assembly size set by the values of the interaction energies.

## Theory of Defect-Induced Size Control

To formalize our proposal, we consider a set of identical hexagonal subunits that can bind to their neighbors through the short-range orientation-dependent interaction rules illustrated in Fig. 2A. The subunits bind through a crystalline or a defect interaction with binding (free) energies $\epsilon_c$ and $\epsilon_d$.

Tiling the plane using subunits with identical orientations gives rise to crystalline interactions everywhere. By contrast, there can be no dense tiling displaying defect interactions everywhere
(see proof in Supplementary Text and Fig. S12). As a result, a large assembly that includes defect interactions must also comprise many additional crystalline interactions.

To determine the most favorable assembly structure, we express the energy of two idealized cases. We thus consider a piece of crystal in the shape of a regular hexagon (which optimizes the crystalline energy), and a six-fold symmetric assembly whose defects impose a vortex-like pattern of subunit orientations (inset of Fig. 2B). The geometric and energetic factors favoring this defect pattern over others are discussed in the Supplementary Text and Fig. S12. To minimize the total energy of a fixed number of subunits, we equivalently minimize the average energy per subunit e. Denoting the assembly radius by r, this quantity reads (see Supplementary Text, Fig. S3 and S4):

$$e^{(\text{crystal})}(r) = \frac{9r^2 + 3r}{3r(r+1) + 1} \epsilon_c \quad (1a)$$

$$e^{(\text{vortex})}(r) = \frac{3\epsilon_c r^2 + (-3\epsilon_c + 4\epsilon_d)r - 2\epsilon_d}{r(r+1)}. \quad (1b)$$

Here the crystalline case includes isolated monomeric subunits as its r = 0 case, while the r = 1 vortex assembly is the largest possible structure that includes only defect interactions, namely a ring-like hexamer. For large (r ≫ 1) assemblies, Eq. (1b) recapitulates the discussion in the caption of Fig. 1 with $\tilde{\epsilon}_c = 3\epsilon_c$, $\tilde{\epsilon}_d = 4\epsilon_d - 6\epsilon_c$ and $\tilde{\epsilon}_p = 6(\epsilon_c - \epsilon_d)$.

The most favorable assembly as a function of the binding energies is shown in the phase diagram of Fig. 2B. Its nature depends only on the ratio and signs of $\epsilon_c$ and $\epsilon_d$. Since the crystal energy of Eq. (1a) is a monotonic function of r, crystalline ground states consist either in isolated monomers or in infinitely large crystals. By contrast, Eq. (1b) can display a minimum at a finite radius $r^*$, yielding finite-size vortex assemblies that dominate the left-hand-side of the phase diagram. For attractive crystalline interactions ($\epsilon_c < 0$), the size of these assemblies diverges as the system approaches the $\epsilon_c = 2\epsilon_d/3$ line. This suggests that any assembly size can be energetically favored under suitable conditions, which could enable size control over arbitrarily large assemblies.

To test the robustness of our conclusions against thermal fluctuations and assess the possible dominance of structures beyond our crystal and ideal vortex assembly, we numerically simulate the self-assembly process. We run a Monte Carlo simulated annealing down to a fixed finite temperature (36), and successfully recover all predicted structures in the expected range of parameters (Fig. 2C). A straightforward finite-temperature extension of our theory moreover accurately predicts the assembly size distribution (Fig. 2D).
<figure><img src="image_3.png" /><figcaption>**Figure 2: Modeling demonstrates size control in randomly assembling subunits.**</figcaption></figure>

(A) We consider hexagonal subunits with anisotropic two-body interactions visualized by *colored patches*. All relative binding orientations not pictured in this panel are strongly penalized in our model. (B) Analytical phase diagram derived from Eqs. (1). *Colored symbols (triangles, square, hexagon, star, circle) link to the nomenclature defined in panel C. Inset: vortex assembly structure considered in our analytical calculations.*

(C) Final structures of simulated assemblies for the parameter values marked by *symbols* in panel B. See Tab. S1 for detailed parameters. Long protrusions can form at the corner of the smallest assemblies, but are increasingly penalized for smaller values of $\epsilon_c/\epsilon_d$ as discussed in the Supplementary Text and Fig. S9.

(D) Distributions of final assembly sizes at the end of the simulations (*symbols*) (36) and comparison with an ideal-gas-of-cluster theory applied to Eq. (1b) (*solid lines*, see Supplementary Text). The moderate deviations observed for the largest assemblies are likely due to equilibrium entropic effects as well as nucleation kinetics (see Fig. S10).
DNA Origami Implementation

We implement our proposed six-fold subunit design using hollow DNA origami nanocylinders (37). Our subunits interact through a six-fold pattern of 16 nucleotide (nt) long single-stranded DNA linkers each comprising a 4 nt oligothymidine anchor domain at the 5' end for flexibility, followed by a 12 nt sticky domain designed to hybridize with its complementary sequence through base-pairing (Fig. 3A). The length and sequence of this sticky domain enable fine control of the hybridization free energy and ensure that it is associated with a melting temperature within the experimentally accessible range 25°C ≤ *T* ≤ 40°C. The linkers are organized in three vertical layers along the height of the cylinder. The middle layer comprises one linker in each of the six binding directions and serves to implement a crystalline interaction ε~c~ ranging from −13 to −9 kcal·mol⁻¹ (approximately −21 to −15k~B~*T* with *T* = 36°C; see Fig. 3B and (36)). We implement each defect interaction using two pairs of linkers from the top and bottom layer. This allows us to probe the range 0.33 ≤ ε~c~/ε~d~ ≤ 0.6 relevant for size-controlled assembly while retaining experimentally accessible melting temperatures.

We assemble our nanocylinders by incubating them at constant temperature *T* = 36°C for 20 hours in a buffer containing [MgCl₂] = 5 mM, resulting in structures highly reminiscent of our theoretical predictions (Fig. 3C and Fig. S14). We further investigate the orientations of individual nanocylinders by labeling them at a single position with a 5 nm DNA-grafted gold nanoparticle, and find that they are consistent with our predicted vortex assemblies (Fig. 3D). We confirm this result by studying the distribution of orientations over a large number of assemblies. We find that adding Gaussian noise with a standard deviation of *σ*~φ~ = 20° to the ideal orientations, to account for effects such as the flexibility of the tether linking the gold nanoparticle to the nanocylinder, yields almost perfect agreement with the data (Fig. 3E, Kolmogorov-Smirnov test p-value *p* = 0.65). Taken together, these results demonstrate that our DNA nanocylinders indeed form our predicted vortex assemblies.

We next quantitatively examine our size control mechanism by studying the distribution of assembly sizes (36). In Fig. 3F we find that it does not markedly depend on the incubation time after 10 hours, suggesting that they are at equilibrium. We additionally vary the MgCl₂ concentration in our sample based on the fact that an increased salt concentration increases non-specific interactions
<figure><img src="image_4.png" /><figcaption>Figure 3: DNA origami nanocylinders form size-controlled assemblies.</figcaption></figure>

(A) OxDNA molecular model of our nanocylinders. Light (dark) blue linkers provide crystalline (defect) interactions. (B) Top row: In a crystalline interaction, complementary crystalline linkers on opposite nanocylinder sides (e.g., 2 and 5) hybridize. Defect linkers remain unhybridized. Bottom row: In a defect interaction (e.g., between sides 2 and 4) two pairs of defect linkers hybridize. Sides 5 and 6 never form defect interactions and do not carry defect linkers. (C) TEM micrographs of the assemblies formed by our nanocylinders illustrate the dependence of their radius on $\epsilon_c/\epsilon_d$. White arrow: a defect-interaction-induced appendage reminiscent of our simulations (Fig. 2C; more examples in Fig. S15) (D) Labeling of side 3 with gold nanoparticles reveals subunit orientations consistent with our vortex assembly design (red arrows). Apparent defect interaction sites are outlined in dark blue ($\epsilon_c/\epsilon_d = 0.5$). (E) The empirical distribution of orientations closely resembles an ideal $r = 3$ vortex assembly structure with noise added ($\epsilon_c/\epsilon_d = 0.5$). (F) The assembly size distribution is robust to variations in incubation time. Here $[MgCl_2] = 5$ mM, $\epsilon_c/\epsilon_d = 0.5$. (G) The assembly size distribution is robust to variations in salt concentration ($\epsilon_c/\epsilon_d = 0.5$). (H) The assembly radius depends on the binding free energies as predicted by theory. See (36) for statistical information.
in DNA origami (Fig. 3G). This does not lead to a significant change of our assembly sizes, indicating that our structures are primarily based on specific hybridization interactions. Finally, Fig. 3H demonstrates that consistent with our theoretical predictions, the ratio of the crystalline to the defect hybridization energies controls the size of the assembly up to large radii, thus demonstrating our proposed size control mechanism.

## Shape Control through Defect Engineering

Beyond the design and control of vortex assemblies, we suggest that defect engineering could be a general and versatile self-assembly strategy. While its full potential remains to be explored, here we demonstrate additional possible designs and point out some of the associated challenges.

In our vortex design, defect interactions cause the subunit orientation to rotate clockwise when crossing a defect line (Fig. 2A). In Fig. 4A we instead consider defect interactions that reverse the arrow orientation. This again favors triangle-shaped crystalline domains, but as shown in the inset of Fig. 4B and Fig. S5 these form an alternating fiber-like pattern with an energy per subunit

$$e^{(\text{fiber})}(r) = \frac{3\epsilon_c r^2 + (4\epsilon_d - 3\epsilon_c)r + 2(\epsilon_c - \epsilon_d)}{r(r+1)} \quad (2)$$

As in the vortex case, this energy can display a minimum for a finite size $r$ of the triangular domains. Figure 4B shows the resulting phase diagram, which we further validate through numerical simulations in Fig. 4C. Finally, we implement this design using our six-fold DNA nanocylinder platform, and indeed observe the formation of fiber-like assemblies (Fig. 4D and S14).

Further reassignments of the defect interactions allow for multiple other size-controlled assemblies. The design of Fig. 4E thus rearranges our triangular crystalline domains into a patterned bulk with twice as many +1 as -2 defects. A variant of the vortex assembly of Fig. 2 based on square subunits displays size control but features flexible defect lines (Fig. 4F). These lines form diagonal paths on a square lattice, implying that their length depends solely on their total horizontal and vertical amplitude. This results in a high degeneracy, whereby all directed defect lines that start and end at the corners of a rectangular assembly yield the same overall energy. When binding energies are low enough to be comparable to the thermal energy $k_B T$, the entropic contribution from this conformational freedom dominates. As shown in Fig. 4F, increasing $k_B T/\epsilon_c$ while keeping $\epsilon_c/\epsilon_d$ fixed thus favors squares over elongated rectangular assemblies (Supplementary Text and
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Defect engineering provides a rich design space.</figcaption></figure>

Here subunits with identical orientations form crystalline interactions with energy $\epsilon_c$ throughout. (A) Defect interactions and conceptual schematic for fiber-like assemblies. (B) Corresponding phase diagram derived from Eqs. (1a) and (2). It differs from Fig. 2B only through the slopes of the diagonal lines. *Inset*: Structure of an ideal fiber. (C) Simulations conducted as in Fig. 2C. (D) TEM micrograph of a fiber formed as in Fig. 3 but with the new defect interactions. (E) Patterned bulk design. (F) Rectangular size-limited design where entropic effects control the aspect ratio of the assembly. (G) A non-size-controlled design. (H) Three-dimensional size control. The subunits are rhombic dodecahedral whose orientation is indicated by a *dark grey arrowhead*. These arrows point away from a hedgehog topological defect at the center of the assembly. See Tab. S1 for simulation parameters and Supplementary Text for the predictions of the optimal assembly sizes $r^*$.
Fig. S6). An alternative square design shown in Fig. 4G does not lead to size control and serves as a cautionary example of the subtleties of defect engineering. In this case the local geometry of the interactions implies a vanishing $\tilde{\epsilon}_p$ in the language of Fig. 1, and thus a function $e(r)$ that does not have a minimum at finite $r$. The condition $\tilde{\epsilon}_p > 0$, while necessary, is not sufficient to enable size control. For instance, a three-fold vortex design gives rise to alternate structures that are more stable than the intended vortex assembly (Fig. S7).

Figure 4H demonstrates defect engineering in three dimensions. The assembly is comprised of six crystalline domains, each in the rough shape of a square-based pyramid. These domains come into contact through planar grain boundaries to form a truncated cube. Individual subunits are rhombic dodecahedra with $D_4$ dihedral symmetry about an axis joining two of their vertices, and can form three non-equivalent defect interactions as illustrated. Numerical simulations demonstrate successful size control through the energy ratio $\epsilon_c/\epsilon_d$ (see also Supplementary Text).

## Discussion

The induction of energetically costly grain boundaries has been used to design materials for centuries (38). Here we propose to harness modern self-assembly techniques to instead make them energetically favorable. We find that this approach results in an energetic situation that resembles classical nucleation theory in reverse. In classical nucleation, an unfavorable boundary energy works alongside a favorable bulk energy to make assemblies of intermediate sizes, the so-called critical nuclei, unstable (39). In our approach, the boundary energy stems not from the outer surface of the assembly, but from the grain boundaries within. Because it is favorable, it produces stable, finite assemblies of arbitrarily large sizes and controllable morphologies in two and three dimensions.

Our strategy plays to the strengths of DNA nanotechnologies developed over the past decades, and could be implemented in other platforms where directional, reversible interactions are controlled quantitatively (40, 41). In contrast with other size control schemes (5), ours can be implemented with a single subunit type and does not require fine control over subunit shape and elasticity. As in other designs however, controlling the morphology of large assemblies becomes increasingly difficult as their size $r$ increases, as the minimum in the specific energy $e(r)$ becomes more and more shallow. Under such conditions the assembly nucleation kinetics can play a non-negligible
role in determining the final assembly size, which could account for the modest discrepancies
between our equilibrium theory and Monte Carlo simulations. This issue is less likely to influence
our DNA origami implementation, which benefits from comparatively long equilibration times.
Instead, the slight shift between the corresponding experimental data and our predictions may be
due to our relatively crude estimate of the subunit binding free energy. This estimate could be
refined by taking into account other (e.g., entropic) effects beyond base pairing.

A crucial requirement for our approach is the ability to favor only a few specific grain boundary
interactions, resulting in topological restrictions that prevent their proliferation throughout the bulk.
In our hexagonal subunit design, this is achieved by singling out two to four favorable interactions
out of twenty-one (42) possible defect-inducing options. The implied number of unexplored choices
hints at an enormous untapped design freedom. This freedom is further broadened when consid-
ering three-dimensional subunits (for instance, rhombic dodecahedra enable ≈ 10<sup>67</sup> distinct defect
interactions patterns – see Supplementary text), the possibility of assigning distinct energies to
different defect interactions, or an independent control of the point energy ε<sub>p</sub> through three-body
or steric interactions between subunits. Exploring this enormous design space and understanding
its structuring principles represents an exciting experimental and theoretical challenge, and could
open new avenues in nanotechnology.
References and Notes

1. S. C. Glotzer, M. J. Solomon, Anisotropy of building blocks and their assembly into complex structures. *Nat. Mater.* **6** (8), 557–562 (2007).

2. A. V. Pinheiro, D. Han, W. M. Shih, H. Yan, Challenges and opportunities for structural DNA nanotechnology. *Nat. Nanotechnol.* **6** (12), 763–772 (2011).

3. D. L. Caspar, A. Klug, Physical principles in the construction of regular viruses, in *Cold Spring Harbor Symposia on Quantitative Biology* (Cold Spring Harbor Laboratory Press), vol. 27 (1962), pp. 1–24.

4. Y. Fu, et al., Single-step rapid assembly of DNA origami nanostructures for addressable nanoscale bioreactors. *J. Am. Chem. Soc.* **135** (2), 696–702 (2013).

5. M. F. Hagan, G. M. Grason, Equilibrium mechanisms of self-limiting assembly. *Rev. Mod. Phys.* **93** (2), 025008 (2021).

6. K. Lund, Y. Liu, S. Lindsay, H. Yan, Self-Assembling a Molecular Pegboard. *J. Am. Chem. Soc.* **127** (50), 17606–17607 (2005).

7. Y. Liu, Y. Ke, H. Yan, Self-Assembly of Symmetric Finite-Size DNA Nanoarrays. *J. Am. Chem. Soc.* **127** (49), 17140–17141 (2005).

8. G. Tikhomirov, P. Petersen, L. Qian, Fractal assembly of micrometre-scale DNA origami arrays with arbitrary patterns. *Nature* **552** (7683), 67–71 (2017).

9. C. G. Evans, E. Winfree, Physical principles for DNA tile self-assembly. *Chem. Soc. Rev.* **46** (12), 3808–3829 (2017).

10. K. F. Wagenbauer, C. Sigl, H. Dietz, Gigadalton-scale shape-programmable DNA assemblies. *Nature* **552** (7683), 78–83 (2017).

11. C. Sigl, et al., Programmable icosahedral shell system for virus trapping. *Nat. Mater.* **20** (9), 1281–1289 (2021).
12. T. E. Videbæk, et al., Economical routes to size-specific assembly of self-closing structures. *Sci. Adv.* **10** (27), eado5979 (2024).

13. D. Hayakawa, et al., Geometrically programmed self-limited assembly of tubules using DNA origami colloids. *Proc. Natl. Acad. Sci. U.S.A.* **119** (43), e2207902119 (2022).

14. C. Karfusehr, M. Eder, F. C. Simmel, Self-assembled cell-scale containers made from DNA origami membranes (2024), doi:10.1101/2024.02.09.579479.

15. M. Lenz, T. A. Witten, Geometrical frustration yields fibre formation in self-assembly. *Nat. Phys.* **13** (11), 1100–1104 (2017).

16. D. M. Hall, M. J. Stevens, G. M. Grason, Building blocks of non-Euclidean ribbons: size-controlled self-assembly via discrete frustrated particles. *Soft Matter* **19** (5), 858–881 (2023).

17. N. W. Hackney, C. Amey, G. M. Grason, Dispersed, condensed and self-limiting states of geometrically frustrated assembly. *arXiv:2303.02121* (2023).

18. J. F. Berengut, et al., Self-limiting polymerization of DNA origami subunits with strain accumulation. *ACS Nano* **14** (12), 17428–17441 (2020).

19. P. W. Rothemund, Folding DNA to create nanoscale shapes and patterns. *Nature* **440** (7082), 297–302 (2006).

20. S. M. Douglas, et al., Self-assembly of DNA into nanoscale three-dimensional shapes. *Nature* **459** (7245), 414–418 (2009).

21. S. Dey, C. Fan, K. V. e. a. Gothelf, DNA origami. *Nat. Rev. Methods Primers* **1** (1) (2021), doi:10.1038/s43586-020-00009-8.

22. B. Wei, M. Dai, P. Yin, Complex shapes self-assembled from single-stranded DNA tiles. *Nature* **485** (7400), 623–6 (2012).

23. Y. Ke, L. L. Ong, W. M. Shih, P. Yin, Three-dimensional structures self-assembled from DNA bricks. *Science* **338** (6111), 1177–83 (2012).
24. S. Y. Park, et al., DNA-programmable nanoparticle crystallization. *Nature* **451** (7178), 553–556 (2008).

25. D. Nykypanchuk, M. M. Maye, D. v. d. Lelie, O. Gang, DNA-guided crystallization of colloidal nanoparticles **451**, 549 – 552 (2008).

26. C. A. Mirkin, S. H. Petrosko, Inspired Beyond Nature: Three Decades of Spherical Nucleic Acids and Colloidal Crystal Engineering with DNA. *ACS Nano* **17** (17), 16291–16307 (2023).

27. P. L. Biancaniello, A. J. Kim, J. C. Crocker, Colloidal interactions and self-assembly using DNA hybridization. *Phys. Rev. Lett.* **94** (5), 058302 (2005).

28. A. McMullen, M. Muñoz Basagoiti, Z. Zeravcic, J. Brujic, Self-assembly of emulsion droplets through programmable folding. *Nature* **610** (7932), 502–506 (2022).

29. S. Sacanna, W. T. Irvine, P. M. Chaikin, D. J. Pine, Lock and key colloids. *Nature* **464** (7288), 575–578 (2010).

30. T. Tigges, A. Walther, Hierarchical Self-Assembly of 3D-Printed Lock-and-Key Colloids through Shape Recognition. *Angew. Chem. Int. Ed.* **55** (37), 11261–5 (2016).

31. M. H. Huntley, A. Murugan, M. P. Brenner, Information capacity of specific interactions. *Proc. Natl. Acad. Sci. U.S.A.* **113** (21), 5841–5846 (2016).

32. W. M. Jacobs, W. B. Rogers, Assembly of Complex Colloidal Systems Using DNA. *Annu. Rev. Condens. Matter Phys.* **16** (1), 443–463 (2025).

33. A. Fotin, et al., Molecular model for a complete clathrin lattice from electron cryomicroscopy. *Nature* **432** (7017), 573–579 (2004).

34. J. SantaLucia, A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. *Proc. Natl. Acad. Sci. U.S.A.* **95** (4), 1460 – 1465 (1998).

35. J. N. Zadeh, et al., NUPACK: Analysis and design of nucleic acid systems. *J. Comput. Chem.* **32** (1), 170–173 (2011).

36. Materials and methods are available as supplementary material.
37. S. F. J. Wickham et al., Complex multicomponent patterns rendered on a 3D DNA-barrel pegboard. *Nat. Commun.* **11** (5768) (2020).

38. A. Williams, *The Sword and the Crucible: A History of the Metallurgy of European Swords up to the 16th Century*, vol. 77 of *History of Warfare* (Brill) (2012).

39. D. W. Oxtoby, Homogeneous nucleation: theory and experiment. *J. Phys. Condens. Matter* **4** (38), 7627 (1992).

40. L. Rossi, et al., Cubic crystals from cubic colloids. *Soft Matter* **7** (9), 4139–4142 (2011), doi:10.1039/c0sm01246g, http://dx.doi.org/10.1039/C0SM01246G.

41. M. Mayarani, J. Laurent, M. Lenz, O. du Roure, J. Heuvingh, Specific yet transient bonds between anisotropic colloids. *arXiv:2408.11569* (2024), https://arxiv.org/abs/2408.11569.

42. L. Koehler, P. Ronceray, M. Lenz, How Do Particles with Complex Interactions Self-Assemble? *Phys. Rev. X* **14**, 041061 (2024).

43. M. E. Fornace, P. N. J., N. A. Pierce, A Unified Dynamic Programming Framework for the Analysis of Interacting Nucleic Acid Strands: Enhanced Models, Scalability, and Speed. *ACS Synth. Biol.* **9** (10), 2665–2678 (2020), doi:10.1021/acssynbio.9b00523.

44. J. N. Zadeh et al., NUPACK: Analysis and design of nucleic acid systems. *J. Comput. Chem.* **32** (1), 170–173 (2011), doi:https://doi.org/10.1002/jcc.21596.

45. D. Doty, B. L. Lee, T. Stérin, scadnano: A Browser-Based, Scriptable Tool for Designing DNA Nanostructures, in *26th International Conference on DNA Computing and Molecular Programming (DNA 26)* (Schloss Dagstuhl – Leibniz-Zentrum für Informatik), vol. 174 of *Leibniz International Proceedings in Informatics (LIPIcs)* (2020), pp. 9:1–9:17, doi:10.4230/LIPIcs.DNA.2020.9.

46. J. Bohlin, M. Matthies, E. Poppleton et al., Design and simulation of DNA, RNA and hybrid protein–nucleic acid nanostructures with oxView. *Nat. Protoc.* **17** (8), 1762–1788 (2022), doi:10.1038/s41596-022-00688-5.
47. E. Poppleton, et al., oxDNA: coarse-grained simulations of nucleic acids made simple. *J. Open Source Softw.* **8** (81), 4693 (2023), doi:10.21105/joss.04693, https://doi.org/10.21105/joss.04693.

48. E. Poppleton, R. Romero, A. Mallya, L. Rovigatti, P. Šulc, OxDNA.org: a public web-server for coarse-grained simulations of DNA and RNA nanostructures. *Nucleic Acids Res.* **49** (W1), W491–W498 (2021), doi:10.1093/nar/gkab324, https://doi.org/10.1093/nar/gkab324.

49. B. Liu, J. Liu, Freezing Directed Construction of Bio/Nano Interfaces: Reagentless Conjugation, Denser Spherical Nucleic Acids, and Better Nanoflares. *J. Am. Chem. Soc.* **139** (28), 9471–9474 (2017), doi:10.1021/jacs.7b04885.

50. S. van der Walt, et al., scikit-image: image processing in Python. *PeerJ* **2**, e453 (2014), doi: 10.7717/peerj.453.

51. N. D. Mermin, The topological theory of defects in ordered media. *Rev. Mod. Phys.* **51** (3), 591–648 (1979).

## Acknowledgments

**Funding:** L.K. was supported by Ecole nationale des ponts et chaussées. P. R. acknowledges funding from the French National Research Agency (ANR-16-CONV-0001), from Excellence Initiative of Aix-Marseille University A*MIDEX and from the European Union (ERC-SuperStoc-101117322). F.C.S. was supported by the TUM Innovation Network "Robotic Intelligence in the Synthesis of Life (RISE)", which is financed through the Excellence Strategy of the German Federal Government. C.K. and F.C.S. acknowledge funding by the Max Planck School Matter to Life, supported by the German Federal Ministry of Education and Research (BMBF) in collaboration with the Max Planck Society. M.L. was supported by Marie Curie Integration Grant No. PCIG12-GA-2012-334053, "Investissements d'Avenir" LabEx PALM (Grant No. ANR-10-LABX-0039-PALM), ANR Grants No. ANR-15-CE13-0004-03, No. ANR-21-CE11-0004-02, No. ANR-22-ERCC-0004-01, and No. ANR-22-CE30-0024-01, as well as ERC Starting Grant No. 677532 and
the Impulscience® program from Fondation Bettencourt-Schueller. M.L.'s group belongs to the
CNRS consortium AQV.

**Author Contributions:** L.K., V.O.R. and A.Z. designed the numerical simulations. L.K. and V.O.R. performed the numerical simulations. L.K. performed the analytical calculations. L.K., P.R. and M.L. designed the concept of defect engineering. L.K., V.O.R., P.R. and M.L. designed the theoretical part of the work with help from A.Z., M.E., C.K., and F.C.S. planned the experiments. M.E. and C.K. conceived the concrete experimental implementation. M.E. performed all experiments with contributions by C.K. for initial experiments. M.E. analysed all experimental data. M.E. wrote all scripts for experimental data analysis, with contributions from C.K., L.K., M.E., V.O.R. C.K., P.R., F.C.S and M.L. wrote the paper.

**Competing Interests:** There are no competing interests to declare.

**Data and Materials Availability:** Github repository to run the 2D simulations: https://github.com/Soft-Biophysics-Group/VortexAssembly.git

Github repository to run the 3D simulations: https://github.com/Soft-Biophysics-Group/frusa_lattice_mc/tree/VortexAssembly3D

**This PDF file includes:**

Materials and Methods

Supplementary Text

Figures S1 to S16

Table S1

**Other Supplementary Materials for this manuscript:**

Data S1 to S4
# Supplementary Materials for
Topological Defect Engineering enables Size and Shape Control
in Self-Assembly

Lara Koehler¹,²†, Markus Eder³†, Vincent Ouazan-Reboul¹, Christoph Karfusehr³,⁴, Andrey Zelenskiy¹, Pierre Ronceray⁵, Friedrich C. Simmel³, Martin Lenz¹,⁶*

¹Université Paris-Saclay, CNRS, LPTMS, F-91405, Orsay, France

²Max Planck Institute for the Physics of Complex Systems, Dresden, Germany

³Department of Bioscience, School of Natural Sciences, Technical University of Munich, Garching, Germany

⁴Max Planck School Matter to Life, Jahnstraße 29, Heidelberg, D-69120, Germany

⁵Aix Marseille Université, CNRS, CINAM, Turing Center for Living Systems, 13288 Marseille, France.

⁶PMMH, CNRS, ESPCI Paris, PSL University, Sorbonne Université, Université Paris-Cité, F-75005, Paris, France.

*Corresponding author. Email: martin.lenz@universite-paris-saclay.fr

†These authors contributed equally to this work.

**This PDF file includes:**

Materials and Methods

Supplementary Text

Figures S1 to S16

Table S1

**Other Supplementary Materials for this manuscript:**

Data S1 to S4
<figure><img src="image_6.png" /><figcaption>**Figure S1:** We measure the assembly size *N* and surface *S* in the simulation and deduce the vortex assembly radius *r* and branch size *b*.</figcaption></figure>

(A) *N* is the number of subunits in the assembly and *S* is the number of subunit edges in contact with an empty site outside of the assembly. (B) The vortex radius *r* is the size of a triangular crystalline domain (in *dark gray*) and the total branch length *b* is the number of subunits within a branch (in *dark red*).

# 1 Materials and Methods

## 1.1 Numerical Materials and Methods

Here we detail the Monte-Carlo simulations presented in Figs. 2 and 4 of the main text. Section 1.1.1 describes the 2D simulations protocol. Section 1.1.2 presents the data analysis procedure used to extract the assembly sizes shown in Fig. 2 of the main text. Section 1.1.3 describes the 3D simulations protocol.

### 1.1.1 Simulation and Annealing Protocol

We perform simulated annealing with 2D lattice subunits interacting as described in Fig. 2A of the main text. We use a triangular lattice of 40 × 40 sites populated with 400 subunits. We perform simulated annealing by linearly increasing the inverse temperature $1/k_B T$ from 0 to 1 in 2000 increments. For each temperature increment, we perform $A = 4 \times 10^7$ Monte Carlo steps (which corresponds to $10^5$ steps per subunit). One Monte Carlo move is attempted at each step and may involve either a single subunit, or several of them. Single-subunit moves consist in either a rotation of a single subunit, or the transfer of a subunit to an unoccupied site, which we refer to as a translation. Each Monte Carlo move is accepted according to the Metropolis criterion.

The energies involved in our simulations are most naturally expressed in units of $k_B T$, while the binding free energies in the experiments are typically listed in kcal/mol. Assuming a temperature of 36°C = 309.15K as during the assembly stage in our experiments, the correspondence is 1kcal/mol = 1.628 $k_B T$/subunit. Therefore, a binding energy of -11.54 kcal/mol, which is the typical binding energy used in the experiments, corresponds to an interaction energy between two subunits of -18.7 $k_B T$. For all simulations in the size-limitation regime, we choose $\epsilon_c$ = -18.7 $k_B T$/subunits = -11.54kcal/mol.

The parameter values used in the 2D simulations simulations are shown in Tab. S1.
<table><caption>**Table S1: Simulation Parameters.**</caption><thead><tr><th>Figure</th><th>Assembly type</th><th>A</th><th># subunits</th><th>ε<sub>c</sub></th><th>ε<sub>d</sub></th><th>ε<sub>c</sub>/ε<sub>d</sub></th><th>ε<sub>f</sub></th></tr></thead><tbody><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>5 × 10<sup>8</sup></td><td>400</td><td>-18.7</td><td>-74.8</td><td>0.25</td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>5 × 10<sup>8</sup></td><td>400</td><td>-18.7</td><td>-45.9</td><td>0.407</td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>5 × 10<sup>8</sup></td><td>400</td><td>-18.7</td><td>-39.0261</td><td>0.479</td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>5 × 10<sup>8</sup></td><td>400</td><td>-18.7</td><td>-34.2305</td><td>0.546</td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-22.44</td><td>5/6</td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>1</td><td>-74.8</td><td></td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>1</td><td></td><td>10</td></tr><tr><td>Fig. 2C</td><td>6-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>1</td><td>1</td><td></td><td>10</td></tr><tr><td>Fig. 4C</td><td>Fiber</td><td>4 × 10<sup>7</sup></td><td>300</td><td>-18.7</td><td>-43.6333</td><td>0.43</td><td>50</td></tr><tr><td>Fig. 4C</td><td>Fiber</td><td>4 × 10<sup>7</sup></td><td>300</td><td>-18.7</td><td>-31.7087</td><td>0.59</td><td>50</td></tr><tr><td>Fig. 4C</td><td>Fiber</td><td>4 × 10<sup>7</sup></td><td>300</td><td>-18.7</td><td>-30.1102</td><td>0.62</td><td>50</td></tr><tr><td>Fig. 4C</td><td>Fiber</td><td>4 × 10<sup>7</sup></td><td>300</td><td>-18.7</td><td>-22.44</td><td>5/6</td><td>50</td></tr><tr><td>Fig. 4C</td><td>Fiber</td><td>4 × 10<sup>7</sup></td><td>300</td><td>-18.7</td><td>1</td><td></td><td>50</td></tr><tr><td>Fig. 4C</td><td>Fiber</td><td>4 × 10<sup>7</sup></td><td>300</td><td>1</td><td>-43.6</td><td></td><td>50</td></tr><tr><td>Fig. S7</td><td>3-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-38.9583</td><td>0.48</td><td>50</td></tr><tr><td>Fig. S7</td><td>3-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-35.9615</td><td>0.52</td><td>50</td></tr><tr><td>Fig. 4F</td><td>4-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-28.9</td><td>0.647</td><td>10</td></tr><tr><td>Fig. 4F</td><td>4-fold vortex</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-6.7</td><td>-10.3545</td><td>0.647</td><td>10</td></tr><tr><td>Fig. 4G</td><td>Non-size-controlled design</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-34.63</td><td>0.54</td><td>10</td></tr><tr><td>Fig. 4G</td><td>Non-size-controlled design</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-31.16</td><td>0.6</td><td>10</td></tr><tr><td>Fig. 4E</td><td>Patterned bulk</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-49.8667</td><td>0.375</td><td>10</td></tr><tr><td>Fig. 4E</td><td>Patterned bulk</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-26.0174</td><td>0.719</td><td>10</td></tr><tr><td>Fig. 4E</td><td>Patterned bulk</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-22.8203</td><td>0.819</td><td>10</td></tr><tr><td>Fig. 4E</td><td>Patterned bulk</td><td>4 × 10<sup>7</sup></td><td>400</td><td>-18.7</td><td>-18.7</td><td>1</td><td>10</td></tr><tr><td>Fig. 4H</td><td>3-dimensional design</td><td>6.64 × 10<sup>8</sup></td><td>1728</td><td>-15</td><td>-40.3169</td><td>0.3720</td><td>100</td></tr><tr><td>Fig. 4H</td><td>3-dimensional design</td><td>6.64 × 10<sup>8</sup></td><td>1728</td><td>-15</td><td>-31.8095</td><td>0.4716</td><td>100</td></tr><tr><td>Fig. 4H</td><td>3-dimensional design</td><td>6.64 × 10<sup>8</sup></td><td>1728</td><td>-15</td><td>-28.7625</td><td>0.5215</td><td>100</td></tr><tr><td>Fig. 4H</td><td>3-dimensional design</td><td>6.64 × 10<sup>8</sup></td><td>1728</td><td>-15</td><td>-25.8621</td><td>0.58</td><td>100</td></tr></tbody></table>

The interaction energies $\epsilon_c$, $\epsilon_d$ and $\epsilon_f$ are expressed in units of $k_B T$. In our simulations we associate a large positive “forbidden interaction” $\epsilon_f$ to all interfaces that are neither of the crystalline nor of the defect type to prevent them from occurring. A is the number of annealing steps per temperature.
### 1.1.2 Measurement of Assembly Radii in Simulations

We now show how we use the simulation results to measure the assembly size distribution. For each assembly within an annealed system, we automatically count its number of subunits N and its number of interfaces between a subunit and an empty site S, as illustrated in Fig. S1. As shown in Fig 2c of the main text, some vortex assemblies display linear protrusions containing extended grain boundaries, which we here refer to as branches. To measure the assembly radius r, we assume that all assemblies have the geometry illustrated in Fig. S1. We use Eq. (S1) to determine r and b, the number of subunits within the branches, from N and S for each assembly in the simulation as

$$N = 3r(r + 1) + b \qquad (S1a)$$

$$S = 6(2r + 1) + 2b. \qquad (S1b)$$

To obtain sufficient statistics on the assembly shapes, we repeat our simulated annealing 20 times. This yields a total of respectively 228, 223, 165 and 103 vortex assemblies, for the ratios $\epsilon_c/\epsilon_d$ listed in Fig. 2D of the main text. We then directly determine the number of assemblies of radius r, m(r), within all those simulations, regardless of the length of their branches b. The averaged radius in the simulations shown in Fig. 3H of the main text is $\sum_k km(k)/\sum_k m(k)$. The measured cumulative probabilities shown in Fig. 2D of the main text are given by

$$P_r(r) = \text{Proba}(\text{radius} \le r) = \frac{\sum_{r' \le r} m(r')}{\sum_{r'} m(r')}, \qquad (S2)$$

where r is an integer in the idealized geometry of Fig. S1, but we sum over all values extracted from the simulations, integer or not.

These measurements allow us to directly assess the impact of the binding energy ratio ($\epsilon_c/\epsilon_d$) on the assembly size. They also allow us to directly compare simulation, experimental and theoretical results in Fig. 3H.

### 1.1.3 Three-dimensional Numerical Simulations

The simulations shown in Fig. 4H are performed on an FCC lattice with 20×20×20 sites populated by 1728 three-dimensional subunits. Each subunit is a rhombic dodecahedron, a shape which on an FCC lattice tiles three-dimensional space similar to hexagons on a triangular lattice. We perform a simulated annealing by decreasing $log_{10}(k_B T)$ from 2 to 0 in 240 constant increments. At each temperature increment we perform $A = 6.64 \times 10^8$ Monte Carlo steps, i.e., 384,000 steps per subunit. Each step consists in attempting one of three single-subunit moves with equal probabilities: moving a subunit to an empty site, rotating a subunit, or moving a subunit to an empty site while rotating it. We run the simulations with a crystalline energy of $\epsilon_c = -15k_B T$, and four different values of the defect energy $\epsilon_d$ as detailed in Tab. S1.

For visualization purposes, we embed the final structure into a cubic domain, as opposed to the parallelepiped-shaped FCC cell which constitutes the actual simulation domain. This conversion is performed by rendering the FCC cell accompanied by its nearest-neighbour periodic copies, and truncating the result to a cube of the same volume as the original FCC cell. One side effect of this visualization trick is that the resulting cubic domain does not fully respect the symmetries of the original cell. As a consequence, not all subunits may be contained in the cube, and several periodic
copies of the same subunit may be included in the visualization box. We also emphasize that the
cubic domains shown in Fig. 4H do not have cubic periodic boundary conditions, but instead those
associated with the FCC unit cell.

The parameter values used in the 3D simulations simulations are shown in Tab. S1.

**1.2 Experimental Materials and Methods**

To experimentally verify our theoretical predictions and numerical findings, we implement a size-
controlled assembly process using DNA origami nanocylinders that interact via base-pairing under
predefined interaction rules. In Sec. 1.2.1, we describe the design of the DNA interaction strands
that implement these rules, as well as the oxDNA simulation used to visualize the structure of
the DNA origami subunit. In Sec. 1.2.2, we then detail the protocols for DNA origami folding,
assembly, and functionalization with gold nanoparticles, which yield the experimental data for
Fig. 3 of the main text. Finally, in Sec. 1.2.3, we outline the procedures for data acquisition and
analysis.

**1.2.1 Design of Interactions and Subunit Simulation**

Here we first explain how the DNA interaction strand sequences for defect and crystalline interac-
tions, as shown in Fig. 3B and listed in data S4, were designed and selected to ensure the precise
interaction energies (and their ratios) required for size-controlled assembly. We then describe how
the oxDNA simulation framework, integrated within oxView, was used to generate the visual
representation of the DNA origami nanocylinder’s mean structure shown in Fig. 3A.

**Interaction Strand Design with NUPACK.** Orthogonal defect and crystalline DNA interaction
strand sets with specified Gibbs free energy differences ($\epsilon_d$ and $\epsilon_c$) were designed using the
NUPACK Python package (43). Initially, exhaustive libraries of all possible sequences of lengths 10,
11, and 12 nucleotides, along with their reverse complements, were generated. For each sequence,
relevant properties such as the Gibbs free energy ($\epsilon$) of the minimum free energy (MFE) proxy
structure and the secondary structure were calculated (model: dna04, temperature: 36 °C, sodium:
50 mM, magnesium: 5 mM, concentration: 10 nM). The sequence library of appropriate length was
then filtered for sequences that exhibit no secondary structure with Gibbs free energy changes
upon binding within a narrow tolerance interval [$\epsilon_d - \delta\epsilon_d$, $\epsilon_d + \delta\epsilon_d$], centered around the intended
value $\epsilon_d$ of the defect interaction. An initial tolerance of $\delta\epsilon_d = 0.05$ kilocal / mol was applied.
If this resulted in too few sequences, $\delta\epsilon_d$ was incrementally increased by 0.1 kilocal / mol until a
sufficient number of candidates (typically a few hundred) was obtained. The same filtering criteria
were subsequently applied to the original sequence library using $\epsilon_c$ to identify potential crystalline
interaction strands. This process yielded sets of defect and crystalline interactions comprising $n_d$
and $n_c$ sequences, respectively.

To ensure intraset orthogonality, these sets were evaluated based on the bound fraction — defined
as the ratio of DNA strands that are correctly bound to their intended reverse complementary partners
to the total number of strands — if all strands in a set were mixed at a concentration of 10 nM each,
with the maximum allowed assembly size set to two. Specifically, sets were deemed orthogonal if
at least 90% of the defect interaction strands and at least 10% of the crystalline interaction strands
bound exclusively to their intended partners. Sets failing to reach these orthogonality thresholds were
excluded from further analysis. After identifying orthogonal defect and crystalline sets, all possible
pairings of one defect set with one crystalline set were generated to form combined sets comprising
$n_d+n_c$ sequences. These new combined sets then underwent an additional intraset orthogonality
check to prevent unintended interactions. Specifically, any set exhibiting binding domains of six or
more consecutive base pairs between strands that are not intended to form complementary pairs
was disregarded as potential candidate for a valid interaction strand set.

From the remaining fully processed interaction strand sets, candidates were randomly selected and manually inspected for potential issues, such as the formation of G-quadruplex structures and long single-base repetitions. Finally, these selected sets were validated using the NUPACK web interface (44). The sequence sets used in the experiments are listed in data S4.

**OxDNA Simulation.** The DNA origami topology and configuration files for the nanocylinder were exported from a scadnano (45) design file (see Fig. S17), imported into oxView (46) and manually pre-arranged, for structure relaxation. Forces used during relaxation were automatically generated in oxView, and a Molecular Dynamics (MD) relaxation step was run on GPUs (ox-serve with nanobase.org webserver; first Bussi-Donadio-Parrinello then Brownian thermostat). Next, a MD simulation was employed (oxDNA2 model (47), 1e9 steps) to simulate the monomer structure. The PDB file for rendering in ChimeraX was generated with the oxDNA analysis tools command line script (48). The monomer structure shown in Fig. 3A of the main text represents a mean structure.

**1.2.2 Experimental Protocols**

This subsection outlines the experimental protocols employed throughout this work. First, we
describe how the DNA origami nanocylinder is folded using an annealing ramp, followed by the
purification steps to remove incomplete or misfolded structures. We then detail the procedure for
assembling properly folded nanocylinders at different experimental conditions (Fig. 3F). Finally, we
explain how gold nanoparticles are functionalized and attached to infer the assemblies’ topological
structure as exemplified in Fig. 3D.

**DNA Origami Folding.** The DNA origami nanocylinders used in this study were originally developed by Wickham et al. (37). Single-stranded p2873 DNA (scaffold strand, for sequence information see data S3) was provided by Prof. Hendrik Dietz’ group (100 nM in ddH₂O). Staple strands were purchased unpurified from Integrated DNA Technologies in 1x TE buffer (10 mM TRIS, 0.1 mM EDTA, pH 8.0) at 200 µM each (for sequence information see data S3 and data S4). Folding reactions were prepared with 50 nM scaffold strand and 200 nM staple strands in 1x FoB18 buffer (5 mM TRIS, 1 mM EDTA, 5 mM NaCl, 18 mM MgCl₂). Folding reactions were annealed in a thermocycler (Eppendorf, Mastercycler nexus GX2). The used protocol was 5 min incubation at 65 °C followed by an annealing ramp of -0.1 °C per 36 s from 60 °C to 30 °C. Folded DNA origami samples were immediately collected for purification.

**DNA Origami Purification.** Samples were purified using an ultrafiltration protocol. A centrifuge (Eppendorf, Centrifuge 5425 R) was pre-heated to 32 °C containing one 2 mL eppi 1x FoB5 washing buffer (5 mM TRIS, 1 mM EDTA, 5 mM NaCl, 5 mM MgCl₂) per sample. Amicon filters (Amicon
Ultra 0.5 ml Ultracel 100k) were placed into 2 mL eppis, loaded with 500 µL pre-heated 1x FoB5 washing buffer and centrifuged (5 min, 10 kilo × g). The supernatant was discarded. Filters were loaded with 60 µL sample solution and 440 µL washing buffer and centrifuged (5 min, 10 kilo × g). The supernatant was discarded. Filters were then loaded with 500 µL washing buffer and centrifuged (5 min, 10 kilo × g). The supernatant was discarded once more. Filters were then loaded with 500 µL of 1x FoBX assembly buffer (5 mM TRIS, 1 mM EDTA, 5 mM NaCl, X=[5 mM, 7.5 mM, 10 mM] MgCl₂) with X being the MgCl₂ concentration used in the experiments as indicated in the main text. The filters were centrifuged (5 min, 10 kilo × g) one last time. Purified and buffer-adjusted samples were extracted from filters with a pipette by repeatedly aspirating and dispensing the solution to dissolve any potential pellets.

**DNA Origami Assembly.** The sample concentration was measured at 260 nm (Implen, NanoPhotometer N50) and normalized to 20 nM with assembly buffer (1x FoB5, 1x FoB7.5 or 1x FoB10). Samples were incubated in a thermocycler at 36 °C for varying amounts of time as indicated in the main text and resulting assemblies are immediately applied to TEM grids.

**Functionalization of Gold Nanoparticles.** To fabricate the DNA-grafted gold nanoparticles (DNA-AuNPs) with a freeze-thaw protocol (49), 110 µL of 90.8 nM AuNPs (Cytodiagnostics, d=5 nm) in 0.01x PBS buffer (0.1 mM) were mixed with 30 µL of HPLC purified 5'-thiol modified 20 nt poly-thymidine DNA strands (Biomers, lyophilized, resuspended to 100 µM in ddH₂O) and 860 µL of ddH₂O to an AuNP:DNA ratio of 1:300. The suspension was stored at -80 °C for 1 hour and subsequently thawed at room temperature. Functionalized nanoparticles were purified by ultrafiltration with an Amicon filter (Amicon Ultra 0.5 ml Ultracel 100k). The filter was pre-washed with 500 µL of ddH₂O, centrifuged (14 kilo × g, 5 min) and the supernatant was discarded. Next, 500 µL of DNA-AuNP sample was loaded, centrifuged (14 kilo × g, 5 min) and supernatant discarded. This step was repeated once. The sample was washed seven times by loading 500 µL of ddH₂O, centrifugating (14 kilo × g, 5 min) and discarding the supernatant. Finally, the filter was placed into a new eppi in inverse orientation and centrifuged (14 kilo × g, 5 min) to retrieve purified DNA-AuNPs in ddH₂O. The concentration of the AuNPs was then measured at 520 nm (Implen, NanoPhotometer N50) and the result of that measurement formed the basis of our dilution of the AuNPs to the standardized concentration used in our labeling experiments.

### 1.2.3 Data Acquisition and Analysis

First, we detail how the data were acquired using a transmission electron microscope to capture images of the assemblies. Next, we explain how vortex assembly sizes were determined from these TEM micrographs using the *scikit-image* Python library, thus enabling the generation of assembly size distributions and the calculation of mean assembly radii (Fig. 3H). Finally, we describe the evaluation method used to determine subunit orientation statistics for Fig. 3E.

**Acquisition of TEM Micrographs.** For TEM imaging, 5 µL of sample solution were incubated on glow-discharged (Electron Microscopy Sciences, K100X; 20 s, 35 mA, negative polarity) formvar carbon Cu400 TEM grids (Science Services) for 5-20 min. The staining solution was prepared by adding 1 µL of 5 M NaOH to 200 µL of a 2 % uranyl formate solution, 10 s vortexing and
subsequent centrifugation (2.5 min, 21 kilo × g). After sample incubation, grids were washed with 5 µL of staining solution and incubated with 20 µL of stain for 25 s. Automated imaging series were conducted with a FEI Tecnai T12 (120 kiloV) equipped with a Tietz TEMCAM-F416 camera, operated with SerialEM. Further sample imaging was conducted with a Philips CM100 (100 kiloV) operated with AmtV600.

**Determination of Assembly Sizes.** TEM image stacks were automatically pre-processed by monomer template matching, assembly recognition and calculation of properties with a custom Python script using the `scikit-image` library (50). The cropped image of a DNA origami monomer was used as a template to match and extract all monomer instances from the first two TEM images in a stack. The average of all instances was calculated and the result used as an optimized template (now adjusted to the imaging conditions of the stack). For every image in a stack, areas occupied by monomers were matched with the optimized template and used as a mask to convert the image to a binary format. Pixels contained in the mask were set to 1, otherwise they were set to 0 (see Fig. S2A). Small gaps between features (pixels set to 1) were closed by using the `skimage.morphology.closing()` method, creating connected areas representing assemblies (see Fig. S2B). Resulting assemblies were labeled with unique identifiers by classifying pixels to be in the same assembly if they neighbor one another and have the same value (see Fig. S2C).

As shown in Fig. S2C, many of the objects in our electron microscopy images are isolated subunits or small assemblies thereof. Each of these objects is assigned a label under the assembly-detection scheme described above. To set them aside from properly assembled vortex assemblies, different assembly properties (bounding box, area, perimeter, and centroid) were extracted for each labelled assembly with the `scikit.image.regionprops()` method. Assemblies were checked for the characteristic central hole, present in vortex assemblies. If hole size, hole eccentricity and monomer neighbor count were correct (a hole must have exactly six neighboring monomers), assemblies were marked as potential candidates and the number of monomers they contained was determined (Fig. S2D). Marked assemblies were individually presented with their extracted properties and their monomer count $N$ manually validated. To obtain a continuous value for the intrinsically discrete assembly size $r$ from $N$, an ideal vortex assembly was assumed and Eq. (S1)A was solved with $b = 0$ to find

$$r = \frac{1}{2} \left( \sqrt{1 + \frac{4N}{3}} - 1 \right), \qquad (S3)$$

which was used in Figs. 3F-H of the main text. Images of vortex assemblies that were identified using this protocol were saved with their corresponding assembly properties and size $r$, then used in the analysis presented in the main text. For Fig. 3G sample sizes were $N_{20h} = 67$, $N_{43h} = 137$ for assemblies at 5mM, $N_{20h} = 131$, $N_{43h} = 278$ at 7.5mM and $N_{20h} = 364$, $N_{43h} = 141$ at 10mM. For Fig. 3H sample sizes were $N_{0.33} = 106$, $N_{0.5} = 562$, $N_{0.55} = 122$ and $N_{0.6} = 39$ for the experiments, and $N_{0.25} = 228$, $N_{0.41} = 223$, $N_{0.48} = 165$ and $N_{0.55} = 103$ for the simulations. Bars in Figs. 3F and G represent Bessel corrected standard deviations.

**Determination of Subunit Orientations.** To evaluate the orientations of subunits labeled with gold nanoparticles (AuNP), assemblies were extracted from TEM image stacks as described above. Using a custom Python script, assembly images were individually presented for manual processing.
<figure><img src="image_7.png" /><figcaption>Figure S2: Automated processing of TEM micrographs with a custom Python script using the *scikit-image* library.</figcaption></figure>

(A) Individual DNA origami nanocylinders are detected by template matching. (B) Closed areas representing assemblies are created by using the `skimage.morphology.closing()` method. (C) Assemblies are labeled with unique identifiers by using the `skimage.measure.label()` method with the 2-connectivity criterion. (D) Filtering for valid vortex assembly candidates is done by checking for the existence and correct size and eccentricity of the characteristic central hole and the presence of exactly six monomers encircling the hole. If all criteria are met, the assembly is saved as candidate for manual validation. The red dot marks the central hole of one such candidate.
Images were mirrored horizontally if the majority of monomers in an assembly exhibited a counter-clockwise direction of rotation. For every AuNP labeled monomer in an assembly, a vector **x**<sub>i</sub> was drawn from the monomer center to the AuNP position. A second vector **y**<sub>i</sub> was drawn center-to-center between the central hole in the assembly and the monomer. After processing all images, vectors **x**<sub>i</sub> and **y**<sub>i</sub> were normalized and their 2D cross product ˆ**x**<sub>i</sub> × ˆ**y**<sub>i</sub> = cos φ was computed. The results are displayed in Fig. 3E of the main text. N = 1375 subunits from three independent experiments were analyzed.
# 2 Supplementary Text

## 2.1 Theory of Defect-Induced Size Control

In this section, we study the stability and size distributions of the assemblies discussed in the main text, leading in particular to the theoretical predictions shown in Fig. 2 and 4. In Sec. 2.1.1, we build the phase diagram of an idealized vortex assembly, including the thermodynamic prediction for the assembly radius. In Sec. 2.1.2 we discuss the stability of alternative assembly geometries including those presented in Fig. 4 of the main text. Finally, in Sec. 2.1.3, we discuss the influence of a finite temperature on the assembly size distribution shown in Fig. 2D of the main text.

### 2.1.1 Phase Diagram of the Vortex Assembly

Here we consider vortex assemblies of different radii $r$ as well as crystals, and determine which structure is most energetically favorable as a function of the interaction energies $\epsilon_c$ and $\epsilon_d$.

**Energy Calculations.** The energy of an assembly is determined by its number of crystalline interactions and defect interactions. We consider perfectly hexagonal assemblies with an integer radius $r$. We compute the energy of a crystalline triangular wedge, corresponding to $1/6^{th}$ of an assembly. We respectively denote by $N_c$ and $N_d$ the total number of crystalline and defect interactions in a triangular wedge. The energy of a triangular wedge then reads

$$E(r) = N_c(r)\epsilon_c + N_d(r)\epsilon_d. \qquad (S4)$$

We determine $N_c(r)$ and $N_d(r)$ from the schematics of Fig. S3 and obtain

$$N_c^{(\text{crystal})}(r) = 3r(r-1)/2 + 2r \qquad (S5a)$$

$$N_d^{(\text{crystal})}(r) = 0 \qquad (S5b)$$

$$N_c^{(\text{vortex})}(r) = 3r(r-1)/2 \qquad (S5c)$$

$$N_d^{(\text{vortex})}(r) = 2r - 1. \qquad (S5d)$$

The total number of subunits in the crystal and vortex assemblies are

$$6N^{(\text{crystal})}(r) = 3r(r+1) + 1 \qquad (S6a)$$

$$6N^{(\text{vortex})}(r) = 3r(r+1). \qquad (S6b)$$

Finally we combine Eqs. (S4-S6) to express energy per subunit in a vortex assembly and a crystal:

$$e^{(\text{crystal})}(r) = \frac{E^{(\text{crystal})}(r)}{N^{(\text{crystal})}(r)} = \frac{(3r^2 + r)\epsilon_c}{r(r+1) + 1/3} \qquad (S7a)$$

$$e^{(\text{vortex})}(r) = \frac{3r^2\epsilon_c + r(-3\epsilon_c + 4\epsilon_d) - 2\epsilon_d}{r(r+1)} \qquad (S7b)$$
<figure><img src="image_8.png" /><figcaption>**Figure S3: Assembly geometries depend on the favored interactions.**</figcaption></figure>

For both (A) crystal and (B) vortex assemblies, we show the assembly schematic. We also show all possible contacts in an interaction map, where the crystalline interaction is colored in green, the defect interaction in blue, and all the other in dark red. For each pair, the orientation of the left and right subunits is read in the boundary and column entries. We compute the energy of the assembly by counting each the occurrence of crystal and defect interactions in a region of the assembly, as well as its total number of subunits.

**Energy Minimization.** To determine size of the most stable assembly we now minimize the energies Eqs. (S7) with respect to $r$. This is equivalent to minimizing the total energy of a large canonical system containing a fixed total number of subunits. Considering relatively large assemblies, we can treat $r$ as a continuous parameter and look for the radii $r^*$ where the following derivatives vanish:

$$ \partial_r e^{(\text{crystal})}(r) = 3\epsilon_c \frac{6r^2 + 6r + 1}{(3r^2 + 3r + 1)^2} \qquad (\text{S8a}) $$

$$ \partial_r e^{(\text{vortex})}(r) = 2 \frac{(3\epsilon_c - 2\epsilon_d)r^2 + \epsilon_d(2r + 1)}{r^2(r + 1)^2} \qquad (\text{S8b}) $$

We find that the energy per crystal subunit is always monotonic. For $\epsilon_c > 0$, the most stable crystal is a monomer ($r^* = 0$, $e^{\text{crystal}} = 0$). For $\epsilon_c < 0$, a bulk crystal is most stable ($r^* = \infty$, $e^{\text{crystal}} = 3\epsilon_c$). By contrast, the energy per subunit of the vortex assembly has a minimum provided that $3\epsilon_c - 2\epsilon_d > 0$, $\epsilon_d < 0$, and $\epsilon_c > \epsilon_d$. It then reads

$$ r^* = \frac{-\epsilon_d + \sqrt{3\epsilon_d(\epsilon_d - \epsilon_c)}}{3\epsilon_c - 2\epsilon_d}, \qquad (\text{S9}) $$

When choosing the parameters to run a numerical simulation, we select an integer target assembly radius $r^*$, and determine the corresponding energy ratio by inverting Eq. (S9):

$$ \frac{\epsilon_c}{\epsilon_d} = \frac{2(r^*)^2 - 2(r^*) - 1}{3(r^*)^2}. \qquad (\text{S10}) $$

**Phase Diagram Boundaries.** For each region of parameter space, we determine whether the vortex assembly or the crystal is more stable, which allows us to construct the phase diagram shown in Fig. 2B of the main text.

* For $\epsilon_c > 0$, $\epsilon_d > 0$, vortex assemblies are unstable: their energy per subunit increases with increasing $r$, so monomers (energy 0) are favored. No self-assembly occurs (light green region in Fig. S4A). This situation is illustrated by the top curve in Fig. S4C.
<figure><img src="image_9.png" /><figcaption>Figure S4: Favorable defects result in vortex assemblies of finite equilibrium sizes for weakly attractive crystalline interactions.</figcaption></figure>

(A) $e^{(\text{vortex})}(r)$ is lower than $e^{(\text{crystal})}(r)$ in the blue regions of the phase diagram. (B) The vortex assembly energy has a minimum (black dot) at finite $r^*$ if $2\epsilon_d/3 < \epsilon_c < 0$. (C) There is no such non-trivial minimum for the crystal energy.

* For $\epsilon_c > 0$, $\epsilon_d < 0$, vortex assemblies are stable with a minimum at $r^* = 1$. They are more stable than crystal assemblies because $e^{(\text{vortex})}(1) = \epsilon_d < 0 < e^{(\text{crystal})}(0)$. Hexamers (vortices of radius 1) are the most stable geometry (dark blue region, Fig. S4A) This situation is illustrated by the top curve in Fig. S4B.

* For $\epsilon_c < 0$, $\epsilon_d < 0$ and $3\epsilon_c - 2\epsilon_d > 0$, we compare $e^{(\text{vortex})}(r^*)$ with $e^{(\text{crystal})}(r = \infty) = 3\epsilon_c$:

$$e^{(\text{vortex})}(r^*) - e^{(\text{crystal})}(r = \infty) = \frac{3(r^*)^2\epsilon_c + r(-3\epsilon_c + 4\epsilon_d) - 2\epsilon_d}{r^*(r^* + 1)} - 3\epsilon_c = \frac{-2r^*(3\epsilon_c - 2\epsilon_d) - 2\epsilon_d}{r^*(r^* + 1)}$$

Using Eq. (S9), we get

$$e^{(\text{vortex})}(r^*) - e^{(\text{crystal})}(r = \infty) = \frac{-2\sqrt{-3\epsilon_d(\epsilon_c - \epsilon_d)}}{r^*(r^* + 1)} < 0.$$

Hence, optimal-size vortex assemblies are always more stable than crystals in this region, and have a finite size (blue region in Fig. S4A). This situation is illustrated by the middle curve in Fig. S4B.

* For $\epsilon_c < 0$, $\epsilon_d < 0$, $3\epsilon_c - 2\epsilon_d < 0$, and $\epsilon_c - \epsilon_d > 0$, both assemblies extend to infinite size. Expanding at large $r$ gives

$$e^{(\text{vortex})}(\infty) - e^{(\text{crystal})}(\infty) = 3\epsilon_c + \frac{1}{r}(-3\epsilon_c + 4\epsilon_d) - 3\epsilon_c - \frac{1}{r}\epsilon_c + O(1/r^2) \\ = \frac{4}{r}(\epsilon_d - \epsilon_c) + O(1/r^2) < 0,$$

so vortices of unlimited size dominate (light blue region in Fig. S4A). This situation is illustrated by the bottom curve in Fig. S4B.

* For $\epsilon_c < 0$, $\epsilon_d < 0$, $\epsilon_c - \epsilon_d < 0$, vortices are infinite but crystals are more stable according to the previous equation (dark green region in Fig. S4). This situation is illustrated by the bottom curve in Fig. S4C.
* For $\epsilon_c < 0$, $\epsilon_d > 0$, $\epsilon_c - \epsilon_d < 0$, vortices are unstable and infinite crystals prevail (dark green region in Fig. S4). This situation is also illustrated by the bottom curve in Fig. S4C.

**Favored Assembly Size.** The most favorable votex radius $r^*$ computed in Eq. (S9) can take non-integer values, which are not compatible with the discrete geometry of Fig. S3. Here we explicitly determine the most favored integer radius of a vortex assembly as a function of $\epsilon_c/\epsilon_d$.

We consider a perfectly hexagonal assembly as in Fig. S3. The assembly size $r^*$ which minimizes the energy satisfies

$$e^{(\text{vortex})}(r^* - 1) > e^{(\text{vortex})}(r^*) \qquad (\text{S11a})$$

$$e^{(\text{vortex})}(r^* + 1) > e^{(\text{vortex})}(r^*), \qquad (\text{S11b})$$

From Eq. (S10), we observe that the non-integer $r^*$ computed above continuously increases with increasing $\epsilon_c/\epsilon_d$. We thus reason that the evolution of an integer $r^*$ proceeds by increments of 1. Determining this evolution than boils down to determining for which values of $\epsilon_c/\epsilon_d$ those increments occur.

We find that $e(r) = e(r + 1)$ for $r = 2/(2 - 3\epsilon_c/\epsilon_d)$. When this quantity is an integer, the preferred integer assembly size increases by one. We deduce the following dependence for the integer radius $r^*$ on $\epsilon_c/\epsilon_d$ that minimizes the assembly energy:

$$r^* = 1 + \text{floor}\left(\frac{2}{2 - 3\epsilon_c/\epsilon_d}\right) \qquad (\text{S12})$$

We use (S12) for the theoretical curve of Fig. 3H of the main text. Inverting this expression now allows us to determine $(\epsilon_c/\epsilon_d)_{\text{lim}}(r^*)$, the energy ratio for which assemblies of size $r^*$ become more favorable than those of size $r^* - 1$.

$$(\epsilon_c/\epsilon_d)_{\text{lim}}(r^*) = \frac{2r^* - 2}{3r^* - 1} \qquad (\text{S13})$$

These ratios determine the boundaries between the stability regions of vortex assemblies of different sizes. We plot them as gray lines in the phase diagram of Fig. 2B of the main text.

### 2.1.2 Generalization of the Defect Engineering Principle

Here we apply the approach of Fig. 2.1.1 to determine the most favorable assembly for the interaction rules discussed in Fig. 4 of the main text, as well as the three-fold design mentioned in relation to it.

**Fiber-like Assembly.** Here we compute the energy per subunit of fiber-like assemblies shown in Fig. 4A-C of the main text and demonstrate that for some parameter values the most stable assembly is a fiber with finite width. In this geometry, the triangular crystalline regions are in contact through defect interactions between subunits with opposite orientations. We derive the energy per subunit of an infinitely long fiber by determining the number of crystalline and defect interactions as a function of the crystalline domain size $r$, as shown in the inset of Fig. S5A, yielding

$$e^{(\text{fiber})}(r, \epsilon_c, \epsilon_d) = \frac{3\epsilon_c r^2 + r(4\epsilon_d - 3\epsilon_c) + 2\epsilon_c - 2\epsilon_d}{r(r+1)} \qquad (\text{S14})$$
<figure><img src="image_10.png" /><figcaption>**Figure S5: Ideal geometry of fibers, vortex assemblies and bulks.**</figcaption></figure>

For each assembly type, we show the assembly schematic, all the possible interactions in a contact map using the color coding of Fig. S3, and an ideal assembly. We count the total number of subunits $N$, the number of crystalline interactions $N_c$ and of defect interactions $N_d$ per the triangular wedge of size $r$ (colored in gray) and use them to compute the energy per subunit $e = (N_c\epsilon_c + N_d\epsilon_d)/N$. (A) fiber-like assembly, (B) patterned bulk (C) four-fold vortex assembly (D) non-size-controlled four-fold design (E) three-fold vortex assembly, (F) width-2 fiber defect.
Similar calculations to those of Sec. 2.1.1 allow to build the phase diagram of Fig. 4B of the main text from Eq. (S14). Fibers of finite width are the most favorable assembly for $2\epsilon_d < 3\epsilon_c < 0$, and their size $r^*$ is given by

$$ \frac{\epsilon_c}{\epsilon_d} = \frac{2(r^*)^2 - 2r^* - 1}{3(r^*)^2 - 2r^* - 1} \qquad (S15) $$

In the simulations in Fig. 4C of the main text, we use energy ratios corresponding to the integer values $r^* = 2, 4$ or $6$ as indicated in the figure.

**Patterned Bulk.** Here we compute the energy per subunit of the patterned bulk shown in Fig. 4E of the main text and demonstrate that for some parameter values the most stable assembly has a finite pattern size. The defect interactions leading to patterned bulks include those involved in the the 6-fold vortex assembly design, with additional interactions to connect several vortices into an infinite two-dimensional assembly. Each triangular wedge now has three defect lines. Using Fig. S5B we find

$$ e^{(\text{patterned})}(r, \epsilon_c, \epsilon_d) = \frac{3\epsilon_c r^2 - r(3\epsilon_c - 6\epsilon_d) - 3\epsilon_d}{r(r+1)} \qquad (S16) $$

This energy has a minimum for a finite $r^*$ if and only if $\epsilon_d < \epsilon_c < 0$, and the value of $r^*$ is given by

$$ \frac{\epsilon_c}{\epsilon_d} = \frac{2(r^*)^2 - 2r^* - 1}{2(r^*)^2} \qquad (S17) $$

In the simulations in Fig. 4C of the main text, we use energy ratios corresponding to the integer values $r^* = 2, 4$ or $6$ as indicated in the figure.

**Four-fold Vortex Assembly.** As stated in the main text, understanding the stability of four-fold vortex assemblies requires a study of both their energy and entropy. Here we first compute the energy per subunit of square vortices, and show that it displays a minimum at a finite size $r$. This energy minimum however remains less favorable than that of very elongated assemblies of width 3, implying that four-fold vortex assemblies should not form at $T = 0$. At higher temperatures, we show that the entropy associated with the many possible defect configurations stabilizes vortex assemblies and makes them the most favorable assemblies in terms of free energy.

We consider square subunits as in Fig. 4G of the main text. Using the schematic and expressions shown in Fig. S5C, we compute the energy per subunit as:

$$ e^{(4, \text{vortex})}(r, \epsilon_c, \epsilon_d) = 2 \frac{\epsilon_c r^2 + r(2\epsilon_d - \epsilon_c)}{(r+1)^2} \qquad (S18) $$

For $\epsilon_d < 3\epsilon_c/2 < 0$, this expression has a minimum $e^{(4, \text{vortex})}(r^*, \epsilon_c, \epsilon_d) = (4\epsilon_d^2 - 4\epsilon_d\epsilon_c + \epsilon_c^2)/[4(\epsilon_d - \epsilon_c)]$ for $r^* = (2\epsilon_d - \epsilon_c)/(2\epsilon_d - 3\epsilon_c)$. This minimum is however always less favorable than the energy per subunit of the width-3 fiber illustrated in Fig. S6A, which reads

$$ e^{(4,3-\text{fiber})}(\epsilon_c, \epsilon_d) = \epsilon_c + \frac{2}{3}\epsilon_d \qquad (S19) $$

This implies that square vortex assemblies are never stable at zero temperature.
<figure><img src="image_11.png" /><figcaption>Figure S6: Increasing temperature stabilizes square 4-fold assemblies over elongated ones</figcaption></figure>

(A) Square vortex assemblies of a given length L and width W may differ by the conformation of their four grain boundaries, the position (ℓ, w) of the central defect (blue cross) and the chirality of the central defect with the corners (left-handed or right-handed). Black arrows illustrates the orientation of the crystalline domain. The defect lines (dotted blue lines) have unconstrained positions. (B) Geometry of an infinitely long assembly of width 3. (C) Free energies per subunit $f^{(4)}$ (solid lines) and $f^{(4,3-fiber)}$ (dotted line) as a function of W for different values of L and the temperature. The legend for both plots is shown in the right panel. The zoom inset in the left panel shows that symmetric assemblies are less favored than assymetric ones.

We now consider rectangular vortex assemblies with length $L \ge 4$ and width $W \ge 4$ (Fig. S6A). The expression of Eq. (S18) generalizes to

$$e^{(4,\text{vortex})}(L, W, \epsilon_c, \epsilon_d) = \frac{(2L + 2W - 4)\epsilon_d + (2LW - 3L - 3W + 4)\epsilon_c}{LW}, \quad (\text{S20})$$

where Eq. (S18) is recovered for $L = W = r + 1$. To compute the conformational entropy associated with the grain boundaries of this assembly, we first consider the grain boundary that runs from the bottom-left corner of the assembly to the central defect as a discrete directed random walk. Each step in the walk is of unit length and proceeds either to the right or upwards. One such path is illustrated in Fig. S6A as a dotted line, with the central defect marked by a cross. The number of paths from (1, 1) to (ℓ − 1, w) respecting those conditions is given by the binomial coefficient $\binom{\ell-2+w-1}{\ell-2}$. We repeat this calculation for all four defect lines, accounting for the 2 possible connectivities between the central defect and the defect lines, and the 2⁴ possible connectivites of the defects to the corners (illustrated in the two examples of Fig. S6A). Finally we sum over all positions of the central defect and find

$$S^{(4,\text{vortex})}(L, W) = k_B \ln \left[ 32 \sum_{\ell=2}^{L-1} \sum_{w=2}^{W-1} \binom{\ell+w-3}{\ell-2} \binom{L-\ell+w-3}{L-\ell-1} \binom{\ell+W-w-3}{\ell-1} \binom{L-\ell+W-3}{L-\ell-2} \right]. \quad (\text{S21})$$

Combining Eqs. (S20) and (S21), we write the free energy per subunit of the square vortex assembly as

$$f^{(4,\text{vortex})}(L, W, \epsilon_c, \epsilon_d) = e^{(4,\text{vortex})}(L, W, \epsilon_c, \epsilon_d) - T \frac{S^{(4,\text{vortex})}(L, W)}{LW} \quad (\text{S22})$$
<figure><img src="image_12.png" /><figcaption>**Figure S7: A three-fold design forms unlimited assemblies and thin defect lines.**</figcaption></figure>

**Left:** Conceptual and subunit design of the three-fold assembly. **Right:** Simulations show the formation of thin defect lines for $\epsilon_c/\epsilon_d < 0.5$, and of unlimited vortices for $\epsilon_c/\epsilon_d > 0.5$. See details of the simulation parameters in Tab. S1.

To determine the relative stability of the square vortex assembly and the width-3 fiber, we show the free energy of $f^{(4,vortex)}/|\epsilon_c|$ for two different values of $k_B T/|\epsilon_c|$ in Fig. S6C. We compare it to the free energy per subunit of the width-3 fiber, which is equal to the expression of $e^{(4,3-fiber)}$ given in Eq. (S19) due to the fact that the fiber has no internal degrees of freedom and thus a vanishing entropy per subunit. This representation makes it apparent that the higher temperature stabilizes the vortex assemblies with respect to the fibers. Furthermore, among the possible vortex assemblies squares (i.e., assemblies with $L = W$) are favored at $k_B T/|\epsilon_c| = 0.15$, while more elongated morphologies have a lower free energy at $k_B T/|\epsilon_c| = 0.05$.

These findings are confirmed by the numerical simulation shown in Fig. 4F of the main text using the same parameters as in Fig. S6C.

**Non-Size-Controlled Design.** Here we show that the four-fold vortex assembly shown in Fig. 4G of the main text does not allow for size control. The schematic of Fig. S5D indeed yields

$$e^{(\text{non-size-controlled})}(r, \epsilon_c, \epsilon_d) = 2 \frac{\epsilon_c r - (\epsilon_d - \epsilon_c)}{r + 1} \quad (\text{S23})$$

This energy is a monotonic function of $r$ and thus cannot display a minimum for a finite and positive $r^*$. Additionally the defect lines associated with this design are straight and thus cannot give rise to the type of entropic effects discussed above.

**Three-fold Vortex Assembly.** Here we study a type of assembly mentioned but not pictured in the main text, namely three-fold vortex assemblies. Similar to the non-size controlled square assemblies above, these do not give rise to size-controlled vortex assemblies. Instead they either form unlimited assemblies or a collection of thin defect lines.

The design considered here involves hexagonal subunits with defect interactions between sub-units whose orientations differ by $2\pi/3$. These allow for the formation of ideal three-fold vortex assemblies, whose energy per subunit reads (Fig. S5E):

$$e^{(3,vortex)}(r, \epsilon_c, \epsilon_d) = \frac{3\epsilon_c r^2 - r(2\epsilon_d - \epsilon_c) - \epsilon_c - \epsilon_d}{r(r+1)} \quad (\text{S24})$$
This energy has a minimum for a finite $r^*$ if and only if $\epsilon_d < 2\epsilon_c < 0$. We however find that in this regime vortices with radii $r \ge 1$ are never as energetically favorable as the two-fiber illustrated in Fig. S5F and whose energy per subunit reads

$$e^{(2-fiber)} = \epsilon_d + \epsilon_c. \qquad (S25)$$

By contrast, in the regime $2\epsilon_c < \epsilon_d < 0$, the most favorable vortex has size $r^* \to \infty$ and is more stable than the width-2 fiber. Consistent with these estimates, the simulations of Fig. S7 show the formation of width-2 fiber for $\epsilon_c/\epsilon_d < 1/2$, and of large vortices for $\epsilon_c/\epsilon_d > 1/2$.

**Three-dimensional Assemblies.** Here we present energetic considerations to rationalize the observation of size control in the three-dimensional design shown in Fig. 4H. We show the full interaction map associated with this design in Fig. S8A. The orientations shown in this interaction map do not represent all 24 possible orientations of a rhombic dodecahedron. Instead it shows the six non-redundant orientations for a subunit with $D_4$ symmetry about an axis that goes through two opposite vertices of the rhombic dodecahedron, as is the case in this design. The geometry of the problem considered here is more complex than those of the two-dimensional systems studied above, and makes an analytical discussion more difficult. We thus instead base our discussion on numerical comparisons between the energies of a range of putative assembly structures.

The Monte-Carlo simulations of Fig. 4H typically produce assemblies with the structure of a truncated cube. We thus adopt this structure as our putative structure. We show our parametrization of the truncated cubes in Fig. S8B. The assembly radius $r$ denotes a length of the side of the cube, expressed in number of subunits. At each corner of the cube, we remove a trirectangular tetrahedron whose right angle coincides with the cube corner, and whose legs of length $l_t$ run along the cube edges. These legs must be shorter than half of the edge length, imposing $1 \le l_t \le \lfloor \frac{r}{2} - 1 \rfloor$. The assembly has six approximately pyramidal crystalline domains as illustrated in the left panel of Fig. 4H.

To generate all valid designs of radius $4 \le r \le 14$, we first compute the intersection of truncated cubes of given radii and truncation lengths with a fully-packed FCC cell. We then orient the subunits in each crystalline domain towards the corresponding aggregate outer face. Due to the mismatch between the cubic geometry of the assembly and the rhombic-dodecahedral shape of the subunits that compose it, two different assemblies are possible for a given $(r, l_t)$. These two possible morphologies differ by a translation of the truncated cube bounding the assembly by the length of one subunit in the direction of any of the cubic faces of the assembly. We take both possible morphologies associated with each $(r, l_t)$ into account.

While this simple geometrical construction recapitulates the overall crystalline domain structure observed in our Monte-Carlo simulations, it does not optimize the specific orientations of its subunits close to the grain boundaries between these domains. To fine-tune this grain boundary structure, we conduct a simulation at $k_B T = 0$ where only subunit rotations are allowed. Our simulation consists in attempting $10^5$ Monte-Carlo moves per subunit. To illustrate the outcome, we show all post-annealing putative structures of radii $4 \le r \le 8$ in Fig. S8C.

For each annealed structure, we extract the total number of crystalline interactions, defect interactions and subunits. We then compute the energy per subunit $e(r)$ as a function of $\epsilon_c$ and $\epsilon_d$ for each structure as in the caption of Fig. S5. This allows us to compare the energies per subunit of all our putative structures as a function of $\epsilon_c/\epsilon_d$ and to determine the most stable one for each
<figure><img src="image_13.png" /><figcaption>**Figure S8: Size limitation in the three-dimensional design.**</figcaption></figure>

(A) Full interaction map of our three-dimensional design. (B) Putative truncated-cubic assembly structure. (C) Putative assemblies for radii $4 \le r \le 8$, and truncations $1 \le l_t \le \lfloor \frac{r}{2} - 1 \rfloor$. We show two structures for each $(r, l_t)$ as discussed in the text. The most favorable structures are outlined by red boxes. (D) Comparison of the optimal energy per subunit (left) and (mean) assembly size (right) obtained on the one hand by comparing ideal putative structures, and on the other through the simulations of Fig. 4D.
value of this energy ratio. We find that the optimal value $r^*$ of the assembly radius depends on $\epsilon_c/\epsilon_d$ (Fig. S8D), while the optimal truncation length is always $\ell_t = 1$ (Fig. S8C, red boxes).

The simulations shown in Fig. 4H of the main text are performed at $\epsilon_c/\epsilon_d$ values where we predict the formation of truncated cubes with assembly radii of 4, 5, and 8, in addition to an extra point at large $\epsilon_c/\epsilon_d$ where we expect an optimal assembly with more subunits than are present in the simulation (see Tab. S1). We find that the energy per subunit in the idealized and annealed structures are in close agreement (Fig. S8D). The latter systematically have larger energy per subunit than the former, probably due to imperfect annealing and/or finite size effects. The corresponding assembly sizes are also in reasonable agreement (Fig. S8E).

### 2.1.3 Effect of Thermal Fluctuations on Vortex Assemblies

The discussion of Sec. 2.1.1 and the resulting phase diagram of Fig. 2B relate to the ground state of our vortex assemblies. Here we assess the effect of thermal fluctuations, which are present in both our experiments and simulations, on the assembly structure. We find that they promote the formation of branches of the type depicted in Fig. S1, but do not significantly affect the assembly radius $r$. The lines in Fig. 2D of the main text are computed using the results of this section.

**Energetic Contribution of Branches.** The simulation results presented in Fig. 2C of the main text show that some vortex assemblies deviate from the perfectly hexagonal structures considered in Sec. 2.1.1. Here we rationalize this observation by showing that forming branches is energetically favorable at low $\epsilon_c/\epsilon_d$.

We consider a hexagonal assembly with branches, and denote by $b$ the total number of subunits within these branches (shown in red in Fig. S1B). The total number of subunits, crystalline interactions and defect interactions in the whole assembly respectively read

$$N(r, b) = 3r(r + 1) + b \qquad (S26a)$$

$$N_c(r, b) = 6 \times [3r(r - 1)/2] + b \qquad (S26b)$$

$$N_d(r, b) = 6 \times (2r - 1) + b. \qquad (S26c)$$

This yields an energy per subunit

$$e^{(\text{branched vortex})}(r, b) = \frac{9r(r-1)\epsilon_c + b\epsilon_c + 6(2r-1)\epsilon_d + b\epsilon_d}{3r(r+1)+b} \qquad (S27)$$

To find the most energetically favorable values of $r$ and $b$ we first seek to minimize this expression with respect to $b$. We find

$$\partial_b e^{(\text{branched vortex})}(r, b) = -\epsilon_d(r-2) \left[ r \left( 1 - 2 \frac{\epsilon_c}{\epsilon_d} \right) - 1 \right] \qquad (S28)$$

This expression is independent of $b$, indicating that the branch length $b^*$ that minimizes $e^{(\text{branched vortex})}$ can only have either one of two values: $b^* = 0$ if $\partial_b e^{(\text{branched vortex})} > 0$ and $b^* = \infty$ if $\partial_b e^{(\text{branched vortex})} < 0$. If $b^* = 0$, $e^{(\text{branched vortex})}(r, 0) = e^{(\text{vortex})}(r)$ and the energy minimum is found for $r = r^*$ [Eq. (S9)]. If $b^* = \infty$, the energy per subunit is independent of $r$ and equal to the energy $e^{(2-fiber)}$ of an infinite fiber of width 2 [Eq. (S25)]. Comparing $e^{(\text{vortex})}(r^*)$ and $e^{(2-fiber)}$,
we find that the former energy per subunit is lower than the latter as long as $\epsilon_c/\epsilon_d > 1/4$. Therefore for $\epsilon_c/\epsilon_d > 1/4$ branch formation is energetically suppressed.

In the case $\epsilon_c/\epsilon_d = 1/4$, $e^{(\text{branched vortex})}$ does not have a single minimum but a line of degenerate minima with $r = 2$ and arbitrary $b$. This implies that branch formation is not energetically suppressed, and thus rationalizes the presence of branches in the leftmost simulations panel of Fig. 2C of the main text. Smaller values of $\epsilon_c/\epsilon_d$ are not well described by our treatment of $r$ as a continuum variable and are not discussed here.

**Assembly Size Distribution.** Thermal fluctuations can result in assemblies with radii other than $r^*$ or the formation of branches even in cases where they are not favored. Here we determine the expected distribution of assembly and branch sizes.

We consider an ideal gas of non-interacting assemblies. We denote by $m(r, b)$ the number of assemblies of radius $r$ with branches containing a total of $b$ subunits. We denote the number of unassembled monomeric subunits as $m(0, 0)$. Due to the lack of interactions between assemblies, the partition function of the total system $Z$ factorizes as

$$Z = \prod_{r,b} \frac{1}{m(r,b)!} [z(r,b)]^{m(r,b)}, \qquad (S29)$$

where the partition function $z(r, b)$ of an individual assembly reads

$$z(r, b) = N_{\text{sites}} \times N_{\text{config}}(b) \times \exp[-\beta N(r, b) e^{(\text{branched vortex})}(r, b)]. \qquad (S30)$$

Here $N_{\text{sites}}$ is the number of accessible sites in the whole system, $N_{\text{config}}(b)$ is the number of possible ways in which a vortex assembly can be decorated by a total of $b$ branch monomers, and $\beta = 1/k_B T$ is the inverse temperature. The energy for a monomer is $e^{(\text{branched vortex})}(0, 0) = 0$. Using the Stirling approximation and introducing the concentration of $(r, b)$-assemblies $c(r, b) = m(r, b)/N_{\text{sites}}$ we obtain a free energy per unit volume

$$f = -\frac{\ln Z}{N_{\text{sites}}\beta} = \sum_{r,b} N(r,b)c(r,b) \left[ \frac{1}{N(r,b)\beta}(\ln c(r,b) - 1) + N_{\text{config}}(b)e^{(\text{branched vortex})}(r,b) \right]. \qquad (S31)$$

We minimize $f$ with respect to the concentrations $c(r, b)$ while holding the total number of subunits constant by setting

$$\frac{\partial}{\partial c(r, b)} \left[ f - \mu \sum_{r,b} N(r, b) c(r, b) \right] = 0 \qquad (S32a)$$

$$\sum_{r,b} N(r,b)c(r,b) = \frac{N_{\text{subunits}}}{N_{\text{sites}}}, \qquad (S32b)$$

where the thus introduced Lagrange multiplier $\mu$ is the subunit chemical potential. This yields

$$c(r, b) = \left[ c_0 N_{\text{config}}(r, b) e^{-\beta e^{(\text{vortex branched})}(r, b)} \right]^{N(r,b)}, \qquad (S33)$$

where $c_0 = c(0, 0) = e^{-\beta\mu}$ denotes the concentration of the monomeric subunits.
<figure><img src="image_14.png" /><figcaption>Figure S9: Distribution of theoretically predicted assembly radius and branch size distributions.</figcaption></figure>

Branches do not form for assemblies with equilibrium radii $r^* > 2$, and their size distribution is determined by entropic effects for $r^* = 2$.

To determine the number of configurations $N_{config}(b)$, we consider all possible ways in which the b undistinguishable monomers can be partitioned into six groups representing the six distinguishable branches of a branchy vortex assembly. To illustrate this choice of a partition, consider a group of b = 5 monomers. These can be distributed by assigning to the ith branch a number of monomer equal to $v_i$, where the vector v can for instance be given by $v = (2, 1, 0, 0, 2, 0)$. We reason that choosing such a vector with $\sum_i v_i = b$ is equivalent to choosing 5 sites out of a set of $b + 5$. Indeed, consider a set of $5 + 5 = 10$ empty squares in our example:

□□□□□□□□□

then choose 5 of them and color them black:

□□■■■■■□□

In this example, the number of white squares between the beginning of the list and the first black square is 2, implying $v_1 = 2$. Counting white squares between the first and the second black square yields $v_2 = 1$, etc. until the whole vector v given above is recovered. There are $\binom{b+5}{5}$ ways of picking the black squares, and therefore

$$N_{config}(b) = \binom{b+5}{5} \qquad (S34)$$

Finally, we determine $c_0$ numerically by solving Eq. (S32b).

We plot the thus determined distributions in Fig. S9 using the parameters of our main text simulations, and find that our formalism predicts very little branching except in the case $r^* = 2$, as well as radii distributions that are highly peaked around the predicted ground-state $r^*$. Finally, we predict the cumulative probability distributions for the assembly radius through

$$P_r(r) = \text{Proba}(\text{radius} \le r) = \frac{N_{\text{sites}}}{N_{\text{subunits}}} \sum_{r'=0}^{r} \sum_{b=0}^{+\infty} c(r', b), \qquad (S35)$$

and compare the result to our simulations by applying Eq. (S2) in Fig. 2D of the main text.
<figure><img src="image_15.png" /><figcaption>Figure S10: Convergence of simulated assembly size toward theoretical prediction with increasing annealing time.</figcaption></figure>

(A) The assembly radius increases with increasing annealing time in large assemblies. (B) Upon sufficient equilibration, the assembly radius distribution defined in Eq. (S2) is well predicted by the theory. (C) A similarly defined branch size distribution also agrees with theory.

**Enhanced Theory-Simulation Agreement via Extended Annealing.** In Fig. 2D and 3H of the main text, we observe a small discrepancy between the predicted and simulated assembly sizes. Here, we show that increasing the annealing time steadily improves agreement with theory, highlighting the reliability of the prediction despite computational constraints. We vary the annealing time $A$ for the binding energy ratios $\epsilon_c/\epsilon_d$ discussed in the main text, and track how the average assembly radius and the distributions of radii and branch sizes evolve. As shown in Fig. S10A, the assembly radius converges quickly toward the predicted equilibrium value $r^*$ for $r^* = 2$ and $r^* = 3$. For larger target sizes ($r^* = 4$ and $r^* = 6$), equilibration requires significantly more computational effort, explaining the small discrepancy observed in Fig. 3H. In panels B and C of Fig. 3, we display the distributions of assembly radii and branch sizes across different annealing protocols, compared to the theoretical predictions derived from Eq. (S35), and confirm that assemblies are not fully equilibrated for $r^* = 4$ and $r^* = 6$. Longer annealing times lead to closer agreement between simulated and theoretical assembly properties. The qualitative agreement of our simulation results with our predictions even in cases of imperfect equilibration highlights the robustness of our defect engineering concept.
## 2.2 Estimate of the Number of Possible Defect-Engineered designs

Here we provide a rationale for the assertion made in the main text that our framework can generate a number of designs of order $10^{67}$. Specifically we show that this number is larger than $8 \times 10^{66}$. Section 2.2.1 introduces our reasoning in the simpler case of hexagonal subunits in two dimensions, and Sec. 2.2.2 derives the estimate for rhombic dodecahedral subunits in three dimensions. We discuss the relevance of these estimates in Sec. 2.2.3.

### 2.2.1 Hexagonal Subunits in two dimensions

**Size of the Interaction Map.** Our reasoning is most easily visualized in the case of two-dimensional hexagonal subunits. Consider the interaction maps shown in Fig. S3. Since each subunit has 6 faces, these maps are $6 \times 6$ matrices whose entries represent the energies associated with each possible face-to-face interaction. In the language of group theory, the set of all orientations of the subunits is obtained by letting an element of the cyclic group of order 6, denoted by $C_6$ in Schönflies notation, act on a reference subunit. This group is of order $n = |C_6| = 6$, resulting in an interaction map of size $n \times n$.

**Number of Independent Interactions.** The interaction between two subunits is unchanged upon an exchange of the two subunits followed by a 180°-rotation of each of the subunits, which is equivalent to subjecting the whole pair of subunits to a 180° solid rotation without unbinding them [Fig. S11A]. As a result, the interaction map is a symmetric matrix and thus has only $n(n+1)/2 = 21$ independent entries. The set of these entries may be visualized as the upper triangular submatrix of the interaction map. Out of these 21 entries, $c = 3$ are dedicated crystalline interactions (green squares in Fig. S3B), leaving $n(n+1)/2 - c = 18$ other interactions. We reason that designing a subunit is equivalent to choosing which of these remaining interactions are repulsive (red squares), and which are defect interactions with energy $\epsilon_d$ (blue squares). Since this boils down to one binary choice per remaining interaction, this reasoning implies that there are $\mathcal{D}_0 = 2^{n(n+1)/2-c}$ possible hexagonal subunit designs.

**Number of Non-Redundant Designs.** Some of the designs enumerated above are redundant. For instance, the hexagonal subunit presented in Fig. S3B where each side $i$ has a crystalline interaction with side $(i + 3 \mod 3)$ and where side 1 has a defect interaction with side 3 and side 2 has a defect interaction with side 6, is equivalent to a subunit with the same crystalline interactions and where side 2 has a defect interaction with side 4 while side 3 has a defect interaction with side 1. Indeed, the latter is identical to the former up to a relabeling of its sides; or equivalently the latter can be obtained from the latter through a rotation by 60°. This multiple counting does not apply to all subunits however: for instance, the rotationally symmetric design of Fig. S3A is counted only once in $\mathcal{D}_0$. A similar reasoning applies to reflections: since two subunits that are mirror images of each other give rise to the same self-assembly behavior (up to a reflection of the whole system), we deduce that chiral designs are counted twice in $\mathcal{D}_0$. By contrast, achiral designs are counted only once. Overall, whether of not a specific subunit design gives rise to multiple counting or not depends on the symmetry group of this subunit.

To provide an estimate of the number of independent subunit designs, we reason that the largest number of such multiple countings is the order of the symmetry group of the hexagon. This group
<figure><img src="image_16.png" /><figcaption>**Figure S11: Symmetries in subunits interactions.**</figcaption></figure>

(A) The energy of an interacting pair of two-dimensional hexagons must be invariant under a rotation by 180° around point O, which exchanges the contacting faces and modifies the subunits' orientations. The red marking is a cue in visualizing the rotation. (B) The energy of an interacting pair of three-dimensional rhombic dodecahedra must be invariant under rotations by 180° around both around the x and y axes. The former flips the orientations of the two contacting faces, while the latter exchanges the faces themselves while modifying the subunits' orientations. The red marking provide visual help in visualizing the rotations.

is the six-fold dihedral group $D_6$ and has $|D_6| = 12$. This upper bound implies a lower bound of the number $\mathcal{D}$ of non-redundant design through

$$ \mathcal{D}^{(\text{hexagons})} \geq \frac{\mathcal{D}_{0}^{(\text{hexagons})}}{|D_6|} = \frac{2^{n(n+1)/2-c}}{|D_6|} = \frac{262\,144}{12} > 21\,845. \quad (\text{S36}) $$

Therefore there are more than 21 845 non-redundant defect engineering designs for hexagonal subunits with crystalline interactions and one defect energy $\epsilon_d$.

### 2.2.2 Rhombic Dodecahedral Subunits in three dimensions

**Size of the Interaction Map.** We now adapt this reasoning to the case of rhombic dodecahedral subunits in three dimensions. The orientations of the rhombic dodecahera are generated by the action of the chiral octahedral symmetry group O with $n = |O| = 24$, implying a most general interaction map with size $n \times n$. We note in passing that the interaction map of the specific rhombic dodecahedral subunits considered in the main text affords a much more compact representation than this general case. Indeed, these special four-fold symmetric rhombic dodecahedral subunits have a rotational symmetry group $C_4(d)$, where d denotes their axis of symmetry, and an associated $n' = |O/C_4(d)| = 6$, consistent with the 6 × 6 interaction map in Fig. S8.

**Number of Independent Interactions.** Going back to the case of generic dodecahedral subunits, we note that the interaction between two subunits lined up on the x axis is invariant both under a rotation of 180° about the x-axis, and under a rotation of 180° about the y-axis [Fig. S11B]. The former rotation maps each of the two interacting faces on itself, but with a reversed orientation. The latter rotation exchanges the contacting faces of the two subunits without modifying their
orientation. The invariance of the interaction under these two transformations dictates the form of the interaction map. We write the 24 × 24 interaction map in such a way that its first 12 coordinates refer to 12 different faces, and that the next 12 recapitulate the same faces with opposite orientations

$$ \text{Interaction map} = \begin{array}{c} \mathbf{I}_{(11)} \quad \mathbf{I}_{(12)} \\ \mathbf{I}_{(21)} \quad \mathbf{I}_{(22)} \end{array} \qquad (S37) $$

The symmetry under rotation around x implies that the interaction map must be invariant under an exchanges of the faces and so that $\mathbf{I}_{(11)} = \mathbf{I}_{(22)}$ and $\mathbf{I}_{(12)} = \mathbf{I}_{(21)}$. The symmetry upon rotation around y functions in a similar fashion as the one discussed in Sec. 2.2.1 and implies a symmetric interaction map: $\mathbf{I}_{(11)} = \mathbf{I}_{(11)}^T$, $\mathbf{I}_{(12)} = \mathbf{I}_{(12)}^T$ and $\mathbf{I}_{(22)} = \mathbf{I}_{(22)}^T$. Taken together, these constraints imply that there are only $2(2+1)/2 \times (n/2)(n/2+1)/2 = 3n(n+2)/8 = 234$ independent entries in the interaction map. Subtracting $c = 6$ crystalline interactions then leaves $3n(n+2)/8 - c = 228$ potential defect interactions, implying $\mathcal{D}_0 = 2^{3n(n+2)/8-c}$ possible rhombic dodecahedral subunit designs.

**Number of Non-Redundant Designs.** As in Sec. 2.2.1, this number $\mathcal{D}_0$ counts equivalent subunit designs multiple times. The largest possible number of multiple countings is given by the order of the full octahedral group $O_h$. As in the previous section we obtain a lower bound of the number $\mathcal{D}$ of possible designs by dividing $\mathcal{D}_0$ by the largest possible number of multiple counts $|O_h| = 48$. This yields

$$ \mathcal{D}^{(\text{rhombic dodecahera})} \ge \frac{\mathcal{D}_{0}^{(\text{rhombic dodecahera})}}{|O_h|} = \frac{2^{3n(n+2)/8-c}}{|O_h|} > 8 \times 10^{66}. \qquad (S38) $$

Therefore there are more than $8 \times 10^{66}$ non-redundant defect engineering designs for rhombic dodecahedral subunits with crystalline interactions and one defect energy $\epsilon_d$.

### 2.2.3 Relevance of these Bounds

The lower bounds computed above overestimate the number of possible designs due to the fact that some of them have additional symmetries (e.g., reflection symmetry) compared to the most generic case. Given the complexity of our interaction maps, we however expect that such cases are a small minority of all possible designs and thus that our lower bounds are quite tight.

We suspect that an overwhelming majority of possible defect engineered designs discussed do not lead to assembly size control or any other useful feature. Nevertheless, the enormous magnitude of $\mathcal{D}$ suggests that setting them aside might still leave a large number of interesting designs to be exploited. We believe that this is enough to motivate a search for a deeper understanding of the principles underlying defect engineering.
## 2.3 Geometrical Considerations

Here we use theory to demonstrate constraints on the spatial arrangement of defects in our assemblies. We specifically prove two key implications of the interaction rules pictured in Fig. 2A of the main text. First, that a large compact assembly must contain a vanishingly small fraction of defect interactions. Second, that alternatives to our six-fold vortex assembly topology are either geometrically forbidden or energetically unstable.

We build our formalism over the first two subsections. In Sec. 2.3.1, we map the discrete description of the orientations of our subunits to an equivalent field formalism. We then express the assemblies' defect energy in this formalism in Sec. 2.3.2. In Sec. 2.3.3, we use our formalism to derive an upper bound for the number of defect interactions in an assembly, which demonstrates our first point. This upper bound forms the basis of the scaling reasoning presented in the caption of Fig. 1 of the main text. Section 2.3.4 then proves our second point by establishing that only some assembly topologies are allowed, and by discussing why vortex assemblies are energetically favored among the permissible topologies.

For the purposes of this section it is convenient to use a different convention for our interaction energies than that described in the main text. The energy used in the main text reads

$$E = N_c \epsilon_c + N_d \epsilon_d, \qquad (S39)$$

where $N_c$ and $N_d$ respectively denote the number of crystalline and defect interfaces, and $\epsilon_c$ and $\epsilon_d$ their energies. We also introduce the number $S$ of interfaces between a subunit and an empty site (represented in orange in all figures). Denoting the total number of subunits in our system as $N$, and recalling that each subunit has six sides, we write the total number of sides in our system as $6N$. That number is moreover also equal to $2N_c + 2N_d + S$. Subtracting the constant $3N$ from Eq. (S39) and using this equality, we obtain the physically equivalent form

$$E = N_d(\epsilon_d - \epsilon_c) - S\frac{\epsilon_c}{2}, \qquad (S40)$$

which boils down to assigning an energy 0 to crystalline interfaces, an energy $\epsilon_d - \epsilon_c$ to defect interfaces and an energy $-\epsilon_c/2$ to interfaces at the surface of the assembly. We adopt this point of view throughout this section.

### 2.3.1 Interpolation of a Field Theory

Here we map the discrete description of our assemblies used in the main text to a continuous field theory. This enable our subsequent use of some convenient tools of vector calculus. Consider a connected assembly of oriented hexagonal subunits on a triangular lattice. In this assembly, interfaces between neighboring subunits are either of the crystalline type or of the defect type illustrated in Fig. 2A of the main text. We denote the length of a subunit edge by $a$.

We define an orientation field $\tilde{\theta}(\mathbf{x})$ at every point $\mathbf{x}$ within the assembly through the procedure illustrated in Fig. S12A. In the procedure, points within a subunit are assigned the orientation of the subunit unless they are within a small boundary region of width $d$ centered around a defect interface. Within this boundary region, the field instead interpolates between the orientations on either side by rotating by an angle $\pi/3$. We consider the $d \ll a$ limit, and within it the detailed form of the interpolating function does not matter for our discussion. We however assume that it is designed in
<figure><img src="image_17.png" /><figcaption>**Figure S12: Geometrical constraints on the presence of defect interactions in an assembly.**</figcaption></figure>

(A) Interpolation of a smooth orientation field $\tilde{\theta}(\mathbf{x})$ (small arrows) from the subunit orientations (large black arrows, with the orientation convention of the main text). The arrows within a neighborhood $d$ of a defect interaction interpolate between the subunit orientations on either side of it (dark blue). (B) Integration domain $\omega$ (dark red) used to relate the discrete and continuum expressions for the defect energy. Here the subunit orientation is pictured with the convention introduced in Sec. 2.3.1. (C) Coarse-graining the corrugated defect boundary of panel B to a straight line produces two possible cases for the crossing of this boundary (dark blue) by a locally straight, oriented curve (dashed red arrow). Cases 1 and 2 respectively correspond to crossing that boundary from right to left, and from left to right. The dashed gray line is the normal to the boundary, and we show the angle $\phi - \theta$ between the orientation of the tangent to the curve and the orientation of the subunits in one of the crystalline domains. (D) Illustration of a convex loop (dashed red line) used to evaluate the winding number of the subunit orientation. Blue lines show the coarse-grained straight grain boundaries. (E) Illustration of the maximum possible increase of the subunit orientation $\theta$ upon one turn around the loop of panel D. (F) Example of maximum possible decrease of the subunit orientation $\theta$ upon one turn around the loop of panel D. (G) Construction of the grain boundary geometry surrounding a topological defect with winding number $-1$. Inset: local structure of the inner boundary of the $-1$ defect. We observe that the tangent vector of every fourth edge of the orange boundary is orthogonal to the subunit orientation, while the others make an angle $\pm\pi/6$ with it. (H) Examples of assemblies with many topological defects. We recall that a contour that encloses several defects has a winding number equal to the sum of the winding numbers of the individual defects (51). From left to right: convex assembly with winding number +1 (the outer boundary encloses five +1 defects and four -1 defects), assembly whose outer boundary is associated with a large positive winding number (the outer boundary encloses five +1 defects and no -1 defects), and assembly with a large negative winding number (the outer boundary encloses no +1 defects and four -1 defects).
such a way that (i) the function $\tilde{\theta}(\mathbf{x})$ is differentiable everywhere within the assembly and (ii) That the field monotonically interpolates between the orientation of two neighboring subunits while choosing the shortest path on the unit circle to do so. For instance, we assume that the interpolation between $\tilde{\theta} = 0$ and $\tilde{\theta} = \pi/3$ proceeds through a monotonic increase over the interval $0 \le \tilde{\theta} \le \pi/3$ as opposed to a decrease over $2\pi \ge \tilde{\theta} \ge \pi/3$.

Within this continuous formalism, we denote the set of points belonging to the assembly (i.e., located within any of the hexagonal subunits of the assembly) as $\Omega$. The orientation within $\Omega$ can be equivalently parametrized by any field $\theta(\mathbf{x}) = \tilde{\theta}(\mathbf{x}) + \theta_0$, where $\theta_0$ is an arbitrary constant. In the following we choose an offset $\theta_0 = -2\pi/3$ compared to the convention of the main text. We define the unitary orientation vector field as

$$ \hat{\mathbf{s}}(\mathbf{x}) = \begin{pmatrix} \cos \theta(\mathbf{x}) \\ \sin \theta(\mathbf{x}) \end{pmatrix}. \qquad (S41) $$

In this new convention, the unit vector field $\hat{\mathbf{s}}$ is tangent to the overall direction of the outer boundary of the ideal vortex assembly shown in the inset of Fig. 2B of the main text.

### 2.3.2 Continuum Expression for the Defect Energy

Here we show that the assembly defect energy takes a convenient form in our field theory. Using our orientation field formalism, we thus propose to write the defect energy of the assembly as

$$ E_d = \frac{2(\epsilon_d - \epsilon_c)}{\sqrt{3}a} \iint_{\Omega} \nabla \times \hat{\mathbf{s}} \, dA, \qquad (S42) $$

where $dA$ denotes the area element. To verify that this expression matches the discrete definition of the main text to lowest order in $d/a$, we first note that the value of its integrand vanishes everywhere except within our boundary regions of thickness $d$. We then compute its value for the representative portion of grain boundary represented in Fig. S12B. We denote the thus-defined domain by $\omega$ and use Stokes' theorem:

$$ E_d = \frac{2(\epsilon_d - \epsilon_c)}{\sqrt{3}a} \iint_{\omega} \nabla \times \hat{\mathbf{s}} \, dA = \frac{2(\epsilon_d - \epsilon_c)}{\sqrt{3}} \oint_{\partial\omega} \hat{\mathbf{s}} \cdot \hat{\mathbf{t}} \, d\ell, \qquad (S43) $$

where $\hat{\mathbf{t}}$ is the tangent vector to domain $\omega$ and $d\ell$ is the associated length element. To lowest order in $d$, the line integral of Eq. (S43) can be approximated by a sum of straight segments which we each label by $i$, yielding

$$ E_d = \frac{2(\epsilon_d - \epsilon_c)}{\sqrt{3}a} \sum_i a \hat{\mathbf{s}}_i \cdot \hat{\mathbf{t}}_i \qquad (S44) $$

The example of Fig. S12B displays the two possible types of defect interfaces associated with our design. It yields $\hat{\mathbf{s}}_1 \cdot \hat{\mathbf{t}}_1 = \hat{\mathbf{s}}_3 \cdot \hat{\mathbf{t}}_3 = 0$ and $\hat{\mathbf{s}}_2 \cdot \hat{\mathbf{t}}_2 = \hat{\mathbf{s}}_4 \cdot \hat{\mathbf{t}}_4 = \sqrt{3}/2$. Combining this result with Eq. (S44) implies that each interface brings a contribution $\epsilon_d - \epsilon_c$ to the defect energy, thus validating the expression proposed in Eq. (S42).
### 2.3.3 Upper Bound on the Number of Defect Interactions

We now use the expression of Eq. (S42) to establish a bound on the assembly defect energy. Since as discussed in Eq. (S40) this energy can alternatively be written as $E_d = N_d(\epsilon_d - \epsilon_c)$, we write

$$N_d = \frac{2}{\sqrt{3}a} \iint_{\Omega} \nabla \times \hat{s} \, dA = \frac{2}{\sqrt{3}a} \oint_{\partial\Omega} \hat{s} \cdot \hat{t} \, d\ell, \quad (S45)$$

where the second equality again uses Stokes' theorem. In Eq. (S45), the line integral runs along the outer perimeter of the whole assembly in the counterclockwise direction, and in the clockwise direction along the surface of any internal holes.

The perimeter of our assembly is comprised of a collection of subunit edges, and the value of the tangent vector $\hat{t}$ thus always matches the discrete orientations of these edges. Let us consider an individual subunit, e.g., in Fig. S12B. It is clear that relative to the orientation vector $\hat{s}$ of this subunit, the orientation of the edge vector $\hat{t}$ can only be $\pm\pi/6$, $\pm\pi/2$ or $\pm 5\pi/6$. We thus always have $\hat{s} \cdot \hat{t} \le (\hat{s} \cdot \hat{t})_{\max} = \sqrt{3}/2$. Using Eq. (S45), this implies

$$0 \le N_d \le \frac{2}{\sqrt{3}a} \oint_{\partial\Omega} (\hat{s} \cdot \hat{t})_{\max} \, d\ell = P. \quad (S46)$$

Here $P$ is the perimeter of the assembly counted in units of the subunit edge length. This perimeter includes the length of the surface of any internal holes. In an assembly where this upper bound would be reached, the defect energy would plays the same role as a negative surface tension (assuming $\epsilon_d < \epsilon_c$).

In the system considered in Figs. 2 and Fig. 3 of the main text, the case that comes closest to the upper bound of Eq. (S46) is that of the ideal vortex assembly. There the dot product $\hat{s} \cdot \hat{t}$ is maximized at every site of the outer perimeter of the assembly, but not in the internal hole, resulting in

$$N_d = P - 12. \quad (S47)$$

This offset between the actual maximal value of $N_d$ and the upper bound of Eq. (S46) results from the microscopic details of our model, and turns out to be crucial in allowing size-controlled assembly. This can be seen in the example of the non-sized-controlled design of Fig. 4H of the main text. There, the maximum number of defect interfaces is exactly equal to $P/2$, and therefore proportional to $P$. The defect energy thus acts exactly as a correction to the assembly surface tension. As a result, the only transition in the assembly ground state morphology is between an absence of growth (for a positive effective surface tension) to unbounded growth (for a negative effective surface tension).

Equation (S46) implies a severe restriction on the number of defect interactions allowed in a large assembly. In an infinite bulk without an extensive number of holes this result imposes a vanishingly small fraction of defect interactions, as stated in the main text. In addition, in assemblies characterized by a single length scale such as the one of Fig. 1A of the main text, the number of defect interactions must at most scale like the lateral size of the assembly rather than its area. This statement forms the basis of the scaling reasoning presented in the caption of Fig. 1A of the main text.
### 2.3.4 Topological and Energetic Constraints Resulting from the Interaction Rules

As implied by its name, a vortex assembly is centered around a single topological defect with winding number +1. One might speculate that the vectorial nature of our order parameter could allow for defects with any other integer winding number. However we do not observe them in simulations, nor do we reproducibly observe single assemblies including more than one vortex. Here we discuss the geometrical and energetic constraints that account for these observations. This additionally allow us to rationalize why our assemblies never take the form of a bulk with an extensive numbers of holes.

**Convex Real-space Loops Can Have Winding Numbers $-1, 0$ or $1$.** We first demonstrate a geometrical limitation imposed by our choice of interactions between subunits. For simplicity we present a coarse-grained argument where the defective grain boundaries between domains with different orientations $\theta$ are straight lines instead of straight-on-average zigzagging lines. As illustrated in Fig. S12C, in that setting the orientations $\theta$ of the subunits on either side of the boundary make angles $\pm\pi/6$ with the normal to the boundary, or $\pm 5\pi/6$ depending on the conventional choice of the direction of this normal.

We consider a differentiable closed loop fully comprised within our assembly (Fig. S12D). We parametrize it by its curvilinear coordinate $u \in [0, \ell]$ where $\ell$ is the total length of the loop. We choose $u$ in such a way that it increases when traveling along the loop in the counterclockwise direction, and denote the angle of its tangent vector $\hat{t}(u)$ with the horizontal axis by $\phi(u)$. We consider a convex curve with a turning number of one, i.e., one whose tangent has a winding number of one (qualitatively, a closed curve that loops only once). This implies

$$ \forall u \in [0, \ell] \quad \frac{d\phi}{du}(u) \ge 0 \quad \text{and} \quad \phi(u) \in [0, 2\pi]. \qquad (S48) $$

We denote by $\theta(u)$ the map that gives the subunit orientation of a point with curvilinear coordinate $u$.

Close to a grain boundary our closed loop is locally a straight line, which we represent by a dashed line in Fig. S12C. Thus any point where the loop intersects the boundary realizes either one of the following scenarios corresponding to the two panels of Fig. S12C:

* **Case 1:** The tangent to the loop is about to poke through the boundary in the same direction as the orientation vectors $\hat{s}$, which is equivalent to stating that $-\frac{\pi}{3} < \phi - \theta < \frac{2\pi}{3}$. Then $\theta(u)$ increases by an amount $\pi/3$ as the loop crosses the boundary.

* **Case 2:** The tangent to the loop is about to poke through the boundary in the opposite direction compared to the orientation vectors $\hat{s}$, which is equivalent to stating that $\frac{2\pi}{3} < \phi - \theta < \frac{5\pi}{3}$. Then $\theta(u)$ decreases by an amount $\pi/3$ as the loop crosses the boundary.

This implies that the angle $\theta$ cannot increase if $\frac{2\pi}{3} \le \phi - \theta \le \frac{5\pi}{3}$, and that it cannot decrease if $-\frac{\pi}{3} \le \phi - \theta \le \frac{2\pi}{3}$. As shown in Fig. S12E-F, these conditions divide the $(\phi, \theta)$ plane in three types of regions: those where only an increase of $\theta$ by increments of $\pi/3$ is possible; those where only an decrease of $\theta$ by decrements of $\pi/3$ is possible; and the boundaries between those two types of regions, where the loop runs parallel to the grain boundary and $\theta$ can neither increase nor decrease.

We point out two trajectories on these schematics. Since $\phi$ is an increasing function of $u$, we parametrize the trajectory along our convex loop by $\phi$ instead of $u$:
* Figure S12E shows a trajectory of maximum possible increase of $\theta$. We start with $\theta \in (-\pi/3, 0)$, which places us in a region where increases are allowed. We immediately increase $\theta$ by three increments of $\pi/3$. The trajectory is then stuck in a no-increase region. We must then wait for $\phi$ to increase to again find ourselves in the region that allows $\theta$ to increase. When that happens, we increase $\theta$ by $\pi/3$. We renew the operation whenever possible until $\phi$ reaches $2\pi$. Clearly this trajectory is constrained to lie strictly below the $\theta = \phi + \pi$ green dashed line. Thus its end point $\theta(\ell)$ lies strictly below the end point of the dashed line, whose vertical coordinate is $3\pi$, i.e., $\theta(\ell) < 3\pi$. Recalling our assumption that $\theta(0) > -\pi/3$, this implies $\theta(\ell) - \theta(0) < 10\pi/3$. In other words, the total increase of $\theta$ over the trajectory is bounded by above by $10\pi/3$.

* Figure S12F shows a trajectory of maximum possible decrease of $\theta$. We start with $\theta \in (4\pi/3, 5\pi/3)$, which immediately allows for three consecutive decreases by $\pi/3$ each. The trajectory is then stuck in a no-decrease region, and must wait until $\phi$ increases by a little less than $\pi$ to undergo three additional consecutive decreases. This scenario then repeats itself once. The whole trajectory must remain above the $\theta = \phi - 11\pi/3$ red dashed line. A similar reasoning to the one detailed above in our discussed of Figure S12e thus implies that the total decrease of $\theta$ over the course of the trajectory is bounded by below by $-10\pi/3$.

This reasoning thus implies

$$ -\frac{10\pi}{3} < \theta(\ell) - \theta(0) < \frac{10\pi}{3}. \qquad (S49) $$

Since the continuity of the orientation field additionally imposes that $\theta(\ell) = \theta(0) \mod 2\pi$, we conclude $\theta$ can only change by $-2\pi$, 0 or $2\pi$. The geometry of our defect interactions thus only allows winding numbers of the order parameter over our convex closed loop with values $-1$, 0 or $+1$.

**Isolated Point-like Defects Can Only Have Winding Number +1.** Consider a large assembly with one or more isolated topological defects, which we take to mean that we can individually enclose each of them in a finite-size convex loop as in Fig. S12D. As a result of the property proved above, the winding number associated with such an isolated defect is either +1 or -1. Here we show that the geometry of our interactions forces the -1 defects to have a finite core size (i.e., prevents them from being point-like), which compromises its local stability.

We recall that the winding number of a defect is defined by considering the change in the angle $\theta(u)$ as $u$ goes around a loop that encloses it. In the case of a counterclockwise loop, a +1 defect is characterized by a $2\pi$-increase of $\theta$, i.e., the vector $\hat{s}$ undergoes one counterclockwise rotation as we go around the loop. Conversely, a -1 defect corresponds to a decrease by $-2\pi$ or equivalently a clockwise rotation of $\hat{s}$. A -1 defect defect must thus be surrounded by crystalline domains with all six possible subunit orientations in the opposite order to that of the orientations found around a vortex assembly, which we illustrate in Fig. S12G. As shown there, the relative orientations of these domains fixes the orientation of the grain boundaries between them according to our interaction rules. This construction makes it evident that lines cannot be extended indefinitely without intersecting their counterparts. Since our interaction rules forbid an interface between two domains whose orientations differ by $2\pi/3$, these intersections must either give rise to a defect or
occur at the outer boundary of the assembly. As shown on the right of Fig. S12G, the orientation of the grain boundaries imposes that the horizontal size of the -1 defect is of the order of the vertical distance between these two defects or outer boundaries; therefore isolated -1 defects are not point-like.

To assess the cost of opening such a defect, we note in the inset of Fig. S12G that a quarter of the edges that constitutes its inner boundary has $\hat{t} \cdot \hat{s} = 0$ while the other three quarters have $\hat{t} \cdot \hat{s} = \sqrt{3}/2$. Equation (S45) thus implies that the total energy per unit length associated with the defect's inner boundary is $3/4(\epsilon_d - \epsilon_c) - \epsilon_c/2 = 3/4\epsilon_d - 5/4\epsilon_c$. In cases where this energy is negative, i.e., whenever $\epsilon_c > 3\epsilon_d/5$, such defects tend to gain energy by extending laterally. As demonstrated in Sec. 2.1.1, this condition is fulfilled in all regimes of the phase diagram of Fig. 2B of the main text where finite vortices form. As a result, an extended isolated -1 defect that forms at the center of a large assembly locally tends to grow until it spans across the whole assembly and splits it in half, thereby deleting itself.

**Convex Assemblies with a Single Vortex are Favored Over Other Defective Assemblies.** Beyond assemblies containing a single isolated topological defect, we now discuss assemblies with multiple defects. This discussion provides some assembly-wide insights on their stability which complement the local arguments discussed in the previous paragraph. For the purpose of this discussion we distinguish between convex and non-convex assemblies, i.e., between assemblies whose outer boundary (colored orange in Fig. S12) is a convex curve and those for which it is not.

We first consider a convex assembly and draw a loop along its outer boundary. Our previous result implies that the net winding number of all defects contained within this loop is -1, 0 or +1. They thus contain as many +1 defects as -1 defects plus or minus one unit. We illustrate one such assembly in Fig. S12H. Let us discuss its defect energy as well as its surface energy while approximating the +1 topological defects as point-like, implying that their contribution to the line integral of Eq. (S45) is negligible. The finite-size -1 defects, which cannot be point-like as previously discussed, contribute to this integral to the rate of 3/4 per unit length. To leading order in the assembly perimeter $P$, the number of defects in the assembly is thus

$$N_d = \left[ (1 - f) + \frac{3}{4}f \right] P + O(P^0) = \left(1 - \frac{f}{4}\right) P + O(P^0), \quad (S50)$$

where $f \in (0, 1)$ is the fraction of the total assembly perimeter associated with an inner boundary. By contrast, breaking up the same number of subunits into a set of vortex assemblies with the same perimeter would result in a larger number of defects $N_d = P + O(P^0)$ [see Eq. (S47)]. The latter option has the same surface energy as the large, multiple-defect assembly, but a more favorable defect energy, and is thus the energetically preferred outcome.

As illustrated in Fig. S12H, non-convex assemblies may display arbitrary large positive or negative winding numbers. Both types of assemblies display the same types of inner boundary which contribute to the integral of Eq. (S45) at an average rate of 3/4, plus additional boundaries which are perpendicular to the subunit orientations and thus do not contribute at all. The basic energetic argument made above for a convex assembly therefore equally applies here, and is compounded by an additional surface energy cost due to the creation of new assembly boundaries that run perpendicular to the local subunit orientation.

Overall, these arguments indicate that multiple-defect assemblies are less energetically favorable
than vortex assemblies. This rationalizes the absence of multiple-defect as well as non-convex assemblies in the results presented in the main text.
<figure><img src="image_18.png" /><figcaption>Figure S13: TEM micrograph of an AuNP labeled vortex assembly.</figcaption></figure>

This micrograph shows the same assembly as Fig. 3D with annotations removed for visibility. Here $\epsilon_c/\epsilon_d = 0.5$.
<figure><img src="image_19.png" /><figcaption>Figure S14: TEM micrographs of size-controlled vortex and fiber assemblies with DNA origami.</figcaption></figure>

(A) Typical field of view of vortex assemblies at $\epsilon_c/\epsilon_d = 0.33$. (B) Typical field of view of vortex assemblies at $\epsilon_c/\epsilon_d = 0.5$. (C) Typical field of view of vortex assemblies at $\epsilon_c/\epsilon_d = 0.6$. (D) Typical field of view of fiber assemblies at $\epsilon_c/\epsilon_d = 0.5$. (E)-(F) Close-ups on standalone fibers under the same conditions. (G)-(I) Close-ups on overlapping and branching fibers under the same conditions.
<figure><img src="image_20.png" /><figcaption>Figure S15: TEM micrographs of vortex assemblies.</figcaption></figure>

These images show defect-interaction-induced appendages reminiscent of those observed in the numerical simulation of Fig. 2C of the main text. Another instance of these appendages is visible on the micrograph of Fig. 3C of the main text.
<figure><img src="image_21.png" /><figcaption>Figure S16: Full histograms of experimental vortex assembly sizes for varying [MgCl₂], assembly times and energy ratios ε_c/ε_d.</figcaption></figure>

(A) Size distributions for isothermally self-assembled samples (ε_c/ε_d = 0.5) for 20 hour and 43 hour respectively. MgCl₂ concentrations of 5 mM (N₂₀h = 67, N₄₃h = 137), 7.5 mM (N₂₀h = 131, N₄₃h = 278) and 10 mM (N₂₀h = 364, N₄₃h = 141) were screened. Errors of mean radii represent Bessel corrected sample standard deviations. (B) Time series of size distribution of isothermally self-assembled samples (ε_c/ε_d = 0.5, [MgCl₂] = 5 mM, N₄h = 113, N₁₂h = 1223, N₂₀h = 562, N₄₃h = 556). A Kolmogorov–Smirnov p-value test reveals the similarity between histograms at different timepoints with p₁₂h→₂₀h = 0.12 and p₂₀h→₄₃h = 0.07, indicating the convergence of the underlying distribution towards a steady state. (C) Size distributions for isothermally self-assembled samples (t = 20 h, [MgCl₂] = 5 mM) for ε_c/ε_d = 0.33, 0.5, 0.55 and 0.6 respectively. (N₀.₃₃ = 106, N₀.₅ = 562, N₀.₅₅ = 122, N₀.₆ = 39).
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Figure S17: Design file of a DNA origami subunit as extracted from the scadnano software.</figcaption></figure>

(A) Spatial arrangement of DNA helices (numbered circles) that form the side wall of the DNA
origami nanocylinder. (B) Main view of the design file showing scaffold routing, staple breaking and
DNA sequences. Positions of single-stranded extensions for defect/crystalline interaction strands
and AuNP linker strands are labeled. The corresponding sequence information can be found in Data
S3 and Data S4.
**Data S1.** scadnano design file for DNAOrigami nanocylinder.

**Data S2.** DNA sequence of p2873 scaffold strand.

**Data S3.** List of mini-scaffold and base staple strand sequences for a monomer of the DNA origami nanocylinder. The prefix *interaction* is used for staple strands acting as base for single-stranded extensions used for modeling defect ($D_YX$) or crystalline ($C_Y$) interactions. Here, X describes the vertical position (top = 1 or bottom = 2) of defect linkers on a side of the DNA origami nanocylinder, whereas Y describes the side itself (see Fig. 3A of the main text).

**Data S4.** List of staple strand sequences for defect ($D_YX$) and crystalline ($C_Y$) interactions used for modelling different energy ratios $\epsilon_c/\epsilon_d$. Notations for interaction staples are the same as in Data S3 with an added $R_Z$. Here, Z represents the energy ratio $\epsilon_c/\epsilon_d$. Subsequences colored in dark and light blue represent defect and crystalline interaction strand extensions (see Fig. 3B of the main text). The rightmost column lists the Gibbs free energy change upon hybridization for each blue subsequence with its reverse complementary sequence.