Phase Behavior of Colloidal Hard Tetragonal Parallelepipeds

Nov 11, 2005 - The impact of particle geometry on the phase behavior of hard colloidal tetragonal parallelepipeds (TPs) was studied by using Monte Car...
4 downloads 0 Views 428KB Size
Download PDF
23008

J. Phys. Chem. B 2005, 109, 23008-23015

Phase Behavior of Colloidal Hard Tetragonal Parallelepipeds (Cuboids): A Monte Carlo Simulation Study Bettina S. John and Fernando A. Escobedo* School of Chemical and Biomolecular Engineering, Cornell UniVersity, Ithaca, New York 14853-5201 ReceiVed: September 12, 2005; In Final Form: October 7, 2005

The impact of particle geometry on the phase behavior of hard colloidal tetragonal parallelepipeds (TPs) was studied by using Monte Carlo simulations in continuum space. TPs or “cuboids” of aspect ratios varying from 0.25 to 8 were simulated by approximating their shapes with multisite objects, i.e., via rigid clusters of hard spheres. Using equation of state curves, order parameters, radial distribution functions, particle distribution functions along three directions, and visual analysis of configurations, an approximate phase diagram for the TPs was mapped out as a function of aspect ratio (r) and volume fraction. For r > 3 and intermediate concentrations, the behavior of the TPs was similar to that of spherocylinders, exhibiting similar liquid crystalline mesophases (e.g., nematic and smectic phases). For r ) 1, a cubatic phase occurs with orientational order along the three axes but with little translational order. For 1 < r < 4, the TPs exhibit a cubatic-like mesophase with a high degree of order along three axes where the major axes of the particles are not all aligned in the same direction. For r < 1, the TPs exhibit a smectic-like phase where the particles have rotational freedom in each layer but form stacks with tetratic order. The equation of state for perfect hard cubes (r ) 1) was also simulated and found to be consistent with that of the rounded-edge r ) 1 TPs, except for its lack of discontinuity at the cubatic-solid transition.

I. Introduction The industrial applications of nanoparticles continue to grow as new discoveries are being made to facilitate their fabrication and to tailor their properties. The realm of application of these particles includes photography, catalysis, biological labeling, photonics, optoelectronics, information storage, surface enhanced Raman scattering, and formation of magnetic ferrofluids.1 Phase stability of a suspension of particles depends not only on their size and chemical composition but also on their shape,2 and this can be critical for the potential applications of the materials.3,4 Particle shape also plays an important role in the kind of liquid crystalline phases formed under lyotropic conditions. Particles with unconventional shapes have recently been synthesized which may exhibit uncharted liquid crystalline behavior.1,5-12 These uncharted phases could have interesting optical and rheological properties which may open up the way to multiple applications.13,14 New developments in synthesis methods have enabled the production of nanoparticles of different shapes. Gold and silver nanocubes have been prepared in solution by Sun and Xia.1 These nanocubes are observed to self-assemble in 2D. Silver nanoprisms have been made from spherical colloids using dual beam illumination by Jin et al.5 Nanorods have been synthesized by different groups using a number of methods.6,7 Tetrapodshaped particles have also been synthesized by Manna et al.7 which could be the precursors for inorganic dendrimers. Monodisperse hexagon Cu2S nanoplates synthesized by Zhang et al.8 were observed to self-assemble into superlattices on evaporation. Nanoacorns were synthesized from anisotropically phase segregated CoPd sulfides by Teranishi et al.9 The potential for fine-tuning nanoparticle properties has increased as a result of the advances in synthesis techniques. * Address correspondence to this author. E-mail: [email protected].

