Topological defects of dipole patchy particles on a spherical surface

We investigate the assembly of the 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_{10}$, 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.


I. INTRODUCTION
Patchy particles are particles of colloidal size and have patches playing 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 the "conventional" materials [1]. Recent developments in synthesis techniques [2][3][4] have made it feasible to fabricate patchy particles of high degrees of freedom [2,[5][6][7][8][9]] and attracted research 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 control defect is important for those tasks. Defects are imperfections and singularities in an ordered structure. From a microscopic point of view, defects can be served 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 can not disappear by continuous deformation of the order parameter. We focus on spherical surface for its substantial applicability such as crystalline membrane, design of crystalline materials, 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 two-dimensional plane is a defect-free triangular lattice containing 6-fold coordinated particles. On a spherical surface, such a translational symmetry is broken, at least twelve 5-fold coordi- * uyen.lieu@aist.go.jp † yoshinaga@tohoku.edu.jp nated 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, the rod-like particles confined to a spherical surface form four +1/2 or two +1 defects [10,14,15] instead of a defect-free nematic phase (longrange 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. Although the appearance of the topological defects is generally accepted as a consequence of the embedded geometry and topology, it remains unclear about: the underlying mechanism on the formation of defects and the dependence on the system size, the influence of the particle-particle and particle-curvature interaction to the defect, and the interplay between positional and orientational order. 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 and behaves somewhat similar to magnetic bead. We dynamically simulate and compare the assembly of such patchy particles confined to a planar geometry and a spherical surface.

Topological defects on curved surface
It is useful to discuss about 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, 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 [16]. If we restrict every face has c vertices, and let N z be the number of vertices which have z connections with the others, then the Euler theorem can be written as (see 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 hence the 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 [17]. When such a particle system is confined to a spherical surface, it is straightforward from eqn (1) that a net charge of z (6 − z) N z = 12 is required. For sphere size below a critical value, the twelve 5-fold coordinated particles are icosahedrally located among the otherwise 6-fold ones owning to the repulsion of the like-sign charge [17]. As the sphere increases, those 5-fold disclinations are screened by additional dislocations which are pairs of 5-and 7-fold coordinated particles [18], noted that it is not limit to the 4-or 8-fold 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 crystal are observed in both experiments and simulation [12,13,19]. 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 [18]. For a square lattice on the sphere, the evidence of how disclinations distribute is still lacking. However, it can be conjectured that by employing a quadrilateral mesh (c = 4) for square lattice, a disclination is now the vertex whose connection differs from four. Eqn (1) becomes z (4 − z) N z = 8, meaning that at least eight disclinations of z = 3 is required for a square lattice on a sphere. The development of dislocations which are now viewed as bound disclination pairs of 3-and 5-fold coordinated particles [20] are 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 [14]. The strength of a defect in this case is characterised via winding number k, defined as how much 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 rotates in one counter-clockwise circuit [21]. On a closed surface, the net strength is equal to Euler characteristic according to Poincaré-Hopf theorem [22]. 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 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,23]. In contrast, in the extreme limit when splay is much softer (or harder) than bending, it is suggest that a +1 defect can split into two +1/2 defect 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-androtate surgery of the sphere with two +1 defects [24]. 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 which may position on the vertices of an anti-cube [14,25] or cube [26,27]. Such difference is perhaps caused by the different types of interaction for simulation and the forms of free energy for continuous description. It is worth noting that the p-atic particles in those mentioned systems has purely rotational degree of freedom, hence the existence of positional order remains unclear.

A. Brownian Dynamics
We employ the Brownian dynamics simulation algorithm for particles in Euclidean space in overdamped limit [28,29]. 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 follows where r(t + ∆t), Ω(t + ∆t) denote the position and orientation of the particle after the time step ∆t; D T , D R are the translational and rotational diffusion coefficients of an isolated particle, respectively; δr, δΩ are the translation and rotation due to thermal fluctuation, satisfying i is independently chosen from a Gaussian distribution with zero mean and unit variance. The 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 [30], in which the tangential parts of the interacting force F and noise δr in eqn (2a) at the point where particle is located are considered; finally an addition 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 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 Y 10 (x) = 3/(4π)p ·x wherep is defined as the orientation of the particle (Fig. 1). This pattern has positive and negative hemispheres similar to the Janus particle [8,31,32]. The interaction potential for a pair of particles i and j comprises an isotropic Week-Chandler-Anderson potential V W CA preventing the overlapping of particle, and an orientation-dependent Morse potential V M : where r ij = r j −r i is the distance vector between particle centre, r = r ij , andr = r ij /r, the unit vectorp i ,p j is the director of particles i, j, respectively, and where ε is the potential well depth, 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) [32]. The anisotropic function F (p i ,p j ,r) depends only on the mutual orientation of the particles and is given as F is normalised so that −1 ≤ F ≤ 1, the negative/positive F indicates attractive/repulsive interaction. The illustration of the pair potential for some given configurations is given in Fig. 2. The most favourable pair is correspondent to the head to tail arrangement. When the head to tail positions are occupied, the secondary stable structure appears in the form of side-by-side anti-parallel alignment. a 2 /D T , respectively. The time step is taken so that under the condition of one unit of force and k B T /ε = 0.1 the particle on average moves 10 −3 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π for the sphere of radius R. The periodic boundary condition in L direction is applied for the planar geometry case. The initial positions and orientations of the particles are randomly distributed. The temperature k B T /ε decreases from 0.5 to 0.05 by intervals of 0.01, and 0.5 × 10 6 simulation steps is performed for each value of k B T /ε. 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 F = 1. Last but not least, the particles are possibly fluctuates between certain thickness 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, there is possibility that metastable states are obtained instead of true ground states. We report a statistical result including all possible states.

