Effect of container shape and walls on solidification of Brownian particles in a narrow system.

We carry out Brownian dynamics simulations and study the ordering of particles under a uniform external force in a narrow system. In our previous studies [M. Sato et al., Phys. Rev. E 87, 032403 (2013); J. Phys. Soc. Jpn. 82, 084804 (2013)], we showed that the ordering of particles depends on the direction of the external force. In the studies, however, the system size and the number of particles are small, so the behaviors we observed corresponds to the motions in the initial stage of crystallization. In this paper, using a longer container and more particles, we investigate how solidification in a narrow system proceeds. We also study the effect of the shape of simulation box on the ordering of particles. When we use a rhomboid as a simulation box, the ratio of the particles with the face-centered cubic structure to those with the hexagonal close-packed structure is larger than that in a cuboid system.


I. INTRODUCTION
Colloidal crystals with a three-dimensional regular structure have received much attention as materials of photonic crystals. Regarding the distance between particles, the colloidal crystals are classified into two types. One is a nonclose-packed crystal, in which the distance between particles is much longer than the diameter of particles, and the other is a close-packed crystal, in which the distance between particles is as large as the diameter of particles. The nonclose-packed colloidal crystals have an attractive feature as photonic crystals that the distance between particles is tunable [1,2]. They are created by use of various methods [3][4][5][6][7][8][9][10], e.g., shear-annealing, controlling the temperature gradient, and controlling the pH gradient.
On the other hand, close-packed colloidal crystals have a possibility to be used as templates for inverse opals with perfect three-dimensional photonic band gaps [11]. Creating large close-packed colloidal crystals is more difficult and requires more time than creating large non-closed packed colloidal crystals. One of the methods used to create closepacked colloidal crystals is the centrifugal method. Recently, Hashimoto et al. [12] succeeded in creating a large colloidal crystal in a container with an inverted triangular bottom using the centrifugal method. They suggested that the formation of the large grains is caused by suppression of nucleation. In a flat-bottomed container, nuclei are formed in many places. Many grains grow simultaneously with competing with each other, so the size of grains becomes small. On the other hand, in an inverted triangular-bottomed container, nuclei are mainly formed around the corner, so a few grains grows largely.
In our previous studies [13,14], keeping the experiment [12] in mind, we carried out simulations of the crystallization of Brownian particles induced by a uniform external force in a cuboidal simulation box. In an early stage of crystallization, the two-dimensional ordering of particles occurs and a triangular lattice is formed on the walls of a simulation box. Then, being affected by the two-dimensional structure on the walls, the three-dimensional ordering in bulk proceeds, and the face-centered cubic (fcc) and hexagonal close-packed (hcp) structures are formed in the solid phase. The two-dimensional structure formed on the walls is the {111} plane of the fcc structure or the {0001} plane of the hcp structure. The planes, however, cannot cross orthogonally. Thus, when we use a cuboidal simulation box, disordered particles are left in the region where the crystals growing from the bottom wall and those growing from a side wall collide with each other. When we change the shape of the container and control the angle between the neighboring walls, the crystal plane formed on the bottom wall and that on the side walls can cross without forming defects in bulk, and the ordering in bulk may become better.
In this paper, carrying out simulations with an external force, we study the ordering of particles in narrow simulation box and investigate the effect of the direction of an external force and the shape of a simulation box on the ordering of particles. In Sec. II, we introduce our model, which is a standard and simple model. In Sec. III, we show the results of our simulation. In our previous studies [13,14], we used a small system, and the formed crystal was thin. In this paper, we carry out simulation using a longer system. We investigate the dependence of the ratio of the number of the fcc-structured particles to that of the hcp-structured particles on the direction of an external force. We show that the number of particles with the fcc structure is higher than those with the hcp structure with a tilted force. We also investigate the effect of the shape of a simulation box on the ratio of the two structures. In Sec. IV, we summarize our results and give brief discussions.