Computer simulations have been extensively used to probe entropy-driven liquid crystalline transitions where the participating particles experience purely repulsive interactions. In such systems, the isotropic fluid typically transforms into an ordered phase at sufficiently high concentrations because the loss of orientational entropy is compensated by the increase in packing entropy associated with the ordered phase. In the following, we briefly review some of the key findings regarding the lyotropic liquid-crystalline behavior of spherocylinders, ellipsoids, platelets, cutspheres, and cylinders. Bolhuis and Frenkel15 mapped out the phase diagram of hard spherocylinders as a function of aspect ratio. Besides isotropic and nematic phases, a smectic and two crystalline phases, AAA stacked and ABC stacked phases, were detected. Camp et al. and Camp and Allen16 studied the phase behavior of hard ellipsoids as a function of aspect ratio. Four homogeneous fluid phases were detected, namely, an isotropic phase, a nematic phase (where the long molecular axes of the molecules are aligned parallel to each other), a discotic phase (where the short molecular axes of the molecules are aligned parallel to each other), and a biaxial phase (with each of the molecular symmetry axes having its own preferred direction of alignment). Eppenga and Frenkel17 found a weak first-order isotropicnematic transition in a simulation study of infinitely thin hard platelets. More recently, Bates and Frenkel18 found a nematiccolumnar transition in a system of hard cylinders in the “thin disk” limit of L/D f 0 by molecular simulation. Veerman and Frenkel19 studied particles shaped as cut spheres and found that the phase behavior strongly depended on the aspect ratio L/D. For L/D ) 0.1, isotropic, nematic, columnar, and solid phases were observed. For L/D ) 0.2, a cubatic phase with cubic orientational order but no translational order was also observed. This cubatic phase, wherein the cut spheres stack up in short columns (the effective L/D of the column being approximately

10.1021/jp0551521 CCC: $30.25 © 2005 American Chemical Society Published on Web 11/11/2005

Phase Behavior of Cuboids 0.8) which are perpendicular to one another, has not been experimentally observed thus far. Simulations of cylinders of aspect ratio L/D ) 0.9 did not reveal a cubatic phase; instead, another ordered phase (in addition to a crystal phase) was observed wherein the orientations of the particles were either along a column or perpendicular to it.20 In a study of small systems of multisite soft-sphere cuboids using Dissipative Particle Dynamics by Elliott and Windle,21 a continuous isotropic-cubatic phase transition was observed. In an MC study of hard multisite cuboids by John et al.,22 the cubatic phase was predicted to exist between the solid and the isotropic phases and to exhibit a first-order isotropization transition. The cubatic phase has also been theoretically predicted in several studies. Blaak and Mulder23 predicted a cubatic phase in a system of cross-like particles (“Onsager crosses”) when the rods were of approximately equal lengths. They used bifurcation analysis and a Gaussian approximation for the nematic orientational distribution function. Frenkel24 predicted the existence of a cubatic phase in a system of particles comprised of three perpendicular rods. Very recently, MartinezRaton25 used fundamental measure functional (FMF) theory to predict the phase behavior of prolate and oblate hard parallelepipeds of restricted orientations (Zwanzig model). Several inhomogeneous phases were observed including the smectic, columnar, discotic smectic, and oriented or plastic solid; however, a cubatic phase was not predicted. The purpose of this study is to investigate the lyotropic phase behavior of tetragonal parallelepipeds (TPs) in continuum space. The term TP is used in this paper to include cubes as well as oblate and prolate parallelepipeds. The terms cuboid and TP are used interchangeably in this paper. The TP shape is used to investigate the role played by nearly flat sides, corners, and edges of the particles on the liquid crystalline phase behavior. A phase diagram is mapped out as a function of the aspect ratio of the TP. Besides the cubatic phase, several other novel mesophases were observed (depending on particle shape and concentration) having varying degrees of orientational and translational order. Further, the effect of “roughness” of the particle surfaces and edges was also studied by contrasting the behavior of multisite TPs and perfect cubes. To the best of our knowledge, Monte Carlo simulations of hard perfect cubes in continuum space have not been carried out before. The method used to find the overlap criterion for cubes can be extended to study hard particles of any geometry (see the Supporting Information). The rest of the paper is organized as follows. Section II explains how the system of TPs was modeled, the details of the Monte Carlo simulation method used, and the thermodynamic and structural properties that were studied. Section III presents and discusses the results of the simulations that were undertaken and compares the results with those of spherocylinders and theoretical predictions. Section IV summarizes the main findings of this work. II. Simulation Methodology The tetragonal parallelepiped (TP) particles (cuboids) are modeled as rigid clusters of tangent spheres fitted in a cuboidal frame. The interior of the TPs is left hollow for computational simplicity. The TPs are hard bodies with fixed distances and angles among their constituent spheres. This model allows easy evaluation of site-site overlap and also provides an approximation to the real particles which may not have the exact shape of a cube or a rectangular TP. In general, a TP as shown in Figure 1 is defined by the relative lengths of its three orthogonal

J. Phys. Chem. B, Vol. 109, No. 48, 2005 23009

Figure 1. Model TP with aspect ratio r ) c/a. The “multisite” TP constructed from spheres used in our simulations is shown on the right.