B. Structure analysis
The coordination number N i of particle i is estimated by combining the conventional Delauney triangulation whose vertices are the particle positions [13], and a distance constraint. The particles further than a certain value are not considered as neighbours due to the relatively short-range interaction (Fig. 2). Such 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,17]: where θ 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 flat space has |ψ 4 | = 1, while the hexagonal one has |ψ 6 | = 1.
Since the head to tail alignment is the most stable configuration (Fig. 2), string-like structures with alternating orientation at their side are expected to occur (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 of the particles i and its coordinated particles 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 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. Figure 3 displays some configurations of a pair of +1/2 defects whose orientation is almost anti-parallel.
Regarding the distance of the defects, aside from the ordinary Euclidean distance, we implement another parameter taking into account the orientation of the defects. 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 of the defect cores, α is angle between the defect direction and the distance vector of the defect cores, and h is the thickness of two consecutive layers. It is shown in Fig. 3 that the number of layer between a pair of +1/2 defects varies although the Euclidean distance is fixed.

A. Particles confined to planar geometry
The assembly of the dipole patchy particles on a flat space is investigated. Figure 4 shows the average properties in terms of bond-orientational order parameter and nematic order parameter, and the typical snapshots at various packing fractions ρ from 0.35 to 0.565. The highly

FIG. 3: Examples of a pair of +1/2 defects with fixed
Euclidean distance between the defect cores. The particle orientations are illustrated by the 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 (a-c) are respectively L = 0, 3, 7.
ordered, almost defect-free structures exhibited via the square lattice ( |ψ 4 | ≈ 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 their sizes, thus the decrease in |ψ 4 | and increase in |ψ 6 | with ρ. We then perform simulations at the packing fraction ρ = 0.48 so that a square lattice structure is formed.

B. Particles confined to spherical surface
Simulations of dipole-like particles confined to a spherical surface are conducted at the conditions identical to the planar geometry case. Figure 5 illustrates the assembly of particles on the sphere at the volume packing for square lattice ρ = 0.48. Different from the almost defectfree highly ordered structures in the flat space, the assembled structure on the sphere includes regions lacking both 4-fold order |ψ 4 | 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 that there are two pairs of +1/2 defects. Figure 6 depicts the positions and characteristics of the four +1/2 defects for the assembly of N = 500 particles. At this size, the two pairs of defects (D 1 , D 1 ) and (D 2 , D 2 ) exhibit symmetric positions via the similar distance and the number of layers between each pair. A great circle can be drawn over the four +1/2 defect while minimising the deviation. The number of layers crossing through the whole great circle can be determined. 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 position of the defect cores has a similar probability. The potential energy of the system is almost identical regardless 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 for various system sizes is presented in Fig. 7. For the small sphere size (N = 144) there is some uneven distribution which is due to the discrete nature of the ratio L pairs /L all . As the sphere size increases, the ratio L pairs /L all is less discrete, but defect is more likely to grow. Such kind of defect affects the position of the four +1/2, for example, the numbers of layers of the pair (D 1 , D 1 ) and (D 2 , D 2 ) are no longer equal. In general, the distribution of L pairs /L all is not significantly different. 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.

Positional order on spherical surface
As mentioned elsewhere, perfect crystals confined to 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 |ψ 4 |. A particle of low |ψ 4 | indicates that either it has a non-square local structure or its coordination number is not equal to four. Therefore, the number of low |ψ 4 | particles is somewhat proportional to the number of disclinations as long as a suitable threshold of |ψ 4 | is chosen. This is similar to the correspondence between the clusters of low |ψ 6 | particles and the scars of 5-, 7-fold coordinated particles for the hexagonal lattice, as shown in Fig. 8a-b. Regarding the square lattice on large sphere size (N = 2000, Fig. 8d-f), the one-dimensional low |ψ 4 | regions are observed as grain boundary scars. There are two grain boundary scars emerging from each +1/2 defects, which makes eight scars in total, then disappearing 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 on the arrangement the scars (Fig. 8a-c) [12]. 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 reduces its size, and become two points if only the particles with |ψ 4 | < 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 lacks of order at the two poles ( Fig. 8l-k). At each pole, the four lowest |ψ 4 | particles create a square. The eight lowest |ψ 4 | points form a square antiprism. This structure is somewhat in agreement with the conjecture of minimum disclinations for square lattice by using Euler theory, that at least eight 3-fold coordinated particles is required on a sphere (see Appendix V A).  Fig. 9. Although the threshold of |ψ 4 | varies, the number of non-square particle n low|ψ4| ∝ N 0.5 , implying that the grain boundary scars increase almost linearly with the radius of the sphere. Such linear dependence has been observed for isotropic particles on a sphere, that the excess disclinations including 5-and 7-fold particles linearly increases with the sphere radius [18,33]. The dependence of the number of low nematic order particles on system size is different from that of |ψ 4 | in terms of scale and trend (Fig. 10). 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 disorder particles below the threshold reduced to a few dozens regardless the size of the sphere. As given in Fig. 8(d-e), there is a strong connection between the low |ψ 4 | 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. In other words, the distinguishing behaviour of the number of non-square particles and the number of disorder nematics to the system size can be originated from the characteristic of the 1-dimensional grain boundary scar and the 0-dimensional +1/2 disclination, plus the interaction between them.  We have performed the assembly of dipole patchy particles on a spherical surface. Defects appear as a requirement by 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 director field is also observed in spherical smectics [34,35], spherical nematics in the extreme elastic constant limit [24], block copolymer assembly [36]. In those 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  applied well to our system. However, we have not found the clear appearance of a two +1 defects in the form of concentric layers of particle orientation; and a deformed form is found instead (see Fig. 3a,6). Such 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 [12]. The grain boundary scars are revealed to be associated with the four +1/2 defects. It is also possible that such mutual interplay between the positional order and the orientational order suppresses the defect-defect interaction, leading to the many position of the four +1/2 defects. Although the 4-fold bond-orientational order parameter ψ 4 is sufficient enough to evaluate the disclination, it is still inadequate to determine the exact 3-fold or 5fold coordinated particles in a square lattice. One should consider further development for the quadrilateral mesh, especially the irregular regions near a defect. This also gives rise to a question on the defect analysis for more complex structures.
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, deformability of the surface, may bring better understanding on the role of the energy-driven factor on the in-plane order. On the materials science aspect, knowing the precise defect structure is important for fabricating building block, for example, the eight scars in the study can be functionalised by chemical linkers, then served as interacting sites.

V. APPENDIX
A. 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, F faces. Euler theorem states that 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 [16]. One may image that such facetted surface is similar to a net embedded on the surface. If we restrict the ring of the net has c vertices, then 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 z (6 − z) N z = 6χ and z (4 − z) N z = 4χ, respectively. It means on a sphere, the minimum number of disclinations is twelve 5-fold coordinated nodes for triangular lattice and eight 3-fold coordinated nodes for quadrilateral one.

CONFLICTS OF INTEREST
There are no conflicts to declare.