II. MODEL
In our simulation, we use a cuboid or a rhomboid as a simulation box. The typical form of simulation box is shown in Fig. 1. We consider hard flat walls in the x and z directions and use a periodic boundary condition in the y direction. The corners of the simulation box are named A, B, C, D, E, F, G, and H (In Fig. 1, the corner E is hidden by particles). We add an external force whose x component is positive, so particles move to the neighboring of the wall ABCD. Thus, we call the wall ABCD as the bottom wall. The two walls parallel to the xy plane, ABFE and DCGH, are called the side walls.
The normal direction of the bottom wall is given by (β,0,1)/ 1 + β 2 and the angle between the lower side wall and the bottom wall is expressed as tan −1 β. The line DC is orthogonal to the line DH. The lengths of the lines DC and DH are L y and L x , respectively. The distance between the two side walls is kept L z irrespective of the tilted angle, so the volume of container is given by L x × L y × L z . We assume that particles act as point masses against the walls, so perfectly elastic collisions occur between a wall and a particle.
The algorithm of our simulation is the same as those in our previous studies [13,14]. For simplicity, hydrodynamic interactions between particles [15][16][17][18][19][20][21][22][23] are neglected in our simulation. The motion of a particle with mass m is given by where r i is the position of the ith particle. The first term in the right-hand side F ext e ext is a uniform external force, where F ext and e ext represent the strength and the direction of the external force, respectively. In the second term, the force F ij represents the interaction from the j th particle. The third term represents the resistance by viscosity, where ζ is the frictional coefficient. The final term F B i , which represents the random force, is isotropic and satisfies F B where T is temperature and k B is the Boltzmann constant.
The interaction force F ij is expressed as −∇U (r ij ), where U (r ij ) is the interaction potential and r ij = |r ij | = |r i − r j | is the distance between the ith and j th particles. We assume that the interaction between two particles is short-range repulsion and use the Weeks-Chandler-Anderson (WCA) potential [24] as U (r ij ). Namely, U (r ij ) is expressed as where σ represents the diameter of particles and r in = 2 1/6 σ .
We assume that the viscosity is high, so the acceleration rate is small and the left hand-side term in Eq. (1) is neglected. The velocity of the ith particle is approximately given by In our simulation, we use σ , ζ σ 2 / , and /σ as the units of length, time, and force, respectively. The normalized difference equation of Eq. (3) is simply expressed as [25] which is the scaled displacement of the ith particle by the random force, satisfies

III. RESULTS OF SIMULATION
In our simulation, the system size L x × L y × L z is 8L × L × L with L = 13.38 and the number of particles N is 10976, so the volume fraction of particles 4πN(σ/2) 3 /3L x L y L z is 0.3. Initially, we put particles in the system at random and move them without an external force for a long time. Then we set the time t to 0 and add an external force. The strength of the external force, the amplitude of random displacement, and the time interval areF ext = 1/ √ 2,R B = 0.1, and t = 1.0 × 10 −4 , respectively. First, we carry out simulations in a cuboidal simulation box. In addition to the simulations with the force normal to the bottom wall, we carry out simulations with a tilted force and investigate the effect of the direction of force on the ordering of particles. Then we use a rhomboidal simulation box and study the effect of the shape of a container on the ordering of particles.

