Topological defects of dipole patchy particles on a spherical surface

Abstract

We investigate the assembly of dipole-like patchy particles confined to a spherical surface by Brownian dynamics simulations. The surface property of the spherical particle is described by the spherical harmonic Y₁₀, and the orientation of the particle is defined as the uniaxial axis. On a flat space, we observe a defect-free square lattice with nematic order. On a spherical surface, defects appear due to the topological constraint. As for the director field, four defects of winding number +1/2 are observed, satisfying the Euler characteristic. We have found many configurations of the four defects lying near a great circle. Regarding the positional order for the square lattice, eight grain boundary scars proliferate linearly with the sphere size. The positions and orientations of the eight grain boundary scars are strongly related to the four +1/2 defect cores.

---

1 Introduction

Patchy particles are particles of colloidal size and have patches acting as specific interactive sites on the particles. Due to the anisotropic interaction, patchy particles are capable of assembling into complex structures whose properties are fundamentally different from “conventional” materials.¹ Recent developments in synthesis techniques^2–4 have made feasible the fabrication of patchy particles of high degrees of freedom^2,5–9 and research has focused on the self-assembly of patchy particles. Practical applications of the particle assembly require knowledge of the principal design of the particle for a specific structure and vice versa, and how to produce high yield target structures. Understanding the principles of controlling defects is important for these tasks. Defects are imperfections and singularities in an ordered structure. From a microscopic point of view, defects can serve as interacting sites for chemical linkers to promote the formation of large scale structures.^10,11 On the macroscopic length-scale, defects are inevitable during crystallisation and affect the overall properties. In particular, using curvature is one method to study defects because the defects are topologically protected, meaning that they cannot disappear upon continuous deformation of the order parameter. We focus on the spherical surface due to its substantial applicability, such as in crystalline membranes, the design of crystalline materials, and the fabrication of patchy colloids.

The in-plane order of a two-dimensional crystal on a curved surface is a much more complicated problem than that on a flat space because the geometrical and topological constraints have to be taken into account. The ground state of isotropic particles on a two-dimensional plane is a defect-free triangular lattice containing 6-fold coordinated particles. On a spherical surface, such a translational symmetry is broken, and at least twelve 5-fold coordinated particles are required to compensate the topology of the sphere. This is similar to the truncated icosahedron pattern of a soccer ball whose twelve pentagons are icosahedrally arranged among the hexagons.^12,13 Aside from affecting the positional order of the particle system, the curved surface also influences the orientational order. For example, rod-like particles confined to a spherical surface form four +1/2 or two +1 defects^10,14–17 instead of a defect-free nematic phase (long-range orientational but no long-range positional order on a flat space).

Patchy particles have both positional and orientational order, which has been studied only separately. For each separate case, although the appearance of the topological defects is generally accepted as a consequence of the embedded geometry and topology, many issues are still recent subjects of study: the underlying mechanism for the formation of defects^12,18 and the dependence on the system size,¹³ the influence of the particle–particle and particle–curvature interaction on the defect. When the interplay between positional and orientational order may appear, our understanding of defects of these two orders remains primitive. In order to address the above problems, we consider the most simple form of anisotropic particle exhibiting both positional and orientational order, that is, a spherical particle with dipole-like patches that behaves somewhat similarly to a magnetic bead. We dynamically simulate and compare the assembly of such patchy particles confined to a planar geometry and a spherical surface.

1.1 Topological defects on a curved surface

It is useful to discuss the basics and the relevant studies on the topological defects of ordered structures embedded in curved surfaces. Suppose a closed surface is facetted and divided into a number of V vertices, E edges, and F faces, Euler theorem states that V − E + F = χ, where the Euler characteristic χ = 2(1 − g), and g is the genus of the closed surface. For instance, V − E + F = 2 is applied to all polyhedra because they are topologically equal to a sphere with χ = 2.¹⁹ If we restrict every face to having c vertices, and let N_z be the number of vertices that have z connections with the others, then the Euler theorem can be written as (see Appendix A.1)

