Impact of Interaction Range and Curvature on Crystal Growth of Particles Confined to Spherical Surfaces

TL;DR

This study uses computer simulations to explore how interaction range and curvature influence crystal growth on spherical surfaces, revealing conditions for ribbon formation and defect incorporation.

Contribution

It demonstrates the formation of ribbons at short interaction ranges and defect-packed structures at longer ranges, highlighting a two-step nucleation process.

Findings

Ribbons form above a critical crystal size.

Short-range interactions favor ribbon growth.

Longer-range interactions promote defect incorporation.

Abstract

When colloidal particles form a crystal phase on a spherical template, their packing is governed by the effective interaction between them and the elastic strain of bending the growing crystal. For example, if growth commences under appropriate conditions, and the circular crystal that forms reaches a critical size, growth continues by incorporation of defects to alleviate elastic strain. Recently it was found experimentally that, if defect formation is somehow not possible, the crystal instead continues growing in ribbons that protrude from the original crystal. Here we report on computer simulations in which we observe both the formation of ribbons at short interaction ranges and packings that incorporate defects if the interaction is longer-ranged. The ribbons only form above some critical crystal size, below which the nucleus is roughly spherically shaped. We find that the scaling of the critical crystal size differs slightly from the one proposed by the Manoharan group, and reason this is because the actual process is a two-step heterogeneous nucleation of ribbons on top of roughly circular crystals.

Figures

Equations

$$U(r)=\left[U_M(r)-U_M(r_c)\right]H(r_c-r),$$

$$U_M(r)=ϵ\left[e^{2\alpha (r-r_0)}-2e^{\alpha (r-r_0)}\right],$$

$$U_0=$$

$$r_{ij}:=$$

$$U_M(r):=$$

$$e^{-\beta (F_1-F_0-C)}=\frac{{⟨f(\beta (U_1-U_0-C))⟩}_0}{{⟨f(\beta (U_0-U_1+C))⟩}_1}$$

$$\Delta G_c=\rho \Delta \mu \pi a^2/4+\gamma \pi a.$$

$$\Delta G_r=\rho \Delta \mu wl+2\gamma (w+l).$$

$$\Delta g_c=\frac{1}{4\pi }\left[\frac{\rho \Delta \mu \pi }{4Y}{\eta }^2+\frac{\pi \gamma }{RY}\eta +\frac{\pi }{24576}{\eta }^6\right].$$

$$\Delta g_r=\frac{1}{4\pi }\left[\frac{\rho \Delta \mu }{Y}\omega \lambda +\frac{2\gamma }{RY}(\omega +\lambda )+\frac{{\omega }^5\lambda }{640(1-{\nu }^2)}\right],$$

$$\Delta g_r=\frac{1}{4\pi }\left[\left(\frac{\rho \Delta \mu }{Y}+\frac{\gamma }{RY}\frac{2}{\omega }+\frac{{\omega }^4}{640(1-{\nu }^2)}\right)\zeta +\frac{\gamma }{RY}2\omega \right].$$

$$4\pi \frac{\partial \Delta g_r}{\partial \omega }=\left[\frac{\rho \Delta \mu }{Y}-\frac{\gamma }{RY}\frac{2}{{\omega }^2}+\frac{{\omega }^3}{160(1-{\nu }^2)}\right]\zeta +2\frac{\gamma }{RY}=0$$

$$\left[\left(\frac{\rho \Delta \mu }{Y}+\frac{\gamma }{RY}\frac{2}{\omega }+\frac{{\omega }^4}{640(1-{\nu }^2)}\right)\right]\zeta <-2\omega \frac{\gamma }{RY}.$$

$$\frac{\rho \Delta \mu }{Y}+\frac{\gamma }{RY}\frac{2}{\omega }+\frac{{\omega }^4}{640(1-{\nu }^2)}+2\omega \frac{\gamma }{RY}\le 0.$$

$$\Delta g_{c+r}=\frac{1}{4\pi }\left[\frac{\rho \Delta \mu }{Y}\left(A_c+A_r\right)+\frac{\gamma }{RY}\left(\sqrt{4\pi A_c}+2\frac{A_r}{\omega }\right)+\frac{A_c^3}{384{\pi }^2}+\frac{{\omega }^4{A}_r}{640(1-{\nu }^2)}\right],$$

$$\mathrm{d}\Delta g_{c+r}=\left(\frac{\partial \Delta g_{c+r}}{\partial A_c}\right)\mathrm{d}A=\frac{\mathrm{d}A}{4\pi R^2}\left[\frac{\rho \Delta \mu }{Y}+\frac{\gamma }{RY}\sqrt{\frac{\pi }{A_c}}+\frac{A_c^2}{128{\pi }^2}\right].$$

$$\mathrm{d}\Delta g_{c+r}=\left(\frac{\partial \Delta g_{c+r}}{\partial A_r}\right)\mathrm{d}A=\frac{\mathrm{d}A}{4\pi R^2}\left[\frac{\rho \Delta \mu }{Y}+\frac{2\gamma }{RY\omega }+\frac{{\omega }^4}{640(1-{\nu }^2)}\right].$$

$$\Delta \Delta g:=$$

$$=$$

Full text

1 Introduction