sides: a, b, and c. In this work, we restrict our attention to TPs with a square cross section, i.e., a ) b, and define the aspect ratio r as r ≡ c/a. The use of spheres to model the TPs becomes computationally intensive when the aspect ratio of the TPs becomes large. We have used Monte Carlo simulations to map out the equation of state for the particles of a given shape. We used a system with N ) 512 TPs in all simulations. Periodic boundary conditions were used to emulate a bulk system. The system was simulated at constant osmotic pressure P* and at a fixed temperature. In this case the temperature is irrelevant because the TPs were hard objects and thus the kinetic energy has no effect on configurational properties. To sample the configuration space the particles were repositioned and reoriented by using translation and rotation moves. The volume of the simulation box was equilibrated at a prescribed pressure by using volume moves.26 In a volume move, the relative centers of mass of the particles remained fixed. Isotropic volume moves (where all box sides change proportionally) were used at low pressures and anisotropic moves (where each box side is allowed to change independently) were used at higher pressures to allow the system to reach its equilibrium particle density. Furthermore, the Parrinello-Rahman method27 was adapted to perform volume moves that allow box shape deformation; this was used with selected oblate TPs and in cubes to ensure that the mesophases obtained with the more conventional volume moves were not being unduly stabilized by tetragonal box shapes. One MC cycle is comprised of one volume move and as many translation and rotation moves as the number of particles. All the moves were accepted with use of the conventional Metropolis criterion.28 The moves were all designed to have an acceptance rate of about 27%. A. Mapping Phase Space for TPs of Varying Aspect Ratio. The aspect ratio r of the TPs was varied from 0.25 to 8. The square cross-section (a ) b) of the TP was 4 × 4 spheres. The ratio c:a ) 1:4, i.e., r ) 0.25, resulted in a square plate. Square layers of spheres were added to the TP to reach the desired aspect ratio. For the limiting cases of r f 0 and ∞, the TPs phase behavior should resemble that of square disks and rods, respectively. The phase space of all the TPs was mapped out through expansion and compression runs in successive simulations where the pressure was gradually decreased or increased, respectively. The equation of state and phases obtained were compared with those of cut spheres (disks) and spherocylinders. The expansion runs from a highly dense system yielded a better, more reproducible description of the phase behavior for most systems. For systems with r e 6.5, we created a densely packed solid phase and expanded it to map out the phase behavior. For the r ) 8 system, we created a densely packed solid as the starting configuration for the expansion simulations and a perfectly aligned nematic phase as the starting configuration for the compression simulations. The equation of state was mapped as pressure vs volume fraction. The volume fraction is defined as the ratio of the total volume occupied by the TPs to the total volume of the simulation box. For simplicity a TP is assumed to occupy a volume ) a × b × c. Pressure is defined

23010 J. Phys. Chem. B, Vol. 109, No. 48, 2005 as Pσ3/, where σ is the diameter of a sphere in the TP, P is the simulation pressure, and  is an arbitrary energy parameter. In principle, a Gibbs-Duhem integration method could be used to map out coexistence lines as a function of r. We did not use this method as it would require (1) a starting coexistence point known with high accuracy and (2) either a molecular dynamics approach15 or a perturbation method16 to simulate the derivative of the system free energy with respect to r. The perturbation method for (2) would be ineffective if perturbations in r involve addition of one 4 × 4 layer of spheres (to all cuboids) because overlaps would be unavoidable at the high phase concentrations of interest. B. Examining the Effect of Surface Roughness. Some of the TPs synthesized experimentally are likely to have rough edges and so the impact of the surface roughness on the phase behavior is important. We probed the effect of surface roughness in a symmetric TP (a ) b ) c). Different degrees of particle roughness were obtained by varying the number of spheres per side of the cube between two (2 × 2 × 2 TPs) and four (4 × 4 × 4 TPs). The equations of state for these systems were compared with that obtained for perfect cubes. The earlier study by John et al. dealt with 3 × 3 × 3 TPs. The equation of state was mapped out in all cases by using both expansion and compression runs. The dimensionless osmotic pressure reported is P* ) Pa3/, where  is an arbitrary energy parameter and a is the side of a TP (a ) b ) c). A system of perfect hard cubes was also simulated at NPT conditions to examine the phase behavior in the limiting case of TPs having sharp edges and no surface roughness. The overlap criterion for two cubes was adapted from the method described by Gottschalk et al.29 (see the Supporting Information). The side of the cube was the same as that of the TPs (a ) b ) c). The same types of moves were used to relax the system and to change the volume as those used for the TPs. C. Characterization of Order and Particle Distribution Function. The order parameter used in this paper was mainly the nematic order parameter Q00. By definition, Q00 ) (1/2)〈3 cos 2θ - 1〉 [≡P2(cos θ)], where θ is the angle between the particle axis and the director. In this study, order parameters with respect to both the long and short axes of the TP were computed. If two or more axes were equivalent the order parameter was evaluated by finding the TP axis (from among the equivalent axes) that aligns best in the direction of interest. This method will be referred to as the “selected particle-axis” nematic order parameter. The director was computed by diagonalizing the ordering tensor comprising the particle axial vectors considered. However, in a few cases, the director was calculated as the average of all particle axes considered. An alternative method to select the director (applied to r ) 1 cuboids primarily) was also implemented as follows. The three axes of all particles were first divided into three sets such that axes in a given set have maximum mutual alignment. The director for each set was then found and one of them was randomly chosen as the director. The nematic and cubatic order parameters (see below) calculated by using the angle subtended by all axes of all cuboids with this director will be referred to as “selected director” order parameters. The cubatic order parameter was calculated for cuboids with r ) 1. It is defined as P4(cos θ) ≡ 〈(35 cos4θ - 30 cos2θ + 3)/8〉. For a system to possess cubatic order, P4(cos θ) should be nonzero and the nematic order parameter P2(cos θ) should be zero.30 The cubatic orientational correlation function was also calculated from gO(r) ) 〈(5∑i,j(u bi(0)‚u bj′(r))4 - 9)/6〉, where b ui

John and Escobedo

