Kibble-Zurek mechanism in curved elastic surface crystals

TL;DR

This paper demonstrates that the Kibble-Zurek mechanism, originally studied in planar systems, also applies to curved elastic surface crystals during stress-induced phase transitions, revealing universal defect formation laws.

Contribution

It extends the validity of Kibble-Zurek scaling laws to non-thermal, first-order transitions in curved, topologically nontrivial geometries like spheres and tori.

Findings

Defect densities follow KZ power laws regardless of surface topology.

Nucleation sequences match experimental observations in spherical crystals.

KZ scaling laws are applicable to a broader class of phase transitions.

Abstract

Topological defects shape the material and transport properties of physical systems. Examples range from vortex lines in quantum superfluids, defect-mediated buckling of graphene, and grain boundaries in ferromagnets and colloidal crystals, to domain structures formed in the early universe. The Kibble-Zurek (KZ) mechanism describes the topological defect formation in continuous non-equilibrium phase transitions with a constant finite quench rate. Universal KZ scaling laws have been verified experimentally and numerically for second-order transitions in planar Euclidean geometries, but their validity for discontinuous first-order transitions in curved and topologically nontrivial systems still poses an open question. Here, we use recent experimentally confirmed theory to investigate topological defect formation in curved elastic surface crystals formed by stress-quenching a bilayer material. Studying both spherical and toroidal crystals, we find that the defect densities follow KZ-type power laws independent of surface geometry and topology. Moreover, the nucleation sequences agree with recent experimental observations for spherical colloidal crystals. These results suggest that KZ scaling laws hold for a much broader class of dynamical phase transitions than previously thought, including non-thermal first-order transitions in non-planar geometries.

Figures

Equations

$$\displaystyle\mathcal{E}=\frac{k}{2}\int_{\omega}d\omega\biggl{[}\gamma_{0}\left(\nabla u\right)^{2}+\frac{1}{12}\left(\triangle u\right)^{2}+au^{2}+\frac{c}{2}u^{4}-\Gamma(u)\biggr{]}$$

$${\tau }_0\frac{\partial }{\partial t}u={\gamma }_0△u-\frac{1}{12}{△}^{2}u-au-c{u}^3-{\delta }_u\Gamma (u)$$

$$a=a_c-\frac{3c}{4}{\Sigma }_e.$$

$${\Sigma }_e(t)=\mu t$$

$$\frac{d}{dt}U\approx \frac{c\mu t}{12{\gamma }_0^4}U.$$

$$\frac{d}{d{t}^{\prime }}U\approx \frac{c{t}^{\prime }}{12{\gamma }_0^4}U.$$

Full text

Introduction. Topological defects are localized perturbations that break the global symmetry of ordered solids or liquids. These defects influence the elastic, magnetic, electronic and optical properties in many natural and man-made systems. Examples range from the domain structures formed in the early universe [1, 2] to liquid crystals [3, 4, 5, 6], grain boundaries in ferromagnets and colloidal crystals [7, 8], and charge transport in graphene [9, 10] and superconductors [11, 12, 13, 14]. Topological defects can be created by varying a control parameter, such as temperature or a magnetic field, rapidly across a phase transition. Understanding the complex nonequilibrium dynamics induced by these nonadiabatic quenches remains an important theoretical challenge. Early progress in the characterization of topological defects was made by Kibble [1] in 1976, while studying domain formation during the rapid cooling of the early universe [15]. About a decade later, Zurek [11] showed how the defect density is related to the quench rate in general second-order phase transitions; this important breakthrough also advanced significantly the understanding of topological defect formation in quantum superfluids and superconductors [16]. Since then, the Kibble-Zurek (KZ) power-law scaling predictions were confirmed experimentally and utilized in a variety of systems, including liquid crystals [5, 17], colloidal monolayers [18], ion crystals [19], Bose-Einstein condensates [20], superfluids [12, 21, 14] and cold atomic clouds [22]. However, since previous theory [1, 11, 23, 24] and experiments focused primarily on continuous second-order transitions in planar Euclidean spaces, relatively little is known about the existence of KZ-type scaling laws for other classes of phase transitions and in more complex geometries. Here, we show that analogous scaling laws hold for nonthermal discontinuous first-order phase transitions on simply and not-simply connected curved surfaces.