Colloidosomes are droplets whose surfaces are densely packed with colloidal particles. [1] Confinement of the colloidal particles to the liquid-liquid interface minimises the contact area between the two liquids. At high surface coverage, the colloids either pack in a disordered, glassy or in an ordered, crystalline fashion. [1, 2, 3, 4] In the former particles are kinetically trapped. In the latter the equilibrium packing of the colloids is determined by the interplay between the curvature of the droplet and the exact nature of the interaction between the colloids. This can give rise to various kinds of defect and defect organisations, producing grain boundary scars [5, 6, 7], pleats [8] and/or growth in ribbon-like shapes that emanate from circular domains without defects at the boundary between the two. [9]

Grain boundary scars and pleats involve point defects in specific arrangements, where a point defect is a particle that has fewer or more than six nearest neighbours. Ribbons form when particles, for whatever reason, refuse to give up local hexagonal order. One such reason can be elastic stress, resulting from a short interaction range. [9, 10] In macroscopic theory, this stress manifests itself in an additional elastic term in the free energy. [11, 12] The experimentally observed ribbon-like crystals of Ref. [9] show an initial, almost isotropic growth that, after a critical size is reached, transitions into the growth of ribbons protruding from the initial crystal. Similar findings are reported in Ref. [10], in which phase field crystal model calculations confirm that growing in ribbons indeed can relieve the elastic stress. Note that the crystal growth is continuous along the interface between the seed and ribbon, i.e., there is no grain boundary between the two because the ribbon preserves the same lattice vectors.

Theoretical arguments for how the largest circular domain should scale with the Young’s modulus and template radius were given by Refs. [9, 10] and are based on finding the optimal size of a circular crystal bent onto a spherical surface, which is, in this topology, bounded by an elastic energy penalty. The theory predicts that ribbons form for a sufficiently short range of interaction, for example of a suitably parametrised Morse potential as employed e.g. Ref. [13], and only then for sufficiently large domain sizes.

We find by means of computer simulation that the range of the interaction potential between the particles indeed dictates whether or not ribbons form. For sufficiently large interaction ranges, we find that ribbon formation is suppressed in favour of incorporating defects into the crystal. For shorter ranges of interaction, the critical size for ribbon formation reduces, which is consistent with Refs [9, 10]. However, the scaling exponent we observe is slightly different, and instead coincides with the scaling of the optimal width of a pure ribbon on a spherical surface proposed in Ref. [14].

The remainder of this article is structured as follows. In Section 2 we discuss the simulation and analysis methods. In Section 3, we discuss our computational findings and show that we indeed observe a critical domain size and extract its scaling relations with respect to the potential range and the radius of curvature. Finally, in Section 4, we underline the most important implications of our findings.

2 Methods

We perform Langevin dynamics (LD) simulations of $N$ particles constrained to a spherical surface with radius $R,$ using the LAMMPS program and a specialised RATTLE algorithm. [15, 16, 17]. The particle number $N$ we couple to the template radius as $Nr_{0}^{2}/16R^{2}=0.4,$ so that if we associate an area $\pi r_{0}^{2}/4$ with each particle of diameter $r_{0},$ we have the same area coverage $\phi=0.4$ for all template radii. Initially the particles are in non-overlapping, random positions on the spherical surface.

We apply a Langevin thermostat with an arbitrary damping time $\tau_{L}$ to the particles to keep the system at a constant temperature $T$ and to make the particles undergo Brownian motion. $\tau_{L}$ is the time it takes for the velocity auto-correlation function of a particle to decay to $1/\mathrm{e}$ of its initial value, and is our reference time unit that in effect measures the ratio of the particle’s mass and the friction constant. For the interaction potential between particle pairs we use a truncated and shifted Morse potential $U(r)$

[TABLE]

where $r$ is the three-dimensional Cartesian inter-particle distance, $r_{c}$ the cut-off distance, $H(r_{c}-r)$ the Heaviside step function that is $1$ if $r_{c}-r>0$ and [math] otherwise and $U_{M}(r)$ is the original Morse potential

[TABLE]

with a well depth $\epsilon,$ equilibrium spacing $r_{0}$ and shape parameter $\alpha.$ Throughout the remainder of this article we use $r_{0}$ as reference length unit and $\epsilon$ as the reference energy unit.

To encourage the growth of a single crystal, we lightly tether the particles to the top of the spherical surface $(x,y,z)=(0,0,R)$ with a harmonic spring with spring constant $\kappa=2.5\epsilon/r_{0}^{2}.$ We initialise the structure with this spring in place for $N_{s}$ time steps at a temperature of $k_{B}T=0.3\epsilon,$ after which we remove the springs. We then equilibrate for another $N_{s}$ steps while linearly ramping down the temperature from $k_{B}T=0.3\epsilon$ to $0.25\epsilon.$ We verified that a longer ramping time had no influence on the formed structures. After this annealing phase we sample for $N_{s}$ time steps at $k_{B}T=0.25\epsilon.$ We analyse snapshots that are $5~{}\tau_{l}$ apart and average over them to get statistics. For $\alpha=40/r_{0}$ a smaller time step of $0.001\tau_{L}$ was required and we set $N_{s}=25\times 10^{6}$ to guarantee proper equilibration and stable dynamics, while for the other values of $\alpha$ we used $N_{s}=250\times 10^{3}$ and a time step size $0.005\tau_{L}.$

The parameter $\alpha$ sets how sharp the potential is peaked around its minimum at $r_{0}$, as the effective spring constant $k=\partial^{2}U(r)/\partial r^{2}|_{r_{0}}=2\alpha^{2}\epsilon.$ See Fig. 1 for a graph comparing the Lennard-Jones and Morse potentials for various values of $\alpha.$ If we ignore the impact of multi-body interactions, the Young’s modulus scales linearly with $k.$ This implies that a larger $\alpha$ leads to a larger Young’s modulus [18]. Therefore, in our computer simulations we can use the parameter $\alpha$ to influence the Young’s modulus of our crystal. By studying the particle vibrations about their lattice site, however, we see that the pair approximation is not very accurate. Hence, we rely on the actual vibrations to estimate the effective spring constant. See Section S1 of the supporting material for a detailed description of the analysis.

