Colloidal interactions and unusual crystallization versus de-mixing of elastic multipoles formed by gold mesoflowers

Colloidal interactions in nematic liquid crystals can be described as interactions between elastic multipoles that depend on particle shape, topology, chirality, boundary conditions and induced topological defects. Here, we describe a nematic colloidal system consisting of mesostructures of gold capable of inducing elastic multipoles of different order. Elastic monopoles are formed by relatively large asymmetric mesoflower particles, for which gravity and elastic torque balancing yields monopole-type interactions. High-order multipoles are instead formed by smaller mesoflowers with a myriad of shapes corresponding to multipoles of different orders, consistent with our computer simulations based on free energy minimization. We reveal unexpected many-body interactions in this colloidal system, ranging from de-mixing of elastic monopoles to a zoo of unusual colloidal crystals formed by high-order multipoles like hexadecapoles. Our findings show that gold mesoflowers may serve as a designer toolkit for engineering colloidal interaction and self-assembly, potentially exceeding that in atomic and molecular systems.

I ntroduced by Einstein within a theoretical framework explaining Brownian motion of tiny particles 1 , colloidal atom paradigm has provided the motivation and means for organizing particles into crystals and other structures, mimicking and even exceeding the diversity of structures in naturally occurring molecular and atomic materials 2 . Long-range elasticity-mediated colloidal interactions between particles 3 in liquid crystal (LC) 4 fluids have enabled a host of anisotropic colloidal self-assemblies and composite materials  . Colloidal inclusions perturb the uniform background of nematic LC's ground-state unidirectional molecular orientations, producing distortions in the molecular ordering described by the coordinate-dependent director field n (r). These director distortions propagate far beyond the physical extent of the particles themselves 3 , though confining surfaces with strong boundary conditions can partially localize and limit the extent of spatial propagation of these distortions 11 . Minimization of director distortions to lower the total free energy cost when colloidal particles are in a close proximity leads to elasticitymediated interactions not present in isotropic host media 3,11 . Under the one-elastic-constant approximation, the governing Euler-Lagrange equation, derived from the minimization of free energy, is of Laplace's type, similar to that of electrostatics, thus allowing for interpreting the nature of long-ranged n(r) distortions based on the multipole expansions 3,11,21 . At small interparticle distances, this multipole description is limited by the presence of topological singularities and non-spherical topographic features of the colloidal particles, whereas surface confinement and boundary conditions on sample surfaces effectively limit this description at large inter-particle distances 11 . The confinement effects effectively screen the long-range nature of interactions and can be accounted for analogously to the screening of electrostatic colloidal interactions in presence of counterions and many other types of screening of physical forces described using the mathematical language of multipoles [27][28][29][30][31] . Deviations of n(r) in opposite directions away from the far-field uniform alignment can be interpreted analogously to opposite charges in electrostatic charge distributions, defining the design principles for achieving diverse types of colloidal interactions and assemblies that mimic the well-understood interactions between electrostatic charge distributions 11 . In addition to theoretical analysis [5][6][7][12][13][14][15]17 , a number of elastic multipoles have been discovered experimentally 3,[8][9][10][11][16][17][18][19][20][21] . Surface anchoring boundary conditions on the particles and size, shape, topology and chirality are all found to be important factors, defining the behavior of nematic LC colloids 11 . Besides, colloidal particles can induce different multipoles depending on the types of defects that occur. For example, colloidal spheres can induce elastic dipoles 3 , quadrupoles 12 , or hexadecapoles 17 depending on whether singular point 3 or "Saturn ring" disclination loop 32,33 or simultaneously both types of defects 17 are formed, respectively. However, such ability of achieving different types of elastic multipoles by colloidal inclusions of the same type is limited, consequently limiting the diversity of assemblies and composite materials that can be achieved. Although two-photon polymerization based fabrication of colloidal particles with complex shapes can allow for designing many different types of elastic LC multipoles 11,16 , it is limited to particles made of polymers and elastomers and cannot be easily scaled, which is a limitation as compared to wet chemical synthesis of colloids.
Here we describe a nematic colloidal system made from mesoflower colloidal particles 34 capable of inducing elastic multipoles of different order. These mesoflowers are mesostructures of gold with highly diverse shapes and with characteristic sharp spikes of sub-micron dimensions 34 . These complex yet diverse particles allow for inducing different elastic multipoles when dispersed in a nematic LC, which for small particles range from dipoles to hexadecapoles, and even higher order multipoles. Larger asymmetric mesoflowers induce elastic monopoles that emerge due to external gravitational torques/forces that are balanced by those originating from the LC's orientational elasticity; this leads to monopole-type interactions with confining substrates and other colloidal objects, as well as to unusual anisotropy in their Brownian motion. Numerical modeling based on Landau-de Gennes free energy minimization, while accounting for the presence of singular defects with corresponding variations of the scalar order parameter in addition to the director configurations, confirms these experimental observations and predicts the existence of other elastic multipoles achieved by systematically varying the distribution of spikes. We reveal changes in physical behavior like pair interactions of colloidal particles stemming from small changes in particle dimensions and shapes, as exemplified by the mesoflower nematic colloidal system. Numerical modeling of interactions in large colloidal systems shows that these effects lead to unusual types of colloidal crystals and transformations between them in the case of high-order multipoles and to de-mixing of elastic monopoles of opposite signs, very differently from electrostatic interactions in atomic systems. In the spirit of the colloidal atom paradigm, our findings reveal that LC colloids have a great potential of not only expanding the length scales of self-assembly from atomic to colloidal scales 1-3 , but also diversifying the forms of colloidal organization by going beyond what is accessible to atomic systems.

