Removing grain boundaries from three-dimensional colloidal crystals using active dopants

Abstract

Using computer simulations we explore how grain boundaries can be removed from three-dimensional colloidal crystals by doping with a small fraction of active colloids. We show that for sufficient self-propulsion, the system is driven into a crystal-fluid coexistence. In this phase separated regime, the active dopants become mobile and spontaneously gather at the grain boundaries. The resulting surface melting and recrystallization of domains result in the motion of the grain boundaries over time and lead to the formation of a large single crystal. However, when the self-propulsion is too low to cause a phase separation, we observe no significant enhancement of grain growth.

---

I. Introduction

One of the key challenges in colloidal self-assembly involves circumventing the long time-scales required for systems to reach their equilibrium states. For example, in colloidal systems where a crystalline phase is thermodynamically stable, the self-assembly of such a phase typically requires crossing a high free-energy barrier. Hence, both nucleation of the equilibrium crystalline phase and the annealing out of crystal defects often occur on prohibitively long time-scales.^1,2

An exciting and promising new avenue to aid colloidal systems in crossing free-energy barriers stems from the recent development of active colloidal particles.^3–9 Such particles actively use energy from their environment to generate directed motion and, as such, systems containing active particles are inherently out of equilibrium. Recent experimental and simulation studies have shown that the dynamics of passive systems can be altered dramatically by the incorporation of active particles.^10–14 Interestingly, even for very low concentrations (<1%), where the self-propelled particles can be viewed as active “dopants”, they can help the underlying passive system to equilibrate by facilitating the crossing of free-energy barriers. Specifically, it has been shown that active dopants can both assist in the nucleation of crystalline domains in high-density hard-sphere glasses,¹⁰ as well as in the removal of grain boundaries in two-dimensional colloidal crystals.¹³

Here we return to the question of grain boundary removal, and explore whether active dopants can also speed up the removal of grain boundaries in three dimensions. In comparison to the 2d case, particles in 3d colloidal crystals are typically more strongly caged, because the coordination number is higher, and as a result the free-energy barriers involved in annealing out grain boundaries and other defects are expected to be significantly higher. Hence, it is interesting to determine whether active dopants can successfully overcome these higher barriers, and thus assist in coarsening a 3d colloidal polycrystal.

It should be noted that large, defect-free colloidal crystals have important photonic uses,¹⁵ and that the growth of such crystals is typically challenging, depending on methods based on, e.g. epitaxial growth from a template,^16–20 the application of external fields,^21–24 or temperature gradients.^25,26 Hence, the possibility of actively removing grain boundaries by the addition of active particles shows great promise in future applications.

II. Method and model

A. Model

We consider a three-dimensional system consisting of N particles that interact via the purely repulsive Weeks–Chandler–Andersen (WCA) potentialwith σ the particle diameter, βε = 40 the energy scale, and β = 1/k_BT, where k_B is the Boltzmann constant and T is the temperature. From the N particles, we identify a subset of particles N_a which we “activate” by applying a constant self-propulsion force f along the self-propulsion axis û_i. We denote the fraction of active particles by α = N_a/N.

We simulate the system using overdamped Brownian dynamics. Thus, the equations of motion for particle i read:

where the translational diffusion coefficient D₀ and the rotational diffusion constant D_r are linked via the Stokes–Einstein relation D_r = 3D₀/σ². Additionally, ξ_i(t) and η_i(t) are stochastic noise terms with zero mean and unit variance. Note that for passive particles f = 0. We measure time in units of the short-time diffusion τ = σ²/D₀.

B. Formation of initial polycrystals

We create starting configurations that are polycrystalline by breaking up the simulation box into polyhedral Voronoi cells which we fill with randomly oriented crystal domains.²⁷ Specifically, we generate a set of randomly placed grain centers, and fill all space closer to a specific center than to any other center with a randomly oriented face-centered cubic crystal domain. In these configurations, it is still possible that two particles are too close to each other. To avoid these “overlaps”, we remove one of the particles if the center of masses are separated by less than the effective hard-sphere diameter, as extracted from the coexistence densities in ref. 28. Finally, the polycrystal is annealed for ∼50τ to ensure that the grain boundary structure is relaxed.

C. Fluid-crystal coexistence simulations