Consider point particles on a surface triangularly facetted (c = 3) by the points, it is well known that on a flat space, the isotropic particles, most of the time, form a triangular lattice resulting in a 6-fold symmetry with z = 6. A disclination in this case refers to a vertex whose z deviates from six and the charge of the disclination is defined as 6 − z.²⁰ When such a particle system is confined to a spherical surface, it is straightforward from eqn (1) that a net charge of is required. For sphere size below a critical value, the twelve 5-fold coordinated particles are icosahedrally located among the otherwise 6-fold ones due to the repulsion of the like-sign charge.²⁰ As the sphere size increases, the 5-fold disclinations are screened by additional dislocations, which are pairs of 5- and 7-fold coordinated particles;²¹ note that the pairs may consist of 4- or 8-fold coordinated particles as long as the sum of charge is neutral. The interaction between the clouds of dislocations eventually leads to the formation of twelve grain boundary scars where each scar consists of pairs of 5- and 7-fold coordinated particles with a net charge of +1. The grain boundary scars in spherical crystals are observed in both experiments and simulations.^12,13,22,23 Due to the formation of the defect scars, the interaction between the disclinations is believed to be screened and less important. The arrangement of these twelve defect scars is unusual, complex and not fully understood, especially in the large sphere limit.²¹ On the ellipsoidal or toroidal surface, the positions of the grain boundary scars are dependent on the local Gaussian curvature of the surface.²⁴ For a square lattice on the sphere, evidence of how disclinations distribute is still lacking. However, it can be conjectured that by employing a quadrilateral mesh (c = 4) for a square lattice, a disclination is now the vertex whose connection differs from four. Eqn (1) becomes , meaning that at least eight disclinations of z = 3 are required for a square lattice on a sphere. The development of dislocations that are now viewed as bound disclination pairs of 3- and 5-fold coordinated particles²⁵ is expected to occur as the system size increases.

The in-plane orientation on a two-dimensional surface is described by a unit vector field tangent to the surface. In general, p-atic vectors are invariant under rotations of 2nπ/p (n is integer) about the surface normal.¹⁴ The strength of a defect in this case is characterised by the winding number k, defined as the amount of rotation of the vector field around a counter-clockwise circuit enclosing the defect core on the order parameter space, k = Δθ/(2π/p), where Δθ is the angle of the vector that rotates in one counter-clockwise circuit.²⁶ On a closed surface, the net strength is equal to the Euler characteristic according to the Poincaré–Hopf theorem.²⁷ For instance, the in-plane order structure of a vector field (p = 1) on a sphere requires at least two +1 defects at the two opposite poles; nematics (p = 2) have either four +1/2 or two +1 defects. The detailed configuration of spherical p-atic order depends on the interaction energy of the system. According to the Frank free energy approach for spherical nematics, in the one elastic constant limit, i.e., splay constant equals bending constant, the ground state exposes four +1/2 defects at the vertices of a tetrahedron inscribed in the sphere.^10,14,28 In contrast, in the extreme limit when splay is much softer (or harder) than bending, some Monte Carlo simulation studies suggest that a +1 defect can split into two +1/2 defects without costing energy. As a result, an infinite number of states of four +1/2 defects lying near a great circle can be obtained by the cut-and-rotate surgery of the sphere with two +1 defects.^16,17 Such behaviours are also confirmed via the Onsager approach minimising the mean-field free energy of hard rods on a sphere.^29–31 It has also been found that the defects of nematic order on a curved surface accumulate at a larger Gaussian curvature.^32,33 Another relevant case is spherical tetratics (p = 4), where the particles are square-shaped or cross-shaped. The low energy states consist of eight +1/4 defects that may position on the vertices of an anti-cube^14,34 or a cube.^35,36 Such a difference is perhaps caused by the different types of interaction for simulation and the forms of free energy for the continuum description. It is worth noting that the p-atic particles in these mentioned systems have a purely rotational degree of freedom, hence the existence of positional order remains unclear.

2 Methods

2.1 Brownian dynamics

We employ the Brownian dynamics simulation algorithm for particles in Euclidean space in the overdamped limit.^37,38 The model is successfully applied for particles confined to a flat space. The translational and rotational motions of the particles confined to a spherical surface are given as followswhere r(t + Δt) and Ω(t + Δt) denote the position and orientation of the particle after the time step Δt; D^T and D^R are the translational and rotational diffusion coefficients of an isolated particle, respectively; δr and δΩ are the translation and rotation due to thermal fluctuation, satisfying and , where each component of δ^G_i is independently chosen from a Gaussian distribution with zero mean and unit variance. Force F and torque T are derived from the pairwise potential. In order to capture the dynamics of particle position on the tangent plane, we apply the algorithm in ref. 39, in which the tangential parts of the interacting force F and noise δr in eqn (2a) at the point where a particle is located are considered; finally harmonic terms F^H = κ(r − R)r/r are added to enforce the confinement after the translation of each time step.

The patchy particle possesses a patterned surface related to its physical or chemical properties, which induces the anisotropic interaction of the particles. Such a pattern can be systematically described by means of spherical harmonics Y_lm. In this study, the dipole-like pattern of a unit particle of radius a = 1 is given as , where p̂ is defined as the orientation of the particle (Fig. 1). This pattern has positive and negative hemispheres similar to the Janus particle.^8,40,41 The interaction potential for a pair of particles i and j comprises an isotropic Week–Chandler–Anderson potential V_WCA preventing the overlapping of the particle, and an orientation-dependent Morse potential V_M:

V = V_WCA(r) − Ξ(p̂ⁱ,p̂^j,r̂)V_M(r) (3)

where r^ij = r^j − rⁱ is the distance vector between particle centre, r = |r^ij|, and r̂ = r^ij/r, the unit vectors p̂ⁱ, p̂^j are the directors of particles i and j, respectively, andwhere ε and ε_M are the depth of the well in the WCA and Morse potential, respectively. We set ε_M = ε in this study. M_d is the Morse potential depth factor (M_d = 2.294a), M_r is the Morse potential range parameter M_r = a, and r_eq is the Morse potential equilibrium position (r_eq = 1.878a).⁴¹ The anisotropic function Ξ(p̂ⁱ,p̂^j,r̂) depends only on the mutual orientation of the particles and is given asΞ is normalised so that −1 ≤ Ξ ≤ 1, and the negative/positive Ξ indicates attractive/repulsive interaction. An illustration of the pair potential for some given configurations is given in Fig. 2. The most favourable pair corresponds to the head to tail arrangement. When the head to tail positions are occupied, the secondary stable structure appears in the form of a side-by-side anti-parallel alignment.

Fig. 1 Illustration of (a) the patterned surface of the Y₁₀ particle and (b) a pair of particles characterised by their orientations and relative positions.

Fig. 2 Pairwise potential as a function of the centre distance for some configurations.

The simulation is conducted in the dimensionless form, where the characteristic length, energy, and time are denoted as a, ε, and a²/D^T, respectively. The time step is taken so that under the condition of one unit of force and k_BT/ε = 0.1, the particle on average moves 10⁻³a in one step. The packing fraction ρ is defined as the ratio of the volume of particles to that of the space confining them, which is L × L × 2a for the planar geometry and 4/3π[(R + a)³ − (R − a)³] for the sphere of radius R. The periodic boundary condition in the L direction is applied for the planar geometry case. The initial positions and orientations of the particles are randomly distributed. The temperature k_BT/ε decreases from 0.5 to 0.05 by intervals of 0.01, and 0.5 × 10⁶ simulation steps are performed for each value of k_BT/ε. In this article, we show the structures at the last time step. Verification of the simulation on repulsive isotropic particles is performed by setting the anisotropic function fixed as Ξ = 1.

Last but not least, the particles possibly fluctuate between certain thicknesses of the spherical layer. This effect, however, is negligible, thanks to the small enough time step and suitable harmonic potential, which is comparable to the particle–particle interaction. The systems in the study contain many interacting particles, and there is a possibility that metastable states are obtained instead of true ground states. We report a statistical result including all possible states.

2.2 Structure analysis

The coordination number N_i of particle i is estimated by combining the conventional Delaunay triangulation whose vertices are the particle positions,¹³ and a distance constraint. The particles further away than a certain value are not considered as neighbours due to the relatively short-range interaction (Fig. 2). Such a distance constraint is based on the first trough of the pairwise distance distribution, which is approximately 2.5a. Then, the local positional order of particle i is calculated by the two-dimensional bond-orientational order parameter ψ_n:^13,20where θ_ij is the angle between particle i and its neighbouring particle j on the tangent plane of particle i if the particles are confined to the spherical surface. |ψ_n| characterises the local degree of the regular n-gon order around a particle; for instance, the perfect square lattice on a flat space has |ψ₄| = 1, while the hexagonal one has |ψ₆| = 1.

Since the head to tail alignment is the most stable configuration (Fig. 2), string-like structures with alternating orientation at their sides are expected to form (see also Fig. 4). Therefore, we evaluate the orientational order of the particles via the nematic order where the head and tail of the vector are treated similarly. The local order parameter tensor Q_i is calculated by averaging the orientation p̂ of particle i and its coordinated particle N_i:

The nematic order parameter s is then three halves the positive eigenvalue of Q_i. The position of the topological defect of the director field is approximated via data of the local nematic order parameter and winding number. The orientation of a defect is simply obtained by averaging the nematic orientation in the defect core region. In particular, for the comet-like +1/2 defect, the defect orientation is defined from the head to the tail of the comet. Fig. 3 displays some configurations of a pair of +1/2 defects whose orientation is almost anti-parallel. Regarding the distance between the defects, aside from the ordinary Euclidean distance, we implement another parameter taking into account the orientation of the defects. As shown in Fig. 3, we may follow the line by connecting the orientation of the particles. Each line forms a layer. The number of layers L between a pair of +1/2 defects intrinsically includes the information of the relative position and orientation of a pair of defects, estimated by L = [d sin α/h], where d is the distance between the defect cores, α is the angle between the defect direction and the distance vector of the defect cores, and h is the thickness of the two consecutive layers. It is shown in Fig. 3 that the number of layers between a pair of +1/2 defects varies although the Euclidean distance is fixed.

