Introduction

"Colloidal molecules” (CMs) refer to clusters of colloidal particles exhibiting specific local angular arrangements that replicate the geometries of molecules1,2,3,4,5,6. These inherent and stable angular orders significantly enhance CMs’ self-assembly capacities. CMs stand out as essential and distinctive building blocks, facilitating the creation of predetermined hierarchical superstructures7,8,9.

Previous studies have utilised solutions to the Thomson problem, which addresses the arrangement of N electrons on a sphere, to achieve specific angular orders and symmetry in CMs10,11,12. Building on this, significant efforts have focused on synthesizing CMs with fixed satellites and precisely controlled bond angles13,14,15. These fixed CMs are crucial for assembling complex crystalline structures that require directional bonding16,17,18.

Moreover, dynamic CMs, featuring configurational flexibility and fluctuations, have garnered increasing attention from physicists and materials scientists19,20,21. Unlike CMs with fixed satellites, dynamic CMs more closely resemble biological structures, such as macromolecules and proteins, where significant configurational fluctuations occur. This flexibility enables them to bind more effectively with irregularly arranged or fluctuating binding sites, leading to more efficient reactions22,23,24. Designing stable and flexible CMs has been achieved through patchy colloids with directional bonding and the application of mobile DNA linkers25,26,27,28,29. Recent studies suggest that dynamic CMs tune their angular orders through a delicate interplay between entropy and electrostatic interactions30. Managing these configurational fluctuations, as seen in biological systems, opens new possibilities for the creation of amorphous and polymorphic materials31,32.

Despite progress in synthesising and assembling CMs, our understanding of their physical characteristics still needs to be improved. For instance, CMs are usually considered symmetric, but how do we adjust the asymmetry of dynamic CMs? On the kinetic side, achieving a CM with a specific angular order is challenging due to multiple and sequential reaction processes among colloidal particles30. However, experimental observation, characterisation, and tuning of CMs’ dynamic structures pose challenges due to both the fast internal and global translational-rotational motion of CMs in the solvent19. These challenges impede a clear understanding of the dynamic features of CMs.

Here, we design a colloidal-emulsion model system that forms micrometre-sized dynamic CMs with tunable electrostatic interactions33,34. By adjusting the ionic strength, we can effectively control the range and the magnitude of electrostatic repulsion between charged satellite particles. This system enables three-dimensional imaging of CMs’ dynamic structures and formation kinetics under confocal microscopy. Through a combination of experiments and numerical simulations, we discover the intrinsic asymmetry in CMs’ dynamic structures. Their angular symmetry evolves continuously from a liquid-like configuration, distinct from a sharp first-order transition. We further propose a method to guide CMs’ ordering kinetics toward the desired target structure by dynamically adjusting the ionic strength in the solvent. This approach offers a solution to the complexities inherent in the formation process.

Results and Discussion

Model system design

Our model CMs consist of satellite PMMA colloidal particles (negatively charged, diameter σs ≈ 1.2 μm, suspended in a density and refrective-index matching cyclohexyl bromide (CHB)-cis-decalin solvent) and central emulsion droplets (negatively charged, diameter σc ranging between 0.8 μm and 1.5 μm, water-glycerol droplet in a CHB-decalin solvent, with sodium dodecyl sulfate (SDS) serving as a surfactant), as illustrated in Fig. 1a. Although gravity could potentially influence the configuration of CMs, we neglect its effect in our model due to the minimal density difference among the PMMA particles, the CHB-decalin solvent, and the emulsions (see Methods for a detailed discussion). The repulsion between satellite particles is modelled by a hard-core Yukawa pair potential, \({U}_{pp}(r)=\frac{{Q}^{2}}{{\sigma }_{{{{\rm{s}}}}}{(1+\frac{\kappa {\sigma }_{{{{\rm{s}}}}}}{2})}^{2}{\epsilon }_{s}}\exp [-\kappa {\sigma }_{{{{\rm{s}}}}}(r/{\sigma }_{{{{\rm{s}}}}}-1)]/(r/{\sigma }_{{{{\rm{s}}}}})\) for r > σs, with κ−1 representing the screening length, and \(\frac{{Q}^{2}}{{\sigma }_{{{{\rm{s}}}}}{(1+\frac{\kappa {\sigma }_{{{{\rm{s}}}}}}{2})}^{2}{\epsilon }_{s}}\), denoted as Ar, being a surface potential parameter. Here, Q denotes the surface charge of the colloidal particle, and ϵs represents the dielectric constant of the solvent. The interaction between particles and the droplet involves a long-range repulsion and a short-range attraction due to charge-image charge attraction and hydrophobic forces33,34. Particles can move freely on the emulsion droplet surface. CMs are formed by mixing colloidal suspensions with emulsions (see Methods). Tetrabutylammonium bromide (TBAB) salt is introduced into the solvent to modify κ−1 and Q. We primarily conducted experiments at two salt concentrations: 1 μM and 5 μM. As the TBAB concentration increases, the surface charge Q rises, the Debye screening length κ−1 decreases, and the surface potential parameter Ar slightly decreases (see Methods for calculation details). By adjusting the stoichiometry ratio to approximately 200:1 between particles and emulsion droplets, we effectively inhibit the formation of large aggregations.

Fig. 1: Dynamic structures of colloidal molecules observed in experiments.
figure 1