Figure 2. Phase diagram of TPs with aspect ratio varying from 0.25 to 8. The dashed lines indicate approximate phase boundaries. P1 is the smectic-like phase found in TPs with r < 1. C1 is the columnar phase for TPs with r < 1. The open squares mark the approximate boundaries of the solid-cubatic coexistence region in r ) 1 cuboids.

and b uj′ are the directors of the two cuboids i and j.30 The directors of the cuboids are the three perpendicular axes of each cuboid. The particle distribution function along the three box axes was obtained by computing the ratio between the average number of TPs whose centers lie in a slab of volume LyLzdx at a distance x from the center of a TP along the X axis and the equivalent number in an ideal gas with the same overall density (Ly is the box length along the Y axis, Lz is the box length along the Z axis, and dx is the slab width). It is an analogue of the radial distribution function in 1D. It can be written as F(x)/F, where F(x) is the number of TPs per unit length between distance x and x + dx from a particle along the X axis and F is the average density. Similar definitions apply for the F(y)/F and F(z)/F functions. III. Results A. Phase Diagram of TPs of Different Aspect Ratio. The phase diagram of the TPs with 0.25 e r e 8 is shown in Figure 2. The two-phase region is shaded and was found from the bounds of the breaks in the equation of state curves; such breaks are the hallmark of first-order transitions. However, several transitions were not first order (e.g., the columnar to smecticlike mesophase for TPs with r < 1) and these are not shown as shaded regions in the phase diagram. The dashed lines in Figure 2 indicate regions where the phase boundary is estimated, but could not be pinpointed from the simulation data (e.g., in cases where continuous transitions may occur). TPs with r ) 6.5 and 8 (Figure 3). On expanding the system from a highly dense and highly ordered solid phase where the particles are arranged in layers, it was observed that a columnar phase was formed in which the TPs had biaxial order (Figure 3a). In such a columnar phase, the particles have high orientational order in all three directions but they can “slide” with respect to one another indicating lack of translational order along the long axis. The smectic phase is obtained upon expansion from the columnar phase. The smectic phase is characterized by translational and orientational order only along the long axis (Figure 3b). The nematic phase is obtained from the expansion run from the smectic phase. The nematic phase has only orientational order along the long axis. Expansion of the nematic phase led to the formation of an isotropic phase. This nematicisotropic transition is marked by a break in the P vs volume

Phase Behavior of Cuboids

J. Phys. Chem. B, Vol. 109, No. 48, 2005 23011

Figure 5. Particle distribution function (top) and snapshots (bottom) of the r ) 2 TPs for the CL mesophase formed at P ) 0.15 (a, b) and the columnar phase formed at P ) 0.4 (c, d). The uniformity of the z-distribution function at P ) 0.4 reflects the lack of layering along the z direction and suggests mutual slipping of the columns.

Figure 3. Snapshots of selected phases for the r ) 8 TPs: (a) columnar phase and (b) smectic phase. The colors are used to help distinguish individual TPs.

Figure 4. Equation of state for the r ) 8 TPs.

fraction curve and hence is first order. The particle distribution functions for r ) 6.5 TPs (not shown) revealed that (i) the isotropic phase did not show any ordering along any of the three axes, (ii) the nematic phase showed weak layering along the x, y, and z axes, (iii) the smectic phase showed strong layering along the z axis only, and (iv) the columnar phase showed appreciable layering along the x and y axes but weak layering along the z axis reflecting the “slippage” of particles along the z axis. Compression runs from the nematic phase yielded the smectic phase without any hysteresis (Figure 4). Such a nematicsmectic transition, however, may only be weakly first order since there was no marked break in the P vs volume fraction curve as observed from Figure 4. Further compression of the smectic phase led to the columnar phase, whose region of stability was consistent with that observed before in the expansion runs (started from the crystal phase). Compression runs from the isotropic phase did not lead to a nematic phase but to a distinct

metastable phase that we named the “parquet” phase. In this phase, the particles are arranged in stacks that were oriented roughly perpendicular to one another. The parquet phase is explained in more detail later in connection with Figure 10b. TPs with r ) 4 and 5. These form a smectic phase upon expansion from the columnar phase. The transition was distinctly first order for r ) 4 and 5. Further expansion of the smectic phase leads to the formation of an isotropic phase via a firstorder transition. A stable nematic phase was not observed. The presence of liquid crystalline phases depends on the molecular shape as well as the aspect ratio. In linear tangent hard sphere chains, the nematic phase first appeared at an aspect ratio of 5.31 The compression runs for this system yielded the parquet phase from the isotropic phase. TPs with 1 < r < 4. These systems underwent a transition at P ) 0.4 and volume fraction 0.9 from the columnar phase to a cubatic-like (CL) mesophase characterized by alignment of the three axes of the particles but not necessarily the same major or minor axes for all particles (Figure 5). The transition was weakly first order as there was a small change in the gradient of the equation of state curve. This CL mesophase melts into the isotropic phase at P ) 0.07 and volume fraction 0.4 upon expansion. Unlike the transition from the columnar phase to the CL mesophase, the transition from the CL mesophase to the isotropic phase was first order as indicated by a break in the equation of state. The particle distribution function revealed layering in two directions in this CL mesophase and a mild preferential ordering along the z axis at intervals equal to the length of the major axis of the TP (Figure 5). Snapshots of the CL mesophase and the columnar phase are also shown in Figure 5. The particle distribution function in the columnar phase shows strong layering along two directions and more uniform distribution along the z direction (Figure 5); this indicates translational disorder in the z direction. The nematic order parameters calculated in the CL mesophase and columnar phase, based on both the long and short axes, indicate orientational correlation along their directions, the order parameter increasing monotonically in the mesophase whereas it remains almost a constant in the columnar phase. For the equivalent short axes, the above