We use direct coexistence simulations to determine the fluid-crystal coexistence densities. More specifically, the coexistence between the crystal and the fluid is studied in an elongated box of L_x = L_y ≈ L_z/4 containing N = 16 000 particles. At the start of the simulation, the entire box is filled with a crystal of uniform density. A fraction α of the particles are then labelled as “active” particles. To speed up equilibration, the active particles are chosen in a narrow slab in the box.³⁴ During the simulation, the crystal partially melts near the region rich in active particles. As a result, the system forms a crystal-fluid coexistence, with the two interfaces perpendicular to ẑ. As the fluid phase has a lower density than the overall system, the density of the crystal region increases. As a result, the crystal becomes compressed along the z-direction introducing undesirable stresses in the crystal. To allow these stresses to relax, we developed a simple algorithm to reshape the box during equilibration. Specifically, we measure the virial stress tensor in the bulk of the crystal region. Note that in this region, the fraction of active particles is negligible, and thus the virial stress is expected to be a good approximation of the true stress. We then average the stress in the crystal region over 0.2τ, and subsequently make a small adjustment to the aspect ratio of the box which reduces the stress, while keeping the total volume constant. We then continue the simulation and repeat this process until the system is in a steady state and the stress in the crystal phase is effectively zero. We then further equilibrate for 100τ with a constant box shape, and finally start measuring the coexistence densities.

D. Cluster algorithm

To distinguish different crystalline domains, we use a cluster algorithm based on the local bond-order parameter q₆.²⁹ As a first step, we identify the neighbours for each particle, defined as all particles located within a distance r_c = 1.5σ of the particle in question. We denote the number of neighbours for particle i as . The bond orientational order parameter of a particle i is then calculated usingwhere Y_6,m(θ,ϕ) are spherical harmonics, θ_ij and ϕ_ij denote the polar and azimuthal angles associated with the center-of-mass distance vector r_ij = r_j − r_i, and m ∈ [−6, 6]. We then quantify the degree of correlation in the environment of two particles i and j using Particles within the same domain exhibit d₆(i,j) ≈ 1 as a result of the translational symmetry within the local crystal lattice. On the other hand, for pairs of particles at the grain boundary d₆(i,j) < 1, as the translational symmetry is broken due to the different crystallographic orientations of the domains.

We define the number of crystalline neighbours ξ(i) of a particle i as the number of neighbours j with d₆(i,j) > 0.95. Subsequently, we consider a particle as being part of the crystal phase when ξ(i) ≥ 9. Finally, we consider two particles to be part of the same domain when they are crystalline neighbours.

III. Results

In ref. 13, we showed that it was possible to anneal out grain boundaries in two-dimensional crystals through doping with active particles. Here, we explore whether this approach is also applicable to three-dimensional systems, where the particles are more strongly caged. To this end, we explore the evolution of a polycrystal both in the case of a purely passive system (α = 0), and when a small fraction (α = 0.01) of the particles is active. A series of snapshots from both scenarios is shown in Fig. 1. Note that in these snapshots, the cluster algorithm described in Section II. D was used to identify the crystalline domains in the system, and that particles identified as disordered are not depicted. As can be seen clearly from the time-series, after t = 350τ, the passive system still exhibits the same grain structure as in the initial configuration: the polycrystal is trapped in a metastable state and does not evolve to its single-crystal equilibrium state [Fig. 1(a–e)]. In contrast, when a small fraction of the particles is active, we observe substantial grain boundary mobility, leading to the rapid coarsening of domains over time [Fig. 1(f–h)]. In this case, when t = 350τ, we are left with a single-domain crystal in coexistence with a fluid. Thus, similar to the scenario in two dimensions, active dopants help the polycrystal evolve to a large-scale defect-free crystal. Moreover, again similar to the 2d case, upon turning off the activity [Fig. 1(i)], the fluid region recrystallizes using the existing crystal as a template. We are left with a single domain crystal with a large stacking fault. Note that this stacking fault is the result of the finite size of our periodic simulation box, as the crystal typically chooses a random orientation incommensurate with the box shape.

