Computer Simulations of Crystal Growth Using a Hard-Sphere Model

Abstract

A review of computer simulation studies on crystal growth in hard-sphere systems is presented. A historical view on the crystallization of hard spheres, including colloidal crystallization, is given in the first section. Crystal phase transition in a system comprising particles without bonding is difficult to understand. In the early days, therefore, many researchers did not accept such crystalline structures as crystals that should be studied in the field of crystal growth. In the last few decades, however, colloidal crystallization has drawn attention because in situ observations of crystallization process has become possible. Next, simulation studies of the crystal/fluid interface of hard spheres are also reviewed. Although colloidal crystallization has now been recognized in the crystal growth field, the stability of the crystal–fluid coexistence state has still not been satisfactorily understood based on a bond-breaking picture, because of an infinite diffuseness of the interfaces in non-bonding systems derived from this picture. Studies of sedimentary colloidal crystallization and colloidal epitaxy using the hard-sphere model are lastly reviewed. An advantage of the colloidal epitaxy is also presented; it is shown that a template not only fixes the crystal growth direction, but also improves the colloidal crystallization. A new technique for reducing defects in colloidal crystals through the gravity effect is also proposed.

1. Introduction

Bonds are commonly formed between various entities in solid materials. In contrast, a class of matter called “soft matter” does not have any bonding. A colloidal system is a typical example of soft matter. The presence of large entities is one of the characteristics of soft matter. Therefore, because of their large size, such entities have slow motion; unlike atomic systems, in situ observation of the crystal growth process at the particle level is possible in colloidal systems.

To get insight into the particle-level mechanism in atomic systems, researchers have to rely on computer simulations. One cannot only follow the rapid motion of atoms in reality and simultaneously look at the atoms in crystal growth processes. These difficulties are not faced with computer simulations such as molecular dynamics (MD) and Monte Carlo (MC) simulations. However, these simulation methods have a limitation in total computation time. For example, a simulation for atomic systems can be conducted for several microseconds at most. In contrast, because colloidal particles have slower motion than atoms, the simulation of colloidal particles can be performed for several days. In other words, simulations corresponding to real situations are possible for soft matter.

This paper is a review on computer simulations using a hard-sphere (HS) model for crystal growth. Information about the structure of the crystal/melt interface at the particle level is necessary to understand phenomena and develop techniques. For example, the relationship between the interface structure and the growth mode, and the colloidal crystallization defect behavior that depends on the interface structure are related to the quality of colloidal crystals. Computer simulations of the HS crystal/fluid interfaces were developed in the 1990s. To some crystal growth researchers, construction of a stable HS crystal/fluid interface in simulations is strikingly similar to the discovery of the crystalline phase in the HS system itself. If the difficulty of arranging HS particles at high number densities could be circumvented, an HS crystal–fluid coexistence state can be simulated starting from a dense HS configuration as an initial state. HS simulations have been further developed in the last few decades. The direction of studies on HS simulations coincides with that of colloidal crystals; a large number of studies have been conducted on controlling defects in colloidal crystals for functionalization.

After briefly introducing various aspects of the HS model, a description of the simulation methods is provided. Simulation of the HS crystal/fluid interface has been provided in Section 3. Simulations of HS systems in a gravitational field have also been discussed in Section 4. Crystal–fluid coexistence state under an external field is also a subject of this review. In the latter part of the review, the effect of different simulation techniques on colloidal crystallization has been extensively investigated.

1.1. Crystalline Phase in a Hard-Sphere System

In 1957, the crystalline phase of the HS system was reported by means of computer simulations [1,2]. Formation of face-centered cubic (fcc) structures by HSs has also been reported [3]. In the early days of the investigation of HS systems, researchers were surprised to see the occurrence of a crystalline phase in a system without any attractive interaction between the particles [4]. That is, it was natural for the researchers of the time to consider bond order as the reason for the occurrence of ordered phases. In general, competition between the configurational entropy and the potential energy of attractive interactions causes phase transition. The free energy of a disordered phase becomes small because of the entropic effect at high temperatures. In contrast, a crystal phase with bonding between atoms becomes energetically stable at low temperatures. This picture leads to the understanding that if there are no interparticle attractions, phase transitions into an ordered phase do not occur. Hence, bond formation can be considered to result from attractive interactions. Crystalline phase transition in the HS system is sometimes termed the Kirkwood–Alder–Wainwright transition (or simply as the Alder transition). Kirkwood was the first to predict this transition [5]. At present, we have an intuitive picture of the Alder transition based on the competition between the two entropic effects (i.e., configurational and vibrational entropies). At low particle number density, the effect of configurational entropy dominates that of vibrational entropy. In contrast, at high density, the vibrational entropy becomes predominant in the ordered phase of the system—even if the configurational entropy is lost, resulting in the total entropy increase. This picture is valid for systems comprising repulsive interactions only; i.e., the Alder transition in a wide sense can be understood in this manner.

Among the researchers who did not believe in the concept of the Alder transition, Onsager was an exception. He had already predicted the existence of a nematic phase in a hard-rod system [6]. He observed that ordering would more easily take place in systems with broken isotropy than in those with spherical particles. Nowadays, computer simulations of hard-rod systems for liquid crystals have been developed as an important branch of the soft matter field. In particular, different conclusions have been reported on the stability of columnar phase several times during the 1980s–1990s, which is both exploratory and interesting (e.g., [7]). The present author believes that the discussion on this development will help in understanding the role of the packing of particles in inducing crystalline order. However, this subject is not a topic of the present review.

Competition between attractive forces and thermal motion is responsible for the gas–liquid phase transition. Therefore, gas and liquid phases are indistinguishable in systems comprising repulsive interactions only, such as supercritical fluids. In terms of the crystal growth, the phase diagram of the HS system drawn on the density axis can be divided into three regions: the fluid region below ${\varphi }_f$, crystal region above ${\varphi }_s$, and a two-phase coexistence region between them; ${\varphi }_f$ and ${\varphi }_s$ are, respectively, the volume fractions of the coexisting fluid and crystal. Here, the density of the HS system is expressed in terms of its volume fraction, $\varphi$ ≡ $(\pi {\sigma }^3/6)(N/V)$, where $\sigma$ is the HS diameter, N is the number of particles, and V is the total volume. The phase diagram in the P–T plane can be divided into two parts by a straight line: the region above $P/T$ = ${(P/T)}_c$ is the crystal region, and the region below it is the fluid region. In 1968, Hoover and Ree applied the single-occupancy cell method and obtained the following: ${\varphi }_f$ = 0.494, ${\varphi }_s$ = 0.545, and ${(P{\sigma }^3/k_{B}T)}_c$ = 11.75 (where $k_B$ is Boltzmann’s constant) [8]. These values were then revised by Davidchack and Laird in 1998 [9]. A direct two-phase coexistence simulation was performed, and the phase transition condition was corrected downward through minimization of stress in the crystal: ${\varphi }_f$ = 0.491, ${\varphi }_s$ = 0.542, and ${(P{\sigma }^3/k_{B}T)}_c$ = 11.55. Recently, the following values have been obtained by constant-pressure MC simulations: ${\varphi }_f$ = 0.492, ${\varphi }_s$ = 0.545, and ${(P{\sigma }^3/k_{B}T)}_c$ = 11.576 [10], while those obtained by MC simulations with a devised umbrella sampling technique are as follows: ${\varphi }_f$ = 0.49188, ${\varphi }_s$ = 0.54312, and ${(P{\sigma }^3/k_{B}T)}_c$ = 11.5727 [11]. It should be noted that the phase transition densities are far lower than the close-pack density (${\varphi }_{cp}=\pi \sqrt{2}/6\cong 0.74$).