a Schematic representation of the experimental system featuring tunable electrostatic interactions. b Confocal microscopy images and reconstructed three-dimensional configurations of N-CMs with varying σc (5 μM TBAB). The scanning speed along Z direction is 7 μm/s (250 nm per slice). c Illustration of the temporal fluctuations of a single CM in experiments, portrayed through real-time trajectories. See also Movies S1-5. d Configurational fluctuations of N-CMs (bottom raw, 1 μM TBAB) in comparison with the corresponding standard Thomson’s solution (top row) for each N. Side and top views are presented, with dots indicating satellite positions after rotation for maximal superimposition. We classify 5-CMs based on their proximity to either a triangular bipyramid or a pyramid configuration.

Full size image

Angular order and dynamic fluctuations of CMs

Our experiments yielded dynamic colloidal molecules (N-CMs) with a satellite particle number N ranging from 2 to 6 (Fig. 1b). At 5 μM TBAB (κσs ≈ 4.5), observed configurations include 2-CMs in a V shape, 3-CMs forming a triangle, 4-CMs adopting a tetrahedral structure, 5-CMs exhibiting both triangular bipyramids and pyramid shapes, and 6-CMs forming an octahedral configuration. Importantly, these N-CMs exhibit stability and flexibility, displaying configurational fluctuations around an equilibrium state characterised by angular orders. Real-time trajectories of single CMs obtained in experiments (Fig. 1c and Movies S1-5) reveal satellite particles undergoing temporal fluctuations and unrestricted movement around the core.

Utilising the Kabsch algorithm for optimal satellite alignment through rotation35,36 (see Methods), we experimentally obtained dynamic structures, including both the average and the fluctuation, by examining multiple CMs with similar core sizes. We emphasize that the system reached equilibrium by the end of our experiment, as the temporal distribution of a single CM was statistically similar to the ensemble distribution obtained from multiple CMs. Interestingly, the average configurations of satellite particles display asymmetry, deviating from the standard symmetric configurations predicted by solutions to the Thomson problem. For instance, the angles in V-shaped configurations may not be the ideal 180 degrees, and triangular configurations are not always equilateral. Intriguingly, despite the triangular bipyramid being the ideal configuration, we also observed the presence of pyramid configurations for 5-CMs.

Figure 1 d shows the distribution of CMs’ configurations with coordination numbers from 2 to 6, derived from snapshots of over 50 independent CMs at equilibrium (1 μM TBAB, κσs ≈ 2.0). For 2-CMs, the angle distribution of the V shape centres around 118 degrees, with a fluctuation of approximately 30 degrees (see also Movie S1). For 3-CMs, the top view reveals an acute triangle as the average configuration instead of an equilateral triangle (Movie S2). For 4-CMs, the tetrahedron is irregular, with two faces nearly perpendicular to each other (Movie S3). For 5-CMs, we classify the entire population into two categories based on their proximity to either the triangular bipyramid or the pyramid (see Methods). Approximately 20% of the observed configurations are closer to the triangular bipyramid, while the majority are closer to the pyramid (Movie S4). The average configuration within each category closely resembles their respective reference. For 6-CMs, four satellites in the middle remain roughly on the same plane below the equator, yet forming distorted pyramids in both the upper and bottom halves (Movie S5).

Our experimental findings apply universally to CMs assembled by charged colloids, e.g., colloidal particles grafted by various types of functional groups, which act as satellites in aqueous solutions or organic solutions5,6,29. In these systems, key experimental parameters, such as ion concentration in the solvent, zeta potential (surface potential), and the size of colloidal particles, play a critical role in quantifying the charged interactions. To systematically and comprehensively understand dynamic CMs’ asymmetry and to assess the universality, we perform canonical Monte Carlo simulations using κσs, Ar, and σc/σs as control parameters (see Methods).

Interestingly, the dynamic CMs’ asymmetry appears to be universal. By manipulating κσs and Ar, we can create CMs with specific dynamic structures. Figure 2a provides guidance for designing asymmetric 2-CMs with an average bond angle, 〈θ〉, ranging from 100∘ to 150∘ (σc/σs = 0.95). Similarly, for 3-CMs, we can adjust their asymmetry, quantified by the solid angle, Ω (the area of the segment of a unit sphere, Ω = 2π in steradian for symmetric 3-CMs), as illustrated in Fig. 2b. Note that the standard deviation of θ (2-CMs) and Ω (3-CMs) decreases consistently with a reduction in CMs’ asymmetry (Figs. S1a and b).

Fig. 2: Structure asymmetry and angular orders of N-CMs as a function of κσs and σc.
figure 2

a Illustration of the asymmetric structure of 2-CMs through the bond angle θ. b Illustration of the asymmetric structure of 3-CMs using the solid angle Ω. c Illustration of the asymmetric structure of 4-CMs represented by δ (the deviation between the average and ideal configurations). Circular symbols represent the average configurations of over 50 independent CMs (σc ≈ σs), obtained from experiments conducted under two different TBAB concentrations: 1 μM (left) and 5 μM (right). The confocal snapshots display representative configurations that closely resemble these average. Standard deviations of θ (2-CMs), Ω (3-CMs), and the magnitude of configurational fluctuation around the average configuration Δ (4- and 6-CMs) are presented in Figs.S1d and e. d and e, Angular order development for 4- and 5-CMs depicted by q3 and q4. Square symbols with error bars represent the average and standard deviation (for σc = 0.95 σs and Ar = 10 kBT). Subpanels illustrate the \(| \Delta \overrightarrow{q}|\)-\({\sigma }_{\overrightarrow{q}}\) master curve for varying Ar (5 to 75 kBT), σc (0.95 to 1.45 σs, incremented by 0.1 σs), and κσs (2 to 30). Note the bifurcation of 5-CMs labelled in blue symbols (κσs = 2, towards triangular bipyramid) and green symbols (κσs = 10, towards pyramid) in the inset of e. The standard triangular bipyramid is used as the reference configuration when calculating \(| \Delta \overrightarrow{q}|\) for 5-CMs.