Fig. 1 Grain growth in a passive system (a–e) vs. a doped system (a and f–i) starting from initial configuration (a), showing significant coarsening of crystal domains in the doped system, compared to its purely passive counterpart. In the lower panel we activate a fraction α = 0.01 of particles for 0 < t/τ < 350 with a self-propulsion of fσ/k_BT = 90. After the activity of the self-propelled particles is switched off we can clearly see that the polycrystal has coarsened. The disordered particles are not depicted. The system contains N = 83 481 particles and the density is ρ_sσ³ = 0.86.

We now turn our attention to the mechanism responsible for the increased mobility of the crystal grain boundaries. To explore this mechanism, we study the dynamics of the active particles during the coarsening process. In Fig. 2 we show the time evolution of a slice of the crystal a short time after the activity is turned on. In this figure, disordered and crystalline particles are coloured grey and red, respectively, and the active particles are highlighted in black. We note that at t = 0, the active particles are randomly distributed throughout the system. However, it is clear from Fig. 2 that by t = 100τ almost all the active particles are located inside the grain boundaries, which have broadened into thick fluid layers. Particles located in these fluid layers exhibit high diffusion rates due to the presence of the active particles, which act like local “stirrers”. As time progresses, the domain in the center gradually melts from the outside inwards, while the adjacent domains grow. This continual surface melting and recrystallization of domains is responsible for the motion of the grain boundaries and leads to the formation of the large single crystal in Fig. 1(h).

Fig. 2 Activity-induced surface melting of polycrystalline domains. We color code disordered and crystalline particles with grey and red, respectively. Active particles are highlighted in black at many times their own diameter. (a–e) Over time we can clearly see that the central grain melts from the surface inwards, until it is completely removed. The active particle fraction is α = 0.01 with a self-propulsion of fσ/k_BT = 90. The system size is N = 83 481 and the density is ρ_sσ³ = 0.86.

The appearance of a substantial fluid region, as in Fig. 1(h), suggests that the presence of active particles shifts the phase boundaries in this system. To confirm this, we studied the coexistence between the crystal and fluid using simulations of N = 16 000 particles in an elongated box as discussed in Section II. C. In Fig. 3(a), we show an example snapshot, where we highlight active particles in red. Clearly, the system depicted in this snapshot has separated into a crystal containing (almost) no active particles and a fluid rich in active particles. To analyse this coexistence in detail we measure the density profile ρ(z) and active particle fraction α(z) along the long axis of the box. Specifically, we divide our system into slabs of length dz = a in the long direction of the box, with a the lattice spacing of the crystal. To calculate the active particle fraction and density, we then determine the mean Voronoi volume per particle in each slab. The resulting profiles, corresponding to the system depicted in Fig. 3(a), are shown in Fig. 3(b). Here, two distinct densities are visible, namely a high-density crystal phase largely free of active particles and a low-density fluid phase rich in active particles. One might further ask whether the active particles accumulate at the fluid-solid interface, as they are known to accumulate at walls.^30,31 However, the smooth decay in active particle fraction (red line in Fig. 3(b)) between the high and low density phases clearly indicates that the active particles do not specifically adsorb to the crystal-fluid interface.

Fig. 3 (a) Direct coexistence between a crystal and a fluid phase at overall density ρ_sσ³ = 0.80 and an active particle fraction α = 0.01 and self-propulsion fσ/k_BT = 50. (b) The corresponding density profile (blue) and the local fraction of active particles (red).

Density profiles, like the one shown in Fig. 3(b) can be used to map out the coexistence regions in our system. The results, for a range of overall system densities ρ_s, are summarized in Fig. 4. From this figure, we first note that for different overall densities the coexistence curves do not collapse in the ρσ³vs. fσ/k_BT representation when the fraction of active particles is kept fixed at α = 0.01. This feature arises from the fact that our system is in fact a binary mixture and hence should not collapse in this representation. We also see from Fig. 4 that for low activity, the system does not phase separate for any of the densities we examined. This was to be expected, since we always chose the overall system density to be higher than the predicted crystal coexistence density for the passive system (ρ_cσ³ = 0.785²⁸) as we were interested in removing grain boundaries from crystals. Importantly, we also see that both the crystal and fluid coexistence densities increase monotonically as a function of the activity of the dopants. Hence, when the activity is sufficiently strong, the melting density is increased beyond the overall system density, and as a result the system partially melts into a dopant-rich fluid and a nearly undoped crystal. We identify this partial melting as the driving mechanism behind the broadening of the grain boundaries observed in e.g.Fig. 2(a). Subsequently, these fluid regions, which display high diffusion rates due to the relatively high concentration of active dopants, enable the fast annealing of the grain boundaries.