Topological defects in curved two-dimensional (2D) crystal structures arise in many biological and physical processes [26], from plant growth [27, 28, 29] and assembly of bacterial cell walls [30, 31], viral capsids [32, 33] and microtubules [34] to the targeted design of carbon nanotube sensors [35] and microlense fabrication [36]. On closed manifolds, Euler’s theorem [37] links the net charge of the topological defects to the genus $g$ of the underlying surface, imposing for example the 12 pentagons on a soccer ball. Depending on the details of the crystallization process (quench rate, geometric constraints, etc.), non-planar 2D crystals typically contain additional excess defects with zero total charge [38, 39]. Recent studies provided important experimental and theoretical insights into the defect statistics and formation dynamics in 2D colloidal crystals assembled on curved liquid-liquid interfaces [40, 41, 25, 42]. Another promising class of experimental systems are curved elastic bilayer systems [43], consisting of a soft substrate and a stiff surface film (Fig. 1A), which can develop hexagonal wrinkling patterns under lateral compression induced by surface swelling [44] or substrate depressurization [45, 46]. Such elastic surface crystals allow the realization of nontrivial shapes of genus $g>0$, such as toroids, which are difficult to achieve in liquids. Moreover, because the transition from the unwrinkled to the hexagonal phase is of first order [46], these soft matter systems also provide an ideal testbed to study generalizations of the KZ scaling laws.

The classical KZ argument for thermally induced second-order transitions builds on the fact that the correlation length $\xi$ and relaxation time $\tau$ diverge as $\xi\propto\vartheta^{-\nu}$ and $\tau\propto\vartheta^{-z\nu}$ (critical slowing down), respectively, when the temperature $\vartheta$ is varied to drive the system continuously from the high-symmetry (e.g. isotropic) phase to the lower-symmetry crystal phase. The critical exponents $\nu$ and $z$ encode the universality class of the transition. For a linear quench $\vartheta=\mu t$ with rate $\mu$, the system will not be able to relax defects during the time interval $|t|\lesssim t_{f}=\tau$, yielding the freeze-out condition $t_{f}\propto(\mu t_{f})^{-z\nu}$ or, equivalently, $t_{f}\propto\mu^{-z\nu/(1+z\nu)}$. The associated correlation length $\xi_{f}=\xi(t_{f})\propto\mu^{-\nu/(1+z\nu)}$ implies the KZ scaling prediction for the defect density at freeze-out, $\rho_{f}\propto\xi_{f}^{-d}\propto\mu^{d\nu/(1+z\nu)}$, where $d$ is the space dimension. The analysis below shows that this argument can be generalized to stress-induced discontinuous pattern formation transitions observed on the 2D surfaces of curved elastic materials. Furthermore, our computations identify a dynamical analogy with recent experimental observations [25] on colloidal crystal formation in curved liquid-liquid interfaces. Altogether, these results suggest that KZ-type scaling laws hold for a much broader class of dynamical phase transitions than previously thought. From a practical perspective, the subsequent analysis offers guidance for how to combine quench dynamics and surface geometry to control both the frequency and localization of topological defects.

I Theory

To investigate topological defect formation in curved geometries, we analyze an experimentally validated continuum model [46] for the surface wrinkling in elastic bilayer materials consisting of a soft core and a stiffer outer shell (Fig. 1A). This generalized Swift-Hohenberg (GSH) theory can reproduce quantitatively the experimentally measured equilibrium phase diagrams [46], but its dynamical implications have not yet been explored.

GSH theory for elastic surface crystals. The GSH equations follow from the nonlinear Koiter shell theory [47] by expanding the elastic energy of film and substrate in the dominant normal displacement field $u$ [46]. Measuring length in units of the film thickness $h$, the surface energy functional of the GSH theory reads [46]

[TABLE]

where $k=E_{f}/(1-\nu^{2})$ for a film with Young’s modulus $E_{f}$ and Poisson ratio $\nu$, and $d\omega$ is the surface element of the undeformed substrate. The nonlinear term $\Gamma(u)=\left[(1-\nu)b^{\alpha\beta}\nabla_{\alpha}u\nabla_{\beta}u+\nu b^{\alpha}_{\alpha}(\nabla u)^{2}\right]u$ represents stretching forces to leading order in the curvature tensor $b^{\alpha\beta}$, with surface gradient $\nabla$ and Laplace-Beltrami operator $\triangle$ (throughout, Greek indices run over the set $\{1,2\}$ and Einstein’s summation convention is used). Taking the variation of $\mathcal{E}$ with respect to $u$ and assuming overdamped dynamics, the surface wrinkling process is described by the GSH equation [46, 39]