Note also that Jackson’s theory of the crystal/melt interface was presented in 1958 [12]. Although his model—the so-called two-level model of the crystal/melt interface—was a mean-field model based on different types of bonds (i.e. solid–solid, melt–melt, and solid–melt bonds), it was successfully applied to various types of materials. This indicates the versatility of the mean-field model. However, many researchers working on crystal growth believe that a bonding picture is essential [13].

1.2. Colloidal Crystals and Effective Hard-Sphere Model

Similar to HS systems, colloidal particles do not form any bonding between the particles. Therefore, crystal growth researchers face a conceptual difficulty when working on colloidal crystals. A colloidal crystal was first reported in 1954 [14]. However, this was not related with the HS phase transition because the colloidal crystals investigated were dried (i.e., they were essentially close-packed crystals). Lux et al. conducted a Bragg diffraction study in the visible light region for colloidal dispersions, the density of which corresponds to that obtained by the Alder transition [15]. A standard theory that is used to describe interparticle interactions between colloidal particles is the DLVO theory [16,17]. In the initial days of research on colloidal crystals, attempts were made to understand colloidal crystallization in the framework of the DLVO theory as a phenomenon in which the particles were trapped in the first minimum of interparticle potential. However, the distance to the first minimum was too small to explain the lattice spacing of colloidal crystals. In 1972, Wadach and Toda successfully explained colloidal crystallization by means of Alder transition [18]. They treated a colloidal sphere of diameter $\sigma$ as an effective HS of diameter ${\sigma }_{eff}$ = $\sigma$ + $\alpha {\kappa }^{-1}$, with $\alpha$ being a factor of order unity, where $\kappa$ was the Debye parameter. Parameter ${\kappa }^{-1}$ is sometimes also referred to as the Debye screening length, which can be regarded as the thickness of the electric double layer around a colloidal particle. In this respect, the HS model was an idealized model for charged colloids. Experimental efforts have been made to realize HS colloids [19,20,21]. Nowadays, colloidal systems exhibiting the HS crystallization behavior are extensively studied using sterically stabilized colloids such as those reported in References [22,23,24,25,26].

A large number of studies on colloidal crystals have been conducted because colloidal crystals—in principle—possess properties similar to those of a photonic crystal. Structures with a specially ordered variation of dielectric constant—where the periodicity of the special order is on the same order as that of the wavelength of light—exhibit a photonic band [27,28,29]. Because the periodicity of colloidal crystals is on the same order as the wavelength of light, a photonic band structure arises which is similar to the band structure of electrons in semiconductors. In contrast with photonic crystals fabricated by nano-manufacturing, colloidal photonic crystals are more cost-effective, because the equipment for colloidal crystallization is cheaper than that for nano-manufacturing. In addition, since colloidal crystallization is a self-assembling phenomenon, it is less time-consuming to form a three-dimensional (3D) structure. However, a shortcoming of colloidal crystallization is the controlling of defects. Colloidal epitaxy is one of techniques that is used to reduce defects [22]. A key idea of the colloidal epitaxy is the uniqueness of stacking sequence in fcc (001) stacking. Growth direction is limited in fcc [001] through the use of a patterned substrate (i.e., a template). A recent development is the use of a template with a defect structure in order to introduce the desired structure in a colloidal crystal [26].

Some studies have also examined the effect of gravity on the defects for defect control. Zhu et al. performed a space shuttle experiment and found that a sediment of HS colloids in microgravity forms a random hexagonal close pack (rhcp) structure [23]. Traversing along the fcc 〈111〉 stacking sequence occurs in a three-fold manner in the regular fcc structure (where the fcc {111} planes can be classified into three with respect to the particle positions parallel to the plane—say, A, B, and C), while it is two-fold in the regular hexagonal close-packed (hcp) structure. That is, ABCABC⋯ stacking sequence corresponds to fcc, and ABAB⋯ corresponds to hcp. In turn, random sequences of A, B, and C correspond to rhcp. Zhu et al. also reported that the effect of gravity reduces the stacking disorder. Small free-energy differences between fcc and hcp crystals account for the difficulty in excluding the stacking disorder. Entropy difference between fcc and hcp crystals has been computed for HS systems, and it was found to be—at most—on the order of ${10}^{-3}{k}_B$ per particle [30,31,32,33,34,35]. In addition, fcc–hcp interfacial free energy was calculated to be on the order of ${10}^{-5}{k}_{B}T$ [33].

2. Simulation Methods

Unlike systems in which particles interact via interparticle potential of a continuous function, some difficulties are faced in hard-particle systems, because the potential function possesses some singular points. In this case, in an MD simulation, the equations of motion cannot be replaced with difference equations when performing the simulation. Differential of potential diverges at the contact of two hard particles. Section 2.1 describes the simulation algorithm for the HS systems. The MC simulation algorithm does not suffer too much from the singularity of potential. In Section 2.2, after a general explanation, the method used for MC simulation of HSs under gravity—on which the present author also worked—is described. Although Brownian dynamics simulation is more suitable for the simulation of suspended particles, this review does not explain it in detail because the present author has not worked on this method. Colloidal particles are subjected to random forces from the particles of dispersion media, in addition to the interparticle force between colloidal particles themselves. Therefore, the Langevin equation governs the time evolution of this system. The random forces should be treated together with the difficulty of potential singularity in the Brownian dynamics simulation.

2.1. MD Method for Hard Spheres in General

Alder and Wainwright presented the algorithm through which they tracked particle coordinates during time evolution by solving collision dynamics between HSs [36]. The review describes the collision dynamics of a system of non-identical HSs; the HS diameter and mass of a particle i are ${\sigma }_i$ and $m_i$, respectively. Let $r_i\left(t\right)$ and $v_i\left(t\right)$ be the position and velocity of particle i at time t, respectively. For all pairs i and j, Equation (1)

$$|r_{ij}+v_{ij}{t}_{ij}|={\sigma }_{ij},$$

(1)

is solved by squaring both its sides to obtain a list {$t_{ij}$} of times after which a collision between particles will occur. Here, ${\sigma }_{ij}$ ≡ $({\sigma }_i+{\sigma }_j)/2$ is the distance between particles i and j at collision. The following equation is obtained

$$t_{ij}=\frac{-b_{ij}-\sqrt{b_{ij}^2-v_{ij}^2(r_{ij}^2-{\sigma }_{ij}^2)}}{v_{ij}^2},$$

(2)