Fig. 4 Coexistence densities as a function of self-propulsion fσ/k_BT, with a fixed active particle fraction α = 0.01 for three different overall system densities ρ_sσ³. Open circles indicate the absence of phase coexistence (pure crystal phase). Closed circles indicate the densities of the coexisting crystal phase, and triangles indicate the coexisting fluid densities. System size: N = 16 000.

Hence, we identify phase separation into a passive crystal and a fluid rich in active particles as the key mechanism behind active grain growth. To further illustrate this we will now return to the study of polycrystalline systems and follow the evolution of grain sizes over time for polycrystals of overall system density ρ_sσ³ = 0.86 upon the activation of a fraction of α = 0.01 of the particles. We now define a time scale t_a associated with the growth of a large domain by calculating the mean time required to form a crystal domain comprising at least a fraction X_c* = 0.4 of all the crystalline particles. This cutoff ensures that the domain has formed at the expense of at least three neighbouring domains. Note that we have checked that the results are robust to variation in X_c*. We plot the time scale of annealing t_a in Fig. 5. At low self-propulsion fσ/k_BT ≲ 30, in the region where no phase separation occurs according to Fig. 4, we observe no significant grain growth and thus t_a is far beyond our simulation time. Only at sufficiently high self-propulsion fσ/k_BT ≥ 30, where the system phase-separates, do we observe a huge increase in grain growth, leading to lower values of t_a. In this regime, we find that increasing the self-propulsion leads to faster annealing of polycrystalline domains. Depending on the self-propulsion f, we also observe that a large fraction of active particles remains kinetically arrested within the crystal phase.

Fig. 5 Time scale associated with formation of a large domain t_a and the fraction of mobile active particles p_m as a function of the self-propulsion f for a density ρ_sσ³ = 0.86.

To quantify the extent of dynamic arrest we calculate the fraction of mobile active particles p_m, defined as particles that have moved more than one lattice spacing Δr ≥ a within a time interval Δt = 2.5τ.³⁵ Clearly, at self-propulsions fσ/k_BT < 30 where we expect no phase separation, the fraction of mobile active particles p_m ≈ 0: all the active particles remain caged inside the domains. At the minimal amount of self-propulsion required for phase separation, fσ/k_BT ≈ 30, we observe a large increase in p_m, as the active particles manage to escape the crystal and create a fluid. However, near this increase only a fraction of active particles become mobile and participate in the active grain growth process. With increasing self-propulsion the fraction of active particles that partake in the process increases until it approaches unity, thus leading to lower values of t_a.

IV. Conclusions

In conclusion, we have explored how grain boundaries can be removed from three-dimensional colloidal crystals by doping with a small fraction of active colloids. We show that at sufficient self-propulsion, the system is driven into a crystal-fluid coexistence. In this case, active dopants become mobile and spontaneously gather at the grain boundaries. The resulting surface melting and recrystallization of domains is responsible for the motion of the grain boundaries and leads to the formation of a large single crystal. However, if the self-propulsion is too low to cause a phase separation, we observe no significant enhancement of grain growth.

Our study thus demonstrates that, similar to the 2d case,¹³ doping with active particles provides an elegant new route to removing grain boundaries in 3d colloidal polycrystals. In both cases, we observe that the active particles are attracted to the grain boundaries. This implies that the active particle autonomously find the location inside the sample where their active properties are most helpful. Therefore, the removal of grain boundaries via active doping is in stark contrast to more conventional methods based on external fields, which act uniformly on all the particles in the sample.^21–24

Interestingly, in the 2d case we observed that the active particles do not interact solely with the grain boundaries, but also exhibit strong interactions with many other types of crystal defects.¹³ This raises the interesting question of what role these other defects, such as for example vacancies and interstitials, will play in this system.

We also show that after a single large domain has formed, the activity can be switched off at any time to recrystallize the fluid regions, using the existing crystal as a template. The required time will depend strongly on the system (and system size) in question. Note that experimentally one may switch off the activity by using active particles that are responsive to certain stimuli (e.g. light-activated colloids) or simply by letting all the fuel be consumed.