Fig. 3 Examples of a pair of +1/2 defects with fixed Euclidean distance between the defect cores. A + 1 defect (a) can split into a pair of + 1/2 defects (b)–(d). The particle orientations are illustrated by the black headless arrows. The defect orientation and position are represented by the red arrow and its tail. The number of layers L between a pair of +1/2 defects from (b)–(d) is, respectively, L = 0, 3, and 7, which can be directly determined by counting the number of layers passing the two cores (the bold dashed line).

3 Results

3.1 Particles confined to planar geometry

The assembly of the dipole patchy particles on a flat space is investigated. Fig. 4 shows the average properties over all particles in terms of the bond-orientational order parameter 〈|ψ_n|〉 and the nematic order parameter 〈s〉, and typical snapshots at various packing fraction ρ from 0.35 to 0.565. The highly ordered, almost defect-free structures exhibited via the square lattice (〈|ψ₄|〉 ≈ 1) and well aligned orientation (〈s〉 ≈ 1) are more frequently obtained at the packing fraction comparable to ρ ≈ 0.48. In these cases, the particle orientations are head to tail, whereas the side by side ones are anti-parallel. At ρ < 0.48, grain boundaries emerge because an excess of voids induces more degrees of freedom of the grains. On the other hand, at dense packing ρ > 0.48, clusters of 5- and/or 6-fold coordinated particles increase in size, thus resulting in a decrease in 〈|ψ₄|〉 and an increase in 〈|ψ₆|〉 with ρ. We then perform simulations at the packing fraction ρ = 0.48 so that a square lattice structure is formed.

Fig. 4 Average properties 〈|ψ₄|〉, 〈|ψ₆|〉, and 〈s〉 of the dipole patchy particles in the planar geometry at various packing fractions. Each snapshot includes the particle positions with the local |ψ₄| (filled circles), and particle orientations with the nematic order parameter s (arrows) in the Voronoi cells. The length of the vector is taken as the particle diameter. The number of particles is N = 500; at least 10 independent runs were conducted at each point.

3.2 Particles confined to the spherical surface

Simulations of dipole-like particles confined to a spherical surface are conducted under conditions identical to the planar geometry case. Fig. 5 illustrates the assembly of particles on the sphere at the volume packing for the square lattice ρ = 0.48. Different from the almost defect-free highly ordered structures in the flat space, the assembled structure on the sphere includes regions lacking both 4-fold order |ψ₄| and nematic order s. There are four low nematic order regions whose winding number is +1/2, thus the Euler characteristic of +2 for the sphere is preserved. The anti-parallel direction of the defect core suggests that there are two pairs of +1/2 defects.

Fig. 5 A self-assembly structure of N = 500 dipole-like particles on a spherical surface with the local (a) 4-fold bond-orientational order parameter |ψ₄| and (b) nematic order parameter s inscribed in the Voronoi cells. The whole sphere consists of a front view and a back view of the two hemispheres connected by an imaginary hinge. The orientations of the particles are shown by vectors of the same size as the particle's diameter. The positions of particles are at the midpoints of the vectors. The four big arrows illustrate the four +1/2 defects. The direction of the two +1/2 defects on the front (or back) is almost anti-parallel.

3.2.1 Orientational order on the spherical surface. Fig. 6 depicts the positions and characteristics of the four +1/2 defects for the assembly of N = 500 particles. A great circle can be drawn over the four +1/2 defects so as to minimise the deviation of geodesic distance from the four defects to a great circle. At this system size, the two pairs of defects (D₁, D_1′) and (D₂, D_2′) are close to a great circle and exhibit symmetric positions via the similar distance and the number of layers between each pair. Since the layers formed by the particle orientations can be observed, we measure the number of layers L_pairs between the pair of defects (D₁, D_1′) or (D₂, D_2′), as well as the number of layers L_all crossing through the whole great circle. The relative distance between the pairs of defects is evaluated as L_pairs/L_all. As shown in Fig. 6, the distribution of the relative distance between the defect cores is approximately a constant function. The potential energy of the system is almost identical regardless of whether the defect cores are closely bound (L_pairs/L_all → 0) or uniformly distributed on a great circle (L_pairs/L_all → 0.5). Similar analysis of various system sizes is presented in Fig. 7. For the small sphere size (N = 144), the distribution does not look uniform. However, this is due to the discrete nature of the ratio L_pairs/L_all. In fact, as the sphere size increases, the distribution of the ratio L_pairs/L_all becomes flatter and independent of the system size. For all the system sizes that we have studied, the energy tends to be independent of the defect position and approaches that of the square lattice in flat space. These findings suggest that the many states of the four +1/2 defects can be applied to a wide range of system size.

