On the influence of a patterned substrate on crystallization in suspensions of hard spheres Sven Dorosz and Tanja Schilling 2 1 Theory of Soft Condensed Matter, Universit´e du 0 2 Luxembourg, L-1511 Luxembourg, Luxembourg n a J Abstract 1 3 We present a computer simulation study on crystal nucleation and growth in supersaturated ] suspensions of mono-disperse hard spheres induced by a triangular lattice substrate. The main t f o result is that compressed substrates are wet by the crystalline phase (the crystalline phase directly s . t a appears without any induction time), while for stretched substrates we observe heterogeneous m nucleation. The shapes of the nucleated crystallites fluctuate strongly. In the case of homogeneous - d n nucleation amorphous precursors have been observed (Phys. Rev. Lett. 105(2):025701 (2010)). For o c heterogeneous nucleation we do not find such precursors. The fluid is directly transformed into [ 2 highly ordered crystallites. v 0 4 PACS numbers: 82.60.Nh, 64.60.Q-,64.60.qe, 64.70.pv,68.55.A- 4 5 . 8 0 1 1 : v i X r a 1 When a supersaturated fluid crystallizes, crystallization is usually induced by the container walls, rather than to proceed from a fluctuation in the bulk of the system. This effect, called heterogeneous nucleation, is of fundamental importance for the kinetics of phase transitions (such as the formation of ice in the supersaturated vapor of clouds), as well as for techno- logical applications, in which the properties of the walls can be designed to influence the properties of the crystals that are formed. In this article we discuss heterogeneous crystal nucleation and growth from the overcompressed fluid of hard spheres. Hard spheres have served successfully as a simple model system for fluids and crystals over the past fifty years. The interaction energy between two hard spheres is either infinite (if they overlap) or zero (if they do not overlap), thus the phase behavior of the model is purely determined by entropy. The simplicity of the potential makes hard spheres particularly suited for computer simulations; and the entropic nature of the phase transition makes them a useful limit case for comparison to other systems, which are governed by an interplay between entropy and enthalpy. Hard spheres are not only of interest to the theoretician, they are also often synthesized on the colloidal scale and used in experiments on fundamental questions of statistical mechanics (see e. g. [1] and references therein). As the topic of our work, crystallization of hard spheres on a substrate, has been studied experimentally [2–7] and theoretically [8–16] before, we briefly lay out in the following, which aspects of this topic have been focused on in the articles cited above. The supersaturated fluid of hard spheres in contact with a planar hard wall has been ad- dressed in computer simulation studies by Dijkstra [9], Auer [8] and Volkov [10]. These studies show that the planar hard wall is wet by the crystalline phase, hence crystallization proceeds layer by layer rather than by the nucleation of crystallites. (For a review on wet- ting and film growth of crystalline phases on structured and unstructured surfaces in various systems, including hard spheres, see the article by Esztermann and Lo¨wen [15].) Also the recent experimental and simulation work by Sandomirski and co-workers [7] dealt with the growth of a crystalline film in contact with a wall. Here the wall was not planar but a 2 fcc layer of spheres. The authors found that the speed of the crystallization front depends non-monotonically on the packing fraction of the fluid and that a depletion zone is present in front of the growing crystal. Heterogeneous nucleation of hardsphere crystals hasmainly been addressed inthe context of template-induced crystallization. Van Blaaderen and co-workers [3–5] showed how to design structured templates to induce the epitaxial growth of large monocrystals and of metastable phases in a sedimenting liquid of hard spheres. Cacciuto and Frenkel studied the effect of finite templates of various sizes and lattice structures on crystallite formation by means of computer simulation [12]. Recently this topic was taken up again and investigated in more detail experimentally and theoretically by the groups of Dijkstra and van Blaaderen [6]. For small two-dimensional seeds of triangular as well as square symmetry they find that nucleation barriers depend on the seed’s symmetry as well as the lattice spacing. This effect is due to defects and changes in crystal morphology that are induced by the seed. Heterogeneousnucleationofhardspheres onaninfinitesubstratehasrecently beenaddressed by Xu and co-workers [11] in a computer simulation study. In this work triangular and square substrates as well as a hcp(1100) pattern were brought in contact with a strongly overcompressed fluid, and the evolution of the density profile perpendicular to the substrate as well as the fraction of crystalline particles were monitored. A metastable bcc-phase that was stabilized for long times was observed. Here we present an extended simulation study of crystallization mechanisms and rates for a fluid of hard spheres brought in contact with a triangular substrate for varying overcom- pression and lattice distortion. To our knowledge there is no systematic study on the effect that distortion of an infinite substrate lattice has on the crystallization mechanism and rate of hard spheres. We would like to close this brief overview by pointing out that there are other useful model systems for crystal nucleation, as e. g. complex plasmas. In contrast to colloidal systems microscopic dynamics in complex plasmas are almost undamped, [17] hence they offer a complementary experimental approach to the topic. 3 I. SETUP OF THE SYSTEM AND SIMULATION DETAILS The simulations were carried out by means of an event driven molecular dynamics program for fixed particle number, volume and energy (for details on event driven MD see refs. [18– 21]). We simulated N = 216,000 hard spheres of diameter σ in contact with a substrate of triangular symmetry formed by N = 4200 spheres of the same diameter σ. The substrate particles were immobile (i.e. they had infinite mass). The simulation box had periodic boundaries in x and y directions. The substrate layers were fixed at z = ±Lz for L = 2 z 30σ...50σ, depending on the overcompression. The initial velocities were drawn from a Gaussian distribution and the initial mean kinetic energy per particle was set to 3 k T. B To monitor crystallinity, we used the local q6q6-bond-order parameter [22, 23], which is defined as follows: For each particle i with n(i) neighbors, the local bond-orientational structure is characterized by n(i) 1 q¯6m(i) := XY6m(~rij) , n(i) j=1 where Y (~r ) are the spherical harmonics with l = 6. ~r is the displacement between 6m ij ij particle i and its neighbor j in a given coordinate frame. A vector ~q (i) is assigned to each 6 particle, the elements m = −6...6 of which are defined as q¯ (i) 6m q (i) := . (1) 6m 6 |q¯ (i)| 1/2 (cid:0)Pm=−6 6m (cid:1) We counted particles as neighbors if their distance satisfied |~r | < 1.4σ. Two neighboring ij particles i and j were regarded as “bonded” within a crystalline region if q~(i)·q~(j) > 0.7. 6 6 We define n (i) as the number of “bonded” neighbors of the ith particle. (In the online b version we use the following colour-coding for the snapshots: if a particle has n > 10, i.e. b an almost perfectly hexagonally ordered surrounding, it is color-coded green, if n > 5 it is b color-coded brown.) We studied various densities between particle number density ρ := Nσ3/V = 1.005 (which corresponds to a volume fraction η = 0.5262) and ρ = Nσ3/V = 1.02 (η = 0.5341). At 4 these densities the chemical potential difference per particle between the metastable fluid and the stable crystalline state is between ∆µ ≃ −0.5 k T and ∆µ ≃ −0.54 k T. The B B overcompressed fluid configurations did not show pre-existing crystallites that might have been created during the preparation process. y 1,02 tilib a ts f o ρ instantaneous growth tim ity 1,015 il s n e d r e b m 1,01 u n nucleation 1,005 1 1,1 1,2 1,3 1,4 lattice constant a [σ] FIG. 1: Representation of all combinations of density ρ and substrate lattice constant a studied in this work. The limit of stability of the homogeneous bulk crystal is indicated by the solid line (green online). At substrate lattice constants smaller than this value (squares) we find complete wetting of the substrate and instantaneous film growth. Systems with a larger substrate lattice constant (circles) exhibit incomplete wetting and heterogeneous nucleation up to a ≤ 1.5σ. Above this stretching, no heterogeneous nucleation event was observed on the scale of the simulation time. Figure 1 shows the densities ρ and substrate lattice constants a (of the fcc-(111) plane) for which we carried out simulations. The lattice constant indicated by the solid line (green online) corresponds to the bulk crystal at the spinodal, i.e. at the density at which the crystal ceases to be metastable with respect to the liquid. We obtained this density by simulation as well as from density functional theory [24]. The corresponding lattice constant is a = 1.15σ (DFT) resp. a = 1.14σ (simulation). One result of our study is that this sp sp 5 line separates the parameter space into regions of different crystallization mechanisms. For a < a , we observed the instantaneous formation of a film, which then grew with time. For sp a > a , the system crystallized via heterogeneous nucleation. The transition between the sp two mechanisms seems to be continuous. For a ≥ 1.5σ no heterogeneous nucleation event was observed on the scale of the simulation time. II. COMPLETE WETTING OF THE SUBSTRATE Forall compressed substrates (a < a )weobserved theformationandgrowthofacrystalline sp film. Typical snapshots are presented in figure 2. (Here, we chose a system at a = 1.1σ, close to a , and a bulk density of ρ = 1.01.) The timescale of the MD simulation is expressed in sp multiples of τ = σ2, with D being the long-time self diffusion coefficient in the bulk fluid 6D obtained in the same MD simulations. In the regime of densities analyzed, the diffusion constant varies by only 5%.) FIG. 2: Snapshots t = τ (left) and t = 100τ after bringing the overcompressed fluid in contact with the substrate, a = 1.1σ (slightly less than a ), ρ = 1.01. Only crystalline sp particles are shown (n > 5) b In order to analyze the crystalline layers quantitatively, we computed the 2-dimensional 6 bond-order parameter ψ for planes perpendicular to thez-direction. (ψ is the 2dequivalent 6 6 of q¯ .) 6 n(i) 1 ψ6(i) := Xei6θij , n(i) j=1 where θ is the angle of the vector ~r and an arbitrary but fixed axis in the plane. We ij ij ∗ impose a cut-off at |~r | < 1.4σ and demand for a crystalline particle that ψ (i)ψ (j) > 0.7 ij 6 6 for six neighbors. To discuss the analysis in detail, we pick three substrate lattice constants a = {1.01σ,1.05σ,1.1σ} at a fixed density ρ = 1.005. 0.008 t = 100 τ t = 50 τ 0.006 t = 1 τ ) z ( ρ y 0.004 t i s n e d 0.002 0 0 5 10 15 distance to the substrate z [σ] FIG. 3: Density profile perpendicular to the substrate for half of the system at different times. ρ = 1.01, a = 1.1σ. A film of layers grows. Figure 3 shows a vertical density profile. As a function of time the layering becomes more pronounced, as seen from the growth of the maxima and the appearance of voids in between thelayers. (Aquantitative analysis ofthegrowthratefordifferent substrate latticeconstants is not reported, because the lateral dimension was too small.) According to these profiles we identify the particles that belong to a given layer and study the hexagonal structure in the 7 plane. The overall defect density η in a given layer n with a total numer of N(n) particles is defined as N(n)−N crys(n) η(n) := , (2) N(n) with N (n) being the number of crystalline particles in layer n. The analysis of the defect crys density isshown infigure4. Wehave also included thetotalnumber of particles N(n) ineach layer n for the three cases of a. The further the substrate is compressed with respect to the equilibrium latticethe largeris thedefect density inthefirst layer. Withlarger distance from the substrate the defect density for all three values of a converges to a substrate independent value. At this point stresses induced by the substrate do not play a role in the growing crystal anymore. Only the tension induced by the shape of the periodic box, which is not commensurate with the equilibrium lattice, matters. 0.5 4200 ) n ( 0.4 N4000 η r y ye t a nsi n l de0.3 s i3800 fect ticle e r d a p 0.2 f 3600 a=1.01σ o a=1.05σ # a=1.10σ 0.1 3400 1 2 3 4 5 1 2 3 4 5 layer n layer n FIG. 4: (left) Defect density as a function of the index of each layer counted from the substrate for three different substrate lattice constants. The data shown has been obtained in the long time limit t > 400 τ and it is averaged over three independent runs each. (right) Number of particles N(n) in each layer n. Figure 5 shows the covering of the substrate for the first three layers after t = 400 τ. There is nopreference offccover hcp. Ananalysisofthesubsequent layers showed thatthestacking is 8 random-hcp. This is in agreement with the small free energy difference of 26±6·10−5k T/σ2 B per particle [25]. Domains of equal structure are much larger for the case a = 1.1σ than for a = 1.01σ, where there are more domain walls. No single crystal phase evolved on the recorded timescales. FIG. 5: Snapshot of the first three layers on top of the substrate for (left) a = 1.01σ and (right) a = 1.10σ. The snapshots correspond to the data analyzed in figure 4. There is no preference of fcc over hcp. III. HETEROGENOUS NUCLEATION NEAR THE SUBSTRATE For the parameter regime 1.15σ ≤ a ≤ 1.4σ we observe the formation of crystallites at the substrate. Figure 6 shows snapshots of typical crystallites at the first nucleation event (figure 6a) and at a much later time (figure 6b). We define the nucleation event as the moment when the first crystalline cluster reaches a size of 100 particles, see figure 6a for a snapshot. In all simulations we observed irreversible growth above this threshold. Below this threshold crystallites appeared and decayed again. Changing this value by ±10 particles does not affect any of the results presented in the 9 (a) t=50τ (b) t=150τ FIG. 6: Snapshots at different times after bringing the overcompressed fluid in contact with the substrate, a = 1.4σ, ρ = 1.01. Crystallite formation at the wall dominates the nucleation process. For clarity, we are not showing the substrate. Figure (a) shows the nucleation event at which the first crystallite reaches 100 solid particles. Figure (b) shows the state of the system at a much later time. following. ρ=1.01 r10000 e ρ=1.0175 ust ρ=1.02 cl n s i e cl arti 1000 p d oli s f o # 1000 100 τ 200 τ 300 τ 400 τ time after critical nucleus is formed FIG. 7: Time evolution of size of the largest cluster for varying density ρ. The data is averaged over 8 independent runs at each given density. 10