A. Solidification in a cuboidal simulation box with a perpendicular force
We first show a typical motion of particles under an external force in Fig. 2. The simulation box is cuboid: the side walls are perpendicular to the bottom wall. The direction of the external force is (1,0,0). Since the x component of the external force is positive, the density of particles near the bottom wall increases with increasing time in an early stage [ Fig. 2(a)]. In the middle stage [ Fig. 2(b)], the region with a high density spreads. In a late stage [ Fig. 2(c)], almost all the particles gather to the side of the bottom wall.
We cannot see the ordering of particles clearly in Fig. 2, so we introduce order parameters Q l (i) and w l (i) [26][27][28][29][30] in order to estimate the local ordering around the ith particle. The definition of Q l (i) is given by where The angles θ ij and φ ij represent the polar and azimuthal angles of r ij , respectively. Y m l (θ ij ,φ ij ) is the spherical harmonics and n n is the number of neighboring particles. Since we use the WCA potential in our simulation, the distance between particles in bulk is about r in if a close-packed structure is formed. However, the distance between particles in solid may be fluctuated by the random displacement. Thus, we regard the j th particle as one of the neighbors of the ith particle when |r ij | is smaller than 1.1r in . The parameter w l (i) is given by where g l is defined as In Eq. (7), the integers m 1 , m 2 , and m 3 run from −l to l while satisfying the condition m 1 + m 2 + m 3 = 0, and the term in parentheses in Eq. (7) is the Wigner 3-j symbol [31].
Since the particles are pressed by an external force, the density of particles in bulk is usually as high as possible. Thus, the expected structures in bulk are the fcc and hcp structures. Using Q 4 and w 4 , we can distinguish the particles with the local fcc structure from those with the local hcp structure. Figure 3 shows the dependence of the number of particles on Q 4 and w 4 , in which we discretize Q 4 and w 4 with the width 2 × 10 −3 and count the number of particles. In Figs. 3(a) and 3(b), there appear two large peaks. In both figures, the positions of the peaks are well separated and the overlap of the two peaks is small. Figure 4 shows the distribution of the density of particles in the Q 4 -w 4 plane, where only the particles having 12 neighbors are taken into account and the density of particles is normalized by the maximum value. The region with high density is colored blue (dark). On the parameter plane, we find that there are two large spots and two small spots. The large spot in the upper left region corresponds to the group of the hcp-structured particles, and the other spots in the lower right region are the groups of the fcc-structured particles. From Fig. 4, we regard a particle as a hcp-structured particle if 0.02 < w 4 < 0.15 and We also introduce another order parameter d l (i,j ) [32,33] in order to judge whether the connection between the ith and j th particles is solidlike or liquidlike. When the ith and j th particles are neighboring particles, i.e., the distance between them is shorter than 1.1r in , the definition of d l (i,j ) is given by where q * l,m (j ) is the complex conjugate of q l,m (j ). If d 6 (i,j ) > d c , the connection between the ith and j th particles is a solidlike one. When the number of the solidlike connection n d is larger than n c , we regard the particle as a solidlike particle. In this paper, the critical value d c = 0.7 and the threshold of the number of solidlike connection n c = 4. Figure 5 show snapshots of particles during solidification, in which the particles are distinguished by the parameter d 6 (i,j ). When n n 4, we regard the particle as a dilute liquidlike particle and color it light gray (lightest). When n n > 4 but n d n c , the particles is a dense liquidlike particle and colored yellow (the second lightest). For red (dark) particles, which represent solidlike particles, n d > n c . In an early stage [ Fig. 5(a)], both dense liquidlike particles and solidlike particles are formed near the bottom. With increasing the time [ Fig. 5(b)], the number of solidlike particles increases. In Figs. 5(a) and 5(b), we find that the interface between the dilute liquid and the dense liquid is perpendicular to the external force, but the interface between the dense liquid and solid is not perpendicular to the external force: Solidification proceeds fast near the side walls z = L z and z = 0. In a late stage [ Fig. 5(c)], almost all the particles are solidlike particles. A few liquidlike particles are left only near the interface.
In order to see the structure of solidified particles in detail, we investigate the spatial distribution of the hcp and fcc structures in bulk. Figure 6 shows the time evolutions of the formation of the hcp and fcc structures, in which the system is seen from the y direction. The criterion to determine the states   ) and blue (the third lightest) particles represent the locally hcp-structured particles and locally fcc-structured particles with 12 neighbors. Red (dark) particles are solidlike particles, but the number of their neighbors is less than 12. We use the same data as those in Fig. 2. TABLE I. Criterion to determine the state of a particle, which is used in Fig. 6. n N and n d are the numbers of neighbors and solidlike connections, respectively. As the threshold of the number of the solidlike connection n c , we use n c = 4. of particles is given in Table I In order to concern that the structures in bulk in a late stage is strongly affected by the side wall, we compare the twodimensional structure in the first layer, which is formed by particles attaching on the wall z = L z , with the three-dimensional ordering in the second layer, which is formed by the particles attaching on the first layer. Since the particles are pressed to the walls by the external force, the triangular lattice is expected to be formed on the side wall z = L z in order to increase the density of particles. Here, we introduce the order parameter ψ(k) [34][35][36], which represents the local two-dimensional ordering on a wall. The parameter ψ(k) is defined as