Fig. 6 Distribution of the four half-strength defects for N = 500. (Top) Some configurations of the assembled particles on the spherical surface. The hemisphere shows the particle orientation with the local order parameter and a pair of defects. The location of the pair (D₁, D_1′) in the front and the pair (D₂, D_2′) in the back with the imaginary great circle is also given. From left to right, the total number of layers of the pairs (D₁, D_1′), (D₂, D_2′) and the number of layers passing through the great circle are (L_pairs, L_all) = (0, 34), (2, 34), (12, 38), and (20, 40). (Bottom) The distribution of the relative distance between the pair of defects shows similar probability. (Middle) The potential energy of spherical structures is almost independent of the defect's core position with the coefficient of variation being 0.2%. The energy is normalised to that of the perfect square lattice on planar geometry. The number of independent configurations is 100.

Fig. 7 Distribution of the relative distance between the pair of +1/2 defects (bar, to the left scale) and the potential energy (circle, to the right scale) at various system sizes. The number of independent runs for N = 144, 1000, 2000, and 4000 is 120, 110, 77, and 42, respectively.

For a large system size, the relaxation time to reach an equilibrium structure becomes longer. Even on the flat surface, defects tend to survive in the large system size. When there are extra defects on a sphere before reaching the equilibrium state, they affect the position of the four +1/2, and as a result, the numbers of layers of the pair (D₁, D_1′) and (D₂, D_2′) are no longer equal. We make sure that the system reaches the equilibrium state so that there is no disclination in the flat space and four +1/2 disclinations on a sphere for the analysis.

3.2.2 Positional order on the spherical surface. As mentioned elsewhere, perfect crystals confined to a spherical surface do not exist; topological defects in the form of disclinations are required due to the topology of the sphere. We determine the disclinations by means of the local |ψ₄|. A particle of low |ψ₄| indicates that either it has a non-square local structure or its coordination number is not equal to four. Therefore, the number of low |ψ₄| particles is somewhat proportional to the number of disclinations as long as a suitable threshold of |ψ₄| is chosen. This is similar to the correspondence between the clusters of low |ψ₆| particles and the scars of 5- and 7-fold coordinated particles for the hexagonal lattice, as shown in Fig. 8a and b. Regarding the square lattice on a large sphere size (N = 2000, Fig. 8d–f), the one-dimensional low |ψ₄| regions are observed as grain boundary scars. There are two grain boundary scars emerging from each +1/2 defect, which makes eight scars in total, which then disappear within the otherwise square lattice particles. The shape of the two scars for each +1/2 defect is always like lines connected by a certain angle. In contrast, the twelve scars of isotropic particles have more degrees of freedom in the arrangement of the scars (Fig. 8a–c).¹² For the small sphere (N = 144), two cases may occur. When the +1/2 defects are clearly observed as presented in Fig. 8g–i, the two scars around the +1/2 defect reduce in size and become two points if only the particles with |ψ₄| < 0.2 are considered. In the case when the +1/2 defects are ill-defined, the particle orientations are well aligned around the sphere's equator while lacking order at the two poles (Fig. 8l–k). At each pole, the four lowest |ψ₄| particles create a square. The eight lowest |ψ₄| points form a square antiprism. This structure is somewhat in agreement with the conjecture of minimum disclinations for the square lattice by using Euler theory, wherein at least eight 3-fold coordinated particles are required on a sphere (see Appendix A.1). In addition, the local coordination number and 3-, 5-, and 7-fold bond-orientational order parameters for the case of N = 2000 particles, as shown in Fig. 8d, are also given in Fig. 9. In contrast with |ψ₄|, the grain boundary scars have higher values of |ψ₃|, |ψ₅|, and |ψ₇|. This is consistent with the coordination number being different from four near the scars. Nevertheless, the positions of the scars are less clear from |ψ₃|, |ψ₅|, and |ψ₇| compared with the scars detected by |ψ₄|.