For our crystals it turns out that effective spring constants $\kappa$ scale significantly less strongly with $\alpha$ than predicted by the harmonic approximation $\kappa\sim\alpha^{2}.$ We obtain scalings between $\kappa\sim\alpha^{1.3}$ to $\kappa\sim\alpha^{1.6},$ with an average exponent of $(1.52\pm 0.08).$ The deviation from the pair potential approximation could be either due to collective effects that are obviously not included in any pair-wise approximation, or due to the non-zero temperature in our simulations, which does influence elastic properties. [19] This is not a complete surprise as hard sphere crystals have an effective spring constant that follows entirely from many-body effects. [20] This shows that it is critical to measure Young’s modulus or at least the effective spring constant of the material at hand rather than relying on the harmonic approximation [18, 10] when performing Langevin dynamics simulations. With the effective spring constant determined, we still need methods to determine the scaling of the chemical potential and the largest circular domain size with $\alpha.$

For larger values of $\alpha$ at which the crystal forms ribbons, most of the direct interaction is governed by nearest-neighbour interactions, as can be inferred from Fig. 1. Because of this, we assume the line tension arises predominantly from the edge particles having fewer than six bonds. In this case, we expect that the line tension and the chemical potential scale similarly with $\alpha,$ and we only have to determine $\Delta\mu.$

To determine the chemical potential, we have to determine the free energy difference between a crystal consisting of $N$ Morse particles and that same crystal in which one of the bulk particles is changed to an ideal gas particle. This is problematic because an ideal gas particle tends to explore all of the spherical surface, while the Morse particle tends to stay near its lattice site. Therefore, we instead determine the free energy difference associated with transforming one Morse particle into an “Einstein particle”, that is, a particle that is tethered to its lattice site by a harmonic spring.

As reference Morse particle we take the particle that has the largest minimum distance to the edge particles. By applying Bennett’s acceptance ratio [21] to the case of $N$ Morse particles and the case of where we transform this reference particle into an Einstein particle, we obtain free energy differences that we can directly relate to the excess chemical potential. In Supporting Material S2 a more detailed description of the procedure is presented. The found scaling is $\Delta\mu\sim\alpha^{(-0.53\pm 0.01)}.$ Note that the free energy difference goes through an appreciable range as a function of $\alpha,$ as we find $\Delta\mu\approx-0.7/\epsilon$ for the largest $\alpha r_{0}=40,$ which is a third of that of $\Delta\mu\approx-2.2/\epsilon$ for the smallest $\alpha r_{0}=4.$ It is thus important to take this scaling into account when assessing the scaling for the largest circular domain size at which the growth transitions to ribbons.

With the scalings for the effective spring constant $\kappa$ and hence Young’s modulus $Y$ and that of the chemical potential $\Delta\mu$ and line tension $\gamma$ determined, we can now determine how the largest circular domain sizes scale as a function of the aforementioned parameters, provided that we can extract what the largest circular domain size is. We describe our approach now in brief.

Finding the largest circular domain size involves identifying the edge of what constitutes a circular domain and extracting the distance from that edge to the centre of the domain, which can be reasonably approximated by finding the particle that is the furthest away from all edge particles. We calculate for all particles the shortest distance to the edge. The particle that is the farthest away from the edge we then consider to be the centre of a circular domain, and its distance to the edge, divided by the template radius $R,$ we take as value for the diameter of the largest circular domain. Note that this particle coincides to that we determine the chemical potential for.

Incidentally, this value also gives a good indication of the transition from incorporated point defects to ribbons, as this number scales differently with the range parameter $\alpha$ when the crystal incorporates point defects, as we shall see in Section 3. This analysis we perform with post-processing scripts that are described in more detail in Supporting Material S3. With all the methods discussed, we now present our results.

3 Results

In Fig. 2 we present typical crystal structures we observe as a function of the range parameter $\alpha.$ The figures show the top half of the spherical template that is covered with particles. For $\alpha r_{0}=4$ the crystal incorporates various defects, predominantly dislocations (Fig. 2(a)). At around $\alpha r_{0}=8$ a transition from incorporated defects to hexagonal packings occurs. From this point on the defects are only located at the edge or around holes (Fig. 2(b)). Increasing $\alpha$ further to $\alpha r_{0}=16$ (Fig. 2(c)) and $\alpha r_{0}=32$ (Fig. 2(d)) leads to a clear formation of ribbon-like structures separated by larger tears. Furthermore the width of the ribbons reduces, which is qualitatively consistent with both the scalings presented in Refs. [9, 10] and Ref. [14].

Note finally that at larger $\alpha$ the effective temperature of the particles appears to have increased slightly, in the sense that more particles appear in a gas-like phase on the side of the template that is not covered by the crystal. The motion of the particles inside the crystal, however, becomes smaller due to the increased sharpness of the potential well depth.

Our crystals are reminiscent of the ones observed in simulations by Cong et al. in Ref. [23]. For $\alpha r_{0}\geq 8$ we observe hexagonal patches separated by tears and holes, i.e., ribbons. Below that $\alpha$ the crystal incorporates pleats and point defects, which are reminiscent of those observed experimentally in Refs. [5, 6, 7, 8]. Hence, as shown experimentally in Ref. [9] and consistent with the calculations of Ref. [10], ribbons only form when the elastic strain, regulated by the range of the potential, prevents the formation of point defects. Interestingly, we do not seem to observe the four-fold branching described in Ref. [10], possibly due to the thermal fluctuations in our simulations, which are absent in a phase-field crystal model.

