Crystalline Order and Topological Charges on Capillary Bridges †

We numerically investigate crystalline order on negative Gaussian curvature capillary bridges. In agreement with the experimental results in [W. we observe for decreasing integrated Gaussian curvature, a sequence of transitions, from no defects to isolated dislocations, pleats, scars and isolated sevenfold disclinations. We especially focus on the dependency of topological charge on the integrated Gaussian curvature, for which we observe, again in agreement with the experimental results, no net disclination for an integrated curvature down to À10, and an approximately linear behavior from there on until the disclinations match the integrated curvature of À12. In contrast to previous studies in which ground states for each geometry are searched for, we here show that the experimental results, which are likely to be in a metastable state, can be best resembled by mimicking the experimental settings and continuously changing the geometry. The obtained configurations are only low energy local minima. The results are computed using a phase field crystal approach on catenoid-like surfaces and are highly sensitive to the initialization. When so materials assemble in the presence of geometric and topological constraints, the regular order favored by local interactions is frustrated, which leads to defect structures in the ground state. The capsid of a spherical virus provides an example. Here proteins form subunits, which arrange in an icosahedral shape, constructed from varying numbers of hex-amers and 12 pentamers. On more general surfaces with spatially varying positive and negative curvature the geometric and topological constraints lead to more complex defect structures. These defects can be chemically functionalized 1,2 and provide the key to self-assembly into complex hierarchical structures with emergent novel macroscopic properties. Geometric and topological constraints thus provide a route to control the self-assembly process and an understanding of the interplay of geometry, topology and defects is therefore of utmost importance, if this route to design novel materials will be explored. See e.g. the themed collection on The geometry and topology of so materials 3 for a broader discussion. Besides this practical motivation, crystalline order on curved surfaces is also a rich academic problem. Various theoretical results are known. Topological constraints are considered in a classic theorem of Euler, which shows for a triangulation of the surface that X i ð6 À iÞv i ¼ 6x, with v i as the number of vertices with i nearest neighbors and x as the Euler characteristic of the surface. Thus for surfaces with the …

When so materials assemble in the presence of geometric and topological constraints, the regular order favored by local interactions is frustrated, which leads to defect structures in the ground state. The capsid of a spherical virus provides an example. Here proteins form subunits, which arrange in an icosahedral shape, constructed from varying numbers of hexamers and 12 pentamers. On more general surfaces with spatially varying positive and negative curvature the geometric and topological constraints lead to more complex defect structures. These defects can be chemically functionalized 1,2 and provide the key to self-assembly into complex hierarchical structures with emergent novel macroscopic properties. Geometric and topological constraints thus provide a route to control the self-assembly process and an understanding of the interplay of geometry, topology and defects is therefore of utmost importance, if this route to design novel materials will be explored. See e.g. the themed collection on The geometry and topology of so materials 3 for a broader discussion.
Besides this practical motivation, crystalline order on curved surfaces is also a rich academic problem. Various theoretical results are known. Topological constraints are considered in a classic theorem of Euler, which shows for a triangulation of the surface that X i ð6 À iÞv i ¼ 6x, with v i as the number of vertices with i nearest neighbors and x as the Euler characteristic of the surface. Thus for surfaces with the topology of a sphere (x ¼ 2) the total charge must be 12. One realization is again the spherical virus, considered as a triangular lattice with sixfold coordination, with 12 vefold disclinations present. However, there are many other possibilities to fulll the theorem of Euler and it is energetics that determines the number, type and arrangement of defects. With each disclination a stretching energy is associated (relative to a perfect triangular lattice in at space) which grows proportional to R 2 , with R as the radius of the sphere. For a xed lattice constant a we have N $ (R/a) 2 , where N is the number of particles. Thus for large N, mechanisms are expected which reduce this extra energy by changing the ground-state conguration. This can be done by surface buckling as considered e.g. in ref. 4 and 5. However, in cases where buckling is limited, the energy is reduced by introducing dislocations and grainboundary scars, see ref. 6. Realizations of water droplets in oil, which are coated with colloidal particles 7 show that these excess dislocations grow linearly with the system size.
The same arguments hold for toroidal crystals and capillary bridges with the topology of a cylinder (x ¼ 0). In these cases, there is no topological need for disclinations and all defects are introduced to relieve the strain induced by the curvature. For toroidal geometries it is found that excess disclinations are energetically favorable for fat torii 8 and that the number of excess disclinations is controlled primarily by the aspect ratio of the torus. 9 Experiments for colloidal particles on oil-glycerol capillary bridges 10 show a more complex behavior. No net topological charge is found on the surface for an integrated Gaussian curvature down to U ¼ À10, with U ¼ 3=p ð G GdG and Gaussian curvature G. For U beyond the threshold of À10, disclinations rapidly ll the surface until 12 disclinations approximately match the integrated curvature of À12. A sequence of transitions is observed, from no defects to isolated dislocations, which organize into pleats, and nally form scars and isolated sevenfold disclinations, see Fig. 1 for an explanation of the various defect types.
If the presence of topological grain-boundary scars for spheres and torii is surprising and unexpected, the occurrence of pleats on capillary bridges is even more surprising. It has been explained in ref. 11 in terms of a curvature screening provided by the stress free boundary.
While the linear increase of excess dislocations and the formation of scars on a sphere and a torus is well understood and reproduced quantitatively by various theoretical approaches and simulation techniques, see e.g., ref. 9 and 12-16 and the review, ref. 17, the behavior on capillary bridges is less understood. In ref. 18 the problem is addressed by mapping the microscopic interacting particle problem to the problem of discrete interacting defects in a continuum elastic background. The critical value for the transition between congurations with and without disclinations or scars is deduced for catenoids. Three different ways are considered. In the rst approach the critical value corresponds to the catenoid for which the integrated Gaussian curvature of the geodesic disc matches the curvature, which would be screened by a sevenfold disclination. Using the denition in ref. 10, in which the integrated Gaussian curvature is dened in units of disclinations as above, the obtained critical waist size c* ¼ 0.85 corresponds to U* ¼ À6.3. The second approach uses an effective screened disclination charge, for which c* ¼ 0.87, corresponding to U* ¼ À6.0. The third approach uses an energetic argument. The transition point for the emergence of a disclination in the interior is found to be c* ¼ 0.8, corresponding to U* ¼ À7.2. All values are signicantly larger than the experimentally obtained U* ¼ À10. The same approach is used in ref. 11 and validated against a numerical energy minimization of a discrete spring model. The results are in qualitative agreement with ref. 10 but no quantitative estimate for U* is given. Simulation results for electrostatically charged particles are obtained for a wide range of system sizes, curvature and interaction potentials, including Yukawa, Coulomb and Lennard-Jones in ref. 19 and 20. The results qualitatively agree with the continuum theory and the experimental data. Depending on U dislocations, pleats, scars and sevenfold disclinations are found, independent of the used potential. The supplemental material in ref. 19 also provides tabulated numbers of positive and negative topological charges and a critical waist size c*, which here correspond to the transition to isolated sevenfold disclinations. It seems to be only weakly dependent on the number of particles, but sensitive to the considered potential. For the Coulomb potential, a critical waist size c* ¼ 0.55 is found, corresponding to U* ¼ À10.01. However, the tabulated charges indicate the presence of charged defects as scars already above that value, in disagreement with the experimental results in ref. 10. The data also strongly depend on the considered constant mean curvature (Delaunay) surface. Within the experiments, the capillary bridge is obtained through a sequence of Delaunay surfaces (unduloids, catenoids and nodoids) each with different geometrical parameters, such as aspect ratio, mean curvature and maximal Gaussian curvature. Such different Delaunay surfaces are considered in ref. 20. However, quantitative results for U* are not given.
Here, instead of a coarse grained elasticity problem for the defects or a discrete particle simulation for the position of N particles, we consider a continuous optimization problem for a number density, with N maxima representing the N particles. The idea is to nd a free energy, which is minimized for a conguration, which corresponds to the minimum of the original problem. We consider the Swi-Hohenberg free energy on the surface G with V G and D G the corresponding surface gradient and surface Laplacian. We dene Thereby r denotes a number density of the particles and f ðrÞ ¼ 1 2 ð1 À 3Þ r 2 þ 1 4 r 4 can have a double-well structure, depending on the parameter 3. Within various parameter regimes, the minimizers of this energy are periodic solutions with a hexagonal structure in a at space, which are geometrically frustrated in the considered setting on G. We here have an interplay between the tendency to form periodic solutions and the geometric frustration which prevents such a periodicity. We consider the H À1 gradient ow of this energy, which can be viewed as an extension of the phase-eld-crystal (PFC) model, we obtain the following system of surface partial differential equations The system is solved towards the steady state and each maximum in the computed number density is identied with a  p i À p j with i-th particle position p i and the distance measured in R 3 , with known results for N [ 12,2790]. However, the approach is also quantitatively related to discrete particle simulations. In ref. 23 the equations are derived from a discrete particle setting via classical dynamic density functional theory, following the derivation in ref. 24 for the at space. Numerical results on a sphere have shown that the obtained minimal energy congurations are insensitive to the specic underlying interparticle potential. These results are in agreement with ref. 19 where ve different potentials are considered and qualitatively the same defect motifs are found.
Both results indicate that the geometric frustration is much stronger than the inuence of the interparticle potential. We therefore stick to the phenomenological model above. Besides this qualitative and quantitative comparison of the phase eld crystal approach with discrete particle simulation results, questions concerning reachable size and efficiency can be asked. The largest problem considered in ref. 23 accounts for N ¼ 99 602 maxima on a sphere and nds in energy E 1 which is only 0.0004% above the lower bound 0.5N 2 À 0.553051N 3/2 from ref. 25. The simulation took about 1.5 days on a single processor with 12 GB RAM. The computational approach certainly is more involved than discrete particle simulations, but allows for simple coupling with other elds, such as surface buckling of the underlying surface, as considered for viral capsided in ref. 26 or uid ow in the surrounding medium, as e.g. considered for bijels 27 and particle stabilized emulsions. 28,29 Different numerical approaches have been considered in ref. 23 to solve the system of surface partial differential equations. We here adapt a parametric nite element setting, which is realized within the simulation toolbox AMDiS. 30 The approach is thereby based on the stable nite element discretization for the PFC model in ref. 31 and is described in detail in ref. 23. The key idea is to use the surface operators on the discrete surface which consists of triangles T. To do the integration on these triangles a parameterization F T :T / T is used, withT the standard element in R 2 . These allow transformation of all integrations onto the standard element using the nite element basis dened also in R 2 . The parameterization F T is given by the coordinates of the surface mesh elements and provides the only difference between solving equations on surfaces and on planar domains. For a surface we have to allow F T : R 2 / R 3 , whereas for a planar domain F T : R 2 / R 2 . With this tiny modication any code to solve partial differential equations on Cartesian grids can be used to solve the same problem on a surface, providing a surface triangulation is given. It is essentially this modication which induces the geometric frustration to the problem.
Our surface is given as an approximation of a catenoid and obtained through a rotation around the x 3 -axis. The parametric equation of the surface is given by with time t and angle f. Only for r 0 ¼ c the approximation is equal to a catenoid, see this approximation, we can scale the surface area by varying the waist radius r 0 , without changing the height s, the vertical curvature radius at the waist c and the integrated Gaussian curvature U. The approach allows for changing the surface to scatter (U˛(À12, 0)) by keeping the surface area A and the outer radius R xed. The surface is discretized using triangular elements with a sufficient mesh size h $ a/10. We use the parameters 3 ¼ 0.4 and the average particle density as r 0 ¼ À0.3, ‡ which correspond to a point in the two-dimensional phase diagram within the hexagonal region.
The boundary conditions are crucial, as they have to be stress free. This is accounted for by specifying Dirichlet conditions using a one-mode approximation, see ref. 32. Choosing the outer radius R in accordance with the distance between neighboring particles allows for a stress free boundary.
We consider four different initializations. We either specify different initial conditions and use (a) random initialization with r ¼ r 0 + h, with white noise h, (b) initialization at the boundary, with the one-mode approximation specied at the boundary, (c) initialization at the waist, with the one-mode approximation specied at the waist or as in ref. 10 start with a cylinder and subsequently change the surface to decrease U. Within this last approach we use a sequence of 42 geometries and compute the steady state on each, with the one from the c. ‡ The considered particle density r is more precisely a non-dimensional density difference with a reference state of an averaged density in the liquid state, see ref. 21. previous surface as the initial condition. An animation is provided in the ESI. † The simulation results are postprocessed in the following way. We identify the maxima of the computed particle density r, interpret them as particle positions and dene their neighborship based on their two-dimensional Voronoi regions. These data are used to evaluate the number of dislocations, scars, pleats, veand sevenfold disclinations, as well as the number of three-and vefold disclinations on the boundary. Within colliding defect structures, clusters are separated preferring dislocations and pleats over disclinations and scars, larger defects over smaller defects and oriented defects over unoriented defects. The visualization is done using the soware Ovito. 33 Fig. 3 shows selected results for the evolving surface. For all initializations, also for the not shown random initialization and the initialization from the boundary or the waist, we observe the whole spectrum of defects: dislocations, pleats, scars, veand sevenfold disclinations and the expected sequence of transitions to dislocations, pleats, scars and isolated sevenfold disclinations. Table 1 summarizes their rst appearance.
The appearance and interactions of defects differ for different initializations. While for random initialization defects are already present for the largest value of U, scars are mainly accompanied by vefold disclinations next to the boundary and sevenfold disclinations appear at the waist, the congurations resulting from initialization at the boundary are characterized by migration of dislocations into the surface and the formation of a second row of dislocations, as well as an increase in size of the occurring pleats. For the initialization from the waist, we observe only oriented dislocation for large U, with the number increasing with decreasing U until the rst defects occur also at the boundary. The appearance of pleats comes together with non-oriented dislocations. As the pleats grow, the dislocations become again oriented. For the evolving geometry, the surface remains defect free up to relatively low values for U, rst defects occur as oriented dislocations at the boundary. Their number stays xed until pleats are formed aer further decreasing U. Here pleats do not grow but fall apart, forming an oriented dislocation in the interior and one at the boundary. Scars and sevenfold disclinations are formed at the same time. In contrast with the size dependent exclusive appearance of disclinations or scars reported in ref. 18, we found both defect types at the same time.
We now concentrate on the rst appearance of scars, vefold or sevenfold disclinations and thus a detached topological charge. For random initializations, scars already appear on surfaces with U ¼ À6. 25. For the initialization at the boundary  and the waist the conguration stays without detached topological charge down to U ¼ À11.14 and U ¼ À11.35, respectively, all in quantitative disagreement with the experimental results in ref. 10. Only the evolving geometry, which resembles the evolution in the experimental setting best, leads to the expected behavior, with no detached topological charge down to À10.44.
To further highlight the quantitative agreement, we compute the detached topological charge q int , as the difference of the number of particles in the interior with seven and ve neighbors, the total topological charge q total , which in addition considers particles at the boundary, which do not have four neighbors and thus should be zero for U˛(À12, 0), and the number of defects N defect , which we scale by the number of particles N to obtain the ratio z ¼ N defect /N. This becomes necessary as N is not xed throughout evolution and can therefore not set to a distinct value. Fig. 4 shows these characteristic quantities as functions of U.
The detached topological charge for the evolving surface simulations shows the expected behavior, with q int ¼ 0 down to U z À10 and an approximately linear increase up to q int z 12 for U z À12. This is in agreement with the experimental results in ref. 10. In agreement with the topological requirement the total charge q total remains zero up to the smallest considered value of U ¼ À11.86 for all initializations, which is just a consistency check of our postprocessing and the number of defects growth for decreasing U. For the evolving surface simulations the congurations stay defect free down to U ¼ À7.75, where defects form suddenly. The number of defects stays almost constant aerwards until a second jump can be seen for U ¼ À10.44. Below that value more and more defects are introduced. This has a clear inuence on the energy, for which we compute again the Coulomb potential p i À p j and a power law potential suggested in p i À p j 3 with i-th particle position p i and the distance measured in R 3 . To eliminate the dependence on N we rescale these energies to obtainẼ 1 ¼ E 1 /N 2 andẼ 3 ¼ E 3 / N 2 . Fig. 5 shows the computed energies as functions of U.
While the different initializations lead to almost identical results for long-range potentialẼ 1 , large differences can be seen for the short range potentialẼ 3 . The curve for the evolving surface resembles nicely the characteristic defects. The energy is increasing down to U ¼ À7.75 with a defect-free conguration. The sudden drop in the energy results from the occurrence of defects, in the form of dislocations at the boundary. The energy again increases until U ¼ À10.44. The small drop here corresponds again to a change in the defect type, which leads to the occurrence of detached topological charge. The energy plots further indicate that the congurations obtained with the evolving surface approach not always lead to the lowest energy. We can assume a highly complex energy landscape, in which all our observed congurations are presumably trapped in local minima. Instead of searching for global minima, as in ref. 19, we concentrate on resampling the experiments in ref. 10, with presumable also only local minima congurations. With the evolving surface approach, which resamples the experimental setting best, quantitative agreement for all considered data q int and z can be achieved. The number of defects is thereby lower than in various putative global minima congurations reported in ref. 19.
The goal to chemically functionalize the defects to control self-assembly into supramolecular structures, requires not only the presence of a certain number and type of defects, but also their position and arrangement to be predictable. The congurations obtained with the evolving surface approach here lead to highly symmetric arrangements, see Fig. 3. The number of defects growth until 12 equally spaced oriented dislocations on the boundary are formed. These defects remain and grow into pleats, without changing their position. This highly symmetric arrangement even remains for smaller U aer the splitting into scars and sevenfold disclinations and only for the lowest values of U the symmetry is lost. To identify this symmetric sequence of transitions as the most favorable path requires computation of the energy barriers Fig. 4 (Left) Detached topological charge q int , (middle) total topological charge q total , and (right) scaled number of defects z functions of U, with green random initialization, black initialization from the boundary, blue initialization from the waist, and purple the evolving surface computation. The left figure shows in addition the experimentally observed linear increase from 0 to 12 for U ¼ À10 to U ¼ À12 in ref. 10. and minimal energy paths, which is currently under investigation using the string method 34 already applied to the phase eld crystal model in ref. 35 and 36. The future goal should be to identify similar symmetric defect arrangements at nanoscales and to further explore the complex interplay of topology, geometry and surface evolution for novel materials design. The considered phase eld crystal approach, even if computationally more expensive than discrete particle simulations, provides a exible tool to consider such an evolution and further coupling with other external elds.