Full size image

For N ≥ 4 CMs, we quantify the asymmetric and fluctuated structure profile by using two parameters, δ, representing the deviation between the average (asymmetric) and standard (symmetric) configuration, and Δ, indicating the magnitude of fluctuations around the average configuration (see Methods for definitions). Figure 1c and Figs. S1c-e display the diagram modifying δ and Δ for 4-CMs and 6-CMs, respectively. Notably, all CMs tend to converge to solutions of Thomson problems as κσs and σc decrease and Ar increases. This behaviour is expected as inter-satellite repulsive interactions become sufficiently long-range to stabilise regular configurations corresponding to solutions to the Thomson problem.

We compare the experimental configurations at different TBAB concentrations with the corresponding simulations (see Methods for technical details on the comparison). As the TBAB concentration decreases from 5 μM (κσs = 4.5, Ar = 4.7 kBT) to 1 μM (κσs = 1.9, Ar = 6.5 kBT), the configurations of 2, 3, and 4-CMs observed in experiments become more regular, though they still exhibit significant asymmetry. This trend is in full agreement with the simulation predictions. Notably, unlike in simulations, where κσ and Ar can be varied independently, in experiments both the Debye screening length and the particle surface charge are functions of the salt concentration. The tunability of the asymmetry strongly depends on the physicochemical properties of the system, including the particle material, solvent choice, and type of salt solute.

Nevertheless, the optimal ranges for asymmetric CM assembly can be universally adapted to other experimental systems. From our simulations, we have identified that to achieve asymmetric CMs, the range of κσs should ideally be between 1 and 20, while the surface potential should range from 10 to 100 kBT. For our system with 1-μm diameter particles, the optimal ζ potential is between 10 and 20 mV, with the TBAB concentration around 1 μM. In an aqueous solvent (where the relative dielectric constant ϵr ≈ 80), given that Ar ∝ ϵrσsζ2 and \(\kappa \propto \sqrt{{n}_{{{{\rm{ion}}}}}/{\epsilon }_{{\mbox{r}}}}\)37, achieving similar conditions would require adjustments such as using a smaller particle size (σs between 200 and 500 nm), lower ζ potential (5 to 10 mV), and an increase in ionic strength (nion between 10 and 100 μM).

The emergence of asymmetric structures can be elucidated through the interplay of entropy and potential energy. Take, for example, the asymmetry in 2-CMs. The potential energy minimum occurs when one particle is positioned with a relative angle to the other of θ = 180∘ for two freely moving and repulsive particles constrained on a spherical surface. However, the entropically favoured relative angle is θ = 90∘. Consequently, the dynamic structures inherently exhibit asymmetry. This simple yet universal principle, previously undisclosed, provides guidance for the design of asymmetric clusters at micro and nano scales.

We further quantified CMs’ angular orders using the Steinhardt bond orientational order parameters, ql (see Methods for definitions). The configuration of 4-CMs gradually converges to a standard tetrahedron with a decrease in κσs (at σc/σs = 0.95), as illustrated by the q3-q4 distributions in Fig. 2d. To characterise the convergence, we defined a vector \(\overrightarrow{q}\)=(q3, q4), with q3 and q4 as its two components, and calculated the deviation from the ideal configuration (\(| \Delta \overrightarrow{q}|\)) and the standard deviation of \(\overrightarrow{q}\) (\({\sigma }_{\overrightarrow{q}}\)). Interestingly, after a stage where the configurations are “liquid-like” (plateau of \({\sigma }_{\overrightarrow{q}}\)), we observed a rapid, continuous improvement of angular orders (\(| \Delta \overrightarrow{q}| \to 0\)), as depicted by the \(| \Delta \overrightarrow{q}|\)-\({\sigma }_{\overrightarrow{q}}\) master curve in the inset of Fig. 2d (for σc/σs ∈ [0.85, 1.45]).

The rapid and continuous enhancement of angular order is marked by a swift convergence towards the ideal configuration. Notably, we emphasise the resemblance of these configurations to those in hydrated ions that have recently been acknowledged for their dynamic structures and fluctuations around perfect angular symmetries38. This similarity hints at the potential for real-time exploration of complex ion-specific assembly behaviours. Similar behaviour is also evident for 3-CMs and 6-CMs (see Fig. S2a-b).

5-CMs exhibit intriguing polymorphic behaviour, distinguishing them from other CMs, as illustrated by the mixing of q3-q4 distributions in Fig. 2e. Close to the “liquid-like stage”, the average configuration resembles an irregular pyramid with significant distortions. Angular order development towards the standard pyramid becomes evident at large κσs, e.g., κσs = 10 (see Fig. S3a), while ordering towards the triangular pyramid is favoured at low κσs, e.g., κσs = 5 (see Fig. S3b). The bifurcation pathways are more clearly visible in the \(| \Delta \overrightarrow{q}|\)-\({\sigma }_{\overrightarrow{q}}\) master curve shown in the inset of Fig. 2e. Note that we use standard triangular bipyramid as the reference configuration when calculating \(| \Delta \overrightarrow{q}|\) for 5-CMs. A systematic calculation of the free-energy difference, ΔF, defined as (FT − FP)/kBT, between the two polymorphs supports the observation that the polymorphic feature is present over a wide range of σc when ΔF ~ 0, as depicted in Fig. S3c.