In Fig. 3 we show an illustrative example of the difference between ribbons and incorporated point defects, obtained by applying the protocol of Section 2 to $N=450$ particles on a template of radius $R=15r_{0}$ for $\alpha r_{0}=4,~{}6$ and $12.$ This example is obtained at an area coverage to $12.5\%$ to illustrate more clearly the different structures. For $\alpha r_{0}=4$ the crystal is roughly spherical and incorporates some defects. For $\alpha r_{0}=6$ we see the onset of a ribbon that is branching out of the initial bulk while there are still defects in the original nucleus as well. Hence, there appears to be a sort of “coexistence” between incorporated defects and ribbon formation. For larger $\alpha r_{0}=8$ we see clear branching, and the defects in the bulk are significantly reduced.

We now analyse the size of the largest circular domain we observe. All simulation results from which the largest circular domain size is extracted below correspond to a higher area coverage of $40\%$ rather than the $12.5\%$ that served as an illustration of the morphology in Fig. 3. First we consider its scaling with the spherical template radius. Figures 4(a) and 4(b) reveal that the critical domain size diameter $a$ scales sublinearly with $R.$ Fitting a power law reveals the scaling exponent to be $(0.83\pm 0.06),$ which is inconsistent with the scaling proposed in [9, 10].

It does match the scaling of the optimal ribbon width proposed in Ref. [14], as well as our own theory based on heterogeneous nucleation of a ribbon-like structure on a circular nucleus, presented in the Supporting Material S4. Note however that the two theories are very similar, predicting either $a\sim R(\rho\Delta\mu/Y)^{1/4}$ or $a\sim R(\gamma/RY)^{1/5}.$ Hence, while our result quantitatively agrees more closely with the two-stage nucleation model, simulations on larger templates are required to give a definitive answer.

The figures also show that the largest circular domain size decreases with increasing $\alpha.$ This is qualitatively consistent with either theory. In Figs. 5 we plot the observed largest domain diameter divided by the scaling argument presented in Refs. [9, 10] (Fig. 3(a)) and the heterogeneous nucleation model of Supporting Material S4 (Fig. 3(b)).

For $\alpha r_{0}\leq 6,$ where the crystal incorporates point defects, the domain size appears to be independent of $\alpha r_{0}.$ This makes sense as the largest circular domain will only depend on the number of particles in the crystal, which, for our simulation setup, scales linearly with $R.$ Since the scaling of the largest domain size with $\alpha$ is clearly not relevant in this regime, we omit the data for $\alpha r_{0}<8$ in the following analysis of the largest circular domain size.

Note that the data presented in Fig. 5(b) collapses onto a single curve, giving a strong indication that in the parameter range we have access to, two-stage nucleation is the main pathway for ribbon formation. For $\alpha r_{0}\geq 8,$ fitting a power law gives us access to the second scaling exponent of $(-0.332\pm 0.005),$ which, in combination with the scaling exponent with $R$ of $(0.83\pm 0.06),$ quantifies our observed scaling of $a\sim R^{(0.83)}\alpha^{-0.33}.$ The heterogeneous nucleation model predicts, with our empirical scaling of $Y$ and $\gamma$ with $\alpha,$ that the second scaling exponent should be $(-0.41\pm 0.02),$ whereas the original scaling proposed in Refs. [9, 10] leads to an exponent of $(-0.51\pm 0.02).$ Again, the heterogeneous nucleation model seems to better fit our observations, although the scaling with the interaction range is certainly less convincing than the scaling with the spherical template radius $R.$

4 Discussion and conclusion

We performed Langevin dynamics simulations of Morse particles of diameter $r_{0}$ on spherical surfaces of varying radius. We can tune the range of interaction between the particles with a parameter $\alpha,$ where larger values of $\alpha$ represent shorter ranges of attraction. We only consider a single interaction strength $\epsilon=4k_{B}T,$ at which we obtained well-equilibrated crystals. Below a critical $\alpha r_{0}\leq 6$ the crystals that form are roughly spherical and incorporate dislocations and point defects. Above this critical value the crystals instead exhibit a smaller, defect-free, circle-like nucleus with protruding, ribbon-like structures. This provides a confirmation of the suggestion that the formation of the ribbons is indeed driven by the elastic instability put forward in Ref. [9].

We quantified the scaling of the size of these circular domains as a function of both the template radius $R$ and the range parameter $\alpha$ and find that the largest circular domain diameter scales as $a\sim R^{(0.83\pm 0.06)}\alpha^{(-0.33\pm 0.01)}.$ This is close to but not consistent with the scaling predicted in Refs. [9, 10], which assume the critical circular domain size follows from the global minimum in the free energy due to the elastic penalty. Rather, we find that the scaling is more consistent with a heterogeneous nucleation model, where the transition is not dictated by the global minimum in the free energy but rather by the fact that for a sufficiently large circular domain, continued growth as a ribbon is energetically more favourable than as a circle. However, we varied both parameters through at most one decade, so further studies on larger crystals are required to determine the full parameter range over which the model is valid.

