Introduction

Photonic and meta- materials provide examples where ordered particles on length scales comparable to electromagnetic wavelengths produce exotic emergent properties1. Colloidal crystallization provides a potential route to self-assemble such materials via processes amenable to scalable manufacturing; however, robust schemes have yet to be identified to obtain the necessary low defect densities. From a broader perspective, obtaining perfect crystals on any length scale remains more art than science (e.g., atoms, molecules, macromolecules)2,3. To design robust crystal growth, recrystallization and annealing schemes to minimize defects, it is necessary to know the basic mechanisms of defect formation and motion.

Although three-dimensional configurations of complex particles are the ultimate goal of self-assembly schemes4,5, here we investigate a relatively simple problem that is still not well understood: how grain boundaries form and move during quasi two-dimensional (2D) crystallization of spherical colloids. Quasi-2D colloidal crystals have been used in studies of melting6, nucleation7, point defect diffusion8 and grain boundary fluctuations9. Other relevant studies include 2D analyses of colloidal crystals to investigate impurity mediated growth10, particle motion within grain boundaries11 and grain boundary pre-melting12. Despite these extensive studies, fundamental understanding of grain boundary formation and motion remains rudimentary13 in all but the most model simulation studies14. 2D crystals are also of interest based on their relevance to thin films15, bubble rafts16 and graphene17. Understanding how grain boundary motion enables relaxation of multi-domain crystals into defect-free crystals is therefore scientifically and technologically interesting.

Results

In this work, we employ real-time microscopy to observe grain boundary formation and motion in a quasi-2D colloidal crystal containing ~200 colloidal particles (Fig. 1a, Supplementary Videos 1–4). Aqueous ~3 micron SiO2 charged colloids crystallize in a quadrupole electrode in MHz AC electric fields (see Methods, Supplementary Information)18, where field-mediated compression of induced dipoles is balanced by the quasi-2D colloid osmotic pressure (i.e., effective hard disk). The electric field amplitude (i.e., applied voltage) is effectively a surrogate for pressure and hence acts as a global thermodynamic variable that determines the relative free energy of all particle configurations for fixed voltage, number and temperature. In the following, we report non-dimensional voltages, where V* = 1 is the voltage required for N particles to produce a hexagonal close packed crystal with hexagonal morphology (as demonstrated in previous work with agreement between microscopy experiments, MC simulations and perturbation theory19, see Supplemental Information for additional details).

Figure 1
figure 1

Reaction coordinates for grain boundary formation and motion can be computed from image analysis of optical microscopy images of electric field mediated colloidal crystallization.

(a) Raw images show representative configurations of 210 ~ 3 μm silica colloids in fluid, bicrystal and single crystal configurations observed over ~7.5 min following a quench (step voltage change) to V* = 0.57 (see main text and Supplementary Methods for definition). Snapshot times are shown by solid black lines in Fig. 1e. Computed reaction coordinates shown by colored particle centers on 8-bit intensity scale for the (b) radius of gyration, , (c) local hexagonal order, , (d) global hexagonal order, and time dependent traces for (e) a single voltage quench and (f) ten consecutive cycles (numbered vertical black lines indicate representative images included in Supplementary Fig. 2).

Full size image

Our previous characterization of electric field mediated colloidal interactions and assembly has yielded kT-scale potentials20,21, feedback control over system size18 and conditions to crystallize N particles19. Based on these findings, we set N = 210 in Fig. 1, which routinely forms bicrystals (i.e., 1 grain boundary between 2 domains) in contrast to single domains in smaller systems and >2 domains in larger systems. As shown in Fig. 1, step-quenches to V* = 0.57 cause an initially dilute fluid phase to first rapidly condense, then form grain boundaries via coalescence of local domains and finally display grain boundary motion as bicrystals relax to single crystals. Reversibility allows repeated quenches between fluid and crystal states to probe the stochastic dynamics of grain boundary formation, diffusion (i.e., random motion) and migration (i.e., drift). We first focus on the V* = 0.57 case, but later investigate step-changes to both lower and higher values of V*.

To interpret and model these measurements, we aim to develop a “low-dimensional model” that quantitatively captures the observed dynamics using “reaction coordinates”22 (rather than enumerating all 2N translational degrees of freedom)23,24. The use of the term “low-dimensional” here does not refer to the Euclidian spatial dimension (i.e., x and y particle center coordinates), which is effectively quasi-two dimensional, but rather, “dimensionality” refers to the number of reaction coordinates necessary to capture the stochastic dynamics of grain boundary formation. It is “low dimensional” because we expect the number of reaction coordinates to be considerably less than the 2 × (N = 210) = 420 dimensions that would be necessary to uniquely specify all possible two-dimensional configurations of 210 particles. Although the term “low-dimensional” may be unfamiliar to some readers, other synonymous terms such as “coarse-grained models” could also be confusing based on an unfortunate overlap of terms relevant to the application of interest in this work.