where the relative position vector and relative velocity vector are expressed as $r_{ij}\left(t\right)$ ≡ $r_i\left(t\right)-r_j\left(t\right)$ and $v_{ij}\left(t\right)$ ≡ $v_i\left(t\right)-v_j\left(t\right)$, respectively, and $b_{ij}$ ≡ $r_{ij}·v_{ij}$. An inappropriate solution to Equation (1), $t_{ij}$ = $[-b_{ij}+\sqrt{b_{ij}^2-v_{ij}^2(r_{ij}^2-{\sigma }_{ij}^2)}]/v_{ij}^2$, is excluded. In addition, positive $b_{ij}$ corresponds to a collision that occurred in the past. After making the list {$t_{ij}$} and sorting it, the time is evolved to $t^{\prime }$ = t + $t_{kl}$, where $t_{kl}$ is the minimum value in the members of {$t_{ij}$}. All positions are updated as $r_i\left(t^{\prime }\right)$ = $r_i\left(t\right)$ + $v_i\left(t\right)t_{kl}$. Velocities remain unchanged, except for particles k and l. From linear and angular momentum conservation at the collision, $v_k$ and $v_l$ are updated as

$$\begin{array}{ccc}\hfill v_k\left(t^{\prime }\right)& =& v_k\left(t\right)-\frac{m_l}{m_k+m_l}\frac{b_{kl}}{{\sigma }_{kl}^2}{r}_{kl}\left(t^{\prime }\right),\hfill \\ \hfill v_l\left(t^{\prime }\right)& =& v_l\left(t\right)+\frac{m_k}{m_k+m_l}\frac{b_{kl}}{{\sigma }_{kl}^2}{r}_{kl}\left(t^{\prime }\right).\hfill \end{array}$$

(3)

The list {$t_{ij}$} is updated by subtracting $t_{kl}$, except for the members involved in the collided particles. The times to collision with particles k and/or i are recalculated using Equation (1). One can solve the time evolution by repeating this procedure.

However, this algorithm cannot be parallelized. Let us make a simple consideration. Consider a case where {$t_{ij}^{\left(1\right)}$} is the list of times to collision for system 1 and {$t_{ij}^{\left(2\right)}$} is the list for system 2; to obtain a list for a system comprising systems 1 and 2, the two lists should be combined. This list can be taken as the list for the composed system if the collision at the interface between the systems is neglected. The average number of collisions doubles until a specific event for the combined system occurs. In other words, the mean time between successive collisions would decrease by half. It can be said that for a large system, a collision would occur at almost any instant. Such a difficulty does not arise for continuous potential, for which the equation of motion is replaced with a difference equation with a non-zero constant time step.

The collision-by-collision method has been developed as an event-driven method (e.g., [37]). At a glance, computational cost for a system with pair-wise interaction is proportional to $N^2$, where N is the number of particles (or system size). However, when the interaction is of short range (e.g., $r_{int}$), the computational cost can be greatly reduced. One may determine interactions within $r_{int}$. A method for this purpose is the linked-cell method [38]. In a usual linked-cell method, the system is divided into different numbers of cells whose size is larger than $r_{int}$. For a particle, the interactions that need to be calculated are limited to within the cell this particle belongs to and its neighboring cells. For the HS system, the cell size is set small so as to include one particle at the most, and in addition to collisions between particles, events through which a particle crosses the cell boundaries are also investigated [39]. This method has been further improved by Isobe [40]. In contrast to Alder and Wainwright’s algorithm, this scheme enables parallelization. An event-driven simulation package has been recently presented in [41].

2.2. MC Method Employed for Hard Spheres under Gravity

MC simulation for the HS system is very simple. If all interparticle distances are larger than or equal to the HS diameter ($r_{ij}$ ≥ ${\sigma }_{ij}$ for non-identical HSs) after an MC move, this MC attempt is accepted. If an overlap between HSs occurs, this attempt is rejected. The present author and coworkers employed the MC method for the simulation of HSs in gravity in order to investigate the sedimentation of HSs on a bottom wall. In this case, in addition to the interparticle interaction, interaction between particles and the wall should be considered. In an MD simulation, particle–wall collision event was also investigated. This is not a very difficult task. However, the MC method was employed because an algorithm to reject a configuration that a particle locates outside the wall was overwhelmingly simple compared to that for investigating collision events between a particle and the wall in the presence of gravity. In the presence of a gravitational field, the effect of gravitational energy should be incorporated. The Metropolis method was employed to investigate the change in the gravitational energy associated with the MC move.

3. Hard-Sphere Crystal/Fluid Interface

It is difficult to construct a crystal/fluid interface in HS systems because of the hard-core singularity. The particle density is relatively high in the simulations of a crystal/fluid interface or an interface between two condensed phases. In the HS system, the volume fraction of HSs in the coexistence region is about 0.5. When crystalline and fluid phases are prepared separately with appropriate coexistence densities and then joined together to create an interface, an overlap would occur between HSs because of the high density. To avoid such an overlap, crystal and fluid blocks should be kept slightly apart from each other or the overlapped particles should be removed. With both methods, the density of the entire system would be less than that of the coexistence one. The method to circumvent this difficulty is described in Section 3.1. Some results of interface simulations are shown in Section 3.2

3.1. Crystal–Fluid Coexistence in Equilibrium

In 1995, the present author and coworkers for the first time successfully performed an MD simulation of crystal–fluid coexistence states in an HS system [42]. It was also in 1995 that Kyrlidis and Brown performed a similar study by MC simulation [43]. To prepare a configuration at a high density, the simulation should be usually started with an fcc configuration, even if one aims at simulating a dense liquid. Overlap between HSs (or contact at a steep repulsive core, in general) occurs for a dense random distribution of particle centers. A crystalline configuration at a density corresponding to a dense liquid is well equilibrated to obtain a random configuration. Equilibration is also key in crystal–fluid coexistence simulations.

A direct method to avoid an overlap between HSs at the interface does not involve separate preparation of crystalline and fluid phases; instead, it involves the preparation of a coexistence state in a system from the first in a course of simulations. Here, one should start with a regular configuration because of the high density. The melting of a regular configuration takes a long time. The present author and coworkers used the following method to circumvent this problem [42]: at first, the system was divided into three parts—the crystal block, melt block, and interfacial region. An fcc crystal of $\varphi$ = ${\varphi }_s$ was located in the crystal block, and a close-packed fcc crystal was present at the center of the melt block. We determined the width of the interfacial region ($d_{int}$) from the following relation:

$$\frac{\pi {\sigma }^3}{6}\frac{N}{(L_{melt}+d_{int}+L_{cryst})A}={\varphi }_{entire},$$

(4)