Fig. 8 Arrangement of positional order of particles on a sphere. (a)–(c) A snapshot of N = 2000 repulsive isotropic particles in terms of |ψ₆|, coordination number, and the icosahedral alignment of the defect scars, which are clusters of 5- and 7-fold coordinated particles. (d)–(i) Snapshot of dipole patchy particles shown in ψ₄, s and the position of the four +1/2 defects for N = 2000 (d)–(f) and N = 144 (g)–(i) particles; the red arrows mark the position and direction of the +1/2 defect. (l)–(k) The antiprism structure for N = 144 particles from different view points.

Fig. 9 Snapshot of dipole patchy particles on the sphere with the local coordination number and |ψ₃|, |ψ₅|, and |ψ₇| for the case N = 2000 shown in Fig. 8d.

Fig. 10 and 11 show the dependence of the number of particles whose local bond-orientational order parameter |ψ₄| and nematic order parameter s are lower than a threshold on the system size N. The relation of the number of low |ψ₄| particles to the total number of particles is displayed in Fig. 10. Irrespective of the threshold of |ψ₄|, the relation between the low |ψ₄| particles and the system size is n_low|ψ4| ∝ N^0.5, implying that the grain boundary scars increase almost linearly with the radius of the sphere. Such a linear dependence has been observed for isotropic particles on a sphere, wherein the excess disclinations including 5- and 7-fold particles linearly increase with the sphere radius.^21,42 The dependence of the number of low nematic order particles on system size is different from that of |ψ₄| in terms of scale and trend (Fig. 11). The number of low s particles increases with the system size for the threshold s < 0.8 and s < 0.7. However, as the threshold is lowered to s < 0.6, the number of low s particles is reduced to a few dozens regardless of the size of the sphere. As given in Fig. 8d and e, there is a strong connection between the low |ψ₄| and low s regions. However, the orientations of the particles at the grain boundary scars are distorted at a different level: the nematic order near the +1/2 defect core is significantly lower than the others (Fig. 8e). In other words, the defects of the orientational order consist of their cores surrounded by fringes on which the nematic order s has an intermediate value. The defects of the positional order at which particles do not form the square lattice appear as 1-dimensional lines, which have uniformly low |ψ₄|. This observation suggests that the nematic order exhibits point (0-dimensional) +1/2 disclinations, but due to the interaction with the positional order, the cores of the +1/2 disclinations may have tails. The number of particles involved in the cores does not grow as the system size increases, whereas the number of particles in the tails may grow similarly to the defects of the positional order. We expect that these tails are unique in our system that has both positional and orientational order.

Fig. 10 Relation of the number of low 4-fold order particles to the total number of particles N at different constraints of |ψ₄| values. The lines are a guide for the eye and have a slope coefficient of 0.5.

Fig. 11 Relation of the number of low nematic order particles to the total number of particles N at different constraints of s values.

4 Discussion and summary

We have performed the assembly of dipole patchy particles on a spherical surface. Defects appear as a requirement due to the topology of the sphere. As for the orientations of the particles, there are many states of the four +1/2 defects near a great circle. The appearance of many stable states of four half-strength defects of the director field is also observed in spherical smectics,^43,44 spherical nematics in the extreme elastic constant limit,^16,17 and block copolymer assembly.⁴⁵ In these studies, the four +1/2 defects can be explained by the cut-and-rotate on a sphere of two +1 defects. This operation seems to be applied well to our system. However, we have not found the clear appearance of two +1 defects in the form of concentric layers of particle orientation (Fig. 3a). Even when the two +1/2 defects are very close and topologically equal to a +1 defect (see Fig. 3b for illustration), the two separated clusters of +1/2 defects are clearly identified (see Fig. 6-upper left for simulation). Such a difference is possibly brought about by either the non-zero temperature simulation or the soft, anisotropic interacting potential used in our study.

Regarding the topological defect of the square lattice order, the emergence of the grain boundary scars and its linear relation with the sphere radius are moderately analogical to that of the hexagonal lattice.¹² The grain boundary scars are revealed to be associated with the four +1/2 defects, rather than sitting at the vertices of a cube. This is expected to be unique in our particle system where both orientational order and positional order are present. In our system, there are two elastic constants associated with the orientational and positional order. They are not independent, but interrelated with each other. For example, when ε_M in eqn (5) vanishes, there is no longer orientational order, but in this case, the particles do not form the square lattice and break the positional order. The elastic constants set the correlation length of the orientational and positional order. Therefore, it would be an interesting future study to extend our system so that the elastic constants of orientational and positional order can be controlled independently.