The kinetics and pathway of CMs’ formation

The formation kinetics of CMs in our system involves multiple and sequential reactions among the core droplet, satellite particles, and free particles, progressing from 1-CM to N-CMs, as schematically illustrated in Fig. 3a. The free energy barrier associated with the addition of satellite particles increases due to particle-CM electrostatic repulsions. In contrast, detachment is irreversible owing to the strong core-satellite attraction. We defined the probability of forming N-CMs at time t as PN(t). The kinetic equation can be written as20:

$$\frac{d{P}_{N}(t)}{dt}={k}_{N-1}{P}_{N-1}(t)-{k}_{N}{P}_{N}(t)$$
(1)

Here, kN represents the transition rate for an N-CM to absorb another satellite: \({k}_{N}={k}_{0}\exp \left(-\frac{\Delta {G}_{N}}{{k}_{{{{\rm{B}}}}}T}\right)\), where k0 is the effective collision frequency between the core and free particles, and ΔGN is the free energy barrier for a free particle to attach to a N-CM. Initially, PN(0) = 1 for N = 0. We assume ΔGN = NΔg, where Δg is a coefficient depending on σc, κσs and Ar. PN(t) can then be solved iteratively:

$${P}_{N}(t)=\left\{\begin{array}{l}\exp (-{k}_{N}t),\quad N=0 \hfill \\ {\sum }_{i=0}^{N}\frac{\mathop{\prod }_{m=0}^{N-1}(-{k}_{m})}{\mathop{\prod }_{n=0,n\ne i}^{N}({k}_{i}-{k}_{n})}\exp (-{k}_{i}t),\quad N \, > \, 0\end{array}\right.$$
(2)
Fig. 3: The kinetics of N-CMs’ growth.
figure 3

a Schematic illustrating the sequential reaction process with respect to N. The forward transition barrier increases linearly with N, while the backward rate is nearly zero. b Evolution of ϕN(t) over time t. Molecular dynamic simulations (circles) align well with the solutions of PN(t). c, The most probable coordination number, Nmax, versus time: dashed line (Δg = 0.1kBT), solid line (Δg = 0.4kBT), and dotted line (Δg = 0.8kBT). Optimisation of N-CMs production within a period is achievable with an appropriate Δg. d, e Distribution of N (indicated by the colour map and symbol sizes) with respect to σc/σs in experiments (in d) and simulations (in e, 35 τc).

Full size image

We also conducted molecular dynamic simulations using HOOMD39. For each core size, we run 200 replicas to investigate the temporal evolution of the fraction of N-CMs within the entire CM population, denoted as ϕN(t) (see Methods). The solutions for typical N < 5 in PN(t) (Δg = 0.4kBT) demonstrate good agreement with simulations, as illustrated in Fig. 3b. Note that κσs and Ar approximately follow the contour lines of δ in Fig. 2c for constant Δg. Deviation occurs for N ≥ 5 when the assumption breaks down due to the core’s accommodation capacity limit. The solutions suggest that optimising the production of N-CMs within a fixed period requires a specific Δg value. For example, setting Δg = 0.4 kBT optimizes 6-CMs within 35 τc (τc ≡ 1/k0), as depicted in Fig. 3c. At various σc/σs, the experimentally observed distributions of N-CMs in Fig. 3d approximately align with simulations in Fig. 3e. The ϕN maxima in experiments are better than in simulations due to the much longer waiting time in experiments. Changing σc can also perform the optimisations, consistent with a previous study20,40.

A kinetic tool box for designing CMs

Efficient production optimisation of N-CMs requires appropriate values of σc and κσs. However, we have also shown the need for sufficiently small σc or κσs to stabilise regular CMs. Thus, efficiently producing regular CMs with maximum yield poses a challenge in pathway design. Inspired by our previous experimental study where colloidal particles crystallize within oil-in-water emulsions34, we observed a decrease in salt concentration within the oil droplets due to ion diffusion into the surrounding water bath. This decrease is substantial in the first few hours after sample preparation, gradually slowing over time. This observation prompted us to explore the possibility that in-situ control over ionic strength could be a means modulate CM assembly kinetics.

Here, we propose a strategy and offer potential design examples. Our strategy involves a ramping protocol, which reserves κσs at a relatively high value with a waiting time optimizing irregular (N − 1)-CMs (e.g., 10 - 30 τc at κσs = 25 for 4-CMs). Subsequently, κσs gradually lowers to a low value (e.g., κσs = 5), ensuring that N-CM possess high angular order and the transition rate from 4-CMs to 5-CMs is not excessively low.

The waiting time and ramping rate are crucial parameters influencing the distribution of N-CMs. In Fig. 4a, b, we present simulation results for three different ramping protocols, analysing the corresponding deviation δ and N over time. In the case of a fast ramping rate and short waiting time (left panels), Nmax remains at 4 for over 20 τc with a yield over 50%. For the same starting point but a lower ramping rate (middle panels), the growth becomes faster, resulting in almost equal yields of 4- and 5-CMs at the end of the simulation. Finally, with a proper waiting time and a slower ramping rate (right panels), the growth is even faster, and 5-CMs dominate at the end of the simulation. Solutions of PN(t) also accurately capture these trends and predict the most probable N over time.

Fig. 4: The dependence of N-CMs’ growth kinetics on ramping protocols derived from simulations.
figure 4

a Three ramping protocols κσs(t) (black dashed lines) and the corresponding configurational fluctuations δ(t)/σS of N-CMs (coloured solid lines) with respect to the scaled elapsed time t/τc. b Probability distributions of N-CMs corresponding to the three different ramping protocols in a. We set Ar = 30 kBT in molecular dynamic simulations. Solid lines represent the most probable N from solutions of PN(t) (Δg = 0.4 kBT). c, Tuning the ionic concentration (blue solid lines) and zeta potential (blue dashed lines) over time in a polystyrene colloidal suspension. Each subpanel corresponds to the respective ramping protocol illustrated in a.

Full size image

These comparisons indicate that the ramping should initiate when (N − M)-CMs are prevalent (M could be 1 or 2), even though their configurations deviate significantly from the standard configuration. Otherwise, the free-energy barrier for another satellite to associate with would be too significant. Following this, an appropriate quench rate towards conditions favouring symmetric N-CMs should be applied to condense the distribution to N with the maximum yield. The precise quantitative relationship between the distribution of N, ramping rates, and waiting time remains further exploration. Our theoretical framework can serve as a valuable tool for guiding the design of such ramping protocols.

To provide a more concrete instance of how the physical parameters in the experimental system should be tuned to align with the proposed ramping protocol, we present an example in Fig. 4c. The system under consideration is a 500 nm polystyrene colloidal suspension in a K4Fe(CN)4 aqueous solution, in which the relationship between zeta potential and ion concentration has been previously measured41. In a dilute system, the collision time τc can be estimated as \({\tau }_{c} \sim {l}_{{{{\rm{free}}}}}^{2}/D\), where lfree is the mean free path and D is the diffusion constant. At a 0.1% volume fraction, the estimated collision time is around 2 minutes, which is significantly longer than the time required for ion equilibration. By mapping the screening length κ and ion concentration through the equation \({\kappa }^{-1}=\sqrt{4\pi {\lambda }_{b}{n}_{{{{\rm{ion}}}}}}\), where λb is the Bjerrum length and nion is the effective ion number density, we convert the ramping protocol in Fig. 4a to the corresponding parameters shown in Fig. 4c. In this annealing protocol, the ion concentration decreases by an order of magnitude, from 0.5 mM to 0.01 mM, over a period of several minutes to hours, resulting in an annealing rate of approximately 1 mM/hr. It is important to note that in practice, the required ramping rate may need to be even slower, since the zeta potential increases from 5 mV to 25 mV as the salt concentration decreases. This increase in zeta potential should raise the energy barrier for satellite attachment and thereby slow the assembly kinetics. This effect should be considered in Eq. (1) for a more accurate solution.

Additionally, to implement this concept experimentally, we propose a two-step approach to control the quench rate of the ionic strengths through the size of oil-in-water emulsion droplets, as illustrated in the flow chart in Fig. S4. Alternatively, dynamic tuning of ionic strengths in experiments can be achieved by adjusting the solvent fraction through dilution or evaporation42,43. We hope this experimental scheme will inspire material engineers to apply the proposed theoretical strategy for controlling the in-situ self-assembly dynamics of CMs.

Summary and Outlook

In this study, we systematically explored the dynamic structures and formation kinetics of CMs by combining 3D real-space imaging and computer simulations. We utilised a colloidal-emulsion model system designed to form micrometre-sized dynamic CMs with tunable electrostatic interactions. The key outcomes of our research unveil the intrinsic asymmetry in the dynamic structures of CMs, with their angular symmetry evolving through continuous and sequential ordering from a liquid-like configuration, contrasting with a sharp first-order transition. Notably, our findings also highlight the polymorphic nature of 5-CMs, introducing a richer diversity of behaviours compared to other CMs.

The emergence of asymmetric dynamic structures arises as a collective mode of motion involving a few or several satellite particles on a curved surface. This asymmetric nature in CMs’ dynamic structures is considered universal for CMs whose angular order results from satellite-satellite repulsive interactions through a Thomson-like mechanism. This behaviour, influenced by thermal fluctuations, reflects a balance between energy and entropy. This asymmetry is more pronounced as the effective satellite-satellite repulsive interaction becomes shorter-ranged. In cases where the angular order is a result of a strong core-satellites anisotropic interaction, such as patchy interactions, the motion of satellite particles becomes localised rather than collective, and the asymmetric nature is anticipated to be less pronounced.

These findings offer valuable insights and potential strategies for researchers engaged in assembling hierarchical structures using flexible colloidal molecules. It is worth noticing that configurational fluctuations play a crucial role in shaping self-assembly dynamics, potentially leading to multi-step solid-solid transitions and bifurcation of the growth44,45,46,47,48,49. Nevertheless, utilising fluctuations to engineer a smoother assembly pathway toward the target structure remains a significant challenge. The insights into configurational distributions arising from thermal fluctuations offer valuable knowledge for designing pathways towards complex hierarchical structures.

On the kinetic side, achieving a CM with a specified angular order is challenging due to the multiple and sequential reaction process among colloidal particles. To overcome this difficulty, we proposed a method that dynamically adjusts the screening of charge-charge interaction, guiding the kinetics of CM ordering toward the target structure. Our kinetic model is versatile and can be applied to systems with similar interaction frameworks. By mapping transition rates to experimental parameters, the model facilitates the identification of the parameter space that yields the highest yield in the minimum time. Furthermore, we highlight that there has been growing interest in utilising non-equilibrium pathways to obtain “shortcuts” towards desired states50,51,52. Our colloidal molecule system provides a complex multi-dimensional parameter space to validate these theoretical models from both computational and experimental perspectives.

Methods

Experiments

Our CMs are composed of satellite PMMA colloidal particles with a diameter of approximately σs ≈ 1.2 μm and a central water(20% glycerol)-in-oil emulsion droplet with a diameter σc ranging from 0.8 μm to 1.6 μm. The formation of these CMs involves mixing PMMA colloidal suspensions with emulsions.

The colloidal suspension comprises polymethyl methacrylate (PMMA) particles and grafted with poly(12-hydroxystearic acid). These particles are dispersed in a density and refrective-index matching oil solvent, a mixture of 70% cyclohexyl bromide (CHB) and 30% cis-decalin (volume ratio), at a 5% volume fraction. The measured particle diameter is σs ≈ 1.2 μm with a size polydispersity of 2.5%. The PMMA particles exhibit a negative charge in the solvent, producing a soft interparticle repulsion. The experimentally determined pair potential, Upp(r), is estimated using the equation:

$${U}_{pp}(r)=\frac{{Q}^{2}}{{\sigma }_{{{{\rm{s}}}}}{(1+\frac{\kappa {\sigma }_{{{{\rm{s}}}}}}{2})}^{2}{\epsilon }_{s}}\exp [-\kappa {\sigma }_{{{{\rm{s}}}}}(r/{\sigma }_{{{{\rm{s}}}}}-1)]/(r/{\sigma }_{{{{\rm{s}}}}}),\,r \, > \, {\sigma }_{{{{\rm{s}}}}}$$
(3)

Here, Q represents the surface charge of the colloidal particle, and ϵs is the dielectric constant of the solvent. The term κ−1 denotes the screening length, and \(\frac{{Q}^{2}}{{\sigma }_{{{{\rm{s}}}}}{(1+\frac{\kappa {\sigma }_{{{{\rm{s}}}}}}{2})}^{2}{\epsilon }_{s}}\) is a surface potential parameter.

We prepared the water(20% glycerol)-in-oil emulsions following the same method outlined in Ref. 33. Tetrabutylammonium bromide (TBAB) salt was introduced to modify the Debye length within the oil mixture (0 - 5 μM). The water phase in the emulsion consisted of a 20% glycerol aqueous solution (in volume), and 10 mM sodium dodecyl sulfate (SDS) was added to prevent particle adhesion onto the water-oil interface. The volume ratio between glycerol and the CHB-decalin solution was maintained at 1:40. By subjecting the sample to 2 minutes of sonication, water droplets ranging in size from hundreds of nanometers to 5 micrometres were achieved. Larger droplets were eliminated through coalescence using a centrifuge (100 g, where g is the gravitational acceleration). It is important to note that the water-glycerol droplets also carry a negative charge in the solvent.

The emulsions were subsequently combined with PMMA colloidal suspensions using a vortex mixer with a stoichiometry ratio of approximately 200:1 between particles and emulsion droplets. The particle-droplet interaction is characterised by a long-range, charge-charge repulsion. Moreover, attractive forces, specifically the charge-image-charge interaction and hydrophobic one, dominate the particle-droplet interaction at short-range33. The resulting mixture was left at room temperature for 24 hours to ensure the formation of a stable suspension of colloidal molecules.

We used a Leica SP8 fast confocal microscope equipped with a 63 × /1.40-NA oil-immersion objective to observe the CMs. PMMA particles, dyed with Nile red, were excited at 552 nm, and their emission was detected between 580 and 700 nm. Water droplets, dyed with fluorescein sodium salt, were excited at 488 nm, with emission detected between 500 and 540 nm. The scanning speed was 28 frames per second, with each frame consisting of 128 pixels  × 128 pixels (voxel size: 134 nm/pixel). The spacing between consecutive frames in the z direction was 250 nm. Since the height of a typical CM is less than 4 μm, fewer than 16 images (taking around 576 milliseconds) were required to scan through the entire structure, enabling us to capture the instantaneous configuration. Given the wide range of droplet sizes in our system, we employed an advanced particle-tracking method to identify these droplets. This algorithm, detailed in Refs. 53,54, is based on the Scale Invariant Feature Transform (SIFT), enabling precise localization of particles from three-dimensional (3D) confocal images without making any assumptions about the target sizes.

The values of Q and κ−1 were estimated from the conductivity σ and zeta potential ζ, following the measurement protocol outlined in Ref. 37. The measured values for (σ,ζ) were (0.2 nS/cm, 14 mV) at 1 μM TBAB and (1.0 nS/cm, 12 mV) at 5 μM TBAB. Using the following relationships: \({\kappa }^{-1}=\sqrt{{\epsilon }_{r}{\epsilon }_{0}{k}_{{{{\rm{B}}}}}T/6\pi \eta a\sigma }\), Q = 2πϵrϵ0σs(1 + κσs/2)ζ, where η is the solvent viscosity (2 mPa ⋅ s), a is the hydrodynamic radius of the solute molecule (0.5 nm), and ϵr is the relative dielectric constant (ϵr ≈ 4), we estimated (κ−1, Q,) to be (630 nm, 45 e) at 1 μM TBAB and (270 nm, 63 e) for 5 μM TBAB. The corresponding values of (κσs,Ar) were approximately (2.0, 6.5 kBT) at 1 μM TBAB and (4.5, 4.7 kBT) at 5 μM TBAB.

Simulations

We employed canonical Monte Carlo simulations to sample the equilibrium distribution of CMs’ configurations. Satellite particles were confined to a spherical surface with a diameter denoted as \(\bar{\sigma }\equiv ({\sigma }_{{{{\rm{s}}}}}+{\sigma }_{{{{\rm{c}}}}})/2\), where σs and σc represent the diameter of satellite particles and core droplets, respectively.

The influence of gravity on our system was found to be negligible. The density of the CHB-decalin solvent (1.18 g/cm3) was closely matched to the average density of CMs. The density difference between the PMMA particles and the solvent was less than 1%, while the difference between the PMMA particles and the water droplets was approximately 3%. The gravitational length h was estimated using the relation: \(\frac{\pi }{6}\Delta \rho g{\sigma }^{3}h \sim {k}_{{{{\rm{B}}}}}T\), which yields h ~ 20 μm, significantly larger than the particle size of 1.2 μm. Moreover, as shown in our SI videos, the CMs rotate freely without any observable accumulation of PMMA particles at the bottom. We also performed a centrifugation experiment at 1000 g for 5 minutes with a density-matched sample, where no noticeable sedimentation or density gradient was observed along the Z direction.

The particle-particle pair potential in simulations, Upp(r), is similar to that in experiments:

$${U}_{pp}(r)={A}_{{{{\rm{r}}}}}\exp [-\kappa {\sigma }_{{{{\rm{s}}}}}(r/{\sigma }_{{{{\rm{s}}}}}-1)]/(r/{\sigma }_{{{{\rm{s}}}}}),\,r \, > \, {\sigma }_{{{{\rm{s}}}}},$$
(4)

where the adjustable parameter Ar represents the surface potential parameter of the satellite particle, and κ corresponds to the inverse of the Debye screening length.

To study the dynamics of growing colloidal clusters, we used the HOOMD package for molecular dynamics simulations. We used the same Upp in experiments and MC simulations. The interaction Upd between a satellite particle and a core droplet comprises a weak, long-range repulsive component and a much stronger attractive component. To speed up the simulation, we neglected the long-range repulsive part. In simulations, we use the following simple, attractive pair potential for Upd(r):

$$U(r)=\frac{{A}_{{{{\rm{a}}}}}}{{r}^{6}},\quad r \, > \, \bar{\sigma }$$
(5)

The parameter Aa is negative, and its magnitude is large enough to prevent the satellite-core detachment. We initiated each simulation with one core droplet and 50 satellite particles in the box. The volume fraction of particles was set to 8 %, an order of magnitude higher than that in experiments. The diameter of satellite particles is 1 μm, while the core size ranges from 0.85 μm to 1.45 μm with an increment of 0.10 μm. The dynamic viscosity was set to be equal to that of water, which is 1.0 mPa ⋅ s. For each core size, we ran 200 replicas to study how the fraction of N-CMs in total CMs, denoted as ϕN, evolves with time. Each simulation replica evolved for a duration of 50 seconds.

Configurations and angular orders of CMs

To determine the maximum superimposition of experimental configurations, we implemented a progressive alignment algorithm. Initially, all CMs were normalised and translated to a common position, aligning their droplet coordinates with the origin. Subsequently, a CM was randomly selected as the initial configuration. To achieve the best alignment, we employed the Kabsch algorithm, which computes an optimal rotation matrix (R) to minimise the root mean squared deviation (RMSD) between the initial and subsequent configurations35. RMSD is given by the following formula:

$${{{\rm{RMSD}}}}=\sqrt{\frac{{\sum }_{i}{\left\vert {{{{\bf{r}}}}}_{i}^{0}-{{{{\bf{r}}}}}_{i}^{1}\right\vert }^{2}}{N}},$$
(6)

where \({{{{\bf{r}}}}}_{i}^{0}\) and \({{{{\bf{r}}}}}_{i}^{1}\) are the satellite particle positions of the initial and subsequent configurations, respectively.

The average configuration was determined by iteratively incorporating experimental configurations of N-CMs into the entire population until convergence was achieved. Subsequently, each new configuration was compared to the previous population average and adjusted to its optimal orientation through rotation. The average configuration was then updated by integrating the contribution of the new configuration. This process was repeated multiple times on our set of experimental configurations until the RMSD values reached convergence. The final configuration distributions were derived from the coordinates obtained during the last iteration.

It is crucial to note that the Kabsch algorithm is specifically applicable to sets of paired points. To establish pairs between a new CM and the iterative configuration, RMSD values were calculated for all possible permutations of its satellites with the iterative configuration. Subsequently, we identified the positions that resulted in the minimal RMSD value.

To quantitatively assess the impact of various simulation parameters on the configuration of N-CMs, we employed two parameters: the deviation between the average and standard configurations, denoted as δ, and the magnitude of fluctuations around the average configuration, denoted as Δ. The parameter δ is defined as follows:

$$\delta=\sqrt{\frac{{\sum }_{i}{\left\vert {\bar{{{{\bf{r}}}}}}_{i}^{e}-{{{{\bf{r}}}}}_{i}^{0}\right\vert }^{2}}{N}}.$$
(7)

Here, \({\bar{{{{\bf{r}}}}}}_{i}^{e}\) represents the mean position of the i-th satellite of N-CM, while \({{{{\bf{r}}}}}_{i}^{0}\) denotes the position of the satellite of the standard configuration. Δ is defined as

$$\Delta=\sqrt{\frac{{\sum }_{i}{\sum }_{j}| {{{{\bf{r}}}}}_{ij}^{e}-{\bar{{{{\bf{r}}}}}}_{i}^{e}{| }^{2}}{MN}}.$$
(8)

Here, \({{{{\bf{r}}}}}_{ij}^{e}\) represents the position of the i-th satellite of the jth N-CMs observed, and M is the total number of N-CMs observed. Figs. S1d and e show the evolution of Δ and δ with respect to κσs for 6-CMs, similar to the trend observed in 4-CMs.

To evaluate the angular symmetry of CMs, we computed the bond orientation order parameters q3 and q4 of CMs and studied their distributions. The bond orientation order parameter ql is defined as

$${q}_{l}=\sqrt{\frac{4\pi }{2l+1}\mathop{\sum }_{m=-l}^{l}| {q}_{lm}{| }^{2}},$$
(9)
$${q}_{lm}=\frac{1}{N}\sum {Y}_{l}^{m}(\theta,\phi ),$$
(10)

where θ and ϕ are the azimuthal and polar angles of satellite particles concerning the core, respectively, and \({Y}_{l}^{m}\) are the spherical harmonics. Figs. S2a-c illustrate the evolution of the bond orientation order parameters for 3-CMs and 6-CMs, mirroring the trend observed in 4-CMs.

To determine whether a 5-CM aligns more closely with a triangular bipyramid or a pyramid configuration, we rotate them to achieve maximal superimposition with standard triangular bipyramid and pyramid configurations having distinct heights. The minimal RMSD of pyramid configurations is then compared with the RMSD of a standard triangular bipyramid. A smaller RMSD indicates proximity to the correspondent category.

Figs. S3a and b present the details of the two bifurcation pathways. At κσs = 10, the average configuration of 5-CMs develops towards the pyramid as the core diameter σc decreases to 0.85 (Fig. S3a). However, when κσs = 5, the development of 5-CMs leans more towards the triangular bipyramid as σc reduces (Fig. S3b). This bifurcation stems from the free-energy difference between the two polymorphs of 5-CMs at different κσs and σc. Due to such a bifurcation, achieving a pure triangular bipyramid for 5-CMs through a pathway of reducing κσs becomes challenging. A more feasible pathway involves reducing σc at κσs < 5.

Comparison between experiments and simulations

To match the experimental conditions with the simulation parameters, we explored a parameter space defined by two fitting parameters: Ar and κσs. The parameter Ar was varied from 5 kBT to 30 kBT in increments of 5 kBT, while κσs ranged from 2 to 16 in steps of 2. To quantify the similarity between experimental and simulated configurations, we defined a parameter δ0 that accounts for both the average and the spread of the configurations. This parameter is defined as:

$${\delta }^{0}=\sqrt{{\delta }_{\,{{\mbox{m}}}}^{2}+{\delta }_{{{\mbox{s}}}\,}^{2}},$$
(11)

where

$${\delta }_{{{{\rm{m}}}}}=\sqrt{\frac{{\sum }_{i}| {\bar{{{{\bf{r}}}}}}_{i}^{e}-{\bar{{{{\bf{r}}}}}}_{i}^{s}{| }^{2}}{N}},$$
(12)

and

$${\delta }_{{{{\rm{s}}}}}=\sqrt{\frac{{\sum }_{i}{(\Delta {{{{\bf{r}}}}}_{i}^{e}-\Delta {{{{\bf{r}}}}}_{i}^{s})}^{2}}{N}}.$$
(13)

Here, δm captures the difference in the mean positions of vertices between experiments and simulations, while δs describes the disparity in positional fluctuations. \({\bar{{{{\bf{r}}}}}}_{i}^{e}\) and \({\bar{{{{\bf{r}}}}}}_{i}^{s}\) represent the mean positions of the i-th vertex of N-CMs in experiments and simulations, respectively. Similarly, \(\Delta {{{{\bf{r}}}}}_{i}^{e}\) and \(\Delta {{{{\bf{r}}}}}_{i}^{s}\) denote the standard deviation of the i-th vertex’s position relative to its mean in both experiments and simulations. Before calculating δ0, we aligned the configurational distributions from experiments and simulations by rotating them to achieve maximal overlap. For each salt concentration, we computed the mean δ0 over a wide range of core particle sizes σc (from 0.85 σs to 1.15 σs) and satellite numbers N (ranging from 2 to 6) to identify the optimal simulation parameters. The optimal parameters were found to be Ar = 7.5 kBT and κσs = 2 for 1 μM TBAB, and Ar = 5 kBT and κσs = 8 for 5 μM TBAB. Under these conditions, δ0 was generally less than 0.1 σs. Although the precise values differ from those estimated from experimental measurements, the trend of asymmetry indicated by our fit aligns with experimental observations. We hypothesize that this deviation may result from the limitations of the Yukawa model at higher salt concentrations, where image charge effects become significant.