where N is the total number of particles, $L_{melt}$ and $L_{crys}$ are, respectively, the lengths (normal to the interface) of the melt and crystal blocks, A is the cross-sectional area of the system (parallel to the interface), and ${\varphi }_{entire}$ is the volume fraction of the entire system, which fell around the center of the ${\varphi }_f$ < $\varphi$ < ${\varphi }_s$ region. That is, $d_{int}$ was calculated before determining the initial configuration of the particles. A snapshot is shown in Figure 1. The close-packed fcc configuration was placed in the melt block in order to obtain a high collision rate between HSs. The simulation was further modified in order to accelerate the randomization of the melt configuration. In Reference [42], the mass of crystal particles was set 1000 times heavier than that of melt particles. Equation (3) shows that heavier particles scarcely move, whereas lighter particles move well. Accordingly, the melt particles adopted a disordered configuration after several MD runs, as shown in Figure 2a. Figure 2 shows snapshots for an N = 3787 system. It took $5\times {10}^5$ total collisions to obtain the configuration shown in Figure 2a. Note that the 3D periodic boundary condition (PBC) was imposed on the configuration. Therefore, the right and left crystals form one block. After the melt block became molten, we made all the masses identical. A snapshot after equilibration of the entire system is shown in Figure 2b. Figure 2b shows a trajectory of the particles after $2\times {10}^6$ total collisions during the simulation of identical masses. This trajectory was constructed by connecting the particle positions observed for over 70 collisions per particle. The numerals indicated on the horizontal axis correspond to the labels of the crystal layers. The mass ratio can also be set to a value other than 1000, such as 100 or 10000. Optimization of the mass ratio would result in acceleration of the equilibration. When the mass ratio was moderate (such as that reported in [42]), the crystal particles showed slight movement. The crystal block will therefore move a little toward equilibrium in an equilibration process of the melt block.

This method was developed by Davidchack and Laird [9]. They initially moved only the melt particles in the equilibration process of the melt block, and then moved only the crystal particles. These processes were repeated. Equilibration of the interface configuration is confirmed by observing the evolution of configuration in repeated processes—in particular, when the change in the switching of the motionless particles was observed. As mentioned in Section 1.1, the revision of the value of the HS phase transition condition was achieved by their well-equilibrated crystal–fluid coexistence state. This direct coexistence method has attracted more attention than the other methods for determining the phase equilibrium condition, as described in Reference [44]. Accurate knowledge of the equilibrium condition would help in determining the crystal growth condition. Distances from phase boundaries to the state point of the crystal growth condition in a phase diagram is a useful piece of information. Direct coexistence simulations also afford atomistic information about the interface, as described in the subsequent subsection.

3.2. Interface Properties

Researchers working on crystal growth classify the interface structure into two types: a sharp interface and a diffuse interface [45]. In the former (such as during crystal growth from vapors), essentially a layer-by-layer growth mode takes place. In the latter (such as for crystal growth from a melt), the crystal growth events occur simultaneously over a number of layers. The experimental observation of an interface between dense phases is difficult because of low accessibility to the interface through a dense phase. Atomistic information from computer simulations complements such an experimental observation.

The information on the diffuseness of an interface is essential to determine whether or not the multilayer mode can occur. Density profiles of the HS crystal/fluid interfaces along the interface normal are shown in Figure 3. The density profile shown in Figure 3b corresponds to the snapshot depicted in Figure 2b. The highlighted numeral “22” is an end position of the crystal block (Figure 1b). The end positions are also highlighted for (100) and (111) interface systems. These interfaces had widths several layers thick and compact interfaces, similar to those in the Lennard–Jones and soft-sphere systems [46,47,48,49,50,51,52,53]. While crystal growth researchers were surprised by the bonding picture (as described in Section 1.1), realization of crystal/fluid interfaces in soft-sphere systems (inverse 12th-power potential [46,50,51] and Weeks–Chandler–Andersen potential [48]) was a convincing proof for the present author that the HS interfaces can be constructed successfully by an MD simulation. If the width of the HS interface was infinitely broad (e.g., inversely proportional to the bond-breaking energy as calculated by Temkin’s multi-layer model [54]; the numerical result shown in References [13,55,56,57] implies this property, and the analytical order parameter profile under continuum approximation [58,59,60] provides evidence), it would become difficult to grow crystals in a two-phase coexistence state via an interface. In other words, when the width of the interface is limited by the system size, the interface structure cannot be controlled through thermodynamic parameters only.

Layer-by-layer investigation of the interfacial region was performed. The trajectories of the particles (not shown in this review) were essentially the same as those reported by Davidchack and Laird [9], except that the system was relatively small. The number of particles used in the simulations reported in Reference [42] were 3888, 3788, and 4230 for (100), (110), and (111) interfaces, respectively. The sizes of the cross-sections parallel to the interface were $9.41\sigma \times 9.41\sigma$, $9.40\sigma \times 9.97\sigma$, and $9.79\sigma \times 9.60\sigma$, respectively. It should be remarked that, with respect to the intra-layer order, several layers were liquid crystalline. That is, despite the definite layered structure, the structure of several layers was far from crystalline with respect to the intra-layer order. Radial distribution function ($g^{\left(2D\right)}\left(r\right)$) was calculated for each layer. The 2D projections of the position vectors parallel to the interface were achieved for each layer and then a histogram of distances between all pairs of projection vectors was obtained. One can find liquid crystalline properties in these profiles. Results for the left interface (for the (100) interface simulation, both the right and left interfaces were almost symmetric) are shown in Figure 4. The second peak at r ≅ $2\sigma$ is indicative of the short-range order of a dense liquid. When two HSs are in contact with each other, the interparticle distance equals the HS diameter (i.e., r = $\sigma$), and when three HSs are in contact with each other on a line, two particles at both ends of the line give r = $2\sigma$. Such configurations often appear for random arrangements of particles at a high density. The peak at r ≅ $2\sigma$ is confirmed in Figure 4a. In Figure 4b–d, a peak at r ≅ $1.6\sigma$ is shown, which reflects the lattice structure. The second-neighbor distance in an fcc (100) plane (a square lattice) is $\sqrt{2}$ × $l_{nn}$, where $l_{nn}$ is the nearest distance that is slightly greater than $\sigma$, because the HS phase transition density is 0.5 in volume fraction. Along with the third peak at r ∼ $2.5\sigma$, a shoulder at r ∼ $2.2\sigma$ can also be seen in Figure 4d; however, the corresponding peak becomes very weak (as shown in Figure 4c). The third peak is attributed to the third-neighbor distance of the square lattice. While fine structures retain highly-ordered lattice, slight disorders such as those caused by thermal vibration broaden the peak, and sometimes it becomes difficult to distinguish shoulders from the main peak. Besides the smectic-A like structure (Figure 4a), a smectic-B like order has also been suggested.

Drift of the interface was significant for both (110) and (111) interfaces (and thus, the right and left interfaces were asymmetric). This review does not discuss radial distributions for (110) and (111) interfaces because of poor statistics. In addition, diffusivity profiles were calculated from mean square displacements. Because of the small size of the system, statistics for mean square displacements could not be improved. Nevertheless, anisotropy in diffusivity was evident for the mean square displacements of particles—in particular, in the left 18th layer corresponding to Figure 4b. The diffusion coefficient in the parallel direction to the interface was greater than that in the normal direction.