Candidate reaction coordinates are computed from particle centers and used to color code images (Figs. 1b–f) including: the radius of gyration, Rg18, to capture condensation from fluid to crystal states, average local hexagonal order, 〈C6〉25, to capture the onset of crystallization and global hexagonal order, ψ626, to capture the degree of polycrystallinity. Rg and 〈C6〉 are normalized by their N particle single crystal values (see Supplementary Methods), so that 〈C6〉 goes from 0–1 for fluids to complete locally ordered states and Rg decreases from arbitrarily high numbers to 1 for complete condensation19. ψ6 is 0 for fluids and 1 for single-domain crystals like 〈C6〉, but in contrast, depends strongly on relative domain size and misorientation (e.g., ψ6 = 0 for bicrystals of identically sized 111 domains with 30° misorientation).

The reaction coordinates trajectories following a single voltage quench (Fig. 1e) show initially decreasing Rg and increasing 〈C6〉 occur in unison (i.e., Rg−1 ≈ 〈C6〉), which demonstrates a close coupling between condensation and local ordering. As Rg and 〈C6〉 plateau, indicating an overall condensed configuration with all particles contained in locally crystalline domains, ψ6 ≈ 0, indicating a bicrystal. For t > 100 s, Rg and 〈C6〉 remain essentially unchanged while ψ6 rises from 0 to ~0.8, which clearly corresponds to grain boundary motion from the bicrystal interior to the periphery where it vanishes (Fig. 1d, Supplementary Video 1).

Ten successive voltage quenches from initial fluid states (Fig. 1f) demonstrate the stochastic nature of the grain boundary dynamics, which is expected from the underlying probabilistic colloidal motion. While Rg and 〈C6〉 reveal condensation and local order emerge in a similar manner for each cycle, the ψ6 trajectories can be categorized into several cases: (1) ψ6 tracks 〈C6〉 indicating the simultaneous emergence of local and global order, (2) ψ6 becomes localized at intermediate values for varying time periods before again increasing, (3) ψ6 becomes arrested for the duration of the observation time and (4) in one case, ψ6 initially increases but then vanishes. The stochastic nature of the observed grain boundary dynamics apparent in an ensemble of trajectories is an important aspect to capture in a quantitative model.

To develop a low-dimensional model of grain boundary formation and motion, it is necessary to determine the number and type of reaction coordinates. For example, it could be speculated that ψ6 is all that is required to track grain boundaries since it visually tracks polycrystallinity in Fig. 1. However, simply tracking ψ6 does not capture how parallel processes of local condensation and crystallization determine the formation and motion of grains of different sizes, shapes and orientations, which ultimately determine the mechanisms of polycrystals relaxing to single crystals. In addition, ψ6 does not uniquely identify some configurations; for example, fluid configurations and maximally misaligned bicrystals both have ψ6 = 0.

To illustrate how reaction coordinate pairs capture additional information, experimental configurations are colored using two coordinates (Figs. 2a, b) and 2D trajectories of (ψ6, 〈C6〉) (Fig. 2c) or (ψ6, Rg) (Fig. 2d) vs. time, where time is indicated by a 256-color scale. We do not plot (Rg, 〈C6〉) since these coordinates are highly correlated in Fig. 1, although such a pair could be useful to distinguish condensed amorphous microstructures (i.e., glasses/gels). Such 2D plots immediately address one issue; plotting ψ6 against either Rg or 〈C6〉 distinguishes fluid and bicrystal configurations (e.g., high Rg, low ψ6 vs. low Rg, low ψ6). These trajectories also show how local ordering (i.e., increasing 〈C6〉) and condensation (i.e., decreasing Rg) influence the emergence of polycrystallinity and subsequent grain boundary motion. Trajectories starting at lower ψ6 tend to become localized at lower ψ6 after condensation, whereas trajectories with initially higher ψ6 tend to rapidly form single crystals.

Figure 2
figure 2

“Two dimensional” trajectories (i.e., two reaction coordinates) capture coalescence of local domains during fast condensation processes to produce bicrystals that relax to single crystals via grain boundary motion over a broad range of timescales.