[TABLE]

where $\gamma_{0}<0$, $c>0$ and $\delta_{u}\Gamma(u)$ denotes the functional derivative of the $\Gamma$-contribution to the energy functional $\mathcal{E}$. Explicit expressions [46] for the coefficients and $\delta_{u}\Gamma(u)$ are summarized in the SI Theory. With no loss of generality, we measure time in units of the damping time $\tau_{0}$ from now on, which is equivalent to setting $\tau_{0}=1$ in Eq. (1).

Linear stability analysis of Eq. (1) implies that wrinkling patterns form via a discontinuous transition (see Fig. 4 in Ref. [46]), when the control parameter $a$ falls below the critical value $a_{c}=3\gamma_{0}^{2}$, corresponding to an increase of the film stress $\sigma$ beyond the critical buckling stress $\sigma_{c}$. Defining the excess film stress $\Sigma_{e}=(\sigma/\sigma_{c})-1$, the bifurcation parameter $a$ is related to $\Sigma_{e}$ by [46]

[TABLE]

In the regime beyond but still close to the wrinkling threshold, $0<\Sigma_{e}\ll 1$, nonlinear stability analysis confirms [46] that the wrinkling solutions adopt a hexagonal pattern (Fig. 1) due to the $\delta_{u}\Gamma(u)$-term, which is breaking the $u\to-u$ symmetry of Eq. (1).

Stress-quenching of elastic surface crystals. To obtain scaling predictions for the topological defect formation in elastic surface crystals, we first solve Eq. (1) numerically for linear stress quenches

[TABLE]

with constant quench rate $\mu$, driving the system from the unwrinkled to the hexagonal phase. In all simulations, material parameters are chosen to match the experimental values reported in Refs. [39, 46] (SI Theory). The three essential differences compared with classical KZ scenarios are that (i) the quench (3) is nonthermal, (ii) the wrinkling phase transition is of first order (Fig. 4b in Ref. [46]), and (iii) the underlying substrate geometries are non-planar.

To break the symmetry of the initially unwrinkled surface $u=0$, a small stationary random field $\epsilon$ with $|\epsilon|\ll 1$ is added to the rhs. of Eq. (1) in simulations (SI Theory). This $\epsilon$-inhomogeneity effectively models initial imperfections in the film displacement, thus mimicking realistic experimental conditions [43]. We simulate Eq. (1) for the linear quench (3) and a given realization of $\epsilon$ using the algorithm described in Ref. [46]. Numerical results presented below are averages over $n$ different realizations of $\epsilon$, with $n$ specified on the corresponding graphs.

II Results

To illustrate the effects of locally varying curvature and surface topology, we compare simulations for spherical and toroidal surfaces.

Surface crystal growth under slow quenching. For slow quasi-adiabatic quenches ($\mu\to 0$), we find that the crystallization process is initiated at isolated nucleation sites and then spreads to cover the entire film (Fig. 1B,D). This behavior is observed for both spheres and tori (Movies 1 and 2). Details of the spreading dynamics and local crystal orientation become evident by reconstructing the corresponding planar crystal patterns (Methods). Starting from a random crystal site and one of its neighbors, we determine the relative positions of all other sites connected to this initial pair, resulting in the planar crystal representations of Fig. 1C,E. The time evolution in these graphs shows how initially separated crystal patches merge to form a single connected crystal covering the entire surface. The radial cuts appearing in Fig. 1C reflect the non-isometric character of the planar representation of spherical crystals, whereas the toroidal crystals unfold nearly isometrically (Fig. 1E). One can see however that, for both geometries, defects form only at the later stages when perfectly hexagonal crystal domains originating from different nucleation sites come in contact with each other. As discussed below, qualitatively similar crystallization sequences were observed in a recent experiment [25] that studied the deposition of charged colloids onto a spherical oil-water interface. To connect with the ideas of Kibble [1] and Zurek [11], we next use the GSH elasticity model to analyze the relation between topological defect formation and quench rate $\mu$, which is difficult to explore in curved colloidal systems due to experimental limitations on the particle deposition rates.