The other crucial parameter is interfacial free energy. The nucleation rate is calculated based on the nucleation work employing the interfacial free energy (e.g., Reference [61]). The equilibrium shape of a crystallite is determined by the orientational dependence of interfacial free energy, also known as the Gibbs–Curie–Wulff rule. Davidchack and Laird developed an interface simulation to calculate the HS crystal/fluid interfacial free energy [62]. In 1986, Broughton and Gilmer computed the interfacial free energy for the Lennard–Jones system by introducing a cleavage potential to calculate the work of formation of the interface [63]. A cleavage wall, instead of the potential, was used to calculate the formation work of the interface. Values published in Reference [62] were ${\gamma }_{100}^{*}$ = 0.62 ± 0.01, ${\gamma }_{110}^{*}$ = 0.64 ± 0.01, and ${\gamma }_{111}^{*}$ = 0.58 ± 0.01 for {100}, {110}, and {111} interfaces, respectively. Here, the asterisk indicates the reduced unit [$k_{B}T/{\sigma }^2$] (i.e., ${\gamma }^{*}$ ≡ ${\sigma }^2\gamma /k_{B}T$). These values were revised in 2005 [64], and recently revised again by Davidchack as ${\gamma }_{100}^{*}$ = 0.5820 ± 0.019, ${\gamma }_{110}^{*}$ = 0.5590 ± 0.020, and ${\gamma }_{111}^{*}$ = 0.5416 ± 0.031, and added ${\gamma }_{120}^{*}$ = 0.5669 ± 0.020 [65]. For comparison with values determined by other methods, see References [11,66].

4. Hard-Sphere System under Gravity