23012 J. Phys. Chem. B, Vol. 109, No. 48, 2005

Figure 6. Particle distribution function (top) and snapshots of the solid phase (bottom) for the r ) 0.5 TPs at P ) 0.8. The snapshots give different views of this phase.

calculation was based on the selected particle-axis nematic order parameter. TPs with r ) 1. The expansion simulations revealed a firstorder transition from a crystal phase to a cubatic phase (for convenience, these results are shown later, in Figure 11). The cubatic phase turns into the isotropic phase on further reducing the pressure at a volume fraction of about 0.5 through a firstorder transition. The equation of state for the compression simulations revealed a first-order isotropic-cubatic transition but did not show any break in the curve at the anticipated pressure for the solid-phase transition. The cubatic phase is characterized by orientational order along the three axes but the packing is not as dense as in the crystal phase. As in ref 22, snapshots of the cubatic phase (not shown) do not reveal much structural difference from those of the solid phase at higher pressures. At the solid-cubatic phase transition, however, the peaks in the particle distribution function of the solid are twice as high as those in the cubatic phase, indicating less translational mobility in the former.22 The particle distribution function (not shown) revealed a high degree of order along three axes in both cubatic and solid phases. The cubatic orientational correlation function gO(r) was close to 1 for distance r > 4σ (results not shown) indicating cubatic order in the range of densities corresponding to the cubatic regime. A Parrinello-Rahman simulation allowing changes in the shape of the box gave similar results. This indicates that the cubatic phase is not an artifact of the tetragonal box. Note that the first-order character of the isotropic-cubatic phase transition has also been confirmed via multicanonical-type simulations designed to pinpoint the coexistence conditions for a system of 3 × 3 × 3 cuboids.32 TPs with r < 1. Typical results are shown for r ) 0.5 in Figures 6-9. On expansion from the crystal phase (Figure 6), the particles formed a columnar mesophase (Figure 7) with order along the two (long) equivalent axes. As the pressure is reduced, this columnar phase is transformed into a smectic-like phase (Figure 8) characterized by the presence of layers along one of the long axes. In this smectic-like phase, the particles are free to rotate within the layers. However, this phase is not strictly a smectic phase as a cross-section of the system reveals that TPs in each layer are not fully uncorrelated but tend to stack together, with the stacks oriented in the same direction (Figure 8), suggestive of tetratic order. Changes in the simulation box shape were allowed in the simulations in the smectic-like and columnar

John and Escobedo

Figure 7. Particle distribution function (top) and snapshots of the mesophase (bottom) formed by the r ) 0.5 TPs at P ) 0.5. The snapshots give the front view (left) and the top view (right) of this phase.

Figure 8. Particle distribution function (top) and snapshots of the smectic-like phase (bottom) for the r ) 0.5 TPs at P ) 0.35. The snapshots give the front view (left) and the top view (right) of this phase.

Figure 9. Equation of state for the r ) 0.5 TPs.

phases. Both phases remained stable. The equation of state for TPs with r ) 0.5 is given in Figure 9. In both the crystal phase (Figure 6) and the columnar phase (Figure 7), the particle

Phase Behavior of Cuboids

Figure 10. (a) Snapshot of the TPs with r ) 0.25 at P ) 1.0 from the compression runs. The TPs pack into stacks reminiscent of the cubatic phase of the cut spheres from ref 19. (b) Snapshot of the parquet phase observed in the expansion run of the r ) 1.5 system of TPs at P ) 0.12.