Results
Experimental generation of elastic multipoles. Gold mesoflowers of size ranging from hundreds of nanometers to micrometers are synthesized by a seed-mediated growth method (see details in "Methods" section) 34 . They are then dispersed in a nematic LC, 4-cyano-4′-pentylbiphenyl (5CB) (Fig. 1). The cetyltrimethylammonium bromide (CTAB) coating on the surface of the particles sets perpendicular boundary conditions for the director n(r), while the sharp spikes sticking out in all directions perturb the uniform far-field alignment n 0 of the LC defined by the rubbing direction of the confining substrates. Although each mesoflower possesses varying number of spikes of different size, the analysis of n(r) around the particles allows for prediction of their colloidal behavior based on electrostatic analogy 11 . In the optical micrographs taken under crossed polarizers with an additional 530 nm phase retardation plate, director distortions manifest themselves as the colored regions different from the background, indicating that n(r) deviates away from n 0 around a particle. The direction of director rotations, or rather the rotation of projection of n(r) to the plane of the sample, in blue (yellow) regions of polarizing optical micrographs is extracted on the basis of addition (subtraction) of phase retardation of the waveplate and the birefringent LC sample with the corresponding director orientation patterns. For our experimental geometry, the blue (yellow) polarized interference colors in optical micrographs reveal positive (negative) x-components of director and clockwise (counterclockwise) rotations of n(r) away from n 0 (see the insets in figures showing the details for particular experiments discussed). The diversity of the mesoflowers leads to a variety of polarized interference color patterns, revealing n(r)-distortions resembling that of elastic multipoles 21,26 , depending on the exact patterns of inter-changing blue and yellow colored regions distributed around the particles (Fig. 1h-k and Supplementary  Fig. 1). In some peculiar cases, the mesoflowers are surrounded by predominately one color (Figs. 1g and 2a, b), implying n(r) rotation to one direction away from n 0 . Such director distortions can be identified as elastic monopoles, although it has long been believed that they should relax to higher order multipoles 11,18,26 because rotations of director in one direction would generate elastic torque that would relax such distortions to minimize free energy. In our system, however, gravity prompts this behavior because the density of gold is much higher than that of the LC and can serve as a source of external torques/forces balancing their elastic counterparts. For sufficiently large particles, gravitational forces and torques can compete with the elastic counterparts and the elastic monopoles can be stabilized rather than relax to higher order multipoles. This behavior is very different from that of conventional colloids in isotropic fluid hosts, where the role of gravity is associated with colloidal particle sedimentation, redistribution along the sample height or destabilization. As an order of magnitude assessment of these unusual gravity effects in LCs, we equate effective gravitational potential of the mesoflower ∝Δρ(4πR 3 /3)gR (here Δρ is the difference between densities of gold 19,320 kg m −3 and the LC 1008 kg m −3 ; g = 9.8 m s −2 is the standard gravity; R is an effective radius of the mesoflower) with the elastic energy / KR (the used LC's average elastic constant K~6.5 pN). We obtain an estimate of effective threshold radius for particles significantly influenced by gravitational effects, R t~2 μm, well within the range of the studied mesoflower dimensions (Fig. 1). Thus, depending also on particle shape, gravitational effects can be a factor for stabilizing different elastic multipoles when the particle size is comparable or larger than R t . For smaller particles, gravity is not strong enough to compete with elastic forces (Fig. 1d, e), but it can serve as a source of symmetrybreaking torques and forces for particles larger than R t . Particles with intermediate dimensions~R t may exhibit monopoles as metastable states due to the interplay between the complex shape of the particle, surrounding director distortions and satellite defects, as well as interactions with the confining substrate. Indeed, when we poke, heat or rotate such particles by laser tweezers, different multipole-like color patterns can appear around the same mesoflower of size~R t ( Fig. 2 and Supplementary Fig. 1), including the monopole-like structures.
The intrinsic viscoelastic anisotropy of the LC host leads to anisotropic Brownian motion of colloidal inclusions 35   further depends on the particle's geometric shape, surface boundary conditions and induced defects 19,36,37 . For mesoflowers with comparable effective dimensions of spikes extending in directions along and perpendicular to n 0 , the angular dependence of the mean square displacement (MSD) of individual particles is directly related to n(r) structures that they induce (Fig. 2). In an isotropic phase of the LC, like colloidal spheres, mesoflowers diffuse in all directions with roughly equal probability during a long period of time 38 , so that the angular dependence of MSD is isotropic. In a nematic phase, even for spherical particles, molecular alignment breaks the symmetry and defines an easy axis for particle diffusion 35 , which yields dumbbell-shaped angular dependencies of MSD, aligned along n 0 and symmetric with respect to it. However, the mesoflower-induced n(r) distortions further enrich this behavior (Fig. 2). For example, the dumbbell-shaped MSD dependence of elastic monopoles has a long axis tilted away from n 0 (Fig. 2a, b), with the tilting direction matching the unidirectionally rotated nearby n(r) orientation and correlating with the elastic monopole sign. In contrast, mesoflowers inducing higher order multipoles exhibit MSD angular dependence symmetric with respect to n 0 (Fig. 2c, d and Supplementary Fig. 2). Interestingly, the diffusion behavior can be altered by switching particle-induced structures between different metastable states associated with reconfiguration of n(r) and elastic multipoles (Fig. 2a, c). Such switching of diffusion anisotropy and medium-mediated long-distance correlation of diffusion anisotropy between particles of the same type, cannot be achieved for mesoflowers or other particles dispersed in an isotropic medium, even for anisotropic particles 38 .
Elasticity-mediated colloidal interactions. Dynamics and interactions of mesoflowers dispersed in the LC are probed under polarizing optical microscopy with a 530 nm retardation plate inserted, so that elastic multipoles can be identified by examining the color patterns. Particles are brought to desired initial positions using laser tweezers and then released. In addition to conventional dipole-dipole ( Supplementary Fig. 3) and quadrupolequadrupole ( Fig. 3) interactions that are common in other nematic colloids 11 , interactions involving elastic monopoles are also observed (Fig. 4a, b). A monopole-like mesoflower surrounded by unidirectionally rotated n(r) (consistent with the blue color in a polarizing micrograph) attracts another mesoflower with a dipolar n(r), surrounded by approximately equal amount of blue and yellow colors within the polarizing micrograph, which is an elastic dipole. The two particles eventually approach each other by sharing blue-colored regions, thus lowering the total energy cost of the ensuing colloidal assembly (Fig. 4a, b). The interaction potential, calculated from the balance of elastic and viscous drag ∝dr c /dt forces, is in the range of hundreds of k B T and its power-law distance dependence ∝−r c −2 is consistent with that of the monopole-dipole interaction (r c is the distance between the centers of interacting colloidal particles). As another example, Fig. 4c, d shows how an assembly of two mesoflowers consisting of one dipole and one monopole attracts another dipolar mesoflower (Fig. 4c, d), with the interactions again consistent with the electrostatic analogy of these nematic colloids. For elastic multipoles of leading orders 2 l and 2 m , the balance of viscous drag and elastic force / 1=r lþmþ2 c yields the anticipated time dependence of inter-particle distance r c t ð Þ ¼ r lþmþ3 0 À Â ðl þ m þ 3Þαt 1=ðlþmþ3Þ , where r 0 is the initial center-to-center distance and α is a fitting parameter, consistent with experiments for pair interactions for all studied multipoles of the same and different orders (Figs. 3 and 4 and Supplementary Fig. 3).

Computer simulations of elastic multipoles and interactions.
Numerical minimization of Landau-de Gennes free energy provides insight into the director configurations around mesoflowers with complex shape and different dimensions (Figs. 5 and 6). These particles induce networks of defect lines meandering on their surfaces, typically along the ridges of spikes, with the longrange distortions of n(r) influenced by the spiky topography of particles but not exactly following it (Figs. 5 and 6). For large and strongly asymmetric particles gravitational torques and forces compete with their elastic counterparts, giving the origin to elastic monopoles ( Fig. 5), whereas smaller particles tend to induce higher order multipoles ( Fig. 6 and Supplementary Fig. 4). For monopole-like particles, the clockwise versus counterclockwise unidirectional rotation of n(r) away from n 0 corresponds to the opposite signs of monopoles, as confirmed by multipole expansion. The x-component of director, n x , which is an effective elastic charge density, has the same sign when plotted on a sphere encompassing the monopole-inducing mesoflowers, though its amplitude is nonuniform, as shown in the Fig. 5a-c, consistent with the leading monopole moment of the elastic charge distribution. Since the twist elastic constant of the LC is the smallest, the bend and splay distortions can also relax through equivalent distortions containing twist, so that the particles orientations can deviate away from n 0 not only in the plane containing g, the gravitational acceleration, but also in the plane orthogonal to g, as experimentally observed in planar LC cells (Figs. 1-4).
Details of the elastic-gravitational torque balance are revealed by simulating an asymmetric particle (similar to low-symmetry particles that tend to induce elastic monopoles in experiments) placed above a wall in xz plane with tangential boundary conditions (Fig. 5). The x-components (see inset in the lower left   Fig. 5d for the definition of the reference system) of the elastic T el x and gravity T g x torques vary with the angle θ between the particle axis (the line connecting the center of the particle core and the tip of one of the spikes) and n 0 . The torques are calculated about the center (see Supplementary Fig. 5), and exhibit different signs within the intervals θ ∈ [−180°, −127°] and θ ∈ [−90°, 0°]. The orientations θ = 0°(−180°) and θ ≃ −127°correspond to the minima of the elastic and gravity energies, respectively. For particles with masses M ( K=g, the elastic energy dominates and the particle adopts one of the two equilibrium orientations either θ eq ≲ 0°, or θ eq ≳ −180°, as shown in Fig. 5f for the case θ eq ≲ 0°. As the particle size increases, the gravity tilts the particle's axis further away from the purely elastic equilibrium θ = 0°towards negative values of θ eq , the equilibrium orientation behaves according to the upper branch of the curve in Fig. 5f. When starting from the second elastic energy minimum at θ = −180°, the gravity tilts the axis towards less negative θ upon increasing the particle size, towards the minimum of the gravitational energy (lower branch of θ eq curve in Fig. 5f). For sizes larger than R t , the particle would adopt its equilibrium orientation in the vicinity of the purely gravitational equilibrium θ ≃ −127°. The equilibrium orientation θ eq in Fig. 5f was obtained for h = 3.5R 0 , where R 0 is the radius of the spherical core of the mesoflower and h is the distance from the particle center to the wall ( Supplementary Fig. 6). As the weight of the particle increases the equilibrium particle-wall separation h eq decreases, as shown in the inset of Fig. 5e (the gravity acts in the directions towards the wall) and depends on balancing of gravitational and the repulsive elastic forces. The latter also depends on the strength of surface anchoring boundary conditions and is calculated here from the free energy profile versus h for the regime of strong boundary conditions corresponding to the experiments. The corresponding free energy results (Fig. 5e) demonstrate particle-wall repulsion, which have similar scaling when using one elastic constant approximation or not, though the elastic torque is almost insensitive to the variation of h for 3.5 ≤ h/R 0 ≤ 20, justifying the approximations used for calculating θ eq (Fig. 5f). The size-dependent levitation of mesoflowers at different h (inset of Fig. 5e) shows how suspensions of mesoflowers in LCs could be potentially used for separating particles of different dimensions, though the range of particle dimensions that can be effectively separated will depend on the density of particles relative to that of the LC host medium. Therefore, this separation method may be limited to highdensity particles, like the ones made from gold that we use in this study.  show color-coded diagrams of n x on the interpolation sphere surrounding the mesoflower when viewed from different directions. The coordinate system is set so that n 0 ||z. c, The xy perspective view of the n x and the color-coded scale of n x . d Elastic (red circles) and gravity (blue line) torques about the center of the particle core ( Supplementary Figs. 5 and 6) versus θ. The gravity acts in the negative y-direction; g is the gravitational acceleration. The mesoflower is placed above a wall with fixed planar boundary conditions along z and distance from the wall h = 3.5 R 0 . Balance of the two torques occurs at an equilibrium orientation θ eq , which is plotted in f versus ΔρR 0 3 g. e, Reduced Landau-de-Gennes free energy F LdG versus h at θ = −40°. Black solid circles correspond to K 11 = K 33 = 2K 22 = 7.8 pN and red solid squares to the oneconstant approximation K 11 ¼ K 33 ¼ K 22 ¼ K; we use T = 298 K and R 0 = 0.2 μm in the minimization. The inset in e shows the probability distribution / expðÀðF LdG þ E g Þ=k B TÞ of the particle-wall separation h for several values of R 0 indicated next to the curves. For the Landau-de Gennes free energy we use F LdG ¼ 5:659 h À1 KR 2 0 (parameters extracted from the black fitting curve in e), and E g is the particle gravitational energy ∝h; average elastic constant K ¼ 6:5 pN, and T = 298 K. g F LdG versus separation distance r c between a pair of like (black solid circles) and opposite (blue solid triangles) monopoles at h = 20R 0 θ 1 = θ 2 = 40°for like monopoles and θ 1 = 40°, θ 2 = −40°for opposite monopoles. Lines in e, g are fitting curves with the function f(x) ∝ 1/x. Director structures around pairs of like-charged h, i and oppositely-charged j, k monopole particles at r c = 5R 0 in h, j and r c = 10R 0 in i, k. The wall-monopole particle repulsive potential ∝ 1/h revealed in Fig. 5e agrees with the predictions of nematostatics 5,6,14 according to which elastic monopoles with opposite signs repel. Indeed, employing the electrostatic analogy, the interaction with the wall can be modeled using an image elastic monopole, yielding the ∝ 1/h potential. The potentials of mean force between two monopole particles with the same as well as opposite signs of monopole moments (Fig. 5g) also agree with nematostatics: 9 like (unlike) elastic monopoles attract (repel) each other with an effective potential ∝ 1/r c . Numerical modeling provides insights that formation of elastic monopoles by mesoflowers of gold is facilitated not only by their relatively large size >R t , but also by symmetry breaking often caused by asymmetric distribution and dimensions of individual spikes ( Supplementary Fig. 5) within the mesoflower particle (both effects enabled by the competition of gravitational forces and torques with their elastic counterparts, which stabilize elastic monopole configurations).
For particles smaller than R t , varying the number and positions of spikes within mesoflowers generates different multipole series of the director distortions ( Fig. 6 and Supplementary Fig. 4). The multipole expansion analysis 21 for the numerically simulated n(r) shows that certain particles induce stable or metastable structures with strongly pronounced multipole moments of different orders, including octupoles (Fig. 6a-c), hexadecapoles (Fig. 6d, e and Supplementary Fig. 4a-c) and even 64-poles (Fig. 6f, g and Supplementary Fig. 4d-f), with the other multipole moments orders of magnitude smaller. The numerically simulated n xdistributions on spheres encompassing these multipolar mesoflowers are consistent with the corresponding charge distributions of high leading order electrostatic multipoles (Fig. 6), though they also contain fine features dictated by detailed geometry of the mesoflowers (compare Fig. 6a, c).
Colloidal crystals and de-mixing of elastic multipoles. A multipolar approximation for an effective interaction between two colloidal particles distance r c apart and with the center-to-center vector forming an angle θ c with the far-field director n 0 reads: 17 where P l (x) is the Legendre polynomial of degree l, Q l characterize the strength of elastic multipole moments, and R eff is the characteristic length scale of the multipole (set in our case by the size of the colloidal particle). The analysis of Eq. (1) shows that a large variety of pair interaction patterns arises for multipoles of different orders (Fig. 7), which can lead to the formation of colloidal crystals and other structures arising from the competition between multipolar elastic and screened electrostatic repulsive interactions 19 . Because our mesoflower colloids allow for the realization of a large variety of multipoles of different orders, one can systematically explore how such colloidal interaction lead to self-organization depending on the order of the leading-order multipole. First, we consider a half-half binary mixture of distinct elastic monopoles confined at a plane coplanar to the far-field director. Contrary to electric charges, similar elastic monopoles attract and dissimilar repel 9,11 , as discussed above for pair interactions. To assess how this behavior impacts collective behavior of many such particles, we assume pairwise additive interaction potential and augment the monopole-monopole elastic interaction with a truncated repulsive Yukawa potential corresponding to the screened electrostatic interactions, which can be tuned by adding xz xz xz n 0 Fig. 6 High-order elastic mesoflower multipoles. a, d, f Pure elastic octupole, hexadecapole and 64-pole with corresponding n(r) in the cross-sectional planes, respectively, depicted on spheres using color-coded diagrams of n x . b, e, g Director structures around a mesoflower with dominant elastic octupole (b), hexadecapole (e), and 64-pole (g) contributions depicted with a perspective views on xz and yz planes. Coordinate system is defined so that n 0 || z. n (r) is shown using rods and defect lines are depicted as red tubes. c Color-coded diagram of n x in the xy, xz, and yz cross sections as well as at the interpolation spheres around the mesoflower with the dominant octupole contribution (b).

counterions through doping LCs with salt and other additives: 39
where r CO is the cutoff distance, A and κ are positive constants, with the latter characterizing electrostatic screening effects due to counterions within the LC 39 . We emphasize that rather strong electro-static repulsion at short particle separations is needed in order to avoid short-range steric, elastic and van der Waals interactions leading to aggregation due to the complex surface geometries of particles. Assuming such screened electro-static repulsions, we perform molecular dynamics simulations using open source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 40 . In this simulation as well as in all the subsequent simulations of self-organizations of higher order elastic multipoles, we neglect possible changes in multipole moments, or generation of additional multipoles, as the interparticle distances and orientations vary. We also neglect the short-range effects that may arise at small inter-particle distances and that possibly cannot be fully described within the multipole expansion approach that we consider. Figure 8 shows a series of snapshots along the system trajectory, where a perfect NaCl-like crystal order (see Fig. 8a) is assumed as the initial condition. In the course of time, the system of oppositely charged elastic monopoles undergo a spatial segregation, Fig. 8f, which is the direct consequence of the fact that like (opposite) monopoles attract (repel) (see also the Supplementary Movie 1). Self-assembly of two-dimensional nematic colloids with the dominant dipolar or quadrupolar elastic interactions has been previously reported 22,41,42 and is consistent with the behavior of colloidal dipoles and quadrupoles formed by gold mesoflowers. Driven by the symmetries of the underlying multipole potentials, the dipole colloids tend to assemble head-to-tail in chains aligned along the far-field director. The chains then assemble in a twodimensional crystalline structure with an anti-ferromagnetic-like alignment of the neighboring chains 22,41 . Quadrupolar colloids reveal two-dimensional crystalline assembly with a rhomboidal unit lattice cell where the colloid-colloid bonds are aligned along the attractive directions of the underlying quadrupole potential 42 . The sectors of attraction and repulsion are oriented with respect to n 0 , which defines the long-range orientations for the crystallographic axes of the ensuing colloidal crystals formed by elastic dipoles and quadrupoles, thus precluding the formation of grain boundaries between them. To the best of our knowledge, crystallization of nematic colloids driven by higher order multipole potentials have not been investigated.
Experimental results [20][21][22]41,42 suggest that the symmetry of the eventual nematic colloids crystals are mainly determined by the distribution of the attractive directions of the corresponding multipole potential. The elastic octupole potential has six directions of attraction forming roughly 60°between each other (Fig. 7a). It is therefore expected that octupole nematic colloids crystalize in two dimensions into a hexagonal lattice, indeed consistent with our LAMMPS-based simulations. The unique alignment of the hexagonal lattice of these elastic-octupolebased crystals with respect to n 0 in this nematic colloidal system precludes the formation of grain boundaries and, because of this, could potentially be used for the formation of large crystal lattices. The case of hexadecapole colloids is fundamentally more interesting, as in this case the potential has eight directions of attractions (Fig. 7b), which permits two types of quadratic lattice arrangements of hexadecapole colloidal particles, with the corresponding unit lattice cells shown in Fig. 9a, b. In order to verify this hypothesis, we performed molecular dynamic simulations (using LAMMPS) of twodimensional hexadecapole colloids, regularized with repulsive Yukawa potential as defined by Eq. (2). We used constant temperature and constant pressure (NPT) ensemble. Initially, the system was subject to a strong isotropic pressure (Fig. 9c, f). In this initial state the colloids formed a glassy structure, with enhanced local hexagonal ordering revealed by the absolute value |q 6 | of the local hexatic bond order parameter (Fig. 9c) and suppressed |q 4 | of the quartic bond order parameter (Fig. 9f). In the course of the simulations, the external isotropic pressure was gradually decreased, which resulted in the spontaneous formation of two types of quadratic colloidal lattices and associated grain boundaries between them (Fig. 9d, e, g, h, and Supplementary Movie 2).
Although the hexadecapole elastic interactions tend to form two-fold degenerate ground state whose lattice units are depicted in Fig. 9a Colored solid lines represent the interaction potential U el ∝ P l+m (cos(θ c )) between two multipoles of the same order, e.g. for l = m. a Interaction potentials between two dipoles, octupoles, and 32-poles (l = 1, 3, and 5, respectively). b Interaction potentials between two quadrupoles, hexadecapoles, and 64-poles (l = 2, 4, and 6, respectively). The double arrow indicates the far-field director n 0 with respect to which θ c is measured. Dashed circles intersecting the potential plots represent equipotential lines of U(θ c ) = 0 for each plot. Multipole pairs mutually repel in the regions where the lines extend beyond the dashed circle, e.g. U el > 0, and attract where the lines lie within the dashed circles, e.g., U el < 0. Such radial directions of interaction are marked by the black arrows. The plots are not presented to scale.
the efficiency of two-dimensional packing at high colloidal concentrations. Colloidal particles are electrostatically charged and mutually repel with screened Coulomb electrostatic potential, where the range and the strength of this potential can be tuned 39 , which provides additional knobs to control the self-assembly of nematic colloids. We have performed illustrative molecular dynamic simulations of self-assembly behavior of hexadecapole colloids with strong, weakly screened repulsive Yukawa potential. The ensuing self-assembled colloidal configurations are presented on Fig. 10 and Supplementary Fig. 7 (see also Supplementary Movies 3-6). At high concentrations the system reveals glassy (Fig. 10a, d) configurations, while at low concentrations particles assemble into repulsion-induced sparse hexagonal lattices (not shown). Surprisingly, at intermediate values of particle concentrations and depending on the strength of the hexadecapole elastic moment, the system can crystalize into a range of rhomboidal lattices, with some examples shown in Fig. 10c, k and Supplementary Fig. 7c, k, and the corresponding unit lattices in Fig. 10q. Importantly, the symmetry of these lattices differs from the ones of hexadecapolar ground states shown in Fig. 9a, b. Additionally, at some number densities we observe coexisting rhomboidal and hexagonal lattices, see Fig. 10c, k and Supplementary Fig. 7c, k. Many interesting energetically comparable colloidal organizations emerge as a result of dense packing and elastic hexadecapolar interactions favoring different lattices. This unusual crystallization behavior is revealed by analyzing bond orientational order parameters within selfassembled crystallites separated by grain boundaries. These examples illustrate the ability of hexadecapolar nematic colloids to self-assemble into a range of two-dimensional low-symmetry crystal structures, with properties that can be designed through controlling the particle shape, concentration as well as internal properties of the nematic host.

Discussion
In this work, we have developed a nematic colloidal system of mesoflowers where complex and diverse shapes of mesostructures induce elastic multipoles of different order, from monopoles to high-order multipoles like 64-poles. By varying the basic building parts of these particles, their colloidal behavior can be effectively controlled. As the shape and size of particles vary, they behave like elastic monopoles that (differently from their electrostatic counterparts) tend to de-mix into the domains of like-charged elastic monopoles, or as high-order elastic multipoles (e.g., hexadecapoles) forming unusual crystals that can be tuned by packing density and strength of electrostatic repulsions. While here we focused only on the analysis of two-dimensional crystal lattices for illustrative purposes, even more complex behavior is expected in three dimensions, where even elastic quadrupoles can form triclinic crystal lattices 39 . As the symmetries of elastic potentials of high-order multipoles, which form the basis of crystals, become intrinsically incompatible with crystalline lattices, various forms of frustration can arise and lead to complex crystallization behavior that can be potentially tuned through changing concentration of counterions and surface charging. In addition to crystals, quasi-crystals and various plastic crystals (with positional order but lacking orientational ordering) could possibly also arise due to diverse multipolar nature of our mesoflower colloids, whereas polydisperse systems are expected to form disordered glassy states. Our findings suggest that these mesostructures may serve as a designer toolkit for engineering pre-defined colloidal interaction and self-assembly. Access to these structures in large quantities through chemical synthesis may facilitate developments of new composite materials fabricated using colloidal self-assembly. On the other hand, the sharp geometric features of the gold mesoflowers within these colloidal particle assemblies may be of interest for plasmonic enhancement and related applications.

Methods
Sample preparation. Gold mesoflowers are synthesized following an aqueous seed-mediated growth procedure 34  accompanied with gentle mixing. The mixture was kept at 80°C for 1 h and cooled down naturally before taking it out for centrifugation at 4000 rpm for 5 min. The residue was collected and washed in deionized water 3 times and then dried in an oven. Liquid crystal suspension of mesoflowers was obtained by sonicating these dried residues within the LC host, after which the suspension was infiltrated between glass substrates with gap defined by monodisperse glass spacers of 10-20 μm in diameter. Prior to being used in cell preparation, the glass substrates were spin-coated with polyimide (PI2555, from HD Microsystem) and rubbed unidirectionally to define the bulk LC alignment. The cell edges were sealed using fastsetting epoxy glue.
Optical video microscopy and laser trapping. An inverted optical microscope (IX81, Olympus) with a charge-coupled device (CCD) camera (Flea, PointGrey) and a holographic laser trapping system operating at 1064 nm was used to take polarizing optical micrographs and videos of mesoflowers in the LC, to probe the stable and metastable director configurations around the particles, as well as to set the initial conditions for studying their pair interactions. The laser was turned off during particle diffusion and interaction when videos were being recorded. To assure good optical imaging resolution and robust laser manipulations, we used an oil-immersion objective lens (UPlanSApo 100×, Olympus) with high numerical aperture NA = 1.40 for both imaging and laser trapping. In order to analyze colloidal interactions and diffusion in directions along and perpendicular to n 0 in the plane of the LC cell, optical videos were analyzed using an image-processing software (ImageJ, freeware from the National Institute of Health) to extract positions of the mesoflowers on a frame-by-frame basis. MSDs that characterize single particle's Brownian motion were calculated by analyzing particle displacements along directions perpendicular and parallel to n 0 (denoted Δr x and Δr z , respectively) based on the particle's positions in each frame of the video. To probe the angular dependence of particle diffusion, a rotation operation was applied Δr 0 where β is the angle of rotation with respect to n 0 and Δr 0 x Δr 0 z is the displacements in the rotated coordinate frames. MSD was then calculated as <Δr x > 2 for β varying from 0°to 360°. Diffusion constants D were extracted from MSDs using the relation <Δr x > 2 = 2Dτ where τ = 1/f and f = 15 Hz is the frame rate of the experimental videos. From Stokes-Einstein relation, the drag coefficients of diffusing particles were calculated as c μ = k B T/D μ (μ = x, z) and the corresponding drag forces as F d = cv, where k B is the Boltzmann constant, T = 300 K is the room temperature and v is the linear speed of the particle. Because of the highly overdamped nature of our colloidal system, the elastic interaction forces could be then found from their balance with F d . The elastic interaction potential was then calculated by integrating the elastic force over distance.
Simulation of director structures and elastic interactions. The generation of elastic monopole and higher-order multipoles in the n(r) field by mesoflowers is also approached computationally, via numerical minimization of the Landau-de Gennes free energy 26 . The Landau-de Gennes free energy functional depends on the tensor order parameter field Q(r) and its spatial derivatives. The functional Red lines indicate the attractive directions from the underlying hexadecapole interaction potential (Fig. 7b). Black lines represent the Cartesian axes with n 0 parallel to the vertical axis. c-h Snapshots of the process of grain boundary formation. Spheres in the top (bottom) rows are colored according to the absolute value of the local quartic order parameter q 4 (j) (hexatic order parameter q 6 (j)); color scale decoding the value for both of these order parameters within 0 to 1 is shown as an inset on the right side of (b). Panels in the same column are taken at the same time. Parameters used are Q 4 ¼ 7 10 À5 ; r CO combines terms which account for variable degree of nematic order, nematic elasticity and biaxiality (at the cores of topological defects), as well as surface anchoring. As such, the Landau-de Gennes approach yields theoretical characterization of n(r) and local changes in the degree of nematic ordering that correspond to global or local minima of the free energy functional 26 . Spiky particles are constructed as a union of a spherical core and a given number of spikes distributed over the surface of the core at predefined locations and orientations (Supplementary Fig. 6). Minimization of the free energy is performed numerically for finite homeotropic surface anchoring at the particle surfaces, and by using variable threedimensional tetrahedral grids 43 . This minimization yields stable or metastable n(r) structures around the particles, which are then compared to the experimentally reconstructed counterparts 16,21 . Nematic configurations around mesoflowers are obtained via numerical minimization of the phenomenological Landau-de Gennes free energy functional 26 where Q ij = Q ji (i, j = 1,..,3) is a traceless tensor order parameter and summation over repeated indices is assumed. In Eq. (4), the parameter a (unlike the constants b and c) is assumed to depend linearly on temperature T: a(T) = a 0 (T − T*), where a 0 is a material dependent constant, and T* is the supercooling temperature of the isotropic phase. Phenomenological parameters L 1 and L 2 are related (via an uniaxial Ansatz for Q ij ) to the Frank-Oseen elastic constants. We describe finite homeotropic anchoring, with the strength coefficient W, at the surface of the particles by using f s ðQ ij Kronecker delta symbol, and ν is the unit outward normal vector to the particle surface 44 ; is the value of the scalar order parameter in the nematic phase, which is thermodynamically favored for 24ac/b 2 < 1. Minimization of the free energy Eq. (4) is then performed numerically by employing adaptive mesh finite elements method as described in more details in ref. 43. This minimization yields theoretical characterization of n(r) and local changes in the degree of order that correspond to global or local minima of the free energy 26 .
Computational geometry of colloidal mesoflowers. We exploit the Open Source Gnu Triangulated Surface (GTS) library 50 to create triangulated surfaces of the mesoflowers. A spiky mesoflower particle is constructed as a union of a spherical core with the radius R 0 and a given number of spike particles ( Supplementary  Fig. 6) decorating the spherical core at predifined locations and orientations. Each spike is chosen as a generalized cylinder with a 5-fold symmetry axis, consistent with the experimental electron microscopy images. The lateral surface of each spike whose symmetry axis coincides with the z-axis of a Carthesian coordinates system (x i , y i , z i ), i = 1, ..., 5, is parametrized as follows: x i s; u ð Þ ¼ A u ð Þ 1 þ 0:3 cos 5s ð Þ ½ cos s ð Þ; y i s; u ð Þ ¼ A u ð Þ 1 þ 0:3 cos 5s ð Þ ½ sin s ð Þ; z i s; u ð Þ ¼ Hu; where 0 ≤ s < 2π, and 0 ≤ u ≤ 1 are parameters; A(u) = A 0 + (A 0 − A 1 )u accounts for the variation of the spike "radius" along the spike symmetry axis with A 0 > A 1 , and H is the height of the spike along its symmetry axis. The triangulated surface of the spike can be manipulated using the functions implemented in the GTS library, including surface translation, rotation, or rescaling. By using these functions together with the GTS function which merges any two surfaces, we generated the spiky mesoflower particles with arbitrary number of spikes placed at predefind relative positions and with predefined orientations. The triangulation of the nematic domain was then performed by using a quality tetrahedral mesh generator 51 , which supports adaptive mesh refinement. Finally, the discretized functional (5) is minimized numerically using INRIA's M1QN3 optimization routine 52 , which implements a limited-memory quasi-Newton method 53 . Molecular dynamics simulations. The self-assembly and formation of colloidal crystals are computer-simulated using an open source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 40 . In the simulations, we use KR eff as the unit energy, and R eff as the unit length. The total number of particles simulated is N~1000, with the system size L x × L y ranging from 26 × 26 to 37 × 37 R 2 eff . The potential energy of the system is assumed to be pair-wise addative, combining both the elastic energy and Yukawa potential: U = U el + U Y . For the computer simulations on the segregation of monopoles (Fig. 8), the system size is 34 × 34 R 2 eff and total number of particles N = 1024 with half of them having positive elastic charges and the other half having negative elastic charges. The interaction potential between particles is taken as where the elastic monopole charges Q a 1 , Q b 1 = ±0.1, and we set The simulation is performed using NVT ensemble.
Temperature T is fixed by applying a Langevin thermostat with k B T KR eff ¼ 10 À3 : For the computer simulations involving hexadecapoles (Figs. 9 and 10, and Supplementary Fig. 7), N = 855 and the effective pair interaction potential is The elastic hexadecapolar moment Q 4 , and the relative Yukawa potential strength A KR 2 eff are varied to obtain different colloidal assembly structures. The simulations are performed in NPT ensemble, using a Nose/Hoover temperature thermostat with k B T KR eff ¼ 10 À3 ; and Nose/Hoover pressure barostat. In order to achieve a spontaneous crystallization an isotropic external pressure is gradually released in the course of the simulations leading to the expansion of the system. The local quartic (hexatic) orientational order parameter q 4 (j) (q 6 (j)), where j is the index of the particle under consideration, is defined as q n j ð Þ ¼ 1 n Σ n k¼1 expðinθ jk Þ; θ jk is the angle between the particles center-to-center vector r jk and x axis, and the sum runs over n nearest neighbors of particle j.

Data availability
The data that support the findings of this study are available from the corresponding author on request.