Earlier in Section 1, it has been mentioned that competition between the potential energy (potential well) and thermal energy gives rise to condensation. Thus, in systems comprising pure repulsive interaction, condensation does not occur. The HS system has only one type of characteristic energy, which is the thermal energy. However, if the system is kept in an external field, another characteristic energy can arise. It is expected that in a colloidal system in gravity (acceleration due to gravity g), a competition between the thermal energy and gravitational energy gives rise to phenomena like phase transitions. Hence, in such cases, it is convenient to introduce a dimensionless parameter $g^{*}$ = $mg\sigma /k_{B}T$. This parameter is sometimes called the “gravitational number” and coincides with the Péclet number, which is defined as the ratio between the drift and diffusion. Parameter $g^{*}$ is also expressed in terms of gravitational length, $l_g$ ≡ $k_{B}T/mg$ (i.e., $g^{*}$ = $\sigma /l_g$. If kinetics is the main concern, the concept of the Péclet number becomes more appropriate. Throughout this paper, we shall use $g^{*}$.

Colloidal suspensions crystallize by sedimentation [67]. This is a typical well-known method of colloidal crystallization, which is as intuitive as colloidal crystallization through the evaporation of dispersion media. Hence, it can be said that densification, in general, is the driving force of the Alder transition. The colloidal crystallization by sedimentation can in turn be understood as a phenomenon caused by the competition between $mg\sigma$ and $k_{B}T$. Because of phase transitions caused by competing factors, it is expected that various phenomena will also occur in the presence of gravity. An instance of such a phenomenon is the effect of gravity that reduces the stacking disorder in crystalline colloidal sediment, as reported by Zhu et al. [23]; this has been discussed at the end of Section 1.2.

4.1. Sedimentation

HS colloidal crystals obtained by settling under normal gravity are a mixture of fcc and hcp [68]. Contrarily, Zhu et al. reported that under microgravity, HS colloidal sediments were of rhcp structure [23]. This discovery cannot be understood intuitively, because the difference in the stacking sequence does not alter the particle number density. A large number of studies have investigated this phenomenon. It was shown by Kegel and Dhont that the rhcp structure, on aging, changes to the faulted-twinned fcc structure under milligravity [24]. The same trend was also obtained by Martelozzo et al. [69]. In contrast, Dolbnya et al. reported a slightly different result: after 1–1.5 years of gravitational annealing, the sediment acquired the fcc structure, including the absence of any stacking disorder [25].

MC simulation study of HSs in gravity were first carried out in 1994 [70]. However, the system sizes were very small, and accordingly, the number of stacking layers was also small. A main concern was the intra-plane order in gravity. Simulation for a large system was performed in 2002 [71]. An intra-plane structure was investigated with respect to the nucleation dynamics at the bottom. It was found that the effect of gravity obstructs the crystal growth from the point of view of nucleation dynamics. This seemingly contradicts the effect of gravity that reduces defects. It was in the mid-2000s that the present author and his coworkers performed computer simulations of HSs to investigate the stacking order in gravity [72,73]. MC simulations with constant volumes were conducted, and the strength of gravity was changed stepwise so as to avoid the trapping of systems in a metastable state such as glassy or polycrystalline states [73]. A series of snapshots is shown in Figure 5. The box length of a system containing N = 1664 particles was $L_x$ = $L_y$ = $6.27\sigma$ and $L_z$ = $49.23\sigma$. Bottom (z = 0) and top (z = $L_z$) walls of this system were hard walls and the PBC was imposed in the x and y directions. After preparing a random configuration, $g^{*}$ was increased up to 1.5 with increments of $\mathrm{\Delta }{g}^{*}$ = 0.1 and then decreased. The system was maintained at each $g^{*}$ for $2\times {10}^5$ MC cycles. Here, an MC cycle (MCC) was defined in a way that an average of one MC particle move per particle was included in one MCC. The maximum displacement of MC moves was fixed at $0.06\sigma$. It is confirmed from Figure 5 that defects in the crystal grown on the bottom wall dramatically vanished at $g^{*}$ ≅ 0.9. Hence, it can be said that we successfully demonstrated that defects are reduced under the effect of gravity. This $g^{*}$ value corresponds to the condition in which the gravitational energy that moves a particle about one interparticle distance becomes comparable to the thermal energy, $k_{B}T$. In other words, the gravitational length is comparable to the interparticle spacing. However, this is a one-particle picture. This gives only a primitive explanation for the gravitational stabilization of the layered structure near the bottom wall. Detailed observation of defect disappearance would deepen the understanding of this mechanism.

A closer look at the evolution of 3D snapshots was undertaken previously [74]. Prior to elaborating on these investigations, it should be noted that the growth direction in Figure 5 is 〈001〉. Unique stacking sequence in this direction is a key idea of colloidal epitaxy [22]. A stacking fault is marked with a red line in Figure 6a–e. Contrary to the intuition that stacking disorder does not change the particle number density, a sink of the center of gravity is accompanied with the shrinkage of the stacking fault, as shown in Figure 6f. A magnification of the snapshot around the lower end of the stacking fault is shown in Figure 7. It should be remembered that a stacking fault is terminated by a Shockley partial dislocation. In addition, it should also be noted that the Shockley partial dislocation terminating an intrinsic stacking fault—such as that shown in Figure 7—accompanies a particle deficiency of 1/3 lattice line. This accounts for the sink of the center of gravity accompanied with the shrinkage of the stacking fault. It was observed that the buoyancy of a partial dislocation acts as a driving force for the disappearance of defects in gravity. In addition to the buoyancy, a stress field yields a driving force that causes the Shockley partial dislocation to move upward along a glide plane [75]. Further, the present author and his coworker calculated a cross term between the elastic field due to a Shockley partial dislocation and that caused by gravity [76]. This effect also acts as a driving force for the upward motion of the Shockley partial dislocation. Although numerical estimation has been carried out using an elastic constant for the HS crystal, this effect is not limited to colloidal crystals. However, the condition of coherent growth (which was reported in Reference [77]) was used.

The glide mechanism of a partial dislocation for the reduction of stacking disorder was revealed in Reference [74]. As mentioned in Section 1.2, the fcc–bcc interfacial free energy was determined by Pronk and Frenkel. Using this value for stacking fault energy (free energy), they estimated the time for stacking fault to be annealed out of a sediment [33]. It was assumed that the melting and recrystallization mechanism was involved here. The estimated time was a few months long, which was fairly consistent with Dolbnya et al.’s result of 1–1.5 years to obtain a perfect fcc stacking. To evaluate the time required for the glide mechanism, the activation barrier for the partial dislocation to move to the next lattice position is needed. While we do not have such a value at present, this mechanism is apparently smooth, because only particles along a dislocation line are involved. Therefore, this mechanism would account for defect reduction in a short period or that occurring immediately after the crystal growth. If a dislocation glide is mediated by a kink motion along the dislocation, the motion becomes smoother. It is possible that defect behaviors with intermediate time scales are related to other mechanisms, such as the climb-up of dislocation and jog-related phenomena.

Immediately after our study [73], a grand canonical MC simulation study of HSs in gravity was performed by Marechal and Dijkstra [78]. They focused on the dynamical aspect of crystallization. They also confirmed the occurrence of crystallization in two successive bottom layers reported by Biben [70]. This phenomenon was later confirmed experimentally [79]. Marechal et al. performed a Brownian dynamics simulation study of the stacking disorder from a dynamical point of view [80]. It was found that fast sedimentation introduces hcp stacking sequence in fcc stacking. Accordingly, a trade-off is suggested: while the gravity effect reduces the stacking disorder, it also accelerates the sedimentation simultaneously. This dilemma has been addressed in detail in Section 4.2. In an earlier study [73], we proposed the use of a centrifuge rotor to control the strength of a gravitational field. Besides examining the stepwise control of a gravitational field in simulations, different types of control are also possible through the variation of the centrifuge rotation velocity.

The results reported in References [73,74] also faced criticism. The driving force of the fcc 〈001〉 growth in these simulations was stress arising from the simulation box of PBC, which indeed was artificial. Our rebuttal was that stresses forcing the fcc 〈001〉 growth could be produced without the PBC. In reality, the fcc 〈001〉 growth has already been realized by utilizing a template in the colloidal epitaxy [22]. This motivated computer simulations on the colloidal epitaxy, as presented in the next subsection.

4.2. Colloidal Epitaxy

In 1997, van Blaaderen et al. developed colloidal epitaxy [22]. As mentioned in Section 1.2, a key idea was the uniqueness of the stacking sequence in fcc 〈001〉. Interfacial free energy of the HS crystal against a flat hard wall is lowest for the fcc {111} face among low-indexed faces such as {110} and {100} [81,82,83]. Therefore, fcc {111} plane faces if one grows an HS crystal on a flat hard wall. This phenomenon has already been reported in an MD simulation study [84], and was also confirmed by a recent MC simulation study [85]. In reality, dense colloidal suspensions exhibit the same trend. Section 1.2 has already explained that it is difficult to avoid the stacking disorder in fcc {111} stacking by using the free energy difference.

In the preceding subsection, we discussed that fcc (001) stacking has many advantages, in addition to its uniqueness. Hence, it is possible to further develop this method. A shortcoming of colloidal epitaxy is that the obtained colloidal crystals are thin. Lin et al. successfully grew thick colloidal crystals with about 30 layers on a square grid template [86]. An advantage of the square grid is that it allows the lattice constant of the bottom layer to be adjusted via on-line matching of the template. Colloidal epitaxy on a square grid template has been reported (e.g., [87,88,89]). The present author employed a square grid as a template. In addition to the advantage just mentioned above, the computational cost to calculate interactions between particles and the pattern on the template is lower for a simple square grid than for an fcc (001) array. Crystallization of HSs on a patterned wall has already been investigated by MC simulations [90,91,92], and has also been recently investigated by MD simulations [93]. Enhancement of crystallization by various patterns has also been shown in the absence of a gravitational field. While gravity should be excluded to investigate the pure substrate effect, interplay between these two effects may yield some benefits. In addition, it is not easy to make the effect of gravity negligible in reality.

As a first stage, we merely replaced the flat bottom wall with a square patterned wall [94,95]. That is, the strength of gravity was increased in a stepwise manner. In Reference [95], the size of a horizontal system was doubled to that of Reference [94]. In simulations on a flat bottom wall (Reference [73]), the stacking sequence was changed to {111} stacking by doubling the size of the horizontal system. In contrast, in the simulation carried out on a square pattern, the stacking direction remained as 〈001〉. The pitch of the square grid was determined so that the lattice constant of the bottom crystal layer became slightly smaller than that for the simulation reported in Reference [73] at $g^{*}$ = 0.9. In References [94,95], our interest was not limited to the detection of the threshold $g^{*}$ value (∼0.9), but was also extended to the phenomena occurring in a relatively wide range around the threshold. Accordingly, this review also examined the defect behavior under conditions when $g^{*}$ is slightly greater than the threshold value. The side-to-side distance between adjacent grooves was $0.338\sigma$. The grid was made of grooves with width $0.707\sigma$. Accordingly, the periodicity of the square grid was $1.478\sigma$. The groove width was selected so that the diagonal length of a square in the cross-section of transverse and longitudinal grooves “coincides” to $\sigma$ (for detail, see note 18 of Reference [96]). Snapshots of the N = 26624 system ($L_x$ = $L_y$ = $25.09\sigma$ and $L_z$ = $200\sigma$) are shown in Figure 8. A two-fold—not three-fold—stacking sequence is an indication of fcc (001) stacking. Thickening of the crystal with increasing $g^{*}$ is also seen. A triangular-shaped defect, as shown in Figure 8, is remarkable. This extended defect was later identified with a stacking fault tetrahedron [97].

After successfully replacing the stress that forced the growth direction to fcc 〈001〉, we examined the sudden switch-on of gravity [101] (some details are given in Section 4.3 because they relate more closely to this section). In the case of a flat bottom wall, this gravitation resulted in polycrystallization [72]. In this respect, it was surprising that in the case of a square grid template, we could successfully avoid the polycrytallization/amorphization without controlling the $g^{*}$ value. This is another advantage of colloidal epitaxy, apart from the uniqueness of the stacking sequence. In other words, crystalline nucleation on a template and successive growth upward can overcome the crystallization inside. Of course, we cannot completely rule out the possibility that the crystal grown from a nucleation inside already becomes large so that it dominates the growth front of the crystal grown from the bottom when the growth front reaches the inside. Hence, it can be concluded that polycrystallization/amorphization can be avoided in colloidal epitaxy under a gravitational condition that the system polycrystallizes if the bottom wall is flat. Such a phenomenon has already been reported experimentally [98].

Here, it is important to note the contrasting wall effect. While the present and preceding subsections focus on the effects that promote the crystal growth in a wider sense, studies on walls that inhibit the crystallization on them are also interesting. Homogeneous nucleation inside a system was investigated using such walls by a Brownian dynamics simulation [99]. Combination of walls that promote and inhibit the crystallization may be used to create complicated structures.

4.3. Gravitational Annealing/Tempering

Annealing is a method of material processing. It is well known that the properties of GaN materials were successfully improved by a thermal annealing [100]; in this example, Mn-doped GaN were treated in ${\mathrm{N}}_2$-ambient above 700 °C for 20 min to lower the resistivity. In contrast, in quenching, the temperature is decreased at a certain rate during the treatment of materials. After an annealing treatment, one can, for example, treat the materials at a temperature higher than the annealing temperature. Through such a tempering treatment, material properties can be modified. The aforementioned simulation under a constant gravitational condition corresponds to annealing. This subsection first briefly describes the gravitational annealing simulations, and then introduces the simulation results of gravitational tempering.

Let us look at some snapshots of simulations under constant gravitational conditions (Reference [101]). Snapshots of the N = 26,624 systems are shown. In this simulation, both the system size and the template were the same as those reported in Reference [95], except for $L_z$ = $1000\sigma$. The results of gravitational annealing are shown in Figure 9 and Figure 10. Figure 9b shows that a crystal first formed at the bottom by sedimentation (Figure 9a) and then grew. A comparison of the $xz$ projections of Figure 9b,c shows an improved crystallinity owing to gravitational annealing. In particular, region $z/\sigma \lesssim 11$ clearly showed better and improved crystallinity than that in the above region. However, one cannot have a simple conclusion on the defect reduction from $yz$ projection. Region $11\lesssim z/\sigma \lesssim 17$ shows increases in the number of defects. On the other hand, the defect around ($x/\sigma$,$y/\sigma$,$z/\sigma$) ∼ (5,7,12) extends in a three-dimensional region. This means that at least 3D analysis is necessary. One cannot judge on the crystallinity, merely based on the extent of defect structures. In addition, a tiny defect structure around ($x/\sigma$,$y/\sigma$,$z/\sigma$) ∼ (−1,0,2) should not be ignored. Hence, it can be said that improvement of crystallinity by gravitational annealing at $g^{*}$ = 1.6 was successful but partial.

Let us look at the snapshots of gravitational annealing at $g^{*}$ = 1.4 (Figure 10). Seemingly, Figure 10b shows that crystallinity of the crystalline layers formed on the bottom (Figure 10a) was already better than that for $g^{*}$ = 1.6. Improved crystallinity can be more clearly seen in Figure 10c than in cases when $g^{*}$ = 1.6. The extent of a less-defective crystalline region is evidently larger for the $g^{*}$ = 1.4 case than for the $g^{*}$ = 1.6 case. In addition, although a tiny defect structure remained for $g^{*}$ = 1.6, no such defects were observed for $g^{*}$ = 1.4. In contrast, the dots around $z/\sigma$ ∼ 7.5 resulting from the overlaid projected points widened. Based on these observations, a deviation of the particle position from its regular lattice point is suggested. Alleviation of gravitational tightness can amplify the lattice vibration of individual particles or the collective motion of particles. The latter may be indistinguishable from lattice defects such as the wandering of lattice lines. However, the overall crystallinity is higher for $g^{*}$ = 1.4 than that for $g^{*}$ = 1.6. Improvement in crystallinity through gravitational annealing has been shown in simulations, although the optimization of annealing under the gravity condition has not yet been completed.

Inspired by these results, we proposed gravitational tempering to reduce defects in colloidal crystals [96]. After growing a colloidal crystal by sedimentation under a rather strong gravitational condition, treatment of the crystal at a moderate gravitational condition may improve the crystallinity. Alleviation of gravitational tightness is key for such a treatment. During crystal growth, defects that hardly move under a gravitational condition become mobile if the gravitation is weakened. However, weakening of the gravitational condition decreases the thermodynamic driving force of crystallization. In addition to the optimization of the tempering gravity, duration of the tempering program should also be optimized. We have successfully demonstrated MC simulations of gravitational tempering to reduce defects [102]. Sedimentation and gravitational tempering were tested at $g^{*}$ = 1.6 and 1.4: $g^{*}$ was kept at 1.6 for the first $2\times {10}^7$ MCCs, then decreased to 1.4 and maintained for $2\times {10}^7$ MCCs, and again 1.6 for the last $2\times {10}^7$ MCCs. The system size and the template were entirely the same as those reported in References [96,101]. Evolution of several parameters and snapshots of a simulation are shown in Figure 11. The result of a single simulation of five runs has been demonstrated here, and a dramatic defect disappearance can be seen (see the animation provided in the Supporting Information at http://pubs.acs.org/doi/suppl/10.1021/cg401884k). Evolution of the center of gravity $z_G$ is plotted in Figure 11a in black. After a steep sink of $z_G$ corresponding to sediment formation, it decreased gradually and became almost constant. Figure 11b is a snapshot of the crystal observed at the end of gravitational annealing after crystal growth by sedimentation at $g^{*}$ = 1.6. Defect structure did not disappear by this gravitational annealing. Just after switching $g^{*}$ to 1.4, $z_G$ promptly increased. After some time, $z_G$ suddenly sank slightly and then became almost constant. Snapshots just before and after this sink are shown in Figure 11c,d, respectively. An extended defect around ($x/\sigma$,$y/\sigma$,$z/\sigma$) ∼ (−10,3,15) disappeared. The slight sink of $z_G$ was caused by this defect disappearance. Hence, it can be concluded that an extended defect that could not disappear through a gravitational annealing treatment did disappear. Defective layers are seen at the top region of $yz$ projections in Figure 11c,d. Shrinking of this defective region was expected by treatment at $g^{*}$ = 1.6 for the last $2\times {10}^7$ MCCs. Unfortunately, this defective region remained unchanged (snapshot is not shown). The numbers of particles of fcc-like, hcp-like, and other crystalline types in the observed region are plotted together in Figure 11. The kink positions in the evolution curves of these numbers coincide with those in $z_G$, except for at the $4\times {10}^7$th MCC. The triangular-shaped extended defect (a stacking fault tetrahedron) disappeared through conversion from hcp or other types of stacking to fcc ones. In Section 4.1, it was mentioned that a Shockley partial dislocation that terminates an intrinsic stacking fault accompanies the particle deficiency of 1/3 lattice line. The particle deficiency at the end of a stacking fault tetrahedron arises from the extended region over a triangle. The particles that underwent conversion between other types of crystals are also distributed over this triangle. Therefore, changes in the numbers and $z_G$ were detectably large. In a previous study [102], we showed another result. Such changes were not detected for the shrinkage of only a few stacked planar defects.

In Section 4.1, we pointed out two aspects of the stacking disorder. It is generally observed that in crystal growth, crystals that grow rapidly contain many defects. A high degree of stacking disorder in rapid sedimentation is natural in this respect. An advantage of colloidal crystallization under a strong gravitational condition is enlargement of the grown crystal, rather than the rapid growth. It has been shown that even if the crystals contain a large number of defects, the defects in large crystals can be reduced by post-growth processing. Another advantage is the possibility of developing a technique of defect control. Hilhorst et al. successfully introduced a desired defect in a colloidal crystal in sedimentation by utilizing a template with a defective lattice structure [26]. For the functionalization of a colloidal crystal such as an optical waveguide, designed defects must be incorporated in the crystal. Such defect engineering can be performed by the sedimentation method using a template with regular lattice by varying the lattice spacing, as pointed out in Reference [97] and as explained as follows. It has been experimentally found that defects in colloidal crystals grown on a template are not triangular, but parallelepiped [103]. The defect structure of partial parallelepiped outer shape has also been reported in Reference [102]; results of two runs of five of the simulations have been reported, and in a run a triangular outer-shaped defect structure was seen as described above, and in the other, the outer shape of the defect structure was partial parallelepiped. As already mentioned, in order to match the template with the lattice spacing of the bottom layer at $g^{*}$ = 1.6, the pitch of the square grid was set slightly smaller than the lattice spacing at $g^{*}$ = 0.9 of Reference [73]. Under the stress provided by this pattern, stacking fault tetrahedra must dominate. By varying the lattice spacing of the template pattern, the dominant defect type may change. The effect of lattice spacing of the template pattern on the defect has already been argued by Jensen et al. [98]. As a recent development, the sedimentation method with colloidal epitaxy has been applied to control the structure of a grain boundary [104].

5. Discussion

A main focus of the present review is the crystal–fluid coexistence in equilibrium, and the other is the crystal growth. The subsequent subsection first introduces an aspect of crystallization included in the equilibrium simulations, and then a feasibility study is presented on gravitational tempering.

5.1. Hard-Sphere Crystallization without External Field

The results of Reference [42] and successive development were discussed in Section 3. However, some of the results were not published because of poor statistics. Results with good statistics were published by other researchers (such as in Reference [9]) before the present author could improve the statistics. Besides the small size of the system, interface drift is another problem. The initial velocities are generally assigned using random numbers with zero average. Even if the averaged velocity over the entire system vanishes, the average over a particular part (e.g., in an interfacial region) cannot be taken as zero. Particular attention must be paid to prevent the interface drift. The present author investigated the crystallization and melting induced by a mass flow in the melt [105]. The driving force of HS crystallization is densification. Therefore, densification induced by a mass flow must cause HS crystallization.

This phenomenon can be understood based on the principle of mass conservation. In a conventional MD simulation, momentum and energy as well as mass are conserved. Treating the materials as a Navier–Stokes–Fourier fluid, the present author derived relations between thermodynamic quantities of two phases coexisting in a non-equilibrium steady-state [106]. Relations for the velocity of the interface, net mass-flow velocity, and thermodynamic variables such as temperature, pressure, and density derived from conservation laws are symmetric with respect to crystallization and melting. In contrast, entropy production changes its sign when we invert the flow velocity; the crystallization and melting are swapped with each other by this operation. This gives rise to a criterion that can be used for determination of the stability of the non-equilibrium steady-state interface. This formulation was extended to the isothermal case [107]. However, related studies have not been performed.

5.2. Processing of Colloidal Crystals Using Centrifugation

The utilization of a centrifuge rotor proposed in Reference [73] has been described in Section 4.1. This section first provides an introductory description on spontaneous sedimentation under normal gravitational condition. Subsequently, we discuss colloidal crystallization by centrifugation sedimentation. In other words, the feasibility of the experimental realization of gravitational tempering is discussed.

Typical values of $g^{*}$ in normal gravity are described prior to determining them using a centrifuge rotor. Van Blaaderen et al. conducted a pioneering work on colloidal epitaxy and reported $g^{*}$ ∼ 1 [22]. Ackerson et al. carried out a light scattering experiment under the gravitational conditions $g^{*}$ ∼ 0.3 and 0.1 [108]. As mentioned in Section 4.1, Kegel and Dohnt conducted a milligravity experiment [24]. Gravity conditions as small as $g^{*}$ on the order of ${10}^{-4}$ have been reported by matching the density of dispersion media to that of a dispersion phase. In a similar experimental study by Martelozzo et al., the gravitational conditions $g^{*}$ ∼ 0.02 and 0.004 were reported [69]. A real-space observation of the sedimentary growth of colloidal crystals was performed by Hoogenboom et al. under the gravitational conditions $g^{*}$ ∼ 1 and 4 [109]. In Dolbnya et al.’s experiment on 1–1.5 year long gravitational annealing, the gravitational conditions were $g^{*}$ ∼ 0.03 [25]. In a real-space observation by Royall et al., the gravitational conditions $g^{*}$ ∼ 1 and 2 were reported [110]. In an experimental study that confirmed the occurrence of crystallization in two successive bottom layers in gravity by Ramsteiner et al., the gravitational condition was $g^{*}$ ∼ 7 [79]. These values indicate that a wide range of the $g^{*}$ value has already been covered. The main issue is in situ tuning of $g^{*}$. Varying $g^{*}$ by changing the density of a dispersion medium and/or mass/diameter of dispersing particles cannot be applied in gravitational tempering because it cannot be carried out in situ. On the contrary, in situ observation is possible by methods such as laser scanning confocal microscopy [22,79,109,110,111].

A number of studies have already been conducted using a centrifuge rotor. Megens et al. investigated colloidal crystals by small-angle X-ray scattering [112]. They centrifuged colloidal suspensions at ∼400 G for more than 5 h to densify them. Here, G stands for normal gravity. In addition to spontaneous sedimentation, centrifugation experiments were also carried out by Ackersonet al. [108]. They used $\sigma$ = 1000 nm and 760 nm particles, and reported that for smaller particles, the system failed to crystallize when the centrifugal force was greater than ∼6 G. Although this threshold value differs by two orders from Megens et al.’s centrifugation condition, it was not an error; in addition to the difference in other conditions, they focused on the kinetic aspect. Although the particle surface was not modified so as to exhibit the HS-like phase behavior, Suzuki et al. grew colloidal crystals by centrifugation at 82, 185, and 329 G [113]. Jensen et al. performed experiments on colloidal epitaxy at 9000 G [98]. These studies indicate that gravity condition can be controlled in centrifugation. In situ control of $g^{*}$ is possible for gravitational tempering programs through centrifugation. A key issue for realizing gravitational tempering has been shown in Suzuki et al.’s gravitational annealing experiment [114]. In this study, we succeeded in removing striations from a silica colloidal crystal grown at 9 G by gravitational annealing at 50 G for 5 days. Photographs before and after gravitational annealing were been presented therein. However, these results were based on ex situ observations. In situ monitoring of disappearing defects is quite difficult during centrifugation. Even more difficult is establishing a setup of an optical system for confocal microscopy in a centrifuge rotor. A plan to rotate a conventional microscope together with a cuvette is under investigation.

6. Summary

We have reviewed here simulation studies of hard-sphere (HS) systems in the context of crystal growth. After a general introduction of models of colloidal particles, simulation methods are discussed. The review consists of two main parts: one on the crystal–fluid coexistence state in equilibrium, and the other the HS systems in gravity. The former part initially describes how the present author constructed the crystal–fluid coexistence state of hard spheres for the first time. Next, some interfacial properties derived from coexistence simulations are introduced. The latter part gives reviews of HS systems in gravity. Firstly, the simulation results of crystallization of hard spheres on a flat wall in gravity are presented. Then, the simulations on a template of a square pattern are described. Finally, the processing of colloidal crystals using gravity is discussed. Gravitational tempering is a new technique proposed by the present author and his coworkers. The review also discusses non-equilibrium phenomena included in equilibrium simulations. It then investigates how mass flow in melt induces crystallization and melting at the interfaces. Feasibility of gravitational tempering has also been argued in the review. The review shows that a wide range of gravitational conditions can even be observed in normal gravity. Gravitational condition can be varied in situ through centrifugation. Difficulty in in situ observation of the defect behavior has also been pointed out.