Geometrically Frustrated Assembly at Finite Temperature: Phase Transitions from Self-Limiting to Bulk States

TL;DR

This study uses Monte Carlo simulations to explore phase transitions in a 2D model of geometrically frustrated assembly, revealing how concentration and temperature influence the shift from self-limiting to bulk states with defect structures.

Contribution

It provides a detailed phase diagram of frustrated assembly, identifying mechanisms of transition driven by concentration and temperature, and highlights the roles of entropy and model parameters.

Findings

Percolation of worm-like domains leads to a heterogeneous network.

Transition thresholds depend on frustration and stiffness ratios.

Self-limiting phase persists at higher frustrations and temperatures.

Abstract

Geometric frustration is recognized to generate complex morphologies in self-assembling particulate and molecular systems. In bulk states, frustrated drives structured arrays of topological defects. In the dilute limit, these systems have been shown to form a novel state of self-limiting assembly, in which the equilibrium size of multi-particle domains are finite and well-defined. In this article, we employ Monte Carlo simulations of a recently developed 2D lattice model of geometrically frustrated assembly~\cite{HackneyPhysRevX.13.041010} to study the phase transitions between the self-limiting and defect bulk phase driven by two distinct mechanisms: (i) increasing concentration and (ii) decreasing temperature or frustration. The first transition is mediated by a concentration-driven percolation transition of self-limiting, worm-like domains into an intermediate heterogeneous network mesophase, which gradually fills in at high concentration to form a quasi-uniform defect bulk state. We find that the percolation threshold is weakly dependent on frustration and shifts to higher concentration as frustration is increased, but depends strongly on the ratio of cohesion to elastic stiffness in the model. The second transition takes place between self-limiting assembly at high-temperature/frustration and phase separation into a condensed bulk state at low temperature/frustration. We consider the competing influences that translational and conformational entropy have on the critical temperature/frustration and show that the self-limiting phase is stabilized at higher frustrations and temperatures than previously expected. Taken together, this understanding of the transition pathways from self-limiting to bulk defect phases of frustrated assembly allows us to map the phase behavior of this 2D minimal model over the full range of concentration.

Equations

$$H=-J\sum _{⟨ij⟩}\cos (\Delta {\theta }_{ij}-A_{ij}){\eta }_i{\eta }_j-K\sum _{⟨ij⟩}{\eta }_i{\eta }_j$$

$$\sum _{\circlearrowright }{A}_{ij}\equiv 2\pi f.$$

$${\nabla }_{\perp }\times \mathbf{A}(\mathbf{x})=2\pi \phi ,$$

$$\Sigma =J+K,$$

$${\ell }_{\mathrm{v}}\approx {\phi }^{-1/2}.$$

$$\mathcal{H}[\mathcal{D}]=\frac{J}{2}{\int }_{\mathcal{D}}{\mathrm{d}}^2\mathbf{x}|\nabla \theta -\mathbf{A}(\mathbf{x}){|}^2+\sigma \mathcal{P}[\mathcal{D}]-{ϵ}_{\mathrm{bulk}}\mathcal{A}[\mathcal{D}]$$

$$\sigma =\frac{\Sigma }{2a}$$

$${\nabla }_{\perp }^2\theta (x,y)={\nabla }_{\perp }\cdot \mathbf{A}$$

$$\hat{n}\cdot {\nabla }_{\perp }\theta {|}_{\partial D}=\hat{n}\cdot \mathbf{A}.$$

$$\mathbf{A}=\pi \phi [y\hat{x}-x\hat{y}]$$

$$ϵ[\mathcal{D}]=\frac{J}{2\mathcal{A}[\mathcal{D}]}{\int }_{\mathcal{D}}{\mathrm{d}}^2\mathbf{x}|\nabla {\theta }_{\ast }-\mathbf{A}(\mathbf{x}){|}^2+\sigma \frac{\mathcal{P}[\mathcal{D}]}{\mathcal{A}[\mathcal{D}]}-{ϵ}_{\mathrm{bulk}}.$$

$${\ell }_{\mathrm{c}}=(\sigma /J{)}^{-1}$$

$${\theta }_{\ast }^{(\mathrm{ribbon})}(\mathbf{x})=\pi \phi xy.$$

$$\Delta ϵ(W)\equiv ϵ(W)+{ϵ}_{\mathrm{bulk}}=\frac{J{\pi }^2{\phi }^2}{6}{W}^2+\frac{2\sigma }{W}.$$

$$W_{*}=\bigg{(}\frac{6}{\pi^{2}}\bigg{)}^{1/3}\ell_{\rm d}$$

$$\ell_{\rm d}=\bigg{(}\frac{\sigma}{J\varphi^{2}}\bigg{)}^{1/3}=\ell^{4/3}_{\rm v}\ell^{-1/3}_{\rm c}$$

$${ϵ}_{\mathrm{sla}}=C_0{J}^{1/3}{\sigma }^{2/3}{\phi }^{2/3}$$

$$s(\mathbf{x})=\sum _{n,m}{s}_{nm}\delta (\mathbf{x}-{\mathbf{x}}_{nm})$$

$$\mathbf{x}_{nm}=d\Big{(}n+\frac{m}{2},\frac{\sqrt{3}}{2}m\Big{)}$$

$$\nabla \times (\nabla \theta )=s(\mathbf{x})$$

$$\frac{J}{2}\int d^{2}{\bf x}\big{|}\nabla\theta-{\bf A}\big{|}^{2}=\frac{JA}{2}\sum_{{\bf G}\neq 0}\frac{\rho_{s}^{2}}{|{\bf G}|^{2}},$$

$$\sum _{\mathbf{G}\ne 0}\frac{1}{|\mathbf{G}{|}^2}\approx {\int }_{1/d}^{{\Lambda }_{\max }}\frac{dq}{q}=\ln (d/a).$$

$${ϵ}_{\mathrm{elas}}(\mathbf{x})\equiv \frac{J}{2}|\nabla \theta (\mathbf{x})-\mathbf{A}(\mathbf{x}){|}^2$$

$$(\mathbf{\hat{n}}\cdot\nabla)\epsilon_{\rm elas}(\mathbf{x})\big{|}_{\partial{\cal C}}=0$$

$$\nabla \theta =\frac{s}{r}\hat{\varphi }.$$

$$\epsilon_{\rm elas}(r)=\frac{J}{2}\Big{(}\frac{s}{r}-\pi\varphi r\Big{)}^{2}.$$

$$\pi R_0^2\phi =s.$$

$$\begin{gathered} \frac{{ϵ}_{\mathrm{defect}}(\bar{r})}{J/R_0^2}=-s^2\frac{\ln \bar{r}}{1-{\bar{r}}^2}-\frac{s}{2}(\phi R_0^2)+\frac{1}{16}(\phi R_0^2{)}^2(1+{\bar{r}}^2) \\ +\frac{2}{1-{\bar{r}}^2}\frac{\bar{r}}{\bar{\lambda }} \end{gathered}$$

$$\bar{\lambda }{s}^2=\frac{1+{\bar{r}}^2}{\frac{1}{\bar{r}}(1-{\bar{r}}^2{)}^2+2\bar{r}\ln \bar{r}-\frac{\bar{r}}{2}(1-\bar{r}{)}^2}$$