Although the 4-fold bond-orientational order parameter ψ₄ is sufficient enough to evaluate the disclination, it is still inadequate to determine the exact 3-fold or 5-fold coordinated particles in a square lattice. This is in contrast to the hexagonal lattice formed by isotropic particles (Fig. 8b). The defect scars in the hexagonal lattice are clearly identified by chains of disclinations with their alternating topological charges. One should consider further development for the quadrilateral mesh, especially the irregular regions near a defect. This also gives rise to a question on defect analysis for more complex structures.

While the simulations in our study focus on patchy particle assembly on a thin shell, it is worth extending the study to the effect of the shell's thickness. For spherical nematics, as the thickness increases, the behaviour of defects becomes more complicated because a three-dimensional problem is involved, and the inhomogenity of thickness may be induced.^28,46 Regarding the effect of surface shape, it is known that the shape has a strong impact on the localisation of defect scars (for isotropic particles)^12,24,47 and disclinations (for nematics).^33,48,49 Investigating such effects for patchy particle systems that have higher degrees of freedom than the aforementioned ones seems promising for theoretical and practical applications.

The complex interaction between the positional order and orientational order of particles on a sphere is of fundamental interest for two-dimensional melting. Tuning the potential, e.g., varying the current short-range interaction to the long-range one, or the type of patchy particle, or deformability of the surface, may give rise to a better understanding on the role of the energy-driven factor in the in-plane order. From the materials science aspect, knowing the precise defect structure is important for fabricating the building block, for example, the eight scars in the study can be functionalised by chemical linkers, which can then serve as interacting sites.

Conflicts of interest

There are no conflicts to declare.

A Appendix

A.1 Derivation of minimum number of positional defects from Euler theorem

Suppose a closed surface is facetted and divided into a number of V vertices, E edges, and F faces. Euler theorem states that

V − E + F = χ (9)

where the Euler characteristic χ = 2(1 − g), and g is the genus or the number of holes of a closed surface, for example χ = 2 for all polyhedra because they are topologically equal to a sphere.¹⁹ One may imagine that such a facetted surface is similar to a net embedded on the surface. If we restrict the ring of the net to having c vertices, then

E = cV/2 (10)

On the other hand, each node or vertex of the net may have z connection with the others i.e., there would be N_z dual polygons with z sides. One may find

After substituting, we obtain

For example, applying the above equation for the triangular lattice (c = 3) and quadrilateral lattice (c = 4) gives and , respectively. This means that on a sphere, the minimum number of disclinations is twelve 5-fold coordinated nodes for the triangular lattice and eight 3-fold coordinated nodes for the quadrilateral one.

Acknowledgements

The authors are grateful to Dr Toru Kajigaya for his helpful discussions. The authors acknowledge the support from JSPS KAKENHI Grant number JP20K14437 to U. T. L., and JP20K03874 and JP20H05259 to N. Y.

References