distribution function along the direction of the two long axes showed sharp peaks at intervals of 4, which correspond to the length of the TP long side (a ) b ) 4σ). In the solid phase (Figure 6), the particle distribution function along the short axis also shows periodic peaks at intervals of c ) 2 (the dimension of the short axis). Such peaks are absent in the smectic-like mesophase (Figure 8) where peaks occur only along one of the long axes of the TP, indicating layering in that direction. Correspondingly, one of the equivalent long axes of the TP shows alignment via the selected particle-axis nematic order parameter, while the remaining particle long axis shows complete disorder. This confirms that the particles are free to rotate within each layer. On compression, TPs with r ) 0.25 packed into small stacks that were roughly perpendicular to one another (Figure 10a). This is reminiscent of the cubatic phase of the cut spheres of aspect ratio 0.2.19 The Parquet Phase. Compression runs for r > 1 TPs starting from an isotropic phase often led to a parquet phase through a first-order transition. The parquet phase is characterized by stacks in which a few TPs are aligned in one direction but these stacks are oriented perpendicular to one another (Figure 10b). Further compression of this phase only resulted in a more dense packing with no phase change. The parquet phase was observed only at one pressure in the expansion runs. It was found in the intermediate region (P ) 0.12) between isotropic and cubatic phases for TPs with r ) 1.5. For such a phase, the “selected director” method described in Section IIC yields a nonzero cubatic order parameter and a vanishing nematic order parameter, indicative of cubatic-like orientational order. The “selected particle-axis” nematic order parameter for this phase with respect to the long axis of the TP was found to be zero indicating that the long axes are randomly oriented. A similar phase has been reported by Schilling and Frenkel as a polycrystalline phase.33 It is likely that in most other cases where the parquet phase forms upon compression of an isotropic fluid, it represents a metastable phase in which the system becomes trapped. It is also likely, however, that such a phase is relevant to experimental work, where trapping in glassy structures is prevalent (note that the moves in our MC simulations approximate the diffusive dynamics in dense systems). To summarize, the phase diagram in Figure 2 shows the results obtained from the expansion runs. The dashed line separating the columnar and solid phases demarks an approximate phase boundary. For r > 4, solid, columnar, nematic, smectic, and isotropic phases were observed. For 1 < r e 4, the only liquid crystalline phases observed were the columnar and the CL mesophase; for r ) 1.5, however, the parquet phase was briefly encountered. TPs with r < 1 gave rise to two additional mesophases: a smectic-like phase (P1) and a columnar phase.

J. Phys. Chem. B, Vol. 109, No. 48, 2005 23013 Comparison between Hard TPs and Spherocylinders. For spherocylinders,15 their aspect ratio is characterize by L/D, where L is the length and D is the diameter of the spherocylinder (so that r ≈ L/D + 1). For L/D e 0.35 ( 0.05, the isotropic fluid freezes to form a plastic crystal (rotator phase). The rotator phase undergoes a first-order transition to an orientationally ordered phase at higher pressures. For TPs with r ) 0.25, a smecticlike mesophase and a columnar phase are observed. For spherocylinders, the smectic phase becomes stable at L/D ) 3.1 and the nematic phase at L/D ) 3.5. For the TPs, the smectic phase becomes stable at r ∼ 4 and the nematic phase at r > 5. For spherocylinders, the nematic-smectic transition and the smectic-solid transition are both first order. Finally, several of the mesophases observed in the TPs are absent in the spherocylinders. Comparison between the Simulated Freely Rotating TPs and TPs with Restricted Orientations (Zwanzig Model). Our results agree qualitatively with the predictions of the FMF theory of Martinez-Raton25 for prolate parallelepipeds with r > 5 wherein the theory also predicts the nematic, smectic, and oriented solid phases. In our simulations, however, a columnar phase is also observed between the smectic and the oriented solid phases. In the region 3 < r < 5, the FMF theory predicts a columnar phase and a discotic smectic phase (besides the isotropic phase). The discotic smectic is a layered phase with the long axes of the parallelepipeds lying within the layers but lacking orientational order in each layer. Although a discotic smectic phase was not detected in the simulations, the CL mesophase was observed wherein either a long or a short axis of the parallelepiped is oriented along the directors of alignment. In such a CL mesophase, significant layering is achieved without constraining the particle long axis to always be parallel to the layer planes. Arguably, the CL mesophase and the discotic smectic are related structures with comparable partial translational order (layering) but with the former having a larger orientational entropy (due to the less constrained orientations of the long particle axes). In the case of 1 < r < 3, the FMF theory predicts that the isotropic phase transitions into a plastic solid at high enough concentrations. A plastic solid has a high degree of translational order but no orientational order. For r ) 1, the theory does not predict the cubatic phase (encountered in simulations) which, in contrast to a plastic solid, has a high degree of orientational order and little translational order. The constrained orientational degrees of freedom in the Zwanzig model assumed by the FMF theory likely preclude the existence of a cubatic-like phase, which would necessarily imply zero orientational entropy. In general, the predictions of the FMF theory are not consistent with our simulation results for oblate parallelepipeds. For 0.4 < r < 1, the theory predicts the existence of the isotropic and plastic solid phases, while for 0.2 < r < 0.4 the theory predicts the appearance of an additional discotic smectic phase. Overall, the differences between the phase behavior predicted by the FMF theory and that observed in simulation can be attributed to the two features of the simulated particles lacking in the Zwanzig TPs: their free rotation and their “bumpy” surfaces and edges. While the former difference is likely the most crucial, the latter must play an increasingly important role in our oblate TPs wherein the number of spheres constituting a particle decreases substantially as r f 0. In the next section, it is shown that at least for r ) 1, even perfectly smooth cubes do not exhibit the phase transitions predicted for the Zwanzig cubes.