Finally, more generally, our findings indicate that, in addition to its well-known influence on the number of local minima and liquid phase stability in free space, [24, 25, 26, 27] the range of attraction also influences the morphology of the formed crystal on a curved surface. This has important consequences for, e.g., the formation of virus capsids, whose building blocks interact through short-ranged interactions [28] but are typically modelled by a long-ranged Lennard-Jones potential [29].

5 Acknowledgements

We are grateful to Mike Hagan, Greg Grason, Jayson Paulose, Vinothan Manoharan and Wouter Ellenbroek for valuable discussions, and Thijs van der Heijden for critical proof-reading. We acknowledge the HFSP for funding under grant RGP0017/2012.

Supporting Material for Impact of Interaction Range and Curvature on Crystal Growth of Particles Confined to Spherical Surfaces

S1 Extracting effective spring constants

In order to estimate the scaling of Young’s modulus, we require the the effective spring constant $\kappa$ that keeps the particles in place in the crystal. For a single particle pair this can be found analytically by harmonically approximating the pair potential around the minimum. For the case of the Morse potential, this leads to $\kappa=2\alpha^{2}\epsilon.$ However, this simple approximation does not take into account collective effects. Therefore, we extract effective spring constants from our simulation data.

From our simulation trajectories, we determine a running average of the particle positions $\mathbf{x}$ by applying an exponentially weighted average to obtain for a time step $t_{n}$ the average $\mathbf{x}_{a}(t_{n})=(1-\xi)\mathbf{x}(t_{n})+\xi\mathbf{x}_{a}(t_{n-1}).$ The parameter $\xi$ controls how fast the past positions are “forgotten”. At each time step, we determine the squared deviation of the particles from the average from the previous frame, $\delta\mathbf{x}^{2}(t_{n}):=(\mathbf{x}(t_{n})-\mathbf{x}_{a}(t_{n-1}))^{2}.$ This serves as an approximation of the squared deviation of the particle from its average lattice site while taking into account that the lattice site might drift in time. Typical squared displacement values we find are of the order of $0.05r_{0}^{2}.$ We empirically tuned $\xi$ to a value of 0.9 at which the fast motions about the average lattice site are suppressed.

This produces for each particle a time trace of the quantity $\delta\mathbf{x}^{2},$ which can be averaged in time to obtain $\left<\delta\mathbf{x}^{2}\right>.$ If we assume particles are bound to their lattice site by a harmonic spring, the effective spring constant then follows from equipartition as $\frac{1}{2}\kappa\left<\delta\mathbf{x}^{2}\right>=k_{B}T.$ In Fig. S1 we present distributions for the squared deviation and effective spring constant for four values of $\alpha.$

In Fig. S2 we present the effective spring constants for all parameter values. Note that for very large $\alpha$ the scaling breaks down. This is because the script used to determine the peak in the histogram simply determines the numerical maximum as an approximation for the mean. This breaks down for noisy and particularly wide distributions, which we obtain for large $\alpha,$ as can be seen from the grey and red distributions in Fig. S1(b).

S2 Determining the chemical potential

To determine the chemical potential of a particle in the crystal, we determine the free energy difference between the original crystal of $N$ Morse particles, and a hybrid Morse-Einstein crystal in which the particle which is most representative of the bulk of the crystal is replaced by an Einstein particle tethered to its initial lattice site with a harmonic spring. To determine the most representative particle, we apply the analysis to extract the largest circular domain size, which we shall discuss later in Appendix S3. This analysis finds the particle that is the furthest away from all edge particles, which is the best available representative of the bulk.

We put the spring constant to $k=50\epsilon/r_{0}^{2},$ at which the Einstein particle and the original Morse particle have comparable root mean square deviations from the average lattice site, calculated following the procedure detailed in Appendix S1. We apply Bennett’s acceptance criterion [21] to the two different cases, whose total potential energies are

[TABLE]

Here, $U_{0}$ is the original truncated, shifted Morse potential $U$ with cut-off distance $r_{c},$ well depth $\epsilon$ and range parameter $\alpha$ acting on all $N$ particles with diameter $r_{0},$ while $U_{1}$ is the original Morse potential acting on $N-1$ particles combined with an “Einstein potential” $U_{E}$ applied to the $N$th particle. $H(x)$ is the Heaviside function with $H(x)=1$ for $x>0$ and [math] otherwise. The Einstein potential is, of course, a harmonic spring with spring constant $\kappa$ that attaches the $N$th particle to its initial position $\mathbf{x}_{N}(0).$ Finally, $r_{ij}$ is the scalar distance between particles $i$ and $j,$ located at positions $\mathbf{x}_{i}$ and $\mathbf{x}_{j}.$

We apply Bennett’s acceptance criterion [21] with energy offset $C=2.4\epsilon$ and sample

[TABLE]

where $\beta$ is the reciprocal thermal energy $1/k_{B}T,$ $f(x)$ is the Fermi function $f(x)=1/(1+\exp(x)),$ and the subscripts [math] and $1$ mean that the average is obtained by sampling the potential energy function $U_{0}$ and $U_{1},$ respectively. The value of $C$ corresponds roughly to the average potential energy difference between the two potential energy functions, $C=\left<U_{1}\right>_{1}-\left<U_{0}\right>_{0},$ which leads to a reasonably good statistics and leads to i. We determined this value for $\alpha=16/r_{0}$ but applied it for all $\alpha.$ For varying $\alpha$ the averages $\left<U_{1}\right>$ and $\left<U_{0}\right>$ do differ significantly, but the method remains usable for this one fixed value of $C.$