1. Q. Chen, S. C. Bae and S. Granick, Nature, 2011, 469, 381–384.
2. Q. Chen, E. Diesel, J. K. Whitmer, S. C. Bae, E. Luijten and S. Granick, J. Am. Chem. Soc., 2011, 133, 7725–7727.
3. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck and D. J. Pine, Nature, 2012, 491, 51–55.
4. R. M. Choueiri, E. Galati, H. Thérien-Aubin, A. Klinkova, E. M. Larin, A. Querejeta-Fernández, L. Han, H. L. Xin, O. Gang, E. B. Zhulina, M. Rubinstein and E. Kumacheva, Nature, 2016, 538, 79–83.
5. F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys., 2010, 132, 184501.
6. F. Romano, E. Sanz, P. Tartaglia and F. Sciortino, J. Phys.: Condens. Matter, 2012, 24, 064113.
7. F. Romano and F. Sciortino, Nat. Commun., 2012, 3, 975.
8. L. Hong, A. Cacciuto, E. Luijten and S. Granick, Nano Lett., 2006, 6, 2510–2514.
9. L. Hong, A. Cacciuto, E. Luijten and S. Granick, Langmuir, 2008, 24, 621–625.
10. D. R. Nelson, Nano Lett., 2002, 2, 1125–1129.
11. Z.-L. Zhang, A. S. Keys, T. Chen and S. C. Glotzer, Langmuir, 2005, 21, 11547–11551.
12. M. J. Bowick and L. Giomi, Adv. Phys., 2009, 58, 449–563.
13. R. E. Guerra, C. P. Kelleher, A. D. Hollingsworth and P. M. Chaikin, Nature, 2018, 554, 346–350.
14. T. C. Lubensky and J. Prost, J. Phys. II, 1992, 2, 371–382.
15. T. Lopez-Leon, V. Koning, K. B. S. Devaiah, V. Vitelli and A. Fernandez-Nieves, Nat. Phys., 2011, 7, 391–394.
16. H. Shin, M. J. Bowick and X. Xing, Phys. Rev. Lett., 2008, 101, 0037802.
17. M. A. Bates, J. Chem. Phys., 2008, 128, 104707.
18. A. M. Turner, V. Vitelli and D. R. Nelson, Rev. Mod. Phys., 2010, 82, 1301–1348.
19. S. Hyde, B. W. Ninham, S. Andersson, K. Larsson, T. Landh, Z. Blum and S. Lidin, The Language of Shape, Elsevier, 1997, pp. 1–42.
20. D. R. Nelson, Defects and Geometry in Condensed Matter Physics, Cambridge University Press, Cambridge, New York, 2002.
21. M. J. Bowick, D. R. Nelson and A. Travesset, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 8738–8751.
22. S. Giarritta, M. Ferrario and P. Giaquinta, Phys. A, 1993, 201, 649–665.
23. T. Einert, P. Lipowsky, J. Schilling, M. J. Bowick and A. R. Bausch, Langmuir, 2005, 21, 12076–12079.
24. F. L. Jiménez, N. Stoop, R. Lagrange, J. Dunkel and P. M. Reis, Phys. Rev. Lett., 2016, 116, 104301.
25. K. J. Strandburg, Rev. Mod. Phys., 1988, 60, 161–207.
26. P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, 1st edn, 1995.
27. R. D. Kamien, Rev. Mod. Phys., 2002, 74, 953–971.
28. V. Vitelli and D. R. Nelson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 021711.
29. W.-Y. Zhang, Y. Jiang and J. Z. Y. Chen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 061710.
30. W.-Y. Zhang, Y. Jiang and J. Z. Y. Chen, Phys. Rev. Lett., 2012, 108, 057801.
31. Q. Liang, S. Ye, P. Zhang and J. Z. Y. Chen, J. Chem. Phys., 2014, 141, 244901.
32. S. Kralj, R. Rosso and E. G. Virga, Soft Matter, 2011, 7, 670–683.
33. V. Koning and V. Vitelli, Fluids, Colloids and Soft Materials, John Wiley & Sons, Ltd, 2016, pp. 369–385.
34. Y. Li, H. Miao, H. Ma and J. Z. Y. Chen, Soft Matter, 2013, 9, 11461.
35. O. V. Manyuhina and M. J. Bowick, Phys. Rev. Lett., 2015, 114, 117801.
36. D. Wang, M. Hermes, R. Kotni, Y. Wu, N. Tasios, Y. Liu, B. de Nijs, E. B. van der Wee, C. B. Murray, M. Dijkstra and A. van Blaaderen, Nat. Commun., 2018, 9, 2228.
37. D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352–1360.
38. E. Dickinson, S. A. Allison and J. A. McCammon, J. Chem. Soc., Faraday Trans. 2, 1985, 81, 591–601.
39. P. Castro-Villarreal, A. Villada-Balbuena, J. M. Méndez-Alcaraz, R. Castañeda-Priego and S. Estrada-Jiménez, J. Chem. Phys., 2014, 140, 214115.
40. M. M. Moghani and B. Khomami, Soft Matter, 2013, 9, 4815.
41. R. A. DeLaCruz-Araujo, D. J. Beltran-Villegas, R. G. Larson and U. M. Córdova-Figueroa, Soft Matter, 2016, 12, 4071–4081.
42. A. R. Bausch, M. J. Bowick, A. Cacciuto, A. D. Dinsmore, M. F. Hsu, D. R. Nelson, M. G. Nikolaides, A. Travesset and D. A. Weitz, Science, 2003, 299, 1716–1718.
43. X. Xing, J. Stat. Phys., 2009, 134, 487–536.
44. F. Serafin, M. J. Bowick and S. R. Nagel, Eur. Phys. J. E: Soft Matter Biol. Phys., 2018, 41, 143.
45. T. L. Chantawansri, A. W. Bosse, A. Hexemer, H. D. Ceniceros, C. J. García-Cervera, E. J. Kramer and G. H. Fredrickson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 031802.
46. D. Seč, T. Lopez-Leon, M. Nobili, C. Blanc, A. Fernandez-Nieves, M. Ravnik and S. Žumer, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 020705.
47. W. T. M. Irvine, V. Vitelli and P. M. Chaikin, Nature, 2010, 468, 947–951.
48. L. Mesarec, W. Góźdź, A. Iglič and S. Kralj, Sci. Rep., 2016, 6, 27117.
49. L. Mesarec, W. Góźdź, A. Iglič, V. Kralj-Iglič, E. G. Virga and S. Kralj, Sci. Rep., 2019, 9, 19742.

---