State of a particle
which represents the sixfold symmetry in the wall around the kth particle. In Eq. (10), n k is the number of kth particle's neighbors on the wall, and the summation is carried out for the particles on the wall. θ km is the angle between the vector (x m − x k ,y m − y k ) and the positive x direction. Figure 7(a) shows the two-dimensional ordering on the wall z = L z , in which we regard the particles with ψ k > 0.9 as FIG. 7. (Color online) Relation between (a) the two-dimensional ordering of particles in the first layer and (b) the three-dimensional ordering in the second layer. We use the same data as those in Fig. 6(c). ordered particles. Note that the meaning of colors in Fig. 7(a) differs from those in Figs. 5, 6, and 7(b). If the ordered particles form a line parallel to the x axis, we color them red (dark gray). With increasing the angle between the line of particles and the edge CD, we gradually change the color from red (dark gray) to yellow (gray). White (light) particles are the particles which do not form a triangular lattice, and green (dark) particles are particles which do not attach to the wall z = L z . In Fig. 7(a), the yellow particles are not seen. Figure 7(b) shows the three-dimensional ordering in the second layer, in which the colors in the second layer are the same as those in Fig. 6 and particles whose y coordinate is smaller than the second layer is green. In the second layer, there are some disordered solidlike particles [red (the darkest) particles] near the bottom wall and the growth interface. Except for the particles, almost all the particles in the layer form a triangular lattice and their structure is the local hcp structure [orange (the second darkest) particles]. Thus, the sheetlike structure, which is formed under a strong influence of the wall z = L z , is the {0001} plane of the hcp structure.
In our simulation, the hcp-and fcc-structured particles are mixed in bulk. We investigate which structure is preferred in bulk. Figure 8 shows the time evolution of numbers of the hcpand fcc-structured particles in bulk. The number of particles with the local hcp structure N hcp and that with the local fcc structure N fcc increase monotonically with increasing time. At each time, N hcp is larger than N fcc . Both of them saturate in a late stage. The saturation of N fcc seems to occur slightly faster than N hcp .   Fig. 6. The external force is (1,0,1)/2 and the other parameters are the same as those in Fig. 2.

B. Solidification in a cuboidal simulation box with a tilted force
As shown in Sec. III A, N hcp is much larger than N fcc when the force perpendicular to the bottom wall. In this section, we change the direction of force and show that the structure in bulk and the ratio N hcp /N fcc change drastically. Figure 9 shows the snapshots of particle's positions with the force (1,0,1)/2. We find that solidification occurs from the walls. In an early stage [ Fig. 9(a) Fig. 9(a) and 9(b)], the interface between the solidlike particles and the dense liquidlike particles, i.e., the interface between the red and yellow (dark and the second lightest) particles, is not perpendicular to the direction of external force. However, the interface between the ordered solid and the disordered solid seems to be approximately perpendicular to the force. In a late stage [ Fig. 9(c)], almost all the particles are solidlike particles. We can see a large domain of the fcc-structured particles. There are a few hcp-structured particles [orange (the fourth lightest) particles], which form lines in Fig. 9(c). At the upper right side, a beltlike region of disordered solid particles are formed toward the lower right side. In the region of the right-hand side of the beltlike region formed by the disordered solidlike particles, the orientation of the domain is changed.
We investigate the structure in the late stage in detail. Figure 10(a) shows the two-dimensional order in the side wall CDHG. The meaning of the colors are same as in Fig. 7(a). The color of almost all the particles is red (dark) in Fig. 7(a). However, the particles which do not have the symmetry of sixfold orientation [white (the lightest) particles] appear at the place marked by arrows in Fig. 10(a). In the right-hand side region from the disordered region, one of the lines of the triangular lattice is tilted from the line CD about 30 • and the color of particles is yellow. Figure 10(b) shows the three-dimensional ordering, which is viewed from the z direction. Note that the meaning of the colors differs from that in Fig. 10(a). Four dense liquidlike particles [yellow (the second lightest) particles] are present at the region marked by arrows. Around the particles, the disordered solidlike particles are present in the lower layers. Figure 10(c) shows particles on the plane BFGC, in which we painted the region with y < l y − σ/2 gray in order to focus on the particles crossing the plane BFGC. In the left-side region from the position marked by an arrow, we can see the arrays of particles, which is roughly parallel to the line RP. The angle θ , which is the angle between the line RP parallel to the array of particles and the line CG, is about 57 • . Taking into account the structure on the wall CGHD and the plane BFGC, we find that the plane that is slightly tilted from one of the {101} planes of the fcc structure appears on the plane BFGC. Thus, if we cut the crystal on the line PR, we see one of the {100} planes of the fcc structure on the cross section. In Fig. 10(c), at the position marked by an arrow, the narrow region of the disordered solidlike particles spreads toward the lower right direction. In the right-hand side of the region, in which the orientation of the triangular lattice on the wall CGHD changes, the particles form a line parallel to the line CG. The difference of the orientation of ordered particles on the wall affects the orientation of the structure in bulk.
Next, we investigate how the number of ordered particles increases with time. Figure 11 shows the time evolutions of N hcp , N fcc , and (N hcp + N fcc ) with the force (1,0,1)/ √ 2. Comparing Fig. 8 with Fig. 11, the number of ordered particles (N hcp + N fcc ) hardly changes. The time evolutions of N hcp and N fcc are also similar to those in Fig. 8. Namely, the numbers increase with increasing time and saturate in a late stage. However, the relation between N fcc and N hcp is opposite to that with the force perpendicular to the bottom wall: N hcp is larger than N fcc in Fig. 8, but N fcc is larger than N hcp in Fig. 11. The change in the relation between N hcp and N fcc by the direction of the external force may be affected by the change in the spatial distribution of ordered particles, so we investigate how the spatial distribution of ordered particles changed by the direction of the external force. Figures 12 and 13 show the spatial distribution of numbers of the hcp-and fcc-structured particles in bulk, in which the numbers are counted in the y direction and averaged over 10 runs. The left ends in the figures correspond to the bottom walls. The number of particles is large in the region with dark color. When the external force is perpendicular to the bottom wall, many hcp-structured particles appear around three walls, x = L x , z = L z , and z = 0, as shown in Fig. 12(a). The width of the region with the hcp structure is wide. The fcc-structured particles mainly appear in the right side in Fig. 12, which is far from the bottom wall. When the external force is tilted from the bottom wall (Fig. 13), the hcp-structured particles appear near the walls, but the width of the structure is thinner than that in Fig. 12.  On the other hand, the width of region with the fcc-structured particles increases.