$${\bar{r}}_{\ast }\simeq \{\begin{array}{ll}\bar{\lambda }& \bar{\lambda }\to 0\\ 1-{\bar{\lambda }}^{-1/3}& \bar{\lambda }\to \infty \end{array}$$

Full text

††preprint: APS/123-QED

I Introduction

Geometric frustration refers to the scenario in which global geometrical constraints on the degrees of freedom of an interacting system obstruct the uniform states of locally preferred order across the entire system [2]. This idea has been used to rationalize the behavior of a wide array of systems ranging from low-temperature magnetism [3, 4] to superconductors in an external field [5]. Further, geometric frustration has been recognized to play an important role in shaping the behavior of many different systems of self-assembling particles in soft matter [6, 7, 8, 9]. Here, deformation of inter-particle arrangement away from the locally preferred state can result in intra-assembly strain gradients that propagate to long range. Examples of such geometrically frustrated assembly include crystallization on curved surfaces [10, 11, 12, 13]; non-tiling polygonal [7, 14, 15] and polyhedral particles [16, 17]; twisted bundles of self-assembling protein filaments [18, 19]; hyperbolic membranes [20, 21, 22, 23, 24]; liquid crystals on curved surfaces [25, 26, 27, 28]; two-dimensional assemblies of bent-core nematic liquid crystals [29, 30]; one-dimensional assemblies of curved colloidal [31, 32] and incommensurate ‘polybrick’ particles [33, 34]; as well as models of viral capsid assembly [35, 36].

Broadly speaking, investigations into the behavior of these systems have focused on the case of either macroscopic bulk or dilute self-assembly. In the first case, the assembling subunits extend to completely fill the available space. Here, the lack of both translational degrees of freedom and free boundaries limit the systems ability to alleviate the accumulation of intra-assembly strain gradients, leading to the proliferation of topological defects throughout the ground state [2, 37, 38]. Conversely, in systems of dilute self-assembly, the added subunit mobility allows for relaxation of the effects of frustration via variation of the size and shape of the individual multi-unit domains. These extra degrees of freedom have been argued to give rise to unique scale-dependent behavior known as self-limiting assembly [9] resulting from the thermodynamic balance between cohesive loss at the boundary and elastic cost of accumulating stress due to frustration. Conceptually, the putative self-limiting state is of particular interest as it implies that frustrated assemblies are able to “measure” and limit their equilibrium dimensions at size scales much larger than the individual subunits or their interaction range. This property is also of potential value for creating synthetic self-assembling materials that can target well-defined mesostructured dimensions that are tunable through programmable misfit of the constituent subunits themselves.

While early theories of self-limitation (e.g. [11, 39]) have been derived from ground-state thermodynamics of continuum elastic models of frustrated assembly, the thermodynamic stability of the self-limiting state at finite temperature has only more recently been studied [1, 34]. Here it was shown that, for fixed and sufficiently dilute concentration, the self-limiting state occurs at intermediate frustration strength. In addition, this state was shown to be flanked by the dispersed state at high frustration and a defect-riddled bulk state at low, but non-zero frustration. In this article, we focus on the possible thermodynamic pathways between the self-limiting and bulk states more broadly, and in the context of a recently developed 2D lattice model of geometrically frustrated assembly [1]. The model is an extension of the uniformly frustrated 2D XY model [38, 40] that incorporates translational degrees of freedom at states of fixed total subunit concentration, which is known to encode the minimal elements of frustration of orientational 2D order on surfaces with fixed and non-zero Gaussian curvature [41, 42]. In particular, it was shown that self-limiting states of this generalized lattice model take the form of finite-width, worm-like domains, while the bulk condensed state takes the form of a sponge-like domain where the “holes” correspond to topologically charged vortices that screen frustration.

While Monte Carlo simulations of the frustrated lattice model have demonstrated the existence of the self-limiting state at finite-temperature and characterized a transition to a bulk-condensed state below a critical frustration strength, two key aspects of the self-limiting to bulk transition remain unclear. First, the existence of a self-limiting to bulk transition with decreasing frustration can be argued purely on the ground state energetics of these states and their distinct frustration dependence. This energetic scaling argument suggests a critical frustration value for the transition proportional to the square of the cohesion to (spin) stiffness ratio, notably independent of temperature. While the critical frustration observed from simulations was shown to increase with cohesion to stiffness, it was also found to be strongly temperature dependent, which implies that entropic factors in the free energy difference between bulk and self-limiting states are non-negligible. At present, it is not clear what is the most significant source of entropy difference between these two states (e.g. translational or conformational fluctuations), how it compares to the relative cohesive and elastic energetic difference, and if and under what conditions this transition ever approaches the purely energetic description. A second, and far more open question, is how the phase behavior of the frustrated lattice model evolves at fixed, intermediate frustration with increasing concentration. While at dilute conditions, the self-limiting state is formed above a (pseudo-)critical aggregation concentration (CAC), it is intuitive to expect a bulk equilibrium state in the limit of large concentration, where the translational degrees of freedom are severely limited and excluded volume effects become strong. It is natural to expect that this high-concentration bulk state will be a space-filling version of the defect-riddled, spongey condensate that forms below the critical value of frustration in dilute conditions (i.e. a semi-regular array of vortices with voided “cores”). However, it is so far unclear how the system proceeds thermodynamically from a state of dispersed, finite-width and variable length aggregates to a solid-like and quasi-uniform array of defect voids. The existence of one or more possible mesophases (e.g. liquid crystalline order) intermediate to self-limiting and bulk states, and the corresponding nature of the thermodynamic transitions between such states, is not known.

In this article, we address these open aspects of the self-limiting to bulk transition in a minimal 2D lattice model of geometrically frustrated assembly [1] using a combination of Monte Carlo simulation and mean-field, continuum theory modeling of its assembly thermodynamics. The scope and key findings of this study are illustrated schematically in the phase diagram shown in Figure 1 along with corresponding cartoons of the different competing phases. First, for the concentration driven transition at fixed intermediate frustration (horizontal arrow in Fig. 1), we show the existence of a phase transition between self-limiting assembly and an intermediate mesophase characterized by percolation of finite-width domains into a heterogeneous “gel-like” network. We show that this percolation transition is weakly dependent on frustration, but strongly dependent on the ratio of cohesion to stiffness, behavior which we compare to an effective bond percolation of self-limiting domains. Upon further increase in concentration, the aggregated and disordered network continues to fill in and continuously evolves to a “quasi-uniform” bulk sponge state that consists of a regular array of similarly-sized, evenly spaced holes surrounding topologically charged vortices. In the limit of full surface coverage, these holes are then uniformly filled in to reach the fully occupied bulk, Abrikosov defect ground state of the uniformly frustrated XY model. We then determine the optimal density of the defect sponge phase, which characterizes the transition between disordered and ordered bulk sponge states, and show that it decreases with frustration and the ratio of inter-particle cohesion to interaction stiffness. Next, we revisit the frustration-driven transition between self-limiting and bulk states at fixed dilute concentrations (vertical arrow in Fig. 1). As indicated in the schematic phase diagram, this transition takes place between a poly-disperse collection of finite sized, worm-like aggregates and a phase separated defect bulk. We analyze distinct microstructural contributions to the relative entropy of self-limiting aggregates and bulk condensates, and propose a mean-field argument to capture the expected leading order contributions from the relatively higher entropy of self-limiting states, and their impact on the relative free energy differences between these states. We show that this model captures the qualitative effect observed in simulation, with the critical frustration decreasing as temperature is increased. We further predict the range of thermodynamic conditions at which the relative thermodynamics of self-limiting and bulk defect states is determined by purely energetic considerations of intra-domain elasticity and cohesion and discuss the implications of this regime for accessible conditions for the state of self-limiting assembly.

The rest of this article is organized as follows. In Section II, we provide a brief review of the 2D lattice model of frustrated assembly, including the microscopic degrees of freedom and key parameters. In Section III, we discuss results from the continuum limit of the two key competing morphologies, and describe the connection with microscopic parameters of the model with key length scales that characterize the mesoscopic effects of frustration. In Section IV.1, we investigate the existence of a frustration dependent percolation transition from self-limiting assembly to a heterogeneous network sponge and examine how this intermediate state approaches the homogeneous sponge at near full concentration. We also compare the internal structure of this sequence of states to that of the phase-separated defect bulk found at weak frustration. In Section V, we revisit the phase separation that occurs at dilute concentration (i.e. below the percolation threshold) and weak frustration; analyzing several different contributions to the total aggregation entropy, including the translational and configurational entropy of worm-like domains. We then describe our mean-field model for the leading temperature-dependence of critical frustration. Finally, in Section VI, we summarize our key results and discuss their implications for the global phase behavior of the frustrated lattice model for the case of fixed, uniform frustration studied here as well as for connections to experimental systems.

II Lattice Model of Geometrically Frustrated Assembly

The model of geometrically frustrated assembly that we will be using throughout this article describes a collection of $N$ subunits occupying a fixed fraction, $\Phi=N/L^{2}$, of sites on an $L\times L$ square lattice with bond spacing $a$. In addition to translational degrees of freedom, each subunit possesses a continuous orientational degree of freedom, $\theta_{i}\in[0,2\pi]$, i.e. an XY-spin or “phase” variable. The energetics of this model can be described by the lattice Hamiltonian [1]:

[TABLE]

where $\langle ij\rangle$ denotes nearest neighbor lattice sites, $\Delta\theta_{ij}=\theta_{i}-\theta_{j}$ defines the nearest neighbor phase difference and $\eta_{i}=0,1$ denotes the occupancy of each lattice site. The first term describes an $XY$-like interaction based off the relative phase difference between subunits and is parameterized by the spin stiffness $J$. The second term describes an Ising-like nearest neighbor interaction parameterized by the cohesive energy $K$. Lastly, the $A_{ij}$ term represents a gauge field along the lattice links that defines the locally preferred nearest neighbor rotation, and crucially introduces the geometric frustration of an aggregated cluster of spins. Following the standard approach for the uniformly frustrated XY model [43, 38], we define the gauge field subject to the constraint

[TABLE]

Here, $f$ is a measure of the microscopic frustration strength and the sum is taken over clockwise oriented minimal lattice plaquettes. This constraint implies that, when $f\neq n$ where $n$ is an integer, propagation of the preferred nearest neighbor rotation, $\Delta\theta_{ij}=A_{ij}$, around the plaquette leads to a net angular mismatch of $2\pi f$, and hence the locally preferred state is incommensurate with the existence of closed loops and the system is said to be geometrically frustrated. The parameter $f\in[0,1/2]$ defines the degree to which the occupied square lattice is frustrated, with $f=0$ corresponding to the unfrustrated case and $f=1/2$ corresponding to the case of maximal frustration [43]. In this article, we focus on the case of fixed and uniform frustration, where the discrete gauge field can be defined by the integral $A_{ij}=\int_{\mathbf{x}_{j}}^{\mathbf{x}_{i}}\rm d\mathbf{x}\cdot\mathbf{A}(\mathbf{x})$ where $\mathbf{A}(\mathbf{x})$ is a two-dimensional vector field defined to have constant curl

[TABLE]

and $\varphi\equiv f/a^{2}$ defines a frustration density. Notably, as one considers occupied domains of increasingly larger size (i.e. upon cohesive assembly into compact clusters), the number of incompatible loops grows as well as the magnitude of phase mismatch around the larger loops (i.e. per Stokes law), implying that the energetic cost of the ground state spin configuration will also grow due to an increasing number and magnitude of local spin misalignment. This accumulation of energetic cost with increasing domain size has been identified as a hallmark behavior of geometrically frustrated assembly [8] and has been shown to give rise to both the self-limiting and defect-riddled bulk states of assembly. In terms of specific physical systems [1], this minimal model most closely describes the frustration of in-plane orientational order (i.e. liquid crystalline) on 2D surfaces of fixed Gaussian curvature [41, 25], proportional to $\varphi$.

The behavior of this model is determined by two microscopic energy scales. The first is defined by the net cohesion of an ideally-aligned bond between subunits

[TABLE]

while the second is simply $J$ itself, which plays an additional role of spin stiffness, penalizing orientational strain gradients. The behavior of this lattice model has been studied by MC simulation in the dilute limit ($\Phi\lesssim 0.2$) and for the case of weak cohesion to stiffness ratios, $\Sigma/J\ll 1$ [1]. In this regime, there exists an apparent binodal in the frustration-concentration plane below which the system phase separates into a topologically defective condensate surrounded by a dispersed monomer gas. Above the binodal (characterized by a critical frustration value), there is a pseudo-critical aggregation transition from a dispersed gas to a state of self-limiting assembly characterized by distribution of uniform finite width, variable length aggregates in coexistence with dispersed monomers. These aggregates exhibit highly anisotropic ribbon-like morphologies with uniform width controlled by the ratio of $f$ and $\Sigma/J$ and is independent of any extensive properties of the system. For a given value of $\Sigma/J$, the critical frustration along the binodal is apparently constant and independent of concentration. As mentioned in the introduction, this motivates the terminology strong versus weak frustration as a means of classifying whether a given set of conditions respectively resides above or below the phase separation binodal.

Conversely, in the fully occupied limit ($\Phi=1$), eq. (1) described uniformly frustrated 2D XY model, which has been used as a statistical framework for studying field-frustrated superconducting arrays [43, 44, 5] and various frustration based glass models [40, 38]. In the fully occupied regime, sufficiently below the Kosterlitz-Thouless temperature of the unfrustrated XY model (i.e. $\beta J\gg 1$), frustration leads to ground states composed of a uniform density of like-signed vortices arranged in a triangular Abrikosov lattice with inter-defect spacing [45]

[TABLE]

In this this article, we apply the Markov chain Monte Carlo algorithm detailed in Ref. [1] to study the phase behavior of frustrated assembly across the full range of concentration. Simulation snapshots of each of the unique states of assembly – self-limiting aggregates, heterogeneous network, quasi-uniform defect sponge and phase separated bulk (sponge) – enumerated in the schematic phase diagram are provided in Fig. 2. Supplementary videos 1-4 show animations of simulated sequences of each of these states.

III Continuum Models of Competing Ground State Morphologies

Here we describe continuum theory predictions for two characteristic morphologies we observed, shown schematically in Fig. 3ab. This analysis is based on the assumption that, for conditions relevant to self-limiting assembly, thermal fluctuations of spin degrees of freedom are relatively weak. As we see below, this is justified by the fact that self-limiting domain sizes larger than a single subunit require $\Sigma/J\ll 1$ while equilibrium assembly necessarily requires sufficiently cohesive interactions $\Sigma\gtrsim k_{B}T$. These conditions imply $J/k_{B}T\gg 1$ and that spin degrees of freedom of bound clusters can be modeled in terms of their ground state configurations.

When the phase strain, $\Delta\theta_{ij}-A_{ij}$, is small, we can take the continuum limit of equation (1) and write a function for the energy of a domain, $\mathcal{D}$, of assembled subunits as:

[TABLE]

where the first term describes the elastic energy of deforming the phase field away from its locally preferred state and the last two terms describe the cohesive energy of a domain with perimeter $\mathcal{P}[\mathcal{D}]$ and area $\mathcal{A}[\mathcal{D}]$. Here, we define the

[TABLE]

as the line energy of the domain boundary (i.e. cohesive cost per unit length) and $\epsilon_{\rm bulk}=2\Sigma/a^{2}$ as the cohesive bulk energy density. Assuming a given fixed domain, $\mathcal{D}$, the Euler-Lagrange equation of for spin degrees of freedom is simply

[TABLE]

with the free boundary conditions

[TABLE]

For simplicity we consider a divergence-free gauge

[TABLE]

for which the right hand side of eq. (8) vanishes, so that the ground states, $\theta_{*}$, are simply harmonic functions satisfying the appropriate phase winding condition at the free boundaries.

For the purposes of understanding the states of optimal assembly, it is necessary to determine the size, shape and topology of domains that optimize the assembly energy density [9],

[TABLE]

Optimal states of assembly of the continuum model derive from optimization of this energy density with respect to domains, and can therefore only depend on two length scales: the length scale set by frustration density, i.e. $\ell_{\rm v}=\varphi^{-1/2}$ and the second determined by the ratio,

[TABLE]

which constitutes a coheso-elastic length scale for the assembly. We now summarize how these length scales determine the structure and energetics of finite-width ribbons and sponge-like, bulk vortex arrays.

III.1 Finite-width ribbon domains

In reference [1], we computed the assembly free energy (per subunit) landscape from the continuum theory of this lattice frustration model for rectangular, defect free domains of arbitrary cross-sectional dimensions. This theory showed that optimal ground-state morphologies break symmetry in highly-anisotropic ribbon-like domains, of finite, self-limiting width and arbitrarily unlimited length. This effect of frustration, which is sometimes referred to as “filamentation” [8], favors elastic misfit gradients distributed across the narrow dimension of the domain and uniform elastic energy along the length of the aggregate and is widely observed in 2D frustrated assemblies [11, 13, 21, 6, 15, 17, 46]. Hence, we only review the energetic selection of the finite width, ribbon-like domains.

We consider a domain extending along the $y$-axis, and finite width, $W$, in the $x$-direction from $x\in[-W/2,+W/2]$, as shown schematically in Fig. 3a. In this geometry, the boundary conditions require that $\partial_{x}\theta=\pi\varphi y$ at $x=\pm W/2$. Additional, uniformity in energy density along the ribbons requires $\partial_{y}(\nabla\theta-\mathbf{A})=0$. These conditions are satisfied by the harmonic function for ground state phase configuration

[TABLE]

In finite width ribbons, ground-state phase-strain, $\nabla\theta-\mathbf{A}=2\pi\varphi x\hat{y}$, is in the longitudinal direction with transverse gradients in magnitude, and the elastic energy density of ribbons grows quadratically with width, proportional to $J\varphi^{2}W^{2}$. The full calculation gives an energy density (relative to cohesive energy $-\epsilon_{\rm bulk}$):

[TABLE]

Minimizing with respect to width gives

[TABLE]

where we have defined the domain scale

[TABLE]

which is itself a combination of length scales set by the frustration density ($\ell_{\rm v}$) and coheso-elastic interplay ($\ell_{\rm c}$). Inserting the selected width into $\Delta\epsilon(W)$, we obtain the ground state cost of finite-width ribbons per unit area

[TABLE]

where $C_{0}=(9\pi^{2}/2)^{1/3}$. Note that, when normalizing this cost relative to the characteristic energy density $J\varphi$, we find $\epsilon_{\rm sla}/(J\varphi)\propto(\ell_{\rm v}/\ell_{\rm c})^{2/3}=(\ell_{\rm d}/\ell_{\rm v})^{2}$, showing that the cost of finite domain formation is controlled by the ratio of the mesoscopic self-limiting width to inter-vortex spacing. Since, as we discuss next, the energy density of the optimal vortex array is roughly $J\varphi$, this suggests that finite-domain formation is energetically favorable compared to defective bulk states when self-limiting domain size is small compared to the characteristic defect spacing. This result is consistent with the intuitive argument put forward in ref. [1].

III.2 Bulk vortex lattice

We next consider the energetics of the defect bulk phase. As previously described by Hackney et. al. [1] and shown in Fig. 2, the bulk condensed state of the lattice frustration model forms a quasi-regular array of vortices with empty (i.e. voided) defect cores, leading to a quasi-uniform sponge morphology. To understand the ground state energetics of this state, we first briefly consider the fully-occupied ($\Phi\to 1$) case, in which defect cores are completely filled in and can be modeled as an infinite plane tiled by a hexagonal Abrikosov lattice of like-charged vortices. Vortices are singular, non-analytic field configurations around which the phase field winds by integer multiples of $2\pi$. We define the areal density of defects as

[TABLE]

where $s_{nm}=\pm 2\pi$ is the topological winding of each elementary vortex and $\mathbf{x}_{nm}$ is the location of the vortex on site $nm$ in the array,

[TABLE]

where $d$ is the nearest neighbor spacing of vortices on a triangular lattice and $n,m\in[-\infty,\infty]$. Following a standard approach [47], we can solve for the phase field produced by the superposition over vortices in the infinite array and compute the elastic energy.

To do so, we note that the intra-aggregate phase field can be related to the areal defect density field via Stoke’s theorem:

[TABLE]

where $\nabla\times{\bf v}=\epsilon_{ij}\partial_{i}v_{j}$ is the 2D curl of vector ${\bf v}$. Following a standard approach [47], we introduce the conjugate field $\tilde{\theta}$ defined by $\partial_{i}\theta=\epsilon_{ij}\partial_{j}\tilde{\theta}$ for which eq. (20) becomes an effective Poisson equation, $\nabla^{2}\tilde{\theta}=s(\mathbf{x})$. From these equations, we can solve for the elastic cost of the defect array and background frustration as

[TABLE]

where {\bf G}=4\pi/(\sqrt{3}d)\big{(}\sqrt{3}h/2,k-h/2\big{)} defines the reciprocal lattice vectors for integer values of $h$ and $k$. Here $\rho_{s}=2s/3d^{2}$ is the areal topological charge density of the array, which must be set to $2\pi\varphi$ to neutralize the far-field elastic cost of the background frustration by removing the divergent ${\bf G}\to 0$ term in the sum above. Assuming that $s=2\pi$, this gives $d=\sqrt{2/(3\varphi)}\propto\ell_{v}$. Notably, the sum over reciprocal vectors in eq. (21) diverges at large $|{\bf G}|$ and must be cutoff at some maximal wave vector $\Lambda_{\rm max}\approx 1/a$, resulting an a logarithmically divergent elastic cost as the microscopic lattice dimension vanishes,

[TABLE]

This illustrates the sensitivity of the bulk state to the core size, $r_{c}$, of the vortices and notably agrees with the observation that, in the bulk phase of the lattice GFA model, the vortex cores adopt a finite radius. Therefore, it is intuitive to understand that, for $\sigma/J\ll 1$, the core size will be much larger than the microscopic cutoff, $a$, due to the competition between the elastic energy, which favors larger cores, and the cohesive energy, which favors smaller cores.

To model the effect of the finite core on the bulk state, we consider a hole of radius $r_{\rm c}$ around each vortex in the assembly. To assess the optimal structure, we consider an approximate unit cell calculation of the phase-strain in which we model the hexagonal cell that surrounds each vortex in the triangular lattice by a circular one (see schematic in Fig. 3b). In an infinite periodic vortex array, the elastic energy density,

[TABLE]

is invariant under symmetries of the arrangement. In particular, normal gradients of the energy density vanish along the cell boundaries, which are lines of mirror symmetry in the triangular lattice, i.e.

[TABLE]

where $\partial C$ denotes the outer boundary of unit cell ${\cal C}$ and $\mathbf{\hat{n}}$ defines its normal. We approximate the unit cell as a circle of outer radius $R_{0}$ surrounding a single vortex at its center, which sits within a voided circular core of radius $r_{c}$ (see Figure 3). To determine the cell dimension, we superpose the phase strain generated by the central vortex

[TABLE]

with the azimuthal gauge field $\mathbf{A}=\pi\varphi~r\hat{\phi}$, to find the axisymmetric elastic energy density

[TABLE]

Applying the cell boundary condition in eq. (24), $\partial_{r}\epsilon_{\rm elas}(R_{0})=0$, results in the following condition of cell size,

[TABLE]

This corresponds to the condition that the defect charges neutralize the background charge of the frustration field and additionally sets the inter-defect spacing $\ell_{\rm v}\approx 2R_{0}\propto\varphi^{-1/2}$.

Given this optimal cell size (i.e. vortex density) we average the free energy density over the occupied annular region $r_{c}<r<R_{0}$ of a single cell, which be written as a function of reduced core size $\bar{r}=r_{c}/R_{0}$,

[TABLE]

which we have written in dimensionless form by introducing the scaled variables $\bar{r}=r_{\rm c}/R_{0}$ and $\bar{\lambda}=(J/\sigma)R_{0}^{-1}\propto\ell_{\rm c}/\ell_{\rm v}\propto\varphi^{1/2}$.

Observing that $\varphi R_{0}^{2}=2s$, we can minimize $\epsilon(\bar{r})$ with respect to $\bar{r}$ to obtain a transcendental equation for the optimal dimensionless hole size $\bar{r}_{*}$:

[TABLE]

which can be approximately satisfied in the limit of small and large $\bar{\lambda}$ as:

[TABLE]

Plugging this back into the expression for $\epsilon(\bar{r})$, we find the limiting form of the optimal sponge energy density:

[TABLE]

Similarly, we can approximate the optimal hole concentration, $\Phi_{\rm hole}=\bar{r}_{*}^{2}$, as:

[TABLE]

In the limit of low frustration, where $\bar{\lambda}\propto\varphi^{1/2}\to 0$, the optimal hole size goes to a constant value of $r_{*}\approx\ell_{\rm c}$ set by the elasto-cohesion scale such that the optimal vortex density increases as the inter-defect spacing diverges in the zero frustration limit. In the large frustration (or low cohesion when $\ell_{\rm c}\gg\ell_{\rm v}$) limit, optimal hole size approaches the inter-vortex spacing, but with a finite inter-vortex “gap” $2(R_{0}-r_{*})\approx R_{0}\bar{\lambda}^{-1/3}\propto\ell_{\rm v}^{2/3}\ell_{\rm c}^{1/3}\sim\varphi^{-1/3}$, which is distinct from (i.e. larger than) the self-limiting domain size scaling.

For low-enough temperature, we assume that the relationship between hole size and concentration derived above will hold for both the uniform vortex sponge at high concentration (when the total density exceeds $1-\Phi_{\rm hole}$) as well as the phase-separated condensates found at weak frustration. Thus, we use this predicted dependence of vortex core size and inter-vortex spacing as a basis for comparing the structural similarity of these two types of bulk states.

IV Concentration Driven Transition from Self-Limiting to Bulk at Strong Frustration

In section I, we suggested that the dilute state of dispersed/self-limiting assembly, found at strong frustration, ultimately evolves into a macroscopically large bulk phase upon increasing concentration via a percolation of finite width domains, first into a heterogneous network of finite width domains and then ultimately into a relatively ordered bulk array of defect-enclosing holes at higher concentration. We investigate this evolution with MC simulations of the lattice model described in section II by running simulations at fixed $f$, $\Sigma/J$, $\beta J$, and increasing concentration, $\Phi$. In particular, we analyze the statistics of clusters of assembled particles, identified as connected graphs of occupied nearest neighbor sites on the square lattice.

IV.1 Percolated networks of finite width domains

The percolation fraction of a given set of conditions is obtained by observing aggregated clusters at regular intervals and counting the fraction of time-steps containing a system spanning aggregate [48]. The percolation fraction as a function of concentration is plotted in Figure 4 for several different values of frustration. Here, we find that the percolation fraction rapidly grows from zero at low concentration to near unity as $\Phi\rightarrow 1$. By defining the percolation threshold, $\Phi_{\rm perc}$, as the concentration at which a spanning cluster is observed half of the time, we quantitatively evaluate this transition. Doing so, we observe $\Phi_{\rm perc}\in[0.4,0.65]$, with the percolation threshold shifting to higher concentration with increased frustration where the self-limiting domains tend to become more narrow in width. This frustration dependence illustrates that the nature of the percolation transition in networks of self-limiting domains is different from the simplest classes of bond or site percolation on a square lattice, which have a fixed percolation threshold of $0.5$ and $\sim 0.59$, respectively [49]. Lastly, we note that the tendency of increased frustration to shift the percolation threshold to higher concentration is consistent with findings of recent study of a similar lattice frustration model that considers phase behavior in the grand canonical ensemble [50].

One likely factor influencing the shift in percolation threshold could be the tendency of frustration to alter the density of branch points in the aggregates, as branchy aggregates are more spatially dense and require more subunits (i.e. higher concentration) to span the lattice [51]. In Figure 5a, we investigate this idea by counting the number of branches per aggregate as a function of morphological skeleton length (i.e. the length of 1D “medial backbone” at the center of the domains) for a range of different concentrations and frustrations. We observe that, in non-percolating aggregates, the number of branches grows roughly linearly with skeleton length, suggesting a uniform branch density below the percolation threshold, or equivalently the mean linear distance between branches is roughly constant and independent of $\Phi$. Furthermore, in Figure 5b, we show that this branch density does indeed increase with frustration, albeit weakly. Going further, we observe that, near and above the percolation threshold, the number of branches grows super-linearly with skeleton length. In this regime, aggregates remain network-like, but realize a branch density that increases with concentration. This suggests a transition in branching morphology of quasi-1D width-limited aggregates that occurs near the onset of percolation.

IV.2 Tree-like versus loopy aggregate topology

Careful analysis of aggregates reveals that this transition to super-linear growth in branch number can be associated with a rapid accumulation of loops of a characteristic dimension set by the preferred vortex size. We characterize this transition in aggregate topology via the fraction of the morphological skeleton contained within tree-like versus cyclic sections. In Figure 6a, we show a sequence of states spanning the transition with the morphological skeleton colored to show the existence of cycles 111When counting looping spans of the graph, we employ an maximal cycle length corresponding to one and a half times the average cycle length measured in the optimal sponge phase corresponding to the energy density minima found in Figure 8, since the optimal vortex size sets a preferred size of loops. Here we see that, at low concentration, the system is dominated by tree-like aggregates with uniform branch density. At high concentration, the system is dominated by defect containing loops in the morphological skeleton. In between these two limits, we observe coexistence of dense clusters of locally loopy (defective) and tree-like (finite-width) aggregates. This crossover from trees to cycles is neatly summarized in Figure 6b, where we show the fraction of morphological skeleton contained within cycles as a function of concentration. Interestingly we find that, at low concentration (but above the pseudo-critical aggregation concentration threshold) the cycle fraction increases as frustration is lowered; this suggests that the aggregates become less tree-like as the weak frustration regime (i.e the phase separation binodal) is approached. Conversely, at high concentration, we find that the cycle fraction increases with frustration, consistent with the increased density (i.e. smaller size) of defect holes. Furthermore, by defining the cycle threshold as the concentration where half of the morphological skeleton is contained within loops that surround defect holes, we can compare the location of this transition in aggregate topology to the formation of a percolating network (see Figure 6c). From this, we see that these two tendencies coincide and that the onset of percolation—which is a global phenomena and is only well defined in the thermodynamic limit—is strongly correlated with the structural evolution of aggregated clusters from predominately tree-like clusters into loopy networks.

IV.3 Evolution from heterogeneous networks to quasi-uniform bulk “sponge”

We next analyze the correlations between defect enclosing loops in the assembled aggregates and their dependence on concentration. Near the onset of percolation, coexistence of the two populations of aggregate morphology leads to a state characterized by dense clusters of defect holes connected by tree-like ribbons—resulting in a bulk phase with a large degree of non-uniformity. This is distinct from the high concentration bulk phase, where the assemblies are more uniformly defect-riddled in terms of hole size and spacing. We describe these distinct structures observed at concentrations near and far above the percolation/cycle threshold as heterogeneous networks and (quasi-)uniform defect sponges, respectively. Additionally, we note the similarity between the structure of the uniform defect sponge and the bulk aggregates observed at weak frustration (i.e. uniformly spaced holes of a given characteristic size) suggesting that the condensed phase represents phase separation between the low concentration monomer gas and the high concentration uniform defect sponge phase.

The distinction between these states can be seen by looking at the average spacing between holes, normalized by the expected center-to-center spacing between defects (i.e. $\ell_{\rm v}\sim\varphi^{-1/2}$) in the fully occupied Abrikosov vortex lattice, shown in Fig. 7. At weak frustration, the normalized hole spacing is a constant of order one across the entire range of concentration (Fig. 7c), consistent with the expected Abrikosov-like structure of the bulk condensates. Above the binodal (Fig. 7a-b), the hole spacing rises from zero to a peak at concentrations near to the percolation threshold. When $\Phi\gg\Phi_{\rm perc}$, the normalized spacing of the strongly frustrated aggregate continuously drops to $\sim 1$. Thus indicating that the structure of the percolated aggregates gradually converges to the Abrikosov-like structure of the phase separated condensates, consistent with the emergence of a space filling network of uniformly sized loops enclosing topological defects. Lastly, we observe a sharp increase in the inter-hole spacing in the limit that $\Phi\sim 1$. This can be easily understood as the limit where the holes start to completely fill in, leading to a divergent spacing between a tiny number of empty cores.

The agreement between the center-to-center hole spacing and the predicted inter-defect spacing suggests that the intra-assembly holes might also exhibit inter-vortex order (e.g. trianglar lattice) at least at sufficiently low temperature. However, despite the regularity in spacing, the observed array of holes don’t display a large amount of orientational order (see Appendix A). This suggests that the temperature range in the vicinity of the phase separation binodal is above the melting point of the Abrikosov defect crystal, at least for the conditions simulated in this study. This defect melting transition has been studied both theoretically [53, 54] and numerically [5, 55] in the (fully occupied) frustrated 2D XY model. Franz and Teitel report a defect melting transition occurring at $T_{\rm m}\simeq 0.007J$ [5], which is much cooler than the lowest temperature of $T=0.014J$ studied here.

In Figure 8a, we show the energy per subunit as a function of increasing concentration. At strong frustration, we see that there is a well defined global minimum at some critical concentration, $\Phi_{\rm min}$. The location of this minimum corresponds to the concentration at which the sponge phase can uniformly achieve the preferred hole size and density over the entire lattice. That is, at this particular density, $\Phi_{\rm min}$, the space-filling vortex sponge achieves both the selected spacing and core size of vortices described in Sec. III.2. Above this point, the holes are forced to fill in to a higher density state leading to additional phase strain around the defect core and an increased energy per subunit. Conversely, at concentrations lower than $\Phi_{\rm min}$, there is an excess of free space that can be filled in for an energetic gain. This excess space results in a larger number of conformational degrees of freedom, leading to an increase in structural heterogeneity. Hence, we use the location of the minimum of energy density as function of concentration as a means of delineating the boundary between the heterogenous network and quasi-uniform vortex sponge.

In section III.2, we introduced the circular cell approximation to describe the energetics of the homogeneous sponge phase. This simplified model predicts that the sponge phase should have an optimal hole density, $\Phi_{\rm hole}$, controlled by the ratio of $\varphi/(\sigma/J)^{2}$ and can be obtained from the transcendental equation for the optimal hole size given in eq. (29). In the limit that $\varphi\ll(\sigma/J)^{2}$ (which corresponds to the zero-temperature definition of weak frustration) we find that $\Phi_{\rm hole}\simeq\varphi/(\sigma/J)^{2}$. Conversely, in the limit that $\varphi\gg(\sigma/J)^{2}$, $\Phi_{\rm hole}$ can be expressed as a polynomial of $\varphi/(\sigma/J)^{2}$ (see eq. (32)). This can be related to the concentration that minimizes the per subunit aggregation energy via the relation:

[TABLE]

In Figure 8b, we test this prediction for the dependence of hole density on frustration and the cohesion-to-stiffness ratio. Comparing our simulation results to the two asymptotic limits, as well as to the full non-linear solution of eq. (29), we find remarkably good agreement with our simplified model of optimal bulk vortex sponges. This comparison confirms the basic trend that $\Phi_{\rm hole}$ increases with frustration and, more specifically, shows that our simulation parameters necessarily fall in the crossover range between the two asymptotic regimes, as this is the regime where self-limiting assembly occurs. Similarly, the dimensionless energy density, $\epsilon/J\varphi$, of the optimal sponge phase can also be obtained from the solution to eq. (29). Comparing this result to our simulation data (Figure 8c) we again find good qualitative agreement between simulation and theory.

Lastly, we observe from Fig. 8a that the energy density landscape near the minimum becomes less convex as frustration is decreased and is essentially flat in the limit of weak frustration. This nearly linear behavior is consistent with the coexistence between a dispersed monomer gas and a single condensed defect bulk observed below the binodal. Conversely, the results described in sec. IV.3 suggest that the coexistence at high frustration is more complex, consisting of phase separation between tree-like self-limiting aggregates of finite width and localized clusters of vortices that can at least partially condense on the self-limiting aggregates themselves.

Throughout this section, we saw that the dispersed/SLA phase found at dilute concentration and strong frustration undergoes a percolation transition to a bulk sponge phase at intermediate concentration. Near the percolation threshold, the internal structure is heterogeneous and marked by the presence of loops and branches of highly variable size. However, the structure becomes more uniform with increasing concentration and gradually approaches the quasi-uniform bulk state. The internal structure of this state is similar to the Abrikosov lattice found in the uniformly frustrated XY model, only with a punctured disc surrounding each defect core, with relatively uniform size for fixed values of $\varphi$ and $\sigma/J$. Additionally, we saw that the structure of this high concentration assembly is locally similar to that of the topologically defective condensate found at weak frustration, confirming that this latter phase represents coexistence between the low concentration dispersed phase and the high concentration homogeneous sponge. In the following section, we turn our attention to this phase separated state and the nature of the transition from self-limiting to bulk condensation at low to intermediate concentration driven by decreasing frustration.

V Frustration Driven Transition from Self-limiting to Bulk Condensation

Below the percolation concentration, assembly of the lattice GFA model exhibits a state of self-limiting aggregation above a critical concentration and, below that, bulk phase separation into vapor and bulk, defect-sponge states. Previously, we found that the critical frustration $\varphi_{c}$ increased with the ratio of cohesion to stiffness ratio as would be expected from zero-temperature arguments [1], but also exhibited a temperature sensitivity that cannot be explained by energetics alone. In this section, we consider thermodynamic effects of finite temperature and their impact on the frustration-driven transition between these states.

The phase boundary separating the self-limiting and condensed states of assembly can be described as the set of conditions where the difference in free energy density between these two states is zero. Thus, we can write down a simple equation for the phase boundary as:

[TABLE]

where $\epsilon_{\rm sla}=C_{0}J^{1/3}\sigma^{2/3}\varphi^{2/3}$ is the excess energy density of the self-limiting domains (see eq. (17)) and $\epsilon_{\rm bulk}\simeq C_{1}J\varphi$ is the excess energy density of the areal (holey) vortex array. We define $\Delta s$ as the entropy density difference between these two states (i.e. the specific entropy change). Here it is important to note that, while the prefactor $C_{1}$ is a constant in the limit of vanishing defect core size (see eq. (21)), it will generically have some dependence on the dimensionless ratio $\varphi/(\sigma/J)^{2}$ which selects the optimal hole size. However, as our circular cell approximation suggests that $C_{1}$ is only weakly dependent on $\varphi/(\sigma/J)^{2}$ (see eq. (31)), we will assume a constant value of $C_{1}\approx 1$. This is in agreement with the results for the dimensionless bulk energy depicted in Figure 8c. When $T=0$, the entropy difference between these states can be neglected and this equation can be solved to find the critical frustration separating these two states:

[TABLE]

with $\varphi>\varphi_{\rm c}$ (i.e. strong frustration) corresponding to the self-limiting state and $\varphi<\varphi_{\rm c}$ (i.e. weak frustration) corresponding to the phase separated state. Note that this criterion for $\varphi/\varphi_{\rm c}$ greater than or smaller than unity is equivalent to consideration of the ratio of optimal vortex spacing to finite domain thickness, $\ell_{\rm v}/\ell_{\rm d}$, originally derived in Reference [1]. Initial comparison of this scaling prediction to simulation results extracted from the lattice GFA model suggested, at best, a substantially lower power-law dependence of $\varphi_{\rm c}$ on $\sigma/J$ and an unpredicted shift with temperature to lower frustration. Although surprising, consideration of this result as the zero-temperature limit of equation (34) suggests the important role that entropy plays in stabilizing the self-limiting phase of assembly relative to an unlimited bulk state. Indeed, this scenario agrees with the result of a recently studied chain model of frustrated assembly [34] that suggests that finite temperature—at least in one dimensions—is a necessary condition for self-limiting assembly. Similarly, it suggests that thermal effects that couple to the relative entropy of assembly, ignored by considerations of the ground state energy alone, are strong under conditions where the lattice GFA model exhibits self-limiting assembly.

While the entropy differences between self-limiting and bulk states of assembly can in principle be computed by thermodynamic integration [56], it not necessarily possible to attribute those quantities to the microscopic origins of their underlying fluctuations. Instead, we analyze several different sources of fluctuations that can distinguish between the overall specific entropy of these two states of assembly. Based on this comparison, we propose a simplistic mean-field type theory for the entropic free energy difference between self-limiting and bulk condensed states, and use it to understand the nature and magnitudes of finite-temperature effects on the critical frustration.

V.1 Conformational entropy

Conformational entropy refers to the entropy associated with fluctuations in aggregate shape and internal order, a property assembled aggregates [57] have in common with macromolecular objects [58]. In the lattice GFA model, aggregates can take on a variety of different conformations that deviate from their locally favored ground state morphology by bending, branching, looping or changing width via capillary (edge) fluctuations. Additionally, the individual subunits in our model have an internal orientational degree of freedom that can fluctuate around its preferred arrangement, which constitute effective “spin wave” fluctuations around the ground state order.

While all of these conformational fluctuations have an associated entropy that contributes to the total of each phase, they do not necessarily distinguish between the two states thermodynamically. In particular, we argue that the specific entropy associated with spin fluctuations and capillary edge fluctuations are not likely to be very different in self-limiting and bulk (vortex sponge) states. In Appendix B we show that phase-fluctuations are largely governed by $\beta J$ and are fairly insensitive to frustration as well as large scale morphological differences between bulk sponge and self-limiting domains. Additionally, in Appendix C we compare the capillary edge fluctuations between the bulk and self limiting states. Since the relatively large core size of the vortex sponge permits a comparable range of edge fluctuations to the self-limiting state, we find that free edge fluctuations of both bulk and self-limiting states are well described by a simple model for capillary fluctuations along a two-dimensional ribbon. The quantitative similarity of these two types of conformational fluctuations suggests that they are unlikely to contribute significantly to the free energy difference between bulk and self-limiting states.

Unlike edge and phase fluctuations, there is a more drastic difference between the 2D paths of the backbones of self-limiting ribbons (i.e. fluctuating “worm-like” ribbons) and bulk sponges (i.e. organized arrays of holes). Here we consider the conformational entropy associated with the backbone of self-limiting domains, akin to the bending fluctuations of a polymer [59]. In the condensed phase, these fluctuations are highly-constrained by the energetic preference for closed loops of a given size and spacing. Therefore, to a first approximation, we estimate that there is no entropy of backbone domain fluctuations in the bulk phase.

We consider a finite-width ($W^{*}$) ribbon domain to have the effective bending elastic free energy of worm-like polymer [60]:

[TABLE]

where $\mathbf{r}(s)$ is the trajectory along the backbone as a function of arc length $s$, $\partial^{2}\mathbf{r}/\partial s^{2}$ is the instantaneous backbone curvature and $B$ is the bending stiffness per unit length. In the context of our model, the bending stiffness of the self-limiting aggregates can be written as:

[TABLE]

where $W_{*}\simeq\ell_{\rm d}=(\sigma/J)^{1/3}\varphi^{-2/3}$ is the preferred self-limiting width. This form derives from considering inverse square radius terms of the excess energy density of finite-width annular domains of the continuum model in the limit of large radius (see Appendix D). Using this, we can evaluate the effect of bending stiffness on conformational fluctuations by calculating the tangent-tangent correlation:

[TABLE]

where $\mathbf{t}=\partial\mathbf{r}/\partial s$ is a vector tangent to the aggregate backbone and $\ell_{\rm p}=\beta B$ defines the persistence length beyond which bend deformations become important. The tangent-tangent correlation measured along the backbone of our self-limited aggregates is given in Figure 9. From this, we observed two different scaling regimes. The first of which is a short range ($s\lesssim 10$) regime where the correlation length has a weaker dependence on frustration than that predicted by the bending stiffness (i.e. $\sim\varphi^{-4/3}$). This short-range behavior is likely the result of both intra-aggregate branching which decreases correlations over a length proportional to the average distance between branches (which range from $~5-6$ according to analysis Fig. 5) as well as fluctuations at the microscale of the square lattice that are not captured by a continuum bending description (see Appendix E). The second is a longer-range regime ($s\gtrsim 10$) where the linear decay of tangent correlations suggests a persistence length that follows the predicted $\sim\varphi^{-4/3}$ scaling. This suggests that, notwithstanding branching along the aggregate backbone, the long range bending fluctuations are well described by an effective semi-flexible ribbon whose backbone orientation fluctuates over length scale $\sim\ell_{\rm p}$ that grows as frustration is reduced and domain width increases.

Evaluating the partition function for a worm-like chain of fixed backbone length $L_{\rm b}$ (see Appendix E) it is straightforward to compute the entropic contribution from bending fluctuations:

[TABLE]

where $\bar{N}=L_{\rm b}/\lambda$ defines the effective number of bendable segments, and $\lambda\approx a$ is a microscopic cutoff for bending fluctuations. Notably, this model predicts that the conformational entropy is extensive in the contour length of a self-limiting ribbon domain. Assuming the ribbons to be roughly rectangular with constant thickness, the number of subunits per ribbon is simply $n=W_{*}L_{\rm b}/a^{2}$ and that the bending entropy of the bulk phase is neglible, we have a specific entropy difference due to conformational fluctuations

[TABLE]

where we are neglecting the much weaker logarithmic dependence encoded in the persistence length of domains.

V.2 Translational entropy

In addition to the conformational entropy, any system of geometrically frustrated assembly will have translational entropy related to the discrete positioning of subunits. Notably, translational entropy in aggregating systems is size-dependent, in effect, because subunits within a single aggregate share a single center-of-mass degree of freedom. This translational entropy is encoded within the aggregate mass distribution, $\phi_{n}$, which denotes the total area fraction of the system contained within aggregates of size $n$. In the context of ideal aggregation theory [9, 61], the mixing entropy of the distribution can be written as

[TABLE]

where $A_{\rm tot}$ is the total system area. In particular, we are concerned with the translational entropy of aggregated structures, i.e. for $n>n_{m}$ where $n_{m}$ is the threshold value of cluster size delineating dispersed from aggregated structures. Throughout this work, we will take $n_{m}=9$, as this aligns with the approximate location of the minima of the aggregate concentration separating the monomer and SLA populations over the range of conditions considered here [1]. Summing over this part of the distribution and dividing by the number of aggregated subunits $\Phi_{\rm agg}A_{\rm tot}/a^{2}$, we have the specific translational entropy of aggregated subunits:

[TABLE]

which is straight forwardly computed from numerical simulations. In Figure 10, we show the specific translational entropy of aggregated subunits, across the condensation transition, for several different values of $\Sigma/J$ 222For this analysis, we use a constant value of aggregation threshold of $n_{m}=??$. From this, we see that, below the critical frustration when assembled subunits have condensed into a single bulk structure, the specific translational entropy of aggregates plummets to negligible values relative to the self-limiting state. This reflects the fact that, in the condensed state, the degrees of freedom of a single collective center of mass are distributed among a macroscopically large number of subunits in the structure, whereas in the self-limiting state, the translational entropy of aggregates is shared by a finite number of subunits, leading to a much larger specific entropy.

We estimate the specific translational entropy of self-limiting aggregates by assuming that the sum over the mass distribution in eq. (42) is dominated by the peak in the distribution, which we have shown previously [1] occurs at mass $n_{*}$ roughly equal to the mass of a square aggregate of ideal self-limiting width, i.e. $n_{*}\approx(W_{*}/a)^{2}$. From this, we estimate the value of the specific translational entropy of self-limiting aggregates as:

[TABLE]

where we used the approximation $\phi_{n_{*}}\approx\Phi_{\rm agg}$. We note that this approximation undercounts the entropy associated with the breadth of the mass-distribution of self-limiting aggregates, in particular the contributions from the exponential tail for which $n\gg n_{*}$. However these large aggregates contribute a relatively small amount of specific entropy due to their large mass and so the effect of their omission is negligible.

Taking the translational entropy of the bulk state to be zero, this gives us the estimate for specific entropy change

[TABLE]

which, neglecting logarithmic factors, also falls off with a power-law with self-limiting size, $W_{*}$, like the conformational entropy estimate.

V.3 Entropic dependence of self-limiting to bulk transition

The prior two sections give estimates of the greater specific entropy of self-limiting aggregates relative to the bulk aggregates. We estimate that both conformational entropy (bending fluctuations) and translational entropy per subunit of self-limiting aggregate decrease with self-limiting width (albeit with different powers of $W_{*}$), in effect due to the larger number of subunits that share the same effective degrees of freedom in larger aggregates. Here, we use these estimates to explore the effect that this excess entropy of the self-limiting aggregates has on the SLA to bulk transition and show that its main role is to depress the value of the critical frustration, $\varphi_{c}$, separating the self-limiting and defect bulk phases.

For simplicity we consider only one of the two entropic contributions and its dependence on frustration and cohesion to stiffness ratio. Based on an estimate of effective bond length $\lambda\approx 10a$ from Fig. 9 in eq. (39) and values of specific entropy in Fig. 10, we expect translational entropy to dominate the specific entropy of self-limiting aggregates for the range of simulation parameters explored here. Thus, using $n_{*}\approx(W_{*}/a)^{2}$ and $W_{*}\sim\varphi^{-2/3}(\sigma/J)^{1/3}$ we estimate a scaling form for the specific entropy change in the self-limiting to bulk transition

[TABLE]

where we neglect the variation of the logarithmic term $\ln(\Phi_{\rm agg}/n_{*})$ near the transition. Substituting this into the difference in specific free energy between self-limiting and condensed states in eq. (34) we have:

[TABLE]

where $C_{\rm T}$ is a numerical prefactor. It is straightforward to solve $\Delta F(\varphi)=0$ for $\varphi_{c}$, and the solutions are plotted in Figure 11, where we convert to dimensionless frustration, $f_{c}$, to compared to simulation results for three values of reduced temperature $k_{B}T/\Sigma$. Here a value of $C_{T}\simeq 45$ was chosen to best match the values extracted from simulation. Notably the curves $f_{c}(\Sigma/J)$ indeed appear to show a roughly linear dependence on the cohesion to stiffness in the range of simulated values, and further shift to increasingly lower values as the reduced temperature $k_{B}T/\Sigma$ is increased. This is consistent with trends observed in the simulations. Physically, the shift to low values $f_{c}(\Sigma/J)$ with increased temperature is a result of the relatively larger specific entropy of the self-limiting aggregates. Indeed the large gap between the ground state $T=0$ value of $f_{c}$ (dashed line) and the values of $k_{B}T/\Sigma$ for the simulations (solid curves) imply that thermal effects are quite strong for simulated assemblies, and that the excess entropy of self-limiting aggregates over bulk states substantially accounts for the broad window of self-limiting assembly exhibited in our simulations.

This mean-field type calculation further allows us to estimate the scale at which thermal effects become strong. It is straightforward to show the perturbative correction to the $T=0$ critical frustration value at low-temperature

[TABLE]

where $\varphi_{\rm c}(T=0)=(C_{0}/C_{1})^{3}(\sigma/J)^{2}$ is the $T=0$ ground-state transition value. We note that our simulations results are carried out for fixed values $k_{B}T/\Sigma$, and hence, increasing cohesion to stiffness ratio (i.e. $\Sigma/J$) corresponds to increasing values of $k_{B}T/J$. In the large temperature limit, we find:

[TABLE]

Hence, we expect a crossover from $\varphi_{\rm c}\propto(\sigma/J)^{2}$ to $\varphi_{\rm c}\propto(\sigma/J)^{1/2}$ around a cohesion scale inversely proportional to the fixed value of reduced temperature $k_{B}T/\Sigma$. From Fig. 11 and the simple mean field estimate, we expect that this crossover occurs in the range of $\Sigma/J\approx 10^{-2}-10^{-1}$. Based on this, we conclude that the apparently linear dependence of critical frustration on $\Sigma/J$ may in fact be a manifestation of a broad crossover between an athermal regime to a strongly thermal regime from low to high values of cohesion to stiffness.

We note by inspection of Fig. 11, that accessing the energetically-dominated athermal regime of the lattice GFA model, where $\varphi_{\rm c}\propto(\sigma/J)^{2}$, would likely require exceptionally low values of $\Sigma/J\lesssim 10^{-3}$ and correspondingly minuscule values of $f$. However, restrictions on the periodic boundary conditions (i.e. where $L$ is an integer multiple of $2/f$) require prohibitively large lattice sizes to access this range of ultra-low frustration.

VI Conclusion and Discussion

In this article, we studied a 2D lattice model of geometrically frustrated assembly to explore the interplay between translational degrees of freedom and frustration over the entire range of subunit concentration. This was done with a particular focus on understanding the various thermodynamic phase transitions separating the dispersed/self-limited aggregate phase from the topologically defective macroscopic bulk phase.

In section IV.1, we studied how, at strong frustration (i.e. above the phase separation binodal), the GFA model evolves from an increasing concentration of shape-fluctuating, finite-width domains above the aggregation transition to a heterogeneous percolated network at intermediate concentrations, and ultimately to a relatively uniform “sponge” of evenly spaced vortices, characterized by large-area cores. We find that the location of this percolation transition shifts to somewhat higher concentration values with increased frustration, an observation that is consistent with both the decreasing overall mass of aggregates and an increase in their branch density with increasing frustration. Additionally, near the percolation threshold, the structure of the sponge phase is quite heterogeneous and is predominantly characterized by large fluctuations in hole size and spacing, which eventually evolves at high concentration to evenly spaced holes of a size and separation that are selected by the energetics of frustration screening (i.e. a “holey” Abrikosov lattice). Notably, in characterizing the loopy versus tree-like structure of the self-limiting domain, we find a very broad crossover, extending from concentrations well below the percolation point up to optimal sponge density. This analysis is summarized in Fig. 12, showing the “short cycle” fraction of aggregated clusters as a function of concentration. This rather gradual crossover between extended finite-width, worm-like domains and defective loops of well-defined size suggests a more complex scenario of coexistence between these morphologies at strong frustration in which defect loops effectively condense in small clusters onto extended tree-like ribbon domains (akin to grapes on the vine). Collectively, these “two-phase” aggregates ultimately mix into a single network, reflecting relatively even proportions of tree-like and “short loop” structures. Increasing concentration generically favors the denser defect-loop morphology, ultimately leading to quasi-uniform sponge.

In section V, we studied contributions to the assembly entropy and how these tend to thermodynamically stabilize self-limiting aggregates relative to the bulk condensed state. We considered simple models of conformational and translational entropy of self-limiting aggregates, which predict that the specific entropy of both effects depends on a reciprocal power of the self-limiting domain size. Based on the scaling of domain size on frustration and cohesion to stiffness ratio in the translational entropy, a mean-field model of the free energy difference between self-limiting and bulk states was used to predict the finite-temperature dependence on the critical frustration, $\varphi_{c}$, separating these states. Comparison between this simple model and simulation results suggest that our GFA simulations are in a strongly thermal regime, in which favorable entropy of dispersed self-limiting aggregates substantially depresses $\varphi_{c}$ below its expected $T=0$ scaling $\varphi_{c}(T=0)\sim(\sigma/J)^{2}$. Consistent with previous observations [1], this model shows that $\varphi_{c}$ decreases with increasing temperature, in effect a thermal melting of condensed, spongey bulk states into a disperse suspension of extended and finite-width aggregates. The apparent linear dependence on cohesion to stiffness observed for $\varphi_{c}$ at fixed temperature is consistent with a broad crossover regime from square dependence to a square-root dependence at low versus high values of $\sigma/J$, although direct numerical evidence of the $T=0$ scaling appears likely out of the accessible range of current simulations.

In summary, this sheds light on key aspects of geometrically frustrated assembly at finite temperature and concentration. Firstly, this shows that, while energetic considerations of elastic cost of frustration and cohesion control the size and structure of the competing states of aggregation, conformational and configurational entropy play a significant, and heretofore underappreciated, role in stabilizing states of self-limiting aggregation. Additionally, this study addresses the fate of self-limiting assembly far above the conditions of ideal aggregation, at concentrations where aggregates are strongly interacting. We demonstrate the concentration bound on the range of conditions that support self-limited assembly, and study how the system evolves from a dispersed solution of aggregates to a relatively ordered and uniform bulk sponge via an intermediate percolating and heterogeneous network mesophase. This scenerio provides a framework for connecting prior studies of bulk and dilute geometrically frustrated assemblies that were previously treated as unrelated. For example, much work has been done studying the bulk defect phase of bend-nematic liquid crystals [30, 63, 64]. Meanwhile, the same mechanism of frustration has been theorized to induce self-limitation under dilute conditions [1, 31, 32]. Thus, the framework outlined here could provide a useful guideline for extending earlier experimental work to the dilute phase, thereby opening up a new avenue for designing novel functional material. Conversely, the pathway could go the other way, helping to extend work done on dilute self-assembly to build complex bulk material from the ground up. Additionally, we note several previous studies of similar phase behavior in the context of colloidal-gels of attractive particles [65]. While many models of these systems assume that the locally aggregated structures form far from equilibrium (i.e. via rapidly quenching [66]), there has been some investigation into the possibility that at least strongly-metastable favored local structures form as a result of geometric frustration of polytetrahedral packing of a small number of sticky spheres [67, 68]. In this context, we anticipate that the lattice GFA model may indeed by a useful analog of such systems in which the frustration may be systematically tuned to adjust the equilibrium dimensions of favored local structures (e.g. domain sizes and inter-defect space). Thus allowing one to probe the statistics, structure and mechanics of the emergent networks.

While the work presented above provides a complete picture of the effect of frustration at finite temperature on assembly over the full range of concentration, several aspects remain to be understood. For one, dispersed monomers and self-limiting aggregates are two types of low concentration solution phases (separated only by a pseudo-critical aggregation crossover) it is far from clear if or how the percolation intersects or merges with the phase separation binodal. Here, there could be a tri-critical point similar to that studied in the Blume-Emory-Griffiths model of phase separation in thin films of super-fluid He3-He4 mixtures [69, 70]. Additionally, the observation that the percolating clusters consist of small, dense regions connected by thin ribbons suggests the possibility of a secondary binodal where the heterogeneous sponge phase is viewed as a coexistence between the homogeneous sponge and self-limiting assembly. Our current speculation (indicated graphically in Fig. 12) is that such a critical point could be hidden within the window between the percolation threshold and the optimal density of the uniform vortex array. However, equilibration of the nearly bulk, yet strongly-fluctuating, configuration of states in this regime is especially challenging.

Beyond the transition between strong- and weak-frustration, the nature of the high frustration (i.e. $f\rightarrow 1/2$) limit is particularly unclear. As frustration increases, the aggregate widths decrease until $W\sim a$. For the parameters used in Figure 12, this happens around $\varphi\simeq 0.06$. Above this, the aggregates are expected to form living polymeric chains of single subunit width. As geometric frustration only acts on closed cycles of connected bonds, these 1D spin chains could avoid the adverse effects of frustration by forming fractal network structures that are self-avoiding and tree-like. Although the concept of a self-limiting aggregate becomes poorly defined in this limit and it is unclear what bearing this regime of phase space has on the behavior of physical assemblies, we speculate that certain conformational statistics, if not also critical phenomena, would likely be described within the universality class of branched lattice animals [71].

Acknowledgements.

The authors are grateful to Jon Machta and Christopher Amey for helpful and detailed discussions about MC simulation approaches and to Kyle Sullivan and Christian Koertje for valuable feedback on this manuscript. This work was supported by the National Science Foundation through Award No. NSF-DMR 2028885 and 2349818. All simulations were performed using the UMass Unity cluster located at the Massachusetts Green High Performance Computing Center.

Appendix A Local Orientational Ordering of Defect Holes

The orientational ordering of the intra-assembly holes can be quantified via the local hexatic order parameter[72]:

[TABLE]

which measures how closely the local arrangement of holes aligns with a perfect triangular lattice. Here $\theta_{nn}$ defines the angle that the vector between the center of two nearest neighboring holes makes with respect to the $\hat{x}$ direction and $N_{nn}$ is the number of nearest neighbors around a given hole. Averaged over every hole in an aggregate, the ensemble average $\langle|\Psi_{6}|\rangle$ assumes a value ranging from 0, corresponding to the total absence of orientational order, to 1, which corresponds to a perfect triangular lattice. This quantity is evaluated over the entire range of concentration for two different values of strong frustration (see Figure 13ab). Here, we observe an initial rise in local orientational order that coincides with the percolation transition. At weak frustration, we observe a similar increase in local orientational order, only here it coincides with the onset of phase separation, which happens at a lower value of concentration than the percolation transition (see Figure 13c). In both cases, the orientational order continues to increase with concentration from $\langle|\Psi_{6}|\rangle\sim 0.4$ to $\sim 0.7$. While this suggests some amount of local orientational ordering, the overall effect is not very strong and there is no discernible differentiation between the heterogeneous, homogeneous and condensed states of assembly under these conditions. This suggests that any amount of apparent ordering is the result of the uniformity in hole size and spacing, rather than any strong hexagonal arrangement.

Appendix B Orientational Entropy

The conformational entropy due to fluctuations in subunit orientation can be calculated within the Einstein crystal approximation [73], which assumes that the lattice link between each pair of neighboring subunits is an independent variable. The validity of this assumption can be immediately seen by looking at the distribution of the gauge invariant phase difference, $\Delta\theta_{ij}-A_{ij}$, along occupied links (see inset of Figure 14) and observing that they are Gaussian in nature. Furthermore, this suggests that we can obtain the entropy of these fluctuations by taking the log of their variance, i.e.

[TABLE]

Noting that the Hamiltonian for this system (eq. (1)) is quadratic in $\Delta\theta-A$ allows us to calculate the variance using the equipartition theorem [74]. Therefore,

[TABLE]

which only depends on the ratio of temperature and spin-stiffness. Comparing this to the orientational entropy measured from a range of simulations sweeping across the transition (see Figure 14), we find good agreement in both the self-limiting and condensed phase. Thus, we see that fluctuations in subunit orientation are largely independent of phase and frustration, suggesting that the change in orientational entropy across the transition should be negligible.

Appendix C Capillary Fluctuations

The aggregates observed in our finite temperature simulations exhibit noticeable capillary fluctuations in aggregate width. The importance of these fluctuations can be evaluated by considering an infinitely long ribbon of varying width, described by $W(z)=W_{*}+\delta W$, where $z$ is the position along the aggregate backbone and $\delta W$ is the deviation of the width away from $W_{*}$. The variation in aggregate width leads to a change in free energy [15]:

[TABLE]

where

[TABLE]

Introducing the Fourier transform width fluctuation:

[TABLE]

we can re-write eq. (52) as:

[TABLE]

The equipartition theorem then tells us that:

[TABLE]

Fourier transforming back to real space, we get the thermal width fluctuation [15]:

[TABLE]

Comparing this to the square of the mean width, $\langle W\rangle^{2}$, we expect that capillary fluctuations will be negligible when:

[TABLE]

The theory discussed above describes harmonic corrections to the mean field theory, which assumes that $\langle W_{*}\rangle$ is independent of temperature and only depends on the dimensionless parameters $f$ and $\Sigma/J$. While this is true in the limit of large $\beta J$, the mean width that we measure via simulation varies with temperature (see Fig 15a). This dependence is largely related to the tendency of increased temperature to drive aggregates towards dissolution, thereby lowering the fraction of subunit mass contained within an aggregate [1]. As this effect is not captured by the mean-field theory, it is difficult to make a direct comparison between theory and simulation. However, the interpretation of the quantity $\langle|\delta W|^{2}\rangle/\langle W_{*}\rangle$ is still the same and our numerical results are at least in qualitative agreement with the theory (see Fig. 15b). Perhaps more presciently, the simulation data shows that—for temperatures near the SLA to bulk transition—the magnitude of width variation is small compared to the mean. This suggests that capillary fluctuations do not provide the dominant correction to the mean field prediction for the critical value of frustration.

Appendix D Bending Energy

In order to find the bending modulus of our self-limiting aggregates, we calculate the curvature dependent term of the excess energy density of an annular domain with finite width $W_{*}$ (see Fig. 16). To start, we introduce the cross-sectional coordinate,

[TABLE]

of our annular domains in terms of the midpoint radius, $R$. This allows us to write down the solution to the Euler-Lagrange equation (eq. (8)) for an annular domain as:

[TABLE]

where $s$ defines the winding around the annular domain. For convenience, we define the gauge field over the annulus as:

[TABLE]

Using these, we can calculate the energy via:

[TABLE]

where, for convenience, we have defined:

[TABLE]

Dividing by the area and minimizing with respect to $s$, we find the optimal phase winding:

[TABLE]

Substituting these back into the Hamiltonian and introducing the parameterizing the infinitesimal arc length as $\rm d\theta=\frac{\rm ds}{R}$, we find the bending energy as a function of curvature:

[TABLE]

This allows us to read off the bending modulus

[TABLE]

given in equation (37).

Appendix E Bending Entropy

In order to calculate the conformational entropy associated with bending, we treat the self-limiting aggregates as 1D chains of $\bar{N}=L_{\rm b}/\lambda$ segments. Here, $L_{\rm b}$ is the total length of the aggregate backbone and $\lambda$ is a microscopic cutoff related to the shortest bendable unit. Thus, we can rewrite the continuum bending energy defined in equation (36) as the discrete sum:

[TABLE]

where $\theta_{i}$ is the angle each segment makes with respect to the $x$-axis. Using this, we can write the partition function as

[TABLE]

Assuming that the orientation of each segment is independent, this can be expressed as a product of $\bar{N}$ Gaussian integrals and evaluated to find:

[TABLE]

where $\ell_{\rm p}=\beta B$ is the persistence length of the chain. From here, it is easy to obtain the free energy:

[TABLE]

This can be differentiated to find the conformational entropy of bending given in equation (39):

[TABLE]