Our results, combined with the recent observation of active particles aiding crystal nucleation in hard-sphere glasses,¹⁰ raise the question whether active dopants can be used in other, hard to equilibrate passive systems. For instance, it will be interesting to discover whether active dopants can help nucleate long sought-after crystals, such as a colloidal diamond phase,³² or the laves phases.³³

Acknowledgements

We acknowledge funding from the Dutch Sector Plan Physics and Chemistry, and funding from a NWO-Veni grant. We thank Vasileios Prymidis, Frank Smallenburg, Simone Dussi, and Siddharth Paliwal for many useful discussions and carefully reading our manuscript.

References

1. P. G. Debenedetti, Metastable liquids: concepts and principles, Princeton University Press, Princeton, 1996.
2. K. Kelton and A. L. Greer, Nucleation in condensed matter: applications in materials and biology, Elsevier, Oxford, 2010.
3. W. F. Paxton, et al. , J. Am. Chem. Soc., 2004, 126, 13424.
4. R. Dreyfus, et al. , Nature, 2005, 437, 862.
5. J. R. Howse, et al. , Phys. Rev. Lett., 2007, 99, 048102.
6. I. Theurkauff, et al. , Phys. Rev. Lett., 2012, 108, 268303.
7. I. Buttinoni, et al. , Phys. Rev. Lett., 2013, 110, 238301.
8. J. Palacci, et al. , Science, 2013, 339, 936.
9. C. Bechinger, et al., 2016, arXiv:1602.00081, arXiv preprint.
10. R. Ni, M. A. C. Stuart, M. Dijkstra and P. G. Bolhuis, Soft Matter, 2014, 10, 6609.
11. F. Kümmel, et al. , Soft Matter, 2015, 11, 6187.
12. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 018301.
13. B. van der Meer, L. Filion and M. Dijkstra, Soft Matter, 2016, 12, 3406.
14. A. Wysocki, R. G. Winkler and G. Gompper, 2016, arXiv:1601.00850, arXiv preprint.
15. V. L. Colvin, MRS Bull., 2001, 26, 637.
16. A. van Blaaderen, R. Ruel and P. Wiltzius, Nature, 1997, 385, 321.
17. K.-h. Lin, et al. , Phys. Rev. Lett., 2000, 85, 1770.
18. P. V. Braun, et al. , Adv. Mater., 2001, 13, 721.
19. M. Allard, E. H. Sargent, P. C. Lewis and E. Kumacheva, Adv. Mater., 2004, 16, 1360.
20. K. Jensen, et al. , Soft Matter, 2013, 9, 320.
21. T. Palberg, W. Mönch, J. Schwarz and P. Leiderer, J. Chem. Phys., 1995, 102, 5082.
22. R. Amos, et al. , Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 61, 2929.
23. S. Gokhale, et al. , Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 20314.
24. P. T. Korda and D. G. Grier, J. Chem. Phys., 2001, 114, 7570.
25. X. Sun, et al. , Langmuir, 2013, 29, 7216.
26. Z. Cheng, W. B. Russel and P. Chaikin, Nature, 1999, 401, 893.
27. J. Schiøtz, T. Vegge, F. Di Tolla and K. W. Jacobsen, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 11971.
28. L. Filion, R. Ni, D. Frenkel and M. Dijkstra, J. Chem. Phys., 2011, 134, 134901.
29. P.-R. Wolde, et al. , Faraday Discuss., 1996, 104, 93.
30. J. Elgeti and G. Gompper, EPL, 2013, 101, 48003.
31. X. Yang, M. L. Manning and M. C. Marchetti, Soft Matter, 2014, 10, 6477.
32. F. Romano, E. Sanz and F. Sciortino, J. Chem. Phys., 2011, 134, 174502.
33. A.-P. Hynninen, et al. , Nat. Mater., 2007, 6, 202.
34. Note that we have also performed simulations where the active particles are initially distributed randomly throughout the box, but over the course of the simulations they are excluded from the crystal region and migrate to the fluid region. Hence, starting in this manner reduces the time needed to reach the steady state.
35. We have checked that our results are robust to variations in both the displacement threshold as well as the time-interval.

---