C. Solidification in a rhomboidal simulation box with a tilted force
In Secs. III A and III B, we use a cuboidal simulation box. Hereafter, we use a rhomboidal simulation box and investigate the relation between the ordering in bulk and the angle between the side wall and the bottom wall. In Fig. 9, where the two walls are perpendicular to each other, the triangular lattice formed on the walls is mainly the {111} crystal plane of the fcc structure. The two {111} planes cannot cross perpendicularly, so the formation of defects in solid cannot be avoid with a cuboidal simulation box. For example, disordered solid and the dense liquidlike particles are left near the bottom wall in Fig. 9(c). However, if we use a rhomboidal simulation box, in which the angle between the bottom wall and the lower side wall is tan −1 (2 √ 2), it may be possible to form the fcc structure whose the fcc structure whose {111} planes appear on the walls and to decrease the number of defects in bulk. Thus, we carry out simulations with a rhomboidal container and concern the effect of the shape of container on the ordering of particles in bulk. Figure 14 show the snapshots of the ordering of particles in a tilted container. The force is given by (1,0, √ 2)/ √ 6, which bisects the angle ∠BCG. Initially, a domain of the hcp-structured particles is formed near the upper corner. In a middle stage [ Fig. 14(b)], however, the domain vanishes and the number of the particles with the local fcc structure increases. A few particles with the local hcp structure are left. In a late stage [ Fig. 14(c)], liquidlike particles are left near the lower side of the bottom wall. Almost the whole region is filled by the fcc-structured particles, but they are separated into some domains by disordered solidlike particles. Thus, contrary to our expectation, the ordering in bulk does not seem to be higher than that in Fig. 9.
In order to show the effect of tilted angle in detail, we need to carry out many runs and see averaged values. Figure 15  show the time evolutions of numbers of the hcp-and fccstructured particles in the tilted container. When we use the tilted simulation box, N fcc is as large as that in the cuboidal simulation box. On the other hand, N hcp in the tilted simulation box is smaller than that in a cuboidal simulation box. The ratio of the number of hcp-structured particles N hcp to the number of ordered particles N hcp + N fcc is about 0.20, which is smaller than that in Fig. 11. Figure 16 shows the spatial distributions of the hcp and fcc structures in bulk. Since the bottom wall is tilted, the upper left region is out of the container. Around the lower left corner, which corresponds to the neighboring of the edge AB, the ordered particles hardly appear. The regions with many hcpstructured particles appear near the bottom wall and the upper side wall. However, comparing Fig. 16(a) with Figs. 12 and 13, we find that the region with many hcp-structured particles becomes narrow and the number of the hcp-structured particles in the region is small. On the other hand, the region with many fcc-structured spreads in the whole of system and the number of fcc-structured particles is larger than that in a rhomboidal simulation box.  We count the number of particles in the y direction. The angle between the bottom wall and the lower side wall is tan −1 2 √ 2 and the direction of the external force is (1,0, √ 2)/ √ 6. The direction of the external force bisects the angle at the upper corner C. The number of ordered particles is large in the dark region. The parameters are the same as those in Fig. 9, and the data are averaged over 10 runs.