23014 J. Phys. Chem. B, Vol. 109, No. 48, 2005

Figure 11. Equation of state curves for multisite r ) 1 cuboids and for perfect cubes. The cuboids have different degrees of coarseness depending on the number of spheres per side used.

B. Effect of Particle Roughness on the Cubatic Phase. For cuboids with r ) 1, the equation of state obtained on varying the number of spheres per side and thus the surface roughness showed a similar trend in all three cases (Figure 11). The isotropic-cubatic transition occurs at higher pressures for the more rounded cuboids. The volume fraction of the cuboidal particles (r ) 1) in Figure 11 was evaluated considering that the effective dimension a of a side of the particle is less than the number of spheres per side because of the roughness of the cuboid facets. This value of a was found by calculating the minimum distance between two cuboids constrained to lie on a flat surface, i.e., for the 4 × 4 × 4 cuboids a ) 3.9.22 The equation of state for perfect cubes is also plotted in Figure 11 where a first-order isotropic-cubatic phase transition is apparent. The pressure required to first form a cubatic phase was the highest for the cuboids with two spheres per side and the lowest for the perfect cubes. This can be explained by the fact that for a given volume fraction more rounded cuboids are more likely to rotate and disorder. The results obtained from the compression and expansion runs of perfect cubes showed very close agreement. No evidence of first- or second-order transition was observed between the cubatic and the crystal phases (Figure 11). Snapshots of the isotropic, cubatic, and solid phases are given in Figure 12, where the mixing of the original color-banded structure of the solid in part d is indicative of translational disorder in part b and (to a lesser extent) in part c. The particle distribution function for the cubatic phase at a particular pressure shows similar ordering along all three axes of the cuboids (results not shown). Figure 13 shows the radial distribution function for the perfect cubes before and after the anticipated cubatic-solid boundary (based on the results of cuboids). Although the radial distribution function at a volume fraction of 0.8 resembles that of a solid phase, it is unclear whether a continuous transition exists to the crystal-like phase. Starting with a well-ordered cubic structure where the rows of cubes are color coded, one can find after long simulations at high densities (in the anticipated solid regime) that some rows appear to “slip” with respect to others. This result contrasts that of the multisite cuboids whose bumpiness could be responsible for inducing crystallization. Groh and Mulder34 found that hard parallel cubes exhibit a solid phase and not a columnar phase. Since at high densities our

John and Escobedo

Figure 12. Snapshots from the expansion run of perfect cubes at different volume fractions: 0.40 (a), 0.52 (b), 0.63 (c), and 0.91 (d). The system is isotropic in part a, cubatic in parts b and c, and solidlike in part d. The cubes are colored based on their positions in a colorbanded solid similar to that shown in part d.

Figure 13. Radial distribution function g(r) for perfect cubes for volume fractions 0.52 (cubatic phase), 0.77, and 0.8. Note that the latter state exhibits a solidlike g(r).

Figure 14. Cubatic and nematic order parameters for the r ) 1 cuboids from configurations collected in expansion runs.

cubes are highly constrained, they would be expected to behave pretty much like parallel cubes and therefore to form a crystal phase. Clearly, further study is needed to clarify the nature of the crystal-like phase exhibited by our system. The snapshots of the cubatic phase indicate orientational order but limited translational order. The cubatic order parameter showed a sharp increase in the cubatic region from the isotropic region where the value was close to zero (Figure 14). The “selected director” nematic order parameter remained close to zero in both regions

Phase Behavior of Cuboids

J. Phys. Chem. B, Vol. 109, No. 48, 2005 23015

(Figure 14). Together, these order parameters provide a signature of the cubatic phase.

This material is available free of charge via the Internet at http:// pubs.acs.org.

IV. Conclusions

References and Notes

The phase diagram for TPs of different aspect ratio was mapped via Monte Carlo simulations. It was found that TPs have a very rich phase diagram with the presence of many unusual mesophases. The TPs of high aspect ratios at low to intermediate densities exhibit a phase behavior similar to that of spherocylinders and to that predicted for Zwanzig TPs. The influence of the edges and corners of the TPs on phase behavior decreases at high aspect ratios and lower concentrations. The mesophases obtained at lower aspect ratios, namely the smecticlike, columnar, and cubatic phases, are particular to our freely rotating TPs as they are not observed in spherocylinders or predicted in Zwanzig TPs. Most of the results were obtained from simulations where the TPs were expanded from the crystal state. Compression from an isotropic phase often gave rise to a parquet phase believed to be metastable in all cases except for r ) 1.5. In examining the effect of particle roughness on the phase behavior of cuboids with r ) 1, it was found that the isotropiccubatic transition for more rounded cuboids required higher pressures compared to more perfect (i.e., with sharper corners and edges) cuboids. The occurrence of the cubatic phase in all these cases provides evidence of its robustness with respect to edge roundedness. The equation of state obtained from the simulation of perfect cubes indicated good agreement with that of the multisite cuboids except for the fact that no discontinuity was observed at the cubatic-solid transition. The layers of cubes at high densities appeared to slip with respect to one another raising questions on the character of the solid phase.