Microscopy images from Figs. 1a with particle centers colored using RGB (Red/Green/Blue) mixing rules for colors represented by two reaction coordinates including (a) ψ6 and C6 and (b) ψ6 and Rg to visualize how global order emerges from local order and during condensation. Ten trajectories following quenches to V* = 0.57 with time represented by a 256-color scale (inset scale bar) for (c) and (d) reaction coordinate pairs with inset plots of single trajectories from Fig. 1e.

Full size image

Although tools exist to identify the minimum dimensionality (i.e., number of reaction coordinates) from observed dynamics (e.g., diffusion mapping27,28,29), such methods are not currently able to predict physically meaningful reaction coordinates22. Ultimately, the number and types of reaction coordinates can be determined empirically by finding what is necessary to produce a quantitative stochastic dynamic model. Because the measured trajectories display both drift and diffusion along reaction coordinates, which appear to be mediated by free energy gradients (i.e., driving forces) and fluctuations/friction (i.e., randomness/resistance), it is assumed that such processes can be captured by a low-dimensional Smoluchowski equation given by30,

where p(x, t) is the probability density of finding the system at coordinate x at time, t, W(x) is the free energy landscape, kT is thermal energy and D(x) is the diffusivity landscape, which can be related in the usual way to mobility, m(x) and friction, ζ(x), landscapes as D(x) = m(x)kT = ζ(x)−1kT. The variable x is a vector of reaction coordinates where the number of coordinates is the model “dimension.” The Smoluchowski equation is a special case of the Fokker-Planck equation that also satisfies the fluctuation-dissipation theorem and leads to Boltzmann sampling at equilibrium (i.e., p(x) = exp[-W(x)/kT])23,30. In more descriptive terms, W(x) is the free energy change associated with moving from one particle configuration to another (or one reaction coordinate to another; (e.g., x1 = (ψ6,1, Rg,1) to x2 = (ψ6,2, Rg,2)) and D(x) captures the associated configuration dependent changes in diffusion and friction (i.e., fluctuations and dissipation).

To obtain W(x) and D(x) in Eq. (1) from trajectories like those in Figs. 1 and 2, we analyzed Brownian Dynamic (BD) simulations that were matched to experiments (by capturing all equilibrium and dynamic properties of the quadrupole experiment on the particle scale19,20,21, see Supplementary Methods). This approach was used because statistics on the particle scale are easily obtained to match experiments and simulations, but BD simulations are better suited to generating large numbers of grain boundary trajectories (e.g., each experimental grain boundary trajectory in Fig. 1 is acquired for ~10 min.).

To provide more details of the matching process, inverse Monte Carlo was used to obtain interaction potentials that capture all equilibrium properties (i.e., radial distribution functions), which were then employed in BD simulations that captured all dynamic properties (i.e., particle scale diffusion, reaction coordinate trajectories). It is important to note that it was necessary to include concentration dependent, but not configuration dependent, hydrodynamic interactions in the BD simulations to produce agreement with experiments. In particular, Stokesian Dynamic simulations31 were used to show that approximate local, concentration dependent hydrodynamic interactions were sufficient to match experiments without the need to consider multi-body, configuration dependent hydrodynamic interactions. This finding indicates local microstructure or morphology do not obviously influence the hydrodynamic interactions involved in the colloidal crystal grain boundary dynamics observed in our work.

Statistical methods reported in the literature23 and further developed by us for application to colloidal assembly32, were used to analyze large numbers of BD simulated trajectories to construct W(x) and D(x) (see more details in Supplementary Methods and our previous work29,32). In brief, the displacement and mean squared displacement of reaction coordinate vs. time trajectories can be used to measure drift and diffusion at each value of x, which ultimately yield W(x) and D(x).To assess the quantitative accuracy of candidate low-dimensional dynamic models, we compared first passage time distributions for ensembles of trajectories between different starting and ending states from particle-scale BD simulations and low-dimensional Langevin dynamic (LDLD) simulations. The LDLD simulations are based on a Langevin equation given as,

where the coefficients are the same as in Eq. (1), Δt is the integration time step and Γ(t) is a noise variable. A successful LDLD model will accurately reproduce first passage time distributions obtained from high 2N-dimensional particle-scale BD simulations and hence the experimental trajectories that were quantitatively matched to the BD simulations.

The only quantitatively accurate model to emerge from candidate 1D and 2D models was one based on (ψ6, Rg) (so higher dimensional models were not considered). Fig. 3a shows the W(x) that quantifies the relative free energy of every configuration and the free energy gradients that drive motion. The global free energy minimum at x = (ψ6, Rg) ≈ (0.8, 1.14) indicates a single domain crystal (with a thin fluid envelope at its periphery) is the thermodynamically favored configuration (image VI in Fig. 3c).

Figure 3
figure 3

Free energy and friction landscapes obtained by fitting trajectories to Smoluchowski equation (equation 1) that quantitatively capture grain boundary diffusion and migration vs. relative domain sizes and misorientation angles.

(a) W(ψ6, Rg)/kT with inset scale bar and two trajectories obtained from the experiments in Figs. 1 and 2. (b) Dψ6ψ6/(kT·s) with inset scale bar with the same trajectories as in (a). (c) Representative configurations from microscopy images and simulated renderings for coordinates marked along trajectories in (a) and (b). Labeled misorientation angles shown by red lines and particle centers are colored according to the RGB composite convention in Fig. 2. First passage time distributions for BD (red) and LDLD (black) trajectories projected onto (d) the ψ6-axis between 0.47–0.66 , 0.38–0.56 , 0.28–0.47 , 0.19–0.38 and (e) Rg-axis between 1.27–1.25 ,1.25–1.23 , 1.23–1.21 , 1.21–1.19 , 1.21–1.19 , 1.19–1.17 and (f) for trajectories between a sink at (ψ6 = 0.8 Rg = 1.14) (i.e., global minimum) and sources at (ψ6 = 0.38 Rg = 1.15) and (ψ6 = 0.65 Rg = 1.16) .

Full size image

The diffusivity landscape, D(x), in Eqs. (1) and (2) is a 2 × 2 tensor; it has diagonal components that capture how friction/mobility in each reaction coordinate mediate drift due to free energy gradients in the same coordinate (Dψ6ψ6, DRgRg) and cross-terms that capture how friction/mobility mediate drift due to free energy gradients in orthogonal coordinates (Dψ6Rg, DRgψ6). Fig. 3b shows Dψ6ψ6/kT, which captures how friction mediates drift and diffusion in ψ6 due to free energy gradients in ψ6. The other components of D(x) are reported in Supplementary Fig. 5. DRgRg shares similar features with Dψ6ψ6. The cross-terms (Dψ6Rg, DRgψ6) indicate a weak coupling between driving forces and drift/diffusion (i.e., friction increases as free energy increases) along orthogonal coordinates for ψ6≳0.5, although this is relatively minor compared to the diagonal terms.

Representative configurations (Fig. 3c) and first passage time distributions (Fig. 3d–f) show the resulting low dimensional ψ6, Rg model quantitatively captures the measured grain boundary dynamics. The agreement between the BD and LDLD simulations is excellent (Fig. 3d–f), demonstrating that Eqs. (1) and (2) with the W(ψ6, Rg) and D(ψ6, Rg) in Fig. 3 provide accurate low dimensional dynamic models of the experiments in Figs. 1 and 2.

Discussion

To aid discussion of how features on W(ψ6, Rg) and D(ψ6, Rg) are connected to microscopic mechanisms, two limiting trajectories are shown on these landscapes; one where a grain boundary forms between two domains and does not move out of the crystal on the ~10 min observation time [trajectory 1 (T1): I–II–IV–V] and one where two domains form but the grain boundary quickly moves to the crystal periphery in <1 min to produce a single domain crystal [trajectory (T2): I–III–VI]. Both trajectories are consistent with expectations by showing drift (i.e., migration) along free energy gradients and superimposed stochastic motion (i.e., diffusion) that is most evident where free energy gradients are minimal. The dramatic difference between these two trajectories is most evident as they approach the global minimum (Fig. 3f) (as quantified between “sources” at III and IV and a “sink” at VI). An order of magnitude difference is observed in the most probable first passage times with ~20 s for III–VI and ~200 s for IV–VI; the latter distribution also shows a much longer asymmetric tail with some trajectories taking >1,000 s to traverse the W(ψ6, Rg) plateau at low Rg.

The T1 trajectory corresponds to rapid condensation along a steep free energy gradient where two locally ordered domains coalesce into a bicrystal with a near maximum 30° misorientation angle. At point IV on T1, the grain boundary randomly diffuses on a free energy plateau with a minimal free-energy gradient (i.e., driving force) to drive migration of the grain boundary to the crystal edge; the trajectory is localized between IV–V for ~10 min. In addition to vanishing free energy gradients, the friction in the vicinity of IV–V is increased ~6× compared to uncondensed states. Multi-body hydrodynamic interactions (i.e., near-field lubrication and far-field flow within the particle structure33) increase particle-scale friction during condensation, which is consistent with the increased friction for trajectories at low Rg on D(ψ6, Rg) (Fig. 3b). Although first passage times vary linearly with frictional changes compared to an exponential dependence for free energy changes34,35, diffusion mediated by friction is the rate determining process in the presence of vanishing free energy gradients (plateau of Fig. 3a).

In contrast to T1, the T2 trajectory is initiated with higher global order before moving down the free energy gradient and rapidly continuing towards the global free energy minimum single crystal. Although T2 passes close to T1, low friction at high Rg allows sufficient diffusion towards higher ψ6, As a result, T2 bypasses the free energy plateau region at low Rg to avoid slow diffusion like T1 and instead shows much faster grain boundary migration. Friction uniformly increases with decreasing Rg (due to hydrodynamic interactions) and has almost no dependence on ψ6, so there is no path of least resistance on D(ψ6, Rg). In short, the fastest trajectories are ones that bypass the free energy plateau region.

The microscopic mechanisms associated with these different trajectories can be understood from the images/renderings (Fig 3c) and the physical meaning of the reaction coordinates. Rg clearly captures condensation as shown by large free-energy gradients on W(ψ6, Rg) and increasing resistance to configurational changes on D(ψ6, Rg) as the result of multi-body hydrodynamic interactions. At this point, we can speculate why 〈C6〉 was not part of a successful dynamic model; it is an indirect measure of condensation and thus not as good as Rg and it is the emergence of global order, captured by ψ6, that is most important to track grain boundaries.

The ψ6 dependence of W(ψ6, Rg) indicates that domains coalescing with minimal misorientation produce higher global order from the outset, which also translates to faster grain boundary migration (via free-energy gradients) from the crystal interior to the periphery. Practically, low misorientation angles produce smaller energy barriers to particle-scale motion (in full 2N-dimensional particle-scale space) within (e.g., string-like motion) and across (e.g., cooperative motion) grain boundaries36. In contrast, domains that coalesce near the maximum 30° misorientation display low initial global order that translates into slow grain boundary diffusion on a W(ψ6, Rg) plateau. Such bicrystals represent an unstable equilibrium where the energy (e.g., energy/atom, interfacial energies)15 of the two sides balance, however, the lower free energy state single crystal (~10 kT) is achieved by fluctuations that eventually allow one grain to increase at the expense of the other grain decreasing.

Because ψ6 does not resolve different combinations of domain size and misorientation, but accurately captures the dynamics, it appears that all such configurations relax in an indistinguishable manner. In particular, greater misorientations between dissimilar sized grains produces relaxation rates equivalent to cases where domains have less misorientation but are of similar size. This finding shows how using ψ6 as a reaction coordinate indicates an aspect of grain boundary motion that would not be easily discovered from tracking particle-scale motion alone. Because grain boundary motion involves many particles rearranging in a cooperative fashion based on relative domain sizes and orientations, it is useful to have a global parameter that captures configurational changes of the entire particle ensemble and therefore naturally captures cooperative phenomena. Although ψ6 should be a good reaction coordinate to monitor relaxation of bicrystals of any size (i.e., distinguish two-domain from single-domain crystals), it may be unsuitable for many-domain crystals because it would remain near zero and miss most of the structure evolution until only a few domains remained. In many-domain crystals, it might be more appropriate to capture the evolution of polycrystallinity using a different reaction coordinate, like average domain size or average misorientation angle across grain boundaries.

The approaches used to measure and model grain boundary motion for V* = 0.57 in Figs. 1,2,3 can be applied at other V*. Using the same BD simulations and non-equilibrium analyses, W(ψ6, Rg) were constructed in the range V* = 0.31–0.69 (Fig. 4a–d) along with representative trajectories and global minimum configurations. At the lowest V*, particles are weakly confined in concentrated fluid configurations without crystal grains or boundaries, which produces a relatively featureless W(ψ6, Rg). At the highest V*, W(ψ6, Rg) is qualitatively similar to Fig. 3a, but the plateau stretches from very low to high ψ6 with an even shallower gradient and deeper global minimum at ψ6 → 1. At intermediate V*, the W(ψ6, Rg) show a continuously shifting global minimum towards lower Rg and higher ψ6 and a stretching plateau corresponding to slower grain boundary migration. D(x) vs. V* are not reported here, but the general trend with increasing voltage is a decreasing magnitude (i.e., decreasing mobility, increasing friction), which is consistent with more condensed configurations hindering particle motion and hence motion along both reaction coordinates. In short, with increasing compression, once grain boundaries form, they experience slower migration and diffusion and the single perfect crystal clearly emerges as the global free energy minimum configuration.

Figure 4
figure 4

Free energy landscapes based on ψ6, Rg reaction coordinate pair capture dynamics at all V* to provide quantitative models of grain boundary formation and motion.

W(ψ6, Rg)/kT at (a) V* = 0.31, (b) V* = 0.44, (c) V* = 0.57 and (d) V* = 0.69 with inset renderings of global minimum configuration and representative trajectories from BD simulations that were matched to the experiments (see Supplementary Methods). The free energy scale is indicated by the inset in (a).

Full size image

The results show voltage cannot simply be increased to increase order, since in general this produces increasingly arrested polycrystalline states. It is interesting to consider how grain boundary formation and motion might be manipulated by “switching” between W(ψ6, Rg) at different V*. If V* is increased very slowly, it would be possible to remain in the global free energy minimum configuration at all voltages; this is the thermodynamic equilibrium limit and a known strategy to make single crystals. However, faster schemes are generally desirable. The optimal control policy37 to achieve a single perfect crystal in minimal time could employ the quantitative non-equilibrium dynamic models in Eqs. (1) and (2) to switch between W(ψ6, Rg) at different V* in an automated, informed manner using feedback control38. In particular, monitoring reaction coordinates in real-time could identify slowly relaxing polycrystalline configurations (i.e., due to vanishing free-energy gradients and high friction) and then V* could be tuned to “land” on another W(ψ6, Rg) at the same coordinates where faster relaxation occurs.

In summary, we report agreement between optical microscopy measurements, Brownian Dynamic simulations and low-dimensional models of stochastic grain boundary formation and motion in quasi-2D colloidal bicrystals. Our results show that two reaction coordinates, one for condensation and one for global order, are sufficient to quantitatively capture first passage times between critical configurations at each applied voltage. Free energy and diffusivity landscapes show that the relative misorientation angles and domain sizes formed during condensation determine the subsequent grain boundary motion. Bicrystals with similar sized domains and a near 30° maximum misorientation angle relax via slow grain boundary diffusion mediated by high friction and vanishing free energy gradients, whereas bicrystals with asymmetrically sized and/or less misoriented domains relax via much faster grain boundary migration due to greater thermodynamic driving forces. By quantifying such dynamics as a function of voltage, ongoing work is developing optimal control algorithms to dynamically tune voltages to avoid kinetic bottlenecks associated with slow grain boundary dynamics.

Future work will extend the modeling approaches developed in this work to understand other mechanisms and defect types involved in colloidal crystal formation/relaxation (e.g., nucleation, growth, point defects, multiple domains, etc.), which will require identification of the minimal number and type of order parameters to accurately capture relevant dynamic processes. Additionally, the role of multi-body, configuration dependent hydrodynamic interactions on diffusivity (friction) landscapes will be explored in attractive particle systems, where changes in near-field interactions during clustering39 and percolation40 have previously been shown to have non-trivial dynamic signatures.

Methods

Coplanar gold thin film quadrupole electrodes were patterned on glass microscope coverslips by spin coating photoresist and physical vapor deposition of a 15 nm chromium layer and a 35 nm gold layer. Nominal 3.13 μm diameter SiO2 colloids with ~50 mV zeta potentials were fractionated in DI water and centrifuged/redispersed five times in 0.1 mM NaOH. PDMS o-rings were coated with vacuum grease and sealed between a coverslip with the patterned quadrupole electrode before it was connected in series with a function generator. Microscopy was performed on an inverted optical microscope with a 63× objective and a 12-bit CCD camera that captured 336 pixel × 256 pixel (81 μm × 62 μm) digital images at rate of 8 frames/s. Video capture and image manipulation were performed using algorithms in MATLAB.

BD simulations in the canonical ensemble were performed for 210 colloidal particles at constant voltage using numerical methods described in previous papers31,32,40,41,42. A 0.1 ms time step was used for at least 2 × 107 steps and reaction coordinates were stored every 1250 steps for subsequent analysis. Particles in simulations were confined within 2D planes. Inverse Monte Carlo methods (including image resolution limiting effects)43,44,45 were used to match measured and simulated radial distribution functions to determine parameters in interactions potentials. The diffusivity was matched by comparing measured and simulated mean square displacements. Parameters used in the BD simulations are reported in the Supplementary Information.