IV. SUMMARY AND DISCUSSION
In this paper, keeping the formation of colloidal crystal caused by a centrifugal force in mind, we carried out Brownian dynamics simulations in narrow systems and studied the ordering of Brownian particles by an external force. Especially, we investigated the effect of the walls, the external force, and the shape of container on the ordering in bulk. When the external force is perpendicular to a wall, the layers of ordered particles parallel to the bottom wall are formed in an early stage. Then, crystallization also starts along the side walls. The growth from the side walls is faster than that from the bottom wall, so that in a late stage, the layers of ordered particles parallel to the side walls become dominant in solid.
During crystallization, the two-dimensional ordering on walls initially occur. Then the ordered particles on walls act as substrate and three-dimensional ordering in bulk proceeds. A triangular lattice formed on walls is one of the {111} planes of the fcc structure or one of the {0001} planes of the hcp structure. Thus, when we use a cuboidal simulation box, the formation of defects in bulk cannot be avoided. We expected that defects in bulk may decrease with a rhomboidal simulation box. In our simulation with rhomboidal simulation box, the ratio of fcc-structured particles to the hcp-structured particles increases, but the decrease of the disordered particles in a rhomboidal simulation box was not clearly observed. The formation of defects probably depend on the strength and the direction of force.
When the direction of force is perpendicular to the bottom wall, the number of the hcp-structured particles is larger than that of the fcc-structured particles. On the other hand, the ratio of the hcp-structured particles to the fcc-structured particles is opposite with a tilted force: The number of the fcc-structured particles is larger than that of the hcp structure. The change in the ratio of two structures due to the direction of the external force may be caused by the difference in the growing interface. In our simulation with the tilted force, as shown in Fig. 9(b), the growing interface of the ordered particles is vicinal face of the {100} plane of the fcc structure. The way of stacking of particles is unique and the frequency of formation of defects is low, so the fcc structure can grow largely. On the other hand, as shown in Fig. 6, the growing interface is not the {100} plane with the perpendicular force. Thus, the two ways of stacking is possible and both the hcp and fcc structures are formed. Another reason is probably the effect of walls. When the direction of force is perpendicular to the bottom wall the number of particles attaching to the wall is larger than those with the tilted force. In our simulation, the hcp-structured particles appear around the walls, so the number of the hcpstructured particles is large with the force perpendicular to the bottom wall. Now, we cannot explain why the hcp structure is formed near the wall and why the fcc structure appears from the solid-liquid interface. Thus, carrying out simulations with various directions and strengths of the external force, we intend to clarify the reason.
Comparing Fig. 6 with Fig. 9, the domain size with the tilted external force is larger than that with the perpendicular external force, which is qualitatively consistent with a experiment [12]. Recently, Suzuki et al. [37] used a thin container tilted from a centrifugal force and succeeded in forming a large colloidal crystal. In the experiment [37], a close packed plane is parallel to a upper side wall, which qualitatively agrees with our simulation. However, when we compare our results to the experiments [12,37] quantitatively, we need to note that the parameters we used are not comparable to the experimental ones. For example, the Peclet number in our simulation is given by F ext σ/k B T = 10/ √ 2, which is much larger than the experimental value [12,13]. Thus, the typical domain size in our simulation is much smaller than the experimental value. In our model, we neglect hydrodynamic effects and the acceleration term, which may affect the crystal structure in sedimentation. Due to the restriction of computational resource of our group, now we cannot get rid of the differences between the experiments [12,13] and our simulation, but we need to improve our model and perform simulation with more comparable value of parameters in future.