Colloidal crystals are often prepared by evaporation from solution, and there is considerable interest to link the processing conditions to the crystal morphology and quality. Here, we study the evaporation-induced assembly of colloidal crystals using massive-scale nonequilibrium molecular dynamics simulations. We apply a recently developed machine-learning technique to characterize the assembling crystal structures with unprecedented microscopic detail. In agreement with previous experiments and simulations, faster evaporation rates lead to earlier onset of crystallization and more disordered surface structures. Surprisingly, we find that collective rearrangements of the bulk crystal during later stages of drying reduce the influence of the initial surface structure, and the final morphology is essentially independent of the evaporation rate. Our structural analysis reveals that the crystallization process is well-described by two time scales, the film drying time and the crystal growth time, with the latter having an unexpected dependence on the evaporation rate due to equilibrium thermodynamic effects at high colloid concentrations. These two time scales may be leveraged to control the relative influence of equilibrium and nonequilibrium growth mechanisms, suggesting a route to rapidly process colloidal crystals while also removing defects. Our analysis additionally reveals that solvent-mediated interactions play a critical role in the crystallization kinetics and that commonly used implicit-solvent models do not faithfully resolve nonequilibrium processes such as drying.
I. INTRODUCTION
Colloidal crystals with long-range order are highly desirable materials for applications such as photonics1 and optics.2 There is particular interest in developing inexpensive, high-throughput fabrication methods for such materials by self-assembly from solution.3 However, colloidal particles typically have weak interactions on the order of thermal energy, and competing crystal structures are stabilized against each other by only small differences in free energy,4–6 resulting in polycrystalline solids.7,8 Selectivity toward crystal structures has been enhanced by designing particle interactions with DNA coatings,9 patches,10,11 depletants,12,13 and shape.14,15 However, even such carefully engineered particles must still dynamically assemble into the desired thermodynamically stable crystal. The nonequilibrium aspects of self-assembly, including the assembly mechanism and processing conditions, therefore present important additional opportunities and challenges to control crystal structures.
One common high-throughput method for fabricating colloidal crystals involves dispersal of particles in a volatile solvent followed by evaporation of the solvent to deposit a crystalline solid onto a horizontal or vertical substrate.3 This assembly process is sensitive to, among many factors, the effective particle interactions16–18 and the evaporation rate as controlled by temperature,19,20 pressure,21–23 relative humidity,24 and solvent properties.25 For slow evaporation onto a horizontal substrate, crystals form in the bulk once the concentration is sufficiently high,26 whereas for fast evaporation, particles accumulate at the drying interface to form well-ordered monolayers27 or multilayers.20 Experiments20 and simulations28 have shown that the best crystal quality at the drying interface is obtained when drying is sufficiently fast to increase the local particle density but slow enough that diffusion anneals defects before additional grains or layers form. Much less is known about how three-dimensional crystals form, although it has been hypothesized that bulk crystal growth proceeds by accumulation and incorporation of particles below the surface monolayer.20
Computer simulations with time-resolved microscopic detail are useful tools to study how evaporation influences crystal formation and growth, which can be challenging to monitor experimentally. For example, computer simulations have shed light on how surfaces influence the nucleation of hard-sphere crystals at equilibrium (the slow drying limit).29,30 Recent simulations by Wang and Brady31 examined the microstructure of drying colloidal films and suggested that crystallinity may exhibit a minimum at intermediate drying rates due to a competition between crystallization at the interface and nucleation in the bulk. These simulations neglected hydrodynamic interactions between colloids, a typical approach16,17,31,32 for improving computational efficiency, but hydrodynamic interactions often play an important role in assembly dynamics, e.g., selecting binary crystal structures for DNA-linked colloids.33 To our knowledge, there has only been one previous simulation study of evaporation-induced colloid crystallization with full hydrodynamic interactions,28 but it was computationally restricted to 4 layers of colloids. The effects of evaporation and hydrodynamics on long-ranged three-dimensional crystal morphology thus remain largely unexplored.
We performed massive-scale nonequilibrium molecular dynamics simulations with an explicit-solvent model to study the evaporation-induced assembly of colloidal crystals from solution onto a horizontal substrate. A machine-learning technique was applied to characterize local structures in the growing crystal including particles at interfaces and defects, which are challenging to identify by traditional methods. Our simulations provide new mechanistic insight into evaporation-induced assembly of bulk crystals, which is found to occur in three steps: (1) accumulation of particles at the drying interface, (2) formation of a sufficiently ordered monolayer with subsequent bulk crystal growth templated by this surface structure, and (3) deformation of the bulk crystal due to confinement effects. We find that the dynamics of this process are well-described by the characteristic time scales of drying in step 1 and crystal growth in step 2. We also demonstrate the importance of the solvent in controlling the crystallization dynamics using complementary implicit-solvent simulations without hydrodynamic interactions. Ultimately, however, we find that the final bulk crystal structure does not depend significantly on the dynamics during drying for the range of parameters considered, as the eventual deformation of the crystal in step 3 disrupts initial surface ordering controlled by the evaporation rate.
The rest of the article is organized as follows. Section II presents the explicit-solvent and implicit-solvent simulation models and the machine-learning method used to characterize the forming crystal structures. Our results are discussed and analyzed in Sec. III, focusing on the crystal growth mechanism, influence of evaporation rate, and role of hydrodynamic interactions. Section IV summarizes our findings and discusses applications and extensions of our work.
II. SIMULATION MODELS AND METHODS
A. Explicit-solvent model
Following the model of Cheng and Grest, nearly hard spherical colloids were dispersed into a liquid solvent of Lennard-Jones particles in coexistence with its vapor.28,34,35 In this article, quantities are reported using the system of units defined by the diameter σ, mass m, and interaction energy ε of the solvent particles. The interaction potentials and parameters for our model are described in the appendix of Ref. 36 and more extensively in the supplementary material of this article. Briefly, the colloids were purely repulsive to each other, representing the excluded volume between particles, but were attracted to the solvent. The simulation box was periodic along the x and y dimensions with edge lengths Lx = 138.4 σ and Ly = 159.9 σ. This size and aspect ratio is roughly commensurate with the FCC (111) crystal face for our model based on the minimum distance between particles measured in the colloid pair distribution function. The finite-size tests of Cheng and Grest suggest that this cross sectional area is sufficiently large to limit simulation artifacts due to the size and geometry of the interface.28 The suspension was supported from below by a flat substrate at z = 0 σ that was attractive for the solvent and repulsive for the colloids. The solvent vapor was additionally confined by a purely repulsive potential.
We dispersed 2052 colloids with diameter d = 10 σ into the liquid solvent in coexistence with its vapor (2.3 × 106 solvent particles in total), giving a liquid film height H0 ≈ 215 σ (Fig. 1). The densities of the coexisting liquid and vapor before dispersing the colloids were 0.66 m/σ3 and 0.045 m/σ3, respectively, at temperature T = ε/kB with kB being Boltzmann’s constant. The initial volume fraction based on the colloid diameter was ϕ0 ≈ 0.24, which is comparable to previous simulation studies28,31 but larger than typical initial concentrations in experiments.20,26 We will show how using a higher initial concentration presents an important opportunity to control the crystal growth mechanism.
Snapshot of colloids (yellow) in the explicit solvent (blue) before the start of drying. The solvent thermostat and evaporation regions are indicated at the bottom and top of the simulation box, respectively, along with the initial film height H0. All snapshots in this article were rendered using Visual Molecular Dynamics (version 1.9.2).37 The explicit solvent is omitted from subsequent snapshots for clarity.
Snapshot of colloids (yellow) in the explicit solvent (blue) before the start of drying. The solvent thermostat and evaporation regions are indicated at the bottom and top of the simulation box, respectively, along with the initial film height H0. All snapshots in this article were rendered using Visual Molecular Dynamics (version 1.9.2).37 The explicit solvent is omitted from subsequent snapshots for clarity.
Evaporation was modeled by deleting randomly selected solvent particles at a constant flux j from the vapor in the top 20 σ of the box.28 Solvent particles within 20 σ of the substrate were weakly coupled to a Langevin thermostat38 at T = ε/kB with friction coefficient 0.1 m/τ where is the unit of time. Both regions are indicated in Fig. 1. We stopped the simulations before the crystal was completely dried (55% of solvent removed) to avoid spurious forces due to the neglect of gravity. All simulations were performed using HOOMD-blue (version 2.1.5) on multiple graphics processing units39–41 with a time step of 0.005 τ.
B. Implicit-solvent model
To test the role of hydrodynamic interactions in the crystallization dynamics, we constructed an implicit-solvent model for the colloids that matched the structural properties and long-time diffusion of the colloids in the explicit solvent, but neglected hydrodynamic interactions between them. The structural effects of the solvent were implicitly included through an effective colloid–colloid pair potential (Fig. S1) that was obtained by minimizing the relative entropy between the explicit- and implicit-solvent models.42 Relative entropy minimization ensures that this pair potential, which approximates the multibody colloid potential of mean force, quantitatively reproduces the colloid–colloid pair correlations.43 Colloid dynamics in the solvent were approximated by Langevin dynamics simulations,38,44 which incorporate the effects of thermal fluctuations and the free-draining Stokes drag on the colloids, but neglect hydrodynamic interactions between colloids.
A reference trajectory was prepared by dispersing 250 colloids into a cubic, periodic simulation box with edge length 100 σ (ϕ0 = 0.13) filled with a Lennard-Jones solvent at the liquid coexistence density and deleting any solvent particles within a distance of 5.5 σ of any colloid. After equilibration, colloid configurations were collected every 5 τ during a 75 000 τ simulation in the canonical ensemble. The effective colloid pair potential was truncated at 18 σ, which was commensurate with the first minimum of the colloid pair distribution function from the reference trajectory (Fig. 2). The potential was modeled using cubic B-splines with a knot density of 2 knots per σ, which ensured that the coarse-grained model was able to resolve all structural features in the pair distribution function. Knot values were then optimized to minimize the relative entropy with the reference trajectory using a combination of conjugate gradient and Hessian descent schemes.42 Since the minimization proceeds by stochastically initiating short trial molecular dynamics trajectories (15 000 τ), we verified that independent optimizations with different knot densities gave nearly identical potentials. The Langevin dynamics friction coefficient was then adjusted to 70 m/τ to match the diffusion coefficient measured from the average colloid mean-squared displacement in the reference trajectory. We validated our model by computing the pair distribution function and equilibrium diffusion coefficient of the colloids in the bulk over a range of concentrations, 0.05 ≤ ϕ0 ≤ 0.26 (Figs. 2 and S2), and found excellent agreement between the implicit- and explicit-solvent simulations.
Radial distribution function g(r) for explicit- (top) and implicit-solvent (bottom) models for colloids at volume fraction ϕ0 in a cubic box with edge length 100 σ at T = ε/kB.
Radial distribution function g(r) for explicit- (top) and implicit-solvent (bottom) models for colloids at volume fraction ϕ0 in a cubic box with edge length 100 σ at T = ε/kB.
During evaporation, the drying liquid-vapor interface was represented by the purely repulsive part of a harmonic potential starting with its minimum at H0 and moving down at a fixed speed.32 This harmonic potential, derived by Pieranski using surface energy arguments,45 effectively models compression due to the downward motion of the interface. A similar approach was recently employed by Wang and Brady, who represented the interface as a moving hard wall.31 The spring constant was chosen to be 5 ε/σ2 so that the colloids could fluctuate near the interface, but all remained below it. We found in preliminary tests that similar particle configurations were obtained with larger spring constants, while smaller spring constants gave an interface potential that did not sufficiently contain the colloids throughout the simulation.
C. Crystal characterization
The crystalline character of particles was assessed throughout the simulations using our recently developed neighborhood graph analysis (NGA) method.36 NGA uses the diffusion map dimensionality reduction technique46–48 to identify collective variables that characterize the local neighborhood topology of each particle (Fig. 3). We first used the averaged local bond order parameter49,50 to filter solid-like particles () for characterization by NGA. Particles were then assigned to structural classes based on the position of their neighborhood graph in the diffusion map compared to nine reference structures (eight of these are shown in Fig. 3). These reference structures were chosen based on high observation frequency (C, G), the appearance of large contiguous clusters (A, B, E, F), or both (D, H).
2D projection of the 3D diffusion map obtained from NGA as described in Ref. 36 and colored as in Ref. 51. The collective variable ψ2 describes the bulk (left) to surface (right) character, while ψ4 describes HCP-like (top) to FCC-like (bottom) character. Selected structures along the HCP-like (A-D) and FCC-like (F-H) branches are illustrated schematically, where colors indicate the 3D stacking layer: blue for bottom, orange for middle, and green for top.
2D projection of the 3D diffusion map obtained from NGA as described in Ref. 36 and colored as in Ref. 51. The collective variable ψ2 describes the bulk (left) to surface (right) character, while ψ4 describes HCP-like (top) to FCC-like (bottom) character. Selected structures along the HCP-like (A-D) and FCC-like (F-H) branches are illustrated schematically, where colors indicate the 3D stacking layer: blue for bottom, orange for middle, and green for top.
The diffusion map coordinates of each particle were projected onto lines between each pair of reference structures (j, k) and the nearest such projection (located a distance dproj from the original point) was selected. The particle’s identity was split between the reference structures according to the fraction (fj, fk) of the distance from the projection to each reference point. Each particle was further weighted by = exp(−dproj/0.05), which is related to the scaling of the Gaussian kernel in the diffusion map, in order to penalize structures that were not in the direct path between reference structures. The k-th channel in Fig. 7 then contains the total weight Wk = Σfk in each frame of the simulation, normalized by ΣWk. Particles were colored in the snapshots by assigning hues to each reference structure, interpolating between them, and desaturating according to the weight by mixing with gray.
III. RESULTS AND DISCUSSION
The explicit-solvent model computationally restricts the colloids and films that can be simulated to the nanoscale. Accordingly, we describe the drying process using dimensionless parameters that can be used to generalize our findings to larger scales. The well-known film Péclet number, Pe = H0/D0, characterizes the drying regime by the ratio of the rate of colloid accumulation at the liquid–vapor interface relative to the rate of diffusion,52–54 where H0 is the initial film height, is the interface speed, and D0 is the colloid diffusion coefficient. When Pe ≪ 1, diffusion is fast compared to accumulation and colloids are expected to be uniformly distributed in the film, eventually leading to nucleation in the bulk.26,31 When Pe ≫ 1, colloids are instead expected to accumulate at the liquid–vapor interface, leading to surface-induced nucleation.20,26,27 We first studied the mechanism of crystal growth for evaporative flux j = 3.6 × 10−4 m/σ2τ, which corresponds to Pe = 60,55 and then varied the flux a factor of two slower (Pe = 30) and faster (Pe = 121). We performed 5 independent simulations for Pe = 30 and 15 independent simulations for Pe = 60 and 121 due to the stochastic nature of assembly pathways.
Colloids were initially distributed uniformly in the film, but accumulated at the interface after drying began (Fig. 4). Heterogeneous nucleation of crystallites at the solvent liquid–vapor interface occurred in all our simulations once the surface layer became sufficiently dense [Fig. 4(a)]. The crystal front then moved down into the bulk liquid at a nearly constant speed [Figs. 4(b) and 4(c)]. Once all particles in the film crystallized [Fig. 4(d)], the crystal continued to compact as the solvent was dried further and the interparticle spacing decreased, eventually reaching a polymorphic close-packed structure [Fig. 4(e)]. The morphology of the bulk crystal was closely coupled with the structure at the liquid–vapor interface. For example, grain boundaries between face-centered cubic (FCC) and hexagonal close-packed (HCP) clusters in the first few layers often propagated to the supporting substrate, up to 15 layers away, in the form of stacking faults. Conversely, the removal of defects during later stages of drying sometimes led to buckling or cracking of the surface [Fig. 4(f)]. Such rearrangements considerably complicated the prediction of the final morphology, and we found no clear correlations between the initial monoloayer [Fig. 4(a)] and final bulk [Fig. 4(f)] morphology.
Snapshots from representative simulation trajectories at Pe = 60 for the explicit-solvent (top) and implicit-solvent (bottom) models. Snapshots in column (a) are at t = 37 500 τ when the uncrystallized film height is Hx (see Fig. 5), and subsequent snapshots [(b)–(f)] are shown every 25 000 τ. Particles are colored according to their positions in the diffusion map presented in Fig. 3. The colors can be roughly mapped onto the following structures: HCP (red), FCC (dark blue), FCC (111) surface (light blue), FCC (100) surface (teal), surface defect (orange/yellow), and weakly ordered surface (yellow/green). White particles correspond to .
Snapshots from representative simulation trajectories at Pe = 60 for the explicit-solvent (top) and implicit-solvent (bottom) models. Snapshots in column (a) are at t = 37 500 τ when the uncrystallized film height is Hx (see Fig. 5), and subsequent snapshots [(b)–(f)] are shown every 25 000 τ. Particles are colored according to their positions in the diffusion map presented in Fig. 3. The colors can be roughly mapped onto the following structures: HCP (red), FCC (dark blue), FCC (111) surface (light blue), FCC (100) surface (teal), surface defect (orange/yellow), and weakly ordered surface (yellow/green). White particles correspond to .
A. Evaporation rate
In order to assess how the evaporation rate affects crystallization kinetics, the position of the crystal interface was tracked using , which we also found to be a good filter for solid-like particles for NGA (Sec. II C). The crystal interface with the disordered colloids in the film was identified as the lowest z coordinate, initially neglecting colloids near the substrate, where the average in the plane exceeded our threshold for solid-like character (). We used a bin size of 10 σ, corresponding to the colloid diameter, to calculate the average and interpolated between bins. This definition of the crystal interface [Fig. 5(a)] is closely related to the total fraction of particles that were considered solid-like [Fig. 5(b)], but appears to lag slightly behind in Fig. 5(a) due to the need to accumulate multiple solid-like particles into a layer to identify an interface. Unfortunately, we did not have sufficient time resolution or statistical sampling in our simulations to study the nucleation process at the interface and so instead focused our analysis on the crystal growth after the first layer formed.
(a) Solvent liquid–vapor interface position (dashed lines) and colloid crystal interface position (solid lines) over the course of the simulations at different evaporation rates. Each line is the average of all trajectories with the standard deviation indicated by the shaded bands. Hx is the height of liquid film below the crystal interface at first bulk crystal appearance. is the solvent evaporation rate and is the crystal front speed (see Table S1). (b) Fraction of solidified particles with .
(a) Solvent liquid–vapor interface position (dashed lines) and colloid crystal interface position (solid lines) over the course of the simulations at different evaporation rates. Each line is the average of all trajectories with the standard deviation indicated by the shaded bands. Hx is the height of liquid film below the crystal interface at first bulk crystal appearance. is the solvent evaporation rate and is the crystal front speed (see Table S1). (b) Fraction of solidified particles with .
The average near the interface trended from initially liquid-like toward crystalline values prior to our identification of the first crystal layer. For fixed initial concentration, the time to form the first sufficiently thick layer for bulk-like growth is expected to be inversely proportional to the evaporation rate (see Table S1 for values). As a consequence, the first crystalline layer formed at roughly the same dried film height. Likewise, the height of the remaining uncrystallized film at the time this layer formed, Hx, was also nearly identical for all three evaporation rates, as marked in Fig. 5(a).
By contrast, we found that the speed of the crystal front was less sensitive to the drying rate. For evaporation-driven growth, the crystal front speed was expected to be proportional to , with ≈ /(1 − ϕ0/ϕx) according to a flux balance56 and ϕx being the local crystal volume fraction. Instead, varied by only a factor of two compared to a four-fold increase in [Fig. 5(a) and Table S1]. In all cases, was greater than . This suggests a two-step mechanism in which accumulation first leads to nucleation at the surface after reaching a sufficiently high density and then templated crystal growth from the surface follows. We propose that this crystal growth results from a combination of accumulation to the surface and additional thermodynamically favored crystallization due to the high initial volume fraction. We confirmed this by drying the crystal for a short time (50 000 τ, drying roughly 10% of the film) to form a surface nucleus and then stopping the evaporation, finding that the crystal continued to grow. This finding suggests that the fluid-like phase was thermodynamically metastable with respect to the crystal phase. Such metastability may be surprising on the basis of the estimated initial volume fractions using the bare colloid diameter; however, the effective volume occupied by a particle is larger than its bare volume due to solvent structuring near the surface (see Fig. 2). This thermodynamically favored growth may be connected to the unexpected annealing behavior reported by Cheng and Grest.28 Unlike assembly onto carefully designed crystal templates,57 our simulations exhibited multiple surface structures forming together, resulting in polymorphic crystal growth.
Previous simulations suggested that there is an optimal evaporative flux to obtain high-quality in-plane crystalline order,28 resulting from a balance between in-plane diffusion and accumulation to the interface. Diffusion at the interface must be sufficiently fast compared to accumulation into the interface to promote growth of large crystal grains. The particle Péclet number, Pep = d/D0, characterizes these processes as the ratio of the rate of incorporation of a particle into the interface by drying compared to the rate at which it diffuses in plane by its own diameter d. When Pep ≫ 1, accumulation is fast compared to diffusion and a more disordered surface structure is expected, whereas for Pep ≪ 1, diffusion dominates and a more ordered surface structure with larger grains should form. Note, however, that letting Pep become arbitrarily small will ultimately result in a less crystalline interface since the film Péclet number Pe = Pep(H/d) must still be sufficiently large to accumulate particles at the interface. These considerations thus suggest the existence of an intermediate optimal evaporative flux, as also previously proposed by Cheng and Grest.28
For our model, the slowest drying rate gives Pep = 1.4, while for the fastest drying rate Pep = 5.6, suggesting that the faster drying rates should result in poorer quality surface structures. Indeed, when we compared surface structures early in the drying process, we found that typically fewer, larger grains were formed at slower evaporation rates, consistent with this model (Figs. S3–S5). However, the early surface structure underwent significant rearrangements at later drying times after a sufficient number of crystal layers formed below. Figure 6 compares the evolution of the surface structure for two independent simulations at Pe = 121 (Pep = 5.6). Initially, one surface is more disordered than the other [Fig. 6(a)]. As the drying progresses, both surfaces undergo rearrangements that increase the surface grain size [Fig. 6(b)] and ultimately both become well-ordered single grains [Fig. 6(c)]. Due to such rearrangements, initial surface order did not correlate with overall crystal quality.
Snapshots of the interface structure for two independent simulations at Pe = 121 for times (a) 25 000 τ, (b) 56 250 τ, and (c) 88 750 τ. Colors are the same as in Fig. 4. Snapshots of all interfaces are provided in Figs. S3–S5 of the supplementary material.
Snapshots of the interface structure for two independent simulations at Pe = 121 for times (a) 25 000 τ, (b) 56 250 τ, and (c) 88 750 τ. Colors are the same as in Fig. 4. Snapshots of all interfaces are provided in Figs. S3–S5 of the supplementary material.
No statistically significant differences in the final crystal structures were seen from multiple independent simulations at the different drying rates. Essentially all characteristic morphologies could be obtained at all rates (see Figs. S6–S8 for snapshots from all 35 simulations). The effects of the evaporation rate are well-described at early times by the film Péclet number Pe, controlling the onset of crystallization and associated with diffusion normal to the interface, and the particle Péclet number Pep, affecting the number of grains that form at the surface and associated with in-plane diffusion. However, deformation of the interface at late stages of drying limits the usefulness of initially high-quality surface structures. This suggests a novel route for rapid processing of colloidal crystals. Starting from a sufficiently high colloid concentration, a period of fast drying is used to accumulate a crystalline layer at the interface. Templated growth from this surface is then driven by both the drying process and equilibrium thermodynamics. Finally, when the crystal front reaches the substrate before all solvent is removed, as in our simulations, additional drying deforms the crystal and removes defects inherent to the fast growth process.
B. Structural evolution
The local descriptors obtained by NGA shed additional insight into the balance of equilibrium and nonequilibrium effects in the crystallization mechanism. From initially liquid-like structures, particles proceeded toward an ordered state through one of two branches corresponding to HCP-like or FCC-like symmetry, beginning with weak hexagonal ordering and adopting increasingly higher coordination numbers until reaching a bulk crystal state. To visualize this progression, we divided each branch into four reference structures with increasing coordination, illustrated in Fig. 3. We then recorded the fraction of solid-like particles closest to each reference for each simulation snapshot, shown in Fig. 7. For the first part of the trajectory, the relevant time scale for comparing the drying rates is the film drying time τd = H0/, which characterizes the formation of a crystalline surface layer. The evolution of structures prior to this event, particularly A and E, is indistinguishable for the three different evaporation rates when time is scaled by τd.
Average fraction of the solid-like particles in each of eight different neighborhood topologies (i.e., local crystal structures) as a function of time. Values are normalized at each point in time. Schematic representations of the corresponding crystal structures are shown in Fig. 3. Vertical white lines indicate transition from τd to τx scaling. Black arrows indicate equivalent film heights, as shown in representative snapshots at right.
Average fraction of the solid-like particles in each of eight different neighborhood topologies (i.e., local crystal structures) as a function of time. Values are normalized at each point in time. Schematic representations of the corresponding crystal structures are shown in Fig. 3. Vertical white lines indicate transition from τd to τx scaling. Black arrows indicate equivalent film heights, as shown in representative snapshots at right.
After the surface layer forms (white lines in Fig. 7), the transitions from surface structures (A, E) toward bulk structures (D, H) are better compared using τx = Hx/, the characteristic time of crystal growth. For example, the transition times between structures C and D are best collapsed between evaporation rates when the time after growth of the surface layer is scaled by τx instead of τd. The inability of τd to correctly describe the time evolution of structures is even more apparent in the sample snapshots shown in the right-hand side of Fig. 7 (see Figs. S9–S11 for additional snapshots), which are taken at equal dried film heights (equal reduced drying times). At this film height, the fastest dried film (Pe = 121) is only roughly 50% crystallized, while the slowest dried film (Pe = 30) is mostly crystallized because it has effectively had more time for growth. We note that there are some dissimilarities in the structural evolution near the end of the trajectories (e.g., different fractions of structure H). These small discrepancies are largely due to interactions with the substrate and subsequent deformation of the crystal, which are not coupled to either time scale. This analysis provides strong indication that in the Pe > 1 regime, changes to the evaporation rate affect the onset of crystallization, but the progression of the crystal morphology is only similar when comparable growth times incorporating both evaporation-induced and equilibrium growth are considered.
C. Influence of solvent
Having assessed the crystal growth mechanism and the effects of the evaporation rate, we now consider the role of the solvent itself in controlling the crystallization. The solvent has both a structural effect on the colloids influencing the preferred particle spacing and a dynamic effect mediating hydrodynamic interactions. The typical approach employed in most prior drying studies16,17,31,32 is to treat the structural effects implicitly while also neglecting hydrodynamic interactions to expedite the simulations and to model larger length and time scales. However, this approximation remains untested for evaporation-induced crystallization, and there has been no direct comparison to assess how solvent-mediated hydrodynamic interactions influence crystal growth. To this end, we constructed an implicit-solvent model matching the effective colloid interactions in the explicit solvent and used Langevin dynamics simulations38,44 to capture the long-time colloid diffusion but ignore hydrodynamic interactions between colloids (Sec. II B). Drying simulations were then conducted using this implicit-solvent model.
We started the implicit-solvent simulations from the same initial colloid configurations as in the explicit-solvent simulations, but found that the crystallization dynamics were significantly different between the two models (Fig. 4) despite good agreement for equilibrium properties. Particles accumulated at the liquid–vapor interface faster for the explicit-solvent model [Fig. 4(a)], resulting in earlier onset of crystal growth [Fig. 4(b)]. Additionally, the explicit-solvent model crystallized exclusively down from the liquid–vapor interface, whereas the implicit-solvent model also crystallized up from the substrate [Fig. 4(c)]. This combined growth led to a more disordered crystal when the two growth fronts met in the middle [Fig. 4(d)]. We noted that the implicit-solvent model did not undergo surface buckling as early as the explicit-solvent model [Fig. 4(e)], which might be expected from our approximation of the liquid–vapor interface as a plane.32 This simplified interface model also neglects capillary interactions between colloids, which may contribute to their in-plane assembly at the interface.58,59 The final morphologies only appeared similar once both models underwent significant rearrangements during compression [Fig. 4(f) and Figs. S12 and S13], consistent with the lack of correlation between the initial and final structure in the explicit-solvent model.
Hydrodynamic interactions in a suspension influence the particle motion and hence distribution in the drying film. Colloid diffusion is driven by forces due to chemical potential gradients.52 Brady has shown in the context of phoretic motion of a colloid that hydrodynamic interactions between particles reduce the velocities that result from these forces.60 Assuming comparable dependence of the chemical potential on concentration in both models, slower particle migration is consistent with a higher colloid density near the interface, resulting in earlier onset of crystallization.
There is also a temperature gradient in the explicit-solvent model due to evaporative cooling that is absent from the Langevin dynamics simulations. At the pseudosteady state, the temperature profile within the drying film is predicted to be linear in distance from the substrate with a slope proportional to the evaporative flux, in good agreement with the simulations. We also found that the temperature gradient was larger for the colloid suspensions than for the neat solvent for Pe = 60. The temperature near the liquid-vapor interface was nearly 25% lower than at the substrate for Pe = 121. The temperature gradient may alter the effective colloid interactions near the interface. Thermophoretic motion may additionally increase the colloid density near the liquid–vapor interface, as was recently demonstrated for a colloid mixture in an explicit solvent.61 The temperature gradient may accordingly increase the driving force for crystallization near the liquid–vapor interface and contribute to faster growth from this surface in the explicit-solvent model.
IV. CONCLUSIONS
We analyzed the mechanism of evaporation-induced colloidal crystallization by extensive nonequilibrium molecular dynamics simulations and characterized the crystal structures using a template-free machine-learning technique. In agreement with previous findings, faster evaporation rates led to earlier onset of crystallization with more disordered surface structures. Surprisingly, the final crystal structure did not depend significantly on the evaporation rate for the range of drying conditions investigated because the initial surface structures deformed to remove defects from the bulk crystal. This finding suggests a key difference between crystals grown by sedimentation, where colloids accumulate onto a rigid substrate, and crystals grown by evaporation, where the interface can deform to relieve defects.
Detailed microscopic structural analysis revealed that the initial regime of the assembly process was well-described by the time scales of film drying and crystal growth. Unexpectedly, the characteristic time scale for crystal growth was influenced by both the evaporation rate and thermodynamically favored crystallization at the studied initial colloid volume fraction. Our results suggest an unconventional route for rapid processing of colloidal crystals: after formation of a surface layer to initiate growth, a fundamentally nonequilibrium process, crystal growth toward the substrate can then be driven by equilibrium thermodynamic considerations if the colloid concentration is sufficiently large. Evaporation both drives crystal growth and removes defects by collective rearrangements as the solvent is removed. The initial colloid volume fraction can then be exploited concurrently with the evaporation rate to manipulate the onset and relative rates of these equilibrium and nonequilibrium growth mechanisms.
We were able to resolve bulk-like crystal structures with many more crystal layers than previous simulations.28 Nonetheless, the present study was effectively restricted to nanometer sized particles in a thin film at moderately high concentrations by the explicit solvent. High film Péclet numbers were achieved by fast evaporation rates, which are estimated to be several orders of magnitude larger than those obtained for typical solvents in experiments, such as water. It would be interesting to extend our analysis to larger particles and films. We expect many of the results obtained here to be applicable when the processing parameters are scaled according to the Péclet number. However, it is likely that a larger initial volume fraction (based on the bare colloid diameter) would be required to achieve a metastable fluid phase with respect to the crystal phase as solvent structural effects should become increasingly negligible with larger colloid diameters. Hydrodynamic interactions are expected to continue to play an important role even for larger colloids and films.
Accessing the size and time scales associated with these larger systems is currently infeasible for the employed explicit-solvent model. The computational challenges posed by the explicit-solvent model would typically be overcome through the use of an implicit-solvent model neglecting hydrodynamic interactions. However, we showed that such a simplified model can yield stark differences in the crystal growth kinetics during evaporation despite good agreement for bulk equilibrium structural and dynamic properties. We conclude that implicit-solvent models lacking hydrodynamic interactions do not faithfully resolve the dynamics of nonequilibrium processes such as drying, in agreement with a recent theoretical analysis of drying colloidal mixtures62 and simulations of drying polymer mixtures.63 An appropriate mesoscale simulation method incorporating hydrodynamic interactions seems the most promising for tackling such problems in the future.
SUPPLEMENTARY MATERIAL
See supplementary material for interaction potentials and parameters for the explicit-solvent model, effective potential for the implicit-solvent model, diffusion coefficients of explicit and implicit models, average evaporation and crystallization speeds, and additional snapshots from all simulated trajectories.
ACKNOWLEDGMENTS
The simulations presented in this article were performed on computational resources supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University. M.P.H. and A.Z.P. acknowledge financial support from the National Science Foundation through Grant Nos. DMR-1420541 (Princeton Center for Complex Materials, a Materials Research Science and Engineering Center) and CBET-1402166. W.F.R. received government support under Contract No. FA950-11-C-0028 and awarded by the Department of Defense, Air Force Office of Scientific Research, the National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a. A.N. was supported by the German Research Foundation under Project No. NI 1487/2-1. T.S. and M.S.S. gratefully acknowledge support from the National Science Foundation through Grant No. CHEM-1300770.
REFERENCES
Initially, D0 = 1.95 × 10−3 σ2/τ and ≈ j/ρ = 5.5 × 10−4 σ/τ with ρ = 0.66 m/σ3 being the solvent liquid density.