KZ-type scaling in spherical surface crystals. We start our numerical scaling analysis by considering the case of globally constant curvature as realized in spheres (Fig. 2; Movie 1). In our simulations of Eq. (1), the quench rate is varied over the range $\mu\in[5\cdot 10^{-8},10^{-4}]$, consistent with the assumption of an overdamped dynamics. Fixing the sphere radius $R/h=80$, we characterize the transition from the unwrinkled to the wrinkled phase in terms of the average squared displacement $\langle u^{2}\rangle$, where the brackets indicate an instantaneous surface average. In the adiabatic limit $\mu\searrow 0$, this order parameter vanishes in the unwrinkled phase, $\langle u^{2}\rangle=0$ for $\Sigma_{e}<0$, and jumps to a finite value $\langle u^{2}\rangle>0$ when the excess film stress $\Sigma_{e}$ crosses zero from below. By contrast, for non-adiabatic quenches, we find that the system remains longer in the unwrinkled phase, before eventually breaking symmetry at some finite positive value $\Sigma_{e}>0$ (Fig. 2D). Such delayed symmetry-breaking is also seen in the classical KZ mechanism, reflecting the critical slowing down in the relaxation dynamics near a thermal second-order phase transition. Yet, the wrinkling transition considered here is neither thermal nor of second order, as hexagonal patterns arise through a subcritical bifurcation [46]. To quantify the scaling behavior for this discontinuous first-order transition, we define the net freeze-out film stress $\Delta\Sigma_{f}=\Sigma_{e,f}-\Sigma_{0}$, with $\Sigma_{e,f}$ being the value at which the pattern amplitude is fully developed (Methods). The constant shift $\Sigma_{0}\approx 0.1$ is required to account for the material imperfections modeled by the $\epsilon$-inhomogeneity, as explained in Refs. [48, 49]. Our simulations confirm power-law scaling $\Delta\Sigma_{f}\propto\mu^{1/2}$ (Fig. 2E), implying that the freeze-out time diverges as $t_{f}\propto\mu^{-1/2}$.

The perhaps most interesting observable is the topological defect density $\rho_{f}$ at freeze-out $t_{f}$. Defects are crystal sites with coordination number $Z\neq 6$ and non-zero topological charge $s=Z-6$. To identify the $\mu$-dependence of $\rho_{f}$ in the GSH theory, we determined the coordination number for each lattice site from the Voronoi cells of the displacement field $u$ (Methods). The resulting Voronoi graphs show that the freeze-out density $\rho_{f}$ increases with the quench rate $\mu$ (Fig. 2A-C). After a quench is completed, defect pairs are expected to annihilate by grain boundary movements. We tested this hypothesis in simulations by stopping the quench at the freeze-out value, $\Sigma(t)=\Sigma_{e,f}$ for $t>t_{f}$, so that the spherical surface crystals could relax to a stress equilibrium. For all considered quench rates, we find that the defect density approaches constant asymptotic values, which converge to the equilibrium value $\rho_{\infty}$ as $\mu\searrow 0$ (Fig. 2F). Note that, although the net toplogical charge is always $-12$ in agreement with Euler’s theorem for hexagonal sphere tilings [37], charge-neutral pairs of penta- and heptagonal defects can reduce the elastic energy [50, 51, 52, 41, 39], resulting in a non-zero defect density $\rho_{\infty}$ even at equilibrium. Defining the relative excess defect density at freeze-out $t_{f}$ by $\Delta\rho_{f}=\rho_{f}-\rho_{\infty}$, our numerical data is consistent with a KZ-type power law $\Delta\rho_{f}\propto\mu^{1/2}$ (Fig. 2G).

KZ-type scaling for toroidal surface crystals. Local curvature variations can influence the equilibrium defect localization in curved elastic crystals [39]. To test if the defect scaling laws are also affected by curvature variations and surface genus, we performed additional parameter scans for tori. The locally varying curvature on a torus determines the strength of the symmetry-breaking term $\delta_{u}\Gamma(u)$ in Eq. (1), implying that a purely hexagonal toroidal surface crystals exists only for sufficiently small aspect ratios $r/R$ [39]. We therefore focus on thin tori with $r/R=0.2$ (Fig. 3A,B; Movie 2). Adopting the same small random $\epsilon$-inhomogeneities as for the sphere simulations (and, hence, the same shift value $\Sigma_{0}$), we observe that the freeze-out stress $\Delta\Sigma_{f}\propto\mu^{1/2}$ and the freeze-out time $t_{f}\propto\mu^{-1/2}$ scale exactly as in the case of constant curvature (Fig. 3C,D). Furthermore, one again finds that faster quenches lead to a higher density of defects at freeze-out (Fig. 3A,B,F). Although Euler’s theorem [37] imposes a vanishing total defect charge $S=0$ on a torus, excess defects of opposite charge tend to aggregate at the inner ($Z=7$) and outer ($Z=5$) equators to lower the elastic energy of the toroidal crystal [39, 53, 54, 38], resulting in a non-zero equilibrium defect density $\rho_{\infty}$. As in the sphere case, we find that the relative excess defect density $\Delta\rho_{f}=\rho_{f}-\rho_{\infty}$ of the toroidal surface crystals follows a square root law $\Delta\rho_{f}\propto\mu^{1/2}$ (Fig. 3F), suggesting that local curvature variations and surface genus do not significantly affect the scaling laws. In the next part, we will rationalize these observations by considering the structure of the amplitude equations [55] for the underlying GSH theory.

III Discussion

KZ-type scaling for nonthermal quenches. As outlined in the introduction, the original KZ scaling relations for continuous second-order transitions can be obtained by analyzing how correlation length and relaxation time diverge as a function of the control parameter as one approaches the critical point from the high-symmetry phase. By contrast, for the nonthermal first-order wrinkling transitions considered here, the correlation length is not defined in the unwrinkled phase, so that the standard KZ arguments do not apply. One can nevertheless explain the numerically determined scaling laws by inspecting the amplitude equations for the GSH model. To this end, we assume approximate hexagonal solutions of the form $u=U(t)\sum_{a=1}^{3}[e^{i\mathbf{k}_{a}\cdot\mathbf{x}}+e^{-i\mathbf{k}_{a}\cdot\mathbf{x}}]$, where $\mathbf{k}_{1}=k_{c}(1,0)$, $\mathbf{k}_{2}=k_{c}(-1/2,\sqrt{3}/2)$, $\mathbf{k}_{3}=k_{c}(-1/2,-\sqrt{3}/2)$ and $k_{c}=\sqrt{6|\gamma_{0}|}$. Inserting this ansatz into the GSH equation (1), one finds to leading order in $U$ and curvature (SI Amplitude Equations)

[TABLE]

Dividing by $\sqrt{\mu}$ and defining a rescaled time $t^{\prime}=\mu^{1/2}t$, we can remove the quench rate dependence at leading order, to obtain

[TABLE]

Non-autonomous equations of this type have been extensively studied in dynamical systems theory and are known to describe delayed bifurcations with some characteristic delay time scale $t^{\prime}_{f}$ [56, 48, 57, 49]. Identifying this characteristic delay with the freeze-out time implies $t_{f}\propto\mu^{-1/2}$, in agreement with our numerical results. It may be worth emphasizing again that this power-law scaling is a direct consequence of the nonthermal stress-quench dynamics. This example corroborates that KZ-type scaling can occur even when the conventional critical slowing down of thermal systems is replaced by other delayed bifurcation mechanisms [58]. Moreover, since the quench enters through a purely local $u$-term in Eq. (1), the scaling law is not affected by curvature variations. Finally, in order to explain the scaling law for the excess defect density $\Delta\rho_{f}$ at freeze-out, it suffices to assume that the mean-squared distance $\delta^{2}$ between defect pairs grows diffusively, consistent with quasi-random migration of defect precursors, so that at freeze-out $\delta^{2}\propto t_{f}\propto\mu^{-1/2}$. For a 2D surface, this then implies that $\rho_{f}\propto\delta^{-2}\propto\mu^{1/2}$, in agreement with the numerical results in Fig. 2G and 3F.

Universal nucleation dynamics in curved surface crystals. Crystal growth in planar Euclidean geometry is well understood. By contrast, the complex interplay between kinetics, substrate curvature and defect formation is still being investigated [25, 26]. The current interest in these topics is in parts driven by recent technological advances in the fabrication of graphene [59] and carbon nanotubes [35] and by the development of modern confocal imaging techniques that make it possible to track micron-sized colloids and cells at high spatio-temporal resolution. For instance, the formation of hexagonal crystal structures similar to those described above can be observed during the early developmental stages in the fruit fly Drosophila melanogaster, when nuclei migrate to accumulate underneath the surface of the ellipsoidal shell that encapsulates the embryo [60]. Similarly, ciliated somatic cells form a spherical crystal on the surface of the colonial alga Volvox carteri [61], with the cells’ arrangement determining the phototactic properties of the organism.

Important insights into the kinetics of crystal growth on curved substrates were obtained recently in a joint experimental and theoretical study [25] on the assembly of charged colloids on spherical liquid-liquid interfaces. The crystal formation dynamics observed in these experiments shares striking similarities with the nucleation sequences of the hexagonal wrinkling patterns shown in Fig. 1C. In both cases, one first observes the formation of several smaller highly regular crystal patches, while the defects form during the later stages when two or more of these patches merge. These kinetic parallels suggest a certain universality in the crystal formation processes on curved surfaces in the slow-quench regime. Extrapolating these similarities to higher quench rates, one may hope that the scaling results identified here translate to other physical and biological systems that develop crystalline structures on their curved surfaces.

IV Conclusions

The above analysis shows that the KZ mechanism for continuous second-order transitions has a direct analog in first-order phase transitions, if the underlying amplitude equations exhibit delayed bifurcations. As a specific example, we have identified the power law scaling relation between topological defects densities and linear quench rates for an experimentally validated generalized Swift-Hohenberg theory [46] describing stress-induced discontinuous surface wrinkling transitions in thin stiff films adhered to a curved soft substrates. With regard to applications, these scaling relations offer concrete guidance for controlling the number of topological defects by tuning the quench rate, extending previous work [38, 39] that showed how substrate geometry can be used for defect localization. The nucleation sequences leading to crystalline hexagonal surface patterns in the thin-film model are qualitatively similar to those reported recently for colloidal crystals on spherical liquid-liquid interfaces [25]. This suggests that elastic surface crystals, which have been realized by substrate depressurization [43] or surface swelling [44], can offer a flexible testbed for exploring generic aspects of crystal growth kinetics and topological defect formation in complex geometries.

Methods

Pattern analysis. Simulations were performed using the algorithm detailed in Ref. [46], adopting parameters for centimetre-sized polydimethylsiloxane-coated elastomer hemispheres [45, 46]. To reconstruct the curved crystal structure and detect topological defects, we first threshold the amplitude field $u$ obtained from simulations to find the center of each crystalline lattice site. Each site is then connected to its nearest neighbors via a Delaunay triangulation, and the hexagonal crystal structure as well as defects are obtained from the dual Voronoi graph. To find the flat crystal representation of Fig. 1C,E, we first construct the Voronoi cells of the crystal. We then randomly select a regular crystal site $s_{0}$ and its six neighbors as center of the planar lattice. Based on their positions relative to $s_{0}$, the neighbors can be assigned to one of the six surrounding unit lattice positions. We then repeat this construction sequentially for each neighbor by aligning their closest neighbors with the existing planar lattice structure, following the procedure used in Ref. [25]. To obtain the time evolution, we first construct the final planar crystal from to the fully crystallized configuration at freeze-out time $t_{f}$. For all earlier times $t<t_{f}$, we use a Hungarian matching algorithm [62] to track identical Voronoi sites and identify them in the fully crystallized planar configuration.

Determination of freeze-out stress. We first determine the stress $\Sigma_{\textrm{max}}$ and corresponding time $t_{\textrm{max}}=\Sigma_{\textrm{max}}/\mu$ where $\langle u^{2}\rangle$ has largest slope. Although $\Sigma_{\textrm{max}}-\Sigma_{0}\propto\mu^{1/2}$ (Fig. S1A,B), crystalline patterns are not yet fully developed at this value of the film stress, rendering an analysis of the defect density scaling practically unfeasible. To obtain robust estimates for the defect densities, we therefore add a time delay to $t_{\textrm{max}}$ to allow the hexagonal patterns to develop completely. According to Eq. (5), the dynamics is rate-independent under the rescaling $t^{\prime}=\mu^{1/2}t$ to first order (Fig. S1C). Therefore, in order to give each quench the same relative amount of time for developing hexagons, we choose a time delay proportional $\mu^{-1/2}$, $t_{f}=t_{\textrm{max}}+\mu^{-1/2}\Delta\tau$, with $\Delta\tau=30$ kept constant throughout our analysis (SI Freeze-out). Accordingly, the freeze-out film stress is defined as $\Sigma_{e,f}=\mu t_{f}=\Sigma_{\textrm{max}}+\mu^{1/2}\Delta\tau$.