(1) Sun, Y.; Xia, Y. Science 2002, 298, 2176. (2) Barnard, A. S.; Zapol, P. J. Chem. Phys. 2004, 121, 4276. (3) Narayanan, R.; El-Sayed, M. A. Nano Letters 2004, 4, 1343. (4) Glotzer, S. C.; Solomon, M. J.; Kotov, N. A. AIChE J. 2004, 50, 2978. (5) Jin, R.; Cao, Y.; Mirkin, C. A.; Kelly, K. L.; Schatz, G. C.; Zheng, J. G. Science 2001, 294, 1901. (6) Puntes, V. F.; Krishnan, K. M.; Alivisatos, A. P. Science 2001, 291, 2115. (7) Manna, L.; Scher, E. C.; Alivisatos, A. P. J. Am. Chem. Soc. 2000, 122, 12700. (8) Zhang, H.; Wu, G.; Chen, X. Langmuir 2005, 21, 4281. (9) Teranishi, T.; Inoue, Y.; Nakaya, M.; Oumi, Y.; Sano, T. J. Am. Chem. Soc. 2004, 126, 9914. (10) Wang, Z. L.; Yin, J. S. Mater. Sci. Eng. A 2000, 286, 39. (11) Ahmadi, T. S.; Wang, Z. L.; Green, T. C.; Henglein, A.; El-Sayed, M. A. Science 1996, 272, 1924. (12) Gou, L.; Murphy, C. J. Nano Letters 2003, 3, 231. (13) Xia, Y.; Gates, B.; Li, Z. AdV. Mater. 2001, 13, 409. (14) Colloids and colloid assemblies; Caruso, F., Ed.; Wiley-VCH: Weinheim, Germany, 2004. (15) Bolhuis, P.; Frenkel, D. J. Chem. Phys. 1997, 106, 666. (16) Camp, P. J.; Mason, C. P.; Allen, M. P.; Khare, A. A.; Kofke, D. A. J. Chem. Phys. 1996, 105, 2837. Camp, P. J.; Allen, M. P. J. Chem. Phys. 1997, 106, 6681. (17) Eppenga, R.; Frenkel, D. Mol. Phys. 1984, 52, 1303. (18) Bates, M. A.; Frenkel, D. Phys. ReV. E 1998, 57, 4824. (19) Veerman, J. A. C.; Frenkel, D. Phys. ReV. A 1992, 45, 5632. (20) Blaak, R.; Frenkel, D.; Mulder, B. M. J. Chem. Phys. 1999, 110, 11652. (21) Elliott, J. A.; Windle, A. H. J. Chem. Phys. 2000, 113, 10367. (22) John, B. S.; Stroock, A.; Escobedo, F. A. J. Chem. Phys. 2004, 120, 9383. (23) Blaak, R.; Mulder, B. Phys. ReV. E 1998, 58, 5873. (24) Frenkel, D. In Liquids, Freezing and the Glass Transition; Hansen, J. P., Levesque, D., Zinn-Justin, J., Eds.; North-Holland: Amsterdam, The Netherlands, 1991; pp 689-762. (25) Martinez-Raton, Y. Phys. ReV. E 2004, 69, 061712. (26) Allen, M. P. Liq. Cryst. 1990, 8, 499. (27) Parrinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196. (28) Frenkel, D.; Smit, B. Understanding Molecular Simulation. From Algorithms to Applications, 2nd ed.; Academic Press: New York, 2002. (29) Gottschalk, S.; Lin, M. C.; Manocha, D. Comput. Graph. 1996, 30, 171-180. (30) Blaak, R. FOM-Institute for Atomic and Molecular Physics, Ph.D. Thesis, Kruislaan, Amsterdam, The Netherlands, 1997. (31) Vega, C.; McBride, C.; MacDowell, L. G. J. Chem. Phys. 2001, 115, 4203. (32) Gospodinov, I. D.; Escobedo, F. A. J. Chem. Phys. 2005, 122, 164103. (33) Schilling, T.; Frenkel, D. J. Phys.: Condens. Matter 2004, 16, S2029. (34) Groh, B.; Mulder, B. J. Chem. Phys. 2001, 114, 3653.

Acknowledgment. The authors are grateful for financial support from the Cornell Center for Materials Research (CCMR) under the SEED program. CCMR is funded by the National Science Foundation. Partial support from an ACS-PRF grant is gratefully acknowledged. The authors would also like to thank Cornell Professors A. Stroock (CBE) for useful discussions and Prof. K. Bala (CS) for the reference on Oriented Bounding Boxes. We are also thankful to S. Fraden (Brandeis University) for useful suggestions and to C. Abreu for help with the visualization software. Supporting Information Available: Additional details regarding the method to test for overlaps between perfect cubes.