Applying Bennett’s acceptance ratio means that we generate trajectories corresponding to both $U_{0}$ and $U_{1},$ and average $f(\beta(U_{1}-U_{0}-C))$ over the trajectory generated by $U_{0},$ while we average $f(\beta(U_{0}-U_{1}+C))$ over the trajectory generated by $U_{1}.$ This can then be converted straightforwardly into an estimate for the free energy difference between the Morse and the Morse-Einstein crystal. The excess chemical potential of the Einstein particle $\mu_{E}^{ex}$ can be calculated analytically from the partition function and is $\mu_{E}^{ex}=k_{B}T\ln(\kappa A/2\pi Nk_{B}T),$ where $\kappa$ is the spring constant, $A$ is the spherical template area, $N$ the number of particles on the template and $k_{B}T$ the thermal energy. Hence, the free energy difference can directly be converted into an excess chemical potential by adding $\mu_{E}^{ex},$ which we present in Fig. S3. The free energy difference appears to scale as $\mu_{E}^{ex}\sim\alpha^{-0.52},$ independent of the spherical template radius.

S3 Extracting the largest domain size

To analyse the ribbon widths a series of steps is required. Each of them is explained in detail here. First, we construct a network of nearest neighbours. Particles are considered nearest neighbours if their distance is less than $1.3r_{0}$ apart. This distance coincides with the minimum after the first peak in the pair correlation function. From this information we can identify “edge particles”, i.e., particles that do not have six nearest neighbours.

With this information, we can compute for all particles that are not edge particles the distance to the nearest edge particles. The largest of these distances we take as an estimate for the largest circular domain size, as it represents the largest possible circle diameter that fits inside the crystal. This assumption breaks down for $\alpha$ sufficiently small to allow for incorporation of defects inside the lattice, rather than at the edge. This is observed for $\alpha\leq 6$ in Fig. 3(a). In Fig. S4 we present snapshots of the different stages of analysis. With the aforementioned steps, we obtain a value for the largest circular domain diameter for each combination of spherical template radius $R$ and potential range parameter $\alpha.$

S4 Classical nucleation theory on curved surfaces

S4.1 Flat plane

Macroscopically, the formation of a two-dimensional crystal on a flat surface is described by Classical Nucleation Theory (CNT). [30] Although this theory is well-known, we introduce it here because we later extend it with elastic terms to describe crystal formation on curved surfaces. Furthermore, it will serve as the basis for our two-stage model for ribbon growth, which we introduce in Section S4.3. In CNT the thermodynamic driving force towards crystal formation is opposed by a surface tension that in two dimensions is a line tension. If we assume the line tension to be invariant to the locally exposed crystal plane, the free energy of a crystal nucleus is given by $\Delta G=N\Delta\mu+\gamma L,$ where $N$ is the number of particles in the crystal, $\Delta\mu<0$ the chemical potential difference between the crystal phase and the surrounding liquid or gas, $L$ the crystal circumference and $\gamma$ the line tension [30, 31].

If we assume a circular crystal, we have $L=\pi a,$ with $a$ the diameter of the crystal so that $\Delta G_{c}=N\Delta\mu+\pi a\gamma.$ Note that actual two-dimensional crystals tend to form hexagonal nuclei instead because the underlying interaction potential encourages hexagonal bond order that induces anisotropic growth. [32, 30, 31] However, both a hexagon and a circle exhibit the same scaling of area and circumference with the number of particles and for the purpose of extracting scaling laws it is a good approximation. Limitless growth of the drop occurs if $\Delta G<0,$ i.e., if $N>-\pi a\gamma/\Delta\mu.$ Of course, $a=2\sqrt{N/\pi\rho}$ depends on $N$ through the equilibrium particle surface density $\rho$ of the crystal. In other words, limitless growth occurs for $N>4\pi\gamma^{2}/(\Delta\mu)^{2}\rho.$ Converting the expression for the free energy to the circle diameter $a$ leads to

[TABLE]

Ribbon-like crystals are better described by a rectangle with a circumference $A=2(w+l)$ and a particle number $N=\rho wl,$ leading to

[TABLE]

Unlike the circle, the ribbon has two parameters $w$ and $l,$ and the optimal free energy has $w=l,$ as this gives the largest area for a given circumference. In the flat plane, a circle-like shape will always be preferred over a ribbon-like shape. This can be seen by comparing Eqs. (7) and (8) under the constraint that they have equal areas, $wl=\pi a^{2}/4.$ Their difference is then simply the difference of the line tension terms, $\gamma(\pi a-2(w+\pi a^{2}/4w)),$ which, if we substitute the optimal $w=\sqrt{\pi}a/2,$ is always negative. Therefore, not surprisingly, in a flat plane, a ribbon is always destabilised with respect to a circle. However, as becomes clear later, this trivial conclusion no longer holds on a curved surface.

To see whether or not a ribbon is more likely to be nucleated than a circle, we consider the kinetic barrier heights. For a circle, the maximum in $\Delta G_{c}$ corresponds to a critical diameter $a=-2\gamma/\rho\Delta\mu$ and is equal to $-\pi\gamma^{2}/\rho\Delta\mu,$ where we recall that $\Delta\mu<0.$ For the ribbon the maximum is located at $l=w=-2\gamma/\rho\Delta\mu$ and has a value of $-4\gamma^{2}/\rho\Delta\mu.$ Therefore, in the flat plane, a ribbon-like crystal is also less likely to nucleate, as evidenced by a higher nucleation barrier in the free energy. This means that in a two-dimensional plane a circular crystal will always be favoured over a ribbon-shaped one, both thermodynamically and kinetically.

S4.2 On a sphere

The above picture changes if we consider crystals constrained to a spherical surface of radius $R.$ In this case, an additional energy enters the free energy that takes into account the elastic cost of bending the crystal to accommodate the curved template. [11, 12] First consider again the circular crystal. To make the analysis more straightforward, we scale the free energy to the spherical surface area, $4\pi R^{2},$ times Young’s modulus $Y.$ This leads to a reduced unit $\eta:=a/R$ and a dimensionless free energy [11, 12, 9, 10]

[TABLE]

Although the elastic term does not influence the kinetic barrier height or critical diameter significantly [33], it does influence the thermodynamic stability of the circular crystal. For very large $Y,$ the elastic term dominates to such an extent that the local minimum in the free energy that occurs for sizes larger than the critical nucleus size becomes larger than $0,$ as is pointed out in Refs. [18, 33, 10]. In that case, $\Delta g_{c}\geq 0,$ and the formation of a circular crystal is thermodynamically suppressed.

Furthermore, for $Y>0,$ there is now an optimal circular diameter that minimises the free energy. Under the assumption $\gamma\ll RY$ this minimum is located at $\eta=(-6144\rho\Delta\mu/Y)^{1/4}:=\eta_{0}.$ In Refs. [9, 10] this value is associated with the critical circular crystal size at which growth transitions to a ribbon, as circular domains larger than this radius have a higher free energy. This scaling we shall refer to as the optimal circle scaling. However, as long as $\Delta g_{c}<0,$ the circle can in principle continue to grow, for example if there is not a sufficient amount of material to nucleate a second, stable circular cluster.

Unlike in the flat plane, on a sphere a ribbon-shaped crystal can be more stable than a circle due to a different scaling of the elastic energy with the crystal dimensions. The dimensionless free energy for the ribbon on a spherical surface is given in terms of dimensionless width $\omega=w/R$ and length $\lambda=l/R$ by [11, 12]

[TABLE]

where $\nu$ is the Poisson ratio of the material. Note that in this expression, $\rho\Delta\mu\lambda+2\gamma/R$ has to be negative, and hence it only holds for $\rho\Delta\mu<-2\gamma/R\lambda=-2\gamma/L.$

The reason a ribbon can be stabilised over a circular crystal on a curved surface is because although the elastic strain scales with $\omega^{5},$ it only scales linearly with $\lambda.$ This leads to a different scaling with the total crystal area, which becomes obvious when the free energy of the ribbon is expressed in terms of the width $\omega$ and dimensionless area $A/R^{2}:=\zeta=\omega\lambda,$

[TABLE]

This free energy can in principle be optimised with respect to $\omega$ to obtain an optimal ribbon width for a given area:

[TABLE]

Because the lowest order term $-2\gamma/RY\zeta\omega^{2}$ is negative, there will always be some $\omega$ that optimises Eq. (11).

For the ribbon to be stable at some point, the total free energy should be negative, i.e., $\Delta g_{r}<0.$ Since $2\omega\gamma/RY>0,$ the ribbon can only be stable when

[TABLE]

Since $0\leq\zeta\leq 1,$ this condition can only be satisfied if at least

[TABLE]

The equality can only be achieved for sufficiently small $\rho\Delta\mu/Y$ and $gRY,$ i.e., at sufficiently large Young’s moduli. When this inequality holds, a ribbon can in principle be indefinitely large and still have a negative free energy, unlike a circular crystal.

However, the kinetic barrier associated with a ribbon is always larger than that of a circle. Under the assumption that the elastic contribution is negligible for crystals smaller than the critical nucleus size, we find that the barrier height, i.e., the maximum in $\Delta g$, is located at $\eta_{b}=-2\gamma/\rho\Delta\mu R$ for the circle and at $\omega_{b}=\eta_{b},\zeta_{b}=\omega^{2}_{b}$ for the ribbon. The corresponding free energies are $\Delta g_{c}^{b}=-\gamma^{2}/4\rho\Delta\mu R^{2}$ and $\Delta g_{r}^{b}=-\gamma^{2}/\pi\rho\Delta\mu R^{2}.$ Hence, the kinetic barrier for formation of a circle is lower than that of a ribbon by a factor of $\pi/4.$

This implies that, while under certain conditions a ribbon-shaped crystal has a lower free energy, it is less likely to be nucleated. A more likely scenario that leads to the formation of ribbon-like structures is therefore a two-stage nucleation, in which the ribbon grows out of a circular nucleus. In such an event, the crystal has the kinetic energy barrier of a circlular shape and initially grows as such. However, at some point, the free energy decrease for continued growth as a circle is smaller than that for growth as a ribbon, and hence, the crystal continues growth into ribbon-like structures. Such a growth pathway is also more consistent with the experimental results [9] and crystal phase field calculations. [10] In Section S4.3 we present a free energy for such a growth type.

S4.3 Two-stage nucleation theory on a sphere

In Section S4 we argue that the free energy of a circular crystal is, under certain conditions, lower than that of a ribbon. Hence, the formation of a pure ribbon-like crystal is not probable, as ribbons smaller than the transition size are destabilised with respect to circular crystals with the same area. This motivates us to investigate the possibility of a two-stage nucleation, meaning that a ribbon-shaped structure grows out of a pre-existing, circular nucleus. For this to happen, the free energy gain due to a small area increment in the shape of a ribbon should be smaller than that due to an equal area increment in the shape of the already existing circle. Both growth modes are sketched in Fig. S5.

From the typical crystal morphologies presented in Refs. [9] and [10], it appears that the ribbons grow out of the initial crystals without distorting the lattice vectors. Therefore, we assume there is no line tension associated with the circle-ribbon interface. We furthermore assume that the line tension associated with the width of the ribbon is of the same as the previously exposed circular rim that is now covered by the ribbon. This means there is no line tension term associated with the width, because it is already accounted for by the original line tension term. With these assumptions, the total free energy of the hybrid crystal becomes

[TABLE]

where $A_{r}$ is the dimensionless ribbon area $\omega\lambda R^{2}$ and $A_{c}$ the circle area $\pi(R\eta)^{2}/4.$ The free energy change of the crystal $\mathrm{d}{\Delta g_{r+c}}$ of growing purely as a circle can be determined by expanding $\Delta g_{c+r}$ in a Taylor series around $A_{c}$ for $A_{r}=0,$ leading to

[TABLE]

To determine the free energy change of growing as a ribbon, we expand $\Delta g_{c+r}$ in terms of $A_{r}$ for fixed $\omega$ and $A_{c},$ leading to

[TABLE]

The width that gives the optimal free energy change follows from optimising $\mathrm{d}{\Delta g_{c+r}}$ with respect to $\omega,$ which leads to $\omega=[320(1-\nu^{2})\gamma/RY]^{1/5}.$

In other words, the free energy change due to continued growth as a circle from an area $A_{c}$ with increment $\Delta A$ is given by $(\partial\Delta g_{c+r}/\partial A_{c})\Delta A$ and the free energy change due to growth as a ribbon out of a circular nucleus is given by $(\partial\Delta g_{c+r}/\partial A_{r})\Delta A.$ Therefore, the difference between Eqs. (13) and (14), $\Delta\Delta g,$ indicates whether continued growth favours a ribbon or a circle,

[TABLE]

Continued growth as a ribbon is preferred under two conditions. First of all, it is necessary that $\Delta\Delta g>0,$ as only then the free energy gain for growing as a ribbon is preferred over growing as a circle. Secondly, it is required that $\Delta g_{c+r}<0,$ as else the formation of any type of crystal at all is destabilised. Finally, although not technically a necessity for the initial onset of ribbon formation, we check whether or not $\partial\Delta g/\partial A_{r}<0.$ If this is the case, the crystal can grow indefinitely as a ribbon, since continued growth as a ribbon will lead to a further decrease in the free energy.

To determine the critical area at which the transition takes place, we first numerically determine for a range of $\gamma/RY$ and $\rho\Delta\mu$ the roots of $\Delta\Delta g,$ i.e., the points that satisfy the first criterion. Then, we determine whether or not $\Delta g_{c+r}<0$ and what the sign of $\partial\Delta g_{c+r}/\partial A_{r}$ is. We find that $\Delta g_{c+r}<0$ is only satistfied at the roots of Eq. (15) for significantly large ratios for $|\rho\Delta\mu R/\gamma|.$ In particular, we find that at a ratio of at least $6$ is required to obtain a transition at an area $A_{c}<0.6$ that satisfies all criteria, corresponding to $\gamma/RY=2\times 10^{-4}$ and $\rho\Delta\mu/Y=12\times 10^{-4}.$ Furthermore, although increasing Young’s modulus leads to a decrease in the critial area at which $\Delta\Delta g=0,$ it also leads to a higher free energy, occasionally such that $\Delta g_{c+r}>0$ at the transition. Hence, although larger Young’s moduli favour ribbon formation, they will only actually form at a proportional decrease in the chemical potential.

To determine how the elastic cost influences the area at which the transition occurs, we determine how the critical circular area scales with Young’s modulus for various ratios of the chemical potential and line tension terms. In figure S6 we plot our findings, which clearly show that the critical area scales as $A_{c}\sim(\gamma/RY)^{2/5}$ and hence the diameter of the largest circular domain size scales as $\eta\sim(\gamma/RY)^{1/5}.$ Note that all points collapse on the same curve given by $A_{c}=23.6594(\gamma/RY)^{0.4}.$ The relative standard errors in the linear regression are $1.3\times 10^{-4}\%$ for the prefactor and $3.5\times 10^{-5}\%$ for the exponent. Both the prefactor and scaling exponent are independent of $\rho\Delta\mu/Y.$

We notice a clear breakdown in this scaling for $\gamma/RY>2\times 10^{-4}.$ In this regime, the line tension is prohibitively large and suppresses ribbon formation. There is another breakdown for sufficiently small $\gamma/RY$ in combination with small $\rho\Delta\mu/Y.$ In this regime, the chemical potential is not negative enough to compensate for the combined costs of bending and the line tension, in other words, there is no negative free energy $\Delta g_{c+r}<0$ anymore.

Hence, our analysis suggests that, if the formation of ribbons indeed proceeds as a two-stage, heterogeneous nucleation of a ribbon-like structure on a pre-existing, circular nucleus, the normalised area of this circular nucleus should scale as $A_{c}=\pi a^{2}/16\pi R^{2}\sim(\gamma/RY)^{2/5},$ where $a$ is the diameter of the circular crystal and $R$ the radius of the template. This is different from the scalings previously proposed in Refs. [9, 10], in which the transition area is associated with the minimum in $\Delta g_{c},$ which, in the limit of negligible line tension, scales as $A_{c}\sim(\rho\Delta\mu/Y)^{1/4}.$

The actual scaling of the transition area can be obtained from computer simulations by determining the largest circular domain size as a function of the spherical template radius, the chemical potential, the line tension and Young’s modulus, as is described in the main text and in the previous sections.