Monte Carlo simulations were used to study the influence of particle aspect ratio on the kinetics and phase behavior of hard gyrobifastigia (GBF). First, the formation of a highly anisotropic nucleus shape in the isotropic-to-crystal transition in regular GBF is explained by the differences in interfacial free energies of various crystal planes and the nucleus geometry predicted by the Wulff construction. GBF-related shapes with various aspect ratios were then studied, mapping their equations of state, determining phase coexistence conditions via interfacial pinning, and computing nucleation free-energy barriers via umbrella sampling using suitable order parameters. Our simulations reveal a reduction of the kinetic barrier for isotropic–crystal transition upon an increase in aspect ratio, and that for highly oblate and prolate aspect ratios, an intermediate nematic phase is stabilized. Our results and observations also support two conjectures for the formation of the crystalline state from the isotropic phase: that low phase free energies at the ordering phase transition correlate with low transition barriers and that the emergence of a mesophase provides a steppingstone that expedites crystallization.

## I. INTRODUCTION

Polyhedral nanoparticles can now be synthesized with an unprecedented control through synthetic routes developed in the last few decades, opening the door for their use as building blocks of new superstructures.^{1–3} Indeed, when concentrated, such particles can form crystalline assemblies that possess unique optical properties and find applications in photonics and plasmonics.^{2,4–6} Computational studies have been useful for predicting the phase behavior of polyhedral particles as a function of their shape,^{7,8} with many such predictions having already been experimentally corroborated.^{9}

Particle aspect ratio (AR) is known to be an important determinant of the phase behavior of colloidal nanoparticles.^{7,8} Generally, higher particle anisotropy favors the formation of lyotropic liquid crystalline phases.^{10} For example, cuboids of very large or very small aspect ratios can stabilize a variety of mesophases, such as nematic, smectic, columnar, and cubatic phases.^{11–14} Low-anisotropy, low-asphericity particles (e.g., cuboctahedra, truncated cubes, rhombic dodecahedra)^{15,16} tend to stabilize rotator mesophases, i.e., solid phases where particles have translation order but limited or no orientational order. For aspherical shapes with low particle anisotropy, such as regular triangular prisms,^{7} octahedra,^{9,17} and gyrobifastigia (GBF),^{7} no mesophases are observed mediating their isotropic and crystalline phases. When it comes to kinetics of disorder to order phase transitions, this latter category of shapes is expected to exhibit the highest nucleation free energy barriers (Δ*G**) at a given supersaturation (Δ*μ*_{od}), as reflected by the difficulty to spontaneously nucleate the crystalline state in molecular simulations. In fact, for high enough barriers, no spontaneous transition to ordered phases is observed when disordered phases are compressed in unbiased simulations (which result instead in dense, kinetically arrested disordered states).^{7} Altogether, the above observations for particles with different types of shape anisotropy and asphericity are in line with the conjecture that the existence of a mesophase can kinetically ease the transition.^{18} This conjecture has been supported by the phase behavior observed in both hard polyhedra^{16,19,20} and soft particles forming block-copolymer type mesophases.^{21}

In this study, we focus on the phase behavior and kinetics of disorder–order phase transition in hard gyrobifastigia (GBF). GBF is chosen for being a particle shape that embodies some unique characteristics and challenges among other faceted particles. GBF is one of a handful of regular convex space-filling polyhedra but is remarkable in that it is quite asymmetric in shape and forms an unusual ABCD lattice (also known as *α*-Sn),^{22} implying a multi-layer level of local cooperativity for crystal nucleation. It is also interesting to note that the GBF polyhedral shape exists in many molecular and solid-state structures.^{23} In a previous study,^{24} we calculated a very large Δ*G** for isotropic–crystal transition (compared to other shapes at a given Δ*μ*_{od}), falling in line with other shapes that do not exhibit a mesophase behavior. This observation was partly explained by the discord between the locally favored structures in the isotropic state and the arrangement of particles in the crystal. Furthermore, nucleus-size pinning simulations, where a nucleus size is maintained and allowed to converge to its equilibrium shape, revealed a highly anisotropic nucleus shape with aspect ratio of approximately two. It is unclear whether this nucleus shape anisotropy plays any role on the ordering transition kinetic mechanism and free-energy barrier. In this study, we investigate those results by first performing direct measurements of the disorder–order interfacial free energy of various crystal planes using the cleaving walls method.^{25,26} We then use these results to perform a Wulff construction to predict a nucleus geometry to corroborate our earlier findings.

The self-assembly of anisotropic particles has been of great interest in the simulation literature.^{27–30} Systematic studies of the phase behavior of colloidal rods (spherocylinders) as a function of their aspect ratio (AR) have revealed the existence of liquid crystalline mesophases for high AR values.^{31} Similar studies for faceted particles reveal a richer mesophase behavior.^{11–14} The kinetics of colloidal rods has been found to be rather nuanced; e.g., short rods (AR = 2) may follow nucleation and growth at moderate supersaturations but will get kinetically arrested with a large number of crystallites at higher supersaturations.^{30} At extremely high supersaturation, the system gets arrested in a glassy state. For longer rods (AR = 3.4), it was reported that the isotropic–smectic transition is suppressed due to spinodal instability.^{30} This motivated us to ponder about the effect that particle faceting could have on the relation between aspect ratio and disorder-to-order kinetics. Accordingly, in this study, we selected GBF-shaped particles of different aspect ratios to investigate the effect of particle anisotropy on the kinetics of isotropic–crystal transition and on the nucleus shape anisotropy (as previously reported for regular GBF).^{24}

## II. METHODS

### A. Model

We restrict our study to the case of athermal systems. Accordingly, any two particles *i* and *j* experience a hard pair-potential given by

The overlap is detected by using the separating axis theorem.^{12,32} A GBF (gyrobifastigium) is composed of two regular triangular prisms (*fastigium* pl. *fastigia* meaning roof) attached at a square base with a twist. The aspect ratio (AR) is defined as

where *h* is the height of GBF and *a* is the side of the square base. For a regular gyrobifastigium, AR = $3$. The particle AR in the study is varied by elongating the triangular faces into isosceles triangles while keeping the square base of attachment unaltered. If the side of the square base is two length units, then the height of the triangular facet is equal to the aspect ratio. We study GBF-like particles with six aspect ratios $\u2208{13,13,1,3,23\u22121,3}$, as shown in Fig. 1.

### B. Monte Carlo simulations

We conducted Metropolis^{33} Monte Carlo (MC) simulations in either the canonical (*NVT*) or the isothermal-isobaric (*NpT*) ensemble as necessary, where *N* is the total number of particles, *V* is the volume of the system, *p* is the pressure, and *T* is the temperature. We use scaled units consistent with our previous studies,^{7} with lengths scaled by the circumradius (*a*_{c}) of the shape. Thus, the dimensionless pressure is given by $p=\beta paac3$, where *p*_{a} is the unscaled pressure and $\beta =1kBT$, where *k*_{B} is Boltzmann’s constant. The chemical potential (*μ*) and free energy (Δ*G*) are scaled by *k*_{B}*T*, and the supersaturation is defined as

where *μ*_{o} and *μ*_{d} are the chemical potentials associated with ordered and disordered phases, respectively. The simulations used periodic boundary conditions to mimic bulk behavior. Each MC cycle included *N* translation, *N* rotation, *N* flip, and two isotropic volume moves (for *NpT* ensemble runs only). Flip moves turn over the particle along the axis perpendicular to the prism base. Unless otherwise stated, all configurations were originally obtained for a cuboidal simulation box with *N* = 1728 particles containing 12 layers of 12 × 12 particles arranged in an ABCD crystal lattice of the GBF honeycomb.

### C. Cleaving wall method

We use the cleaving walls method^{24,25} implemented for our hard particle Monte Carlo simulations to directly calculate *γ* for various crystal planes. More details on the method implementation and the initial setup are described in a previous study.^{34} The method essentially calculates the reversible work done per unit area of interface created through the following three steps:

Cleaving: Configurations are simulated for both isotropic and crystalline phases at coexistence conditions, set up such that the boxes have identical cross-sectional dimensions along a desired plane (called the cleaving plane). For the crystal phase, the configurations are generated such that the cleaving plane is in perfect alignment with the crystal plane of interest. A cleaving wall potential is then applied, which gradually moves the particles to create a gap along the cleaving plane wider than twice the circumradius of the particle shape, such that hard particles do not interact across the gap.

Transposition: The two cleaved systems are then transposed into one system such that the isotropic and crystalline phases are placed across the gap. Since particles do not interact across the gap, this step has no work contribution.

Merging: The cleaving potential is gradually removed to allow interactions between particles from the two phases, thus creating an interfacial system with two interfaces.

In this study, we used as cleaving potential flat walls that interact with particle centroids; these walls start being positioned at the midplane and are then gradually moved to create an excluded gap along the cleaving plane. We measured the pressure on the wall using virtual perturbations of the order of 10^{−3} in reduced length units in the wall positions. The pressure was integrated with respect to the wall position to obtain the reversible work imparted on the system in steps (i) and (iii). This reversible work done per unit area of the interface created is the interfacial free energy *γ*. We performed all calculations at the disorder–order coexistence pressures. For all configurations, we first expanded from the tessellated crystal with appropriate crystal orientation at high pressure of *p* = 32 to the coexistence pressure *p*_{ic,GBF} = 10.8. This resulting configuration was then melted to the isotropic phase at *p* = 2 and recompressed to the coexistence pressure with the cleaving plane cross section constrained throughout.

### D. Coexistence pressures and driving forces

The isotropic–crystal coexistence pressure for regular GBFs under scaled units is *p*_{ic,GBF} = 10.8 as calculated in a previous study.^{34} For non-regular GBFs, we calculate the isotropic–crystal coexistence pressures using the interfacial pinning method.^{35,36} Driving forces Δ*μ*_{od} were calculated using thermodynamic integration from the coexistence conditions.^{37} Absolute chemical potentials *μ* at coexistence were calculated via thermodynamic integration over density from the ideal gas limit as described elsewhere in the literature,^{37,38} with system-specific details provided in the supplementary material.

### E. Order parameters

Steinhardt *et al.*^{39} order parameters are used to capture both local translational and orientational order. We use *q*_{6} to quantify the translational order, whose implementation requires a suitable choice of the nearest-neighbor cutoff distance *r*_{c}. More details on our use of this order parameter are provided in previous studies.^{16,24} *q*_{6} values for any two particles (*i* and *j*) can be used to evaluate translational correlation $dqi,j$,

where the asterisk (*) denotes the complex conjugate.

We use the *P*_{2} orientational order parameter^{40} to measure the orientational correlation between two particles (*i* and *j*),

where *θ*_{ij} is the angle between the particles’ z-axis unit vectors, z-axis being the long axis of the regular GBF.

The definition of nucleus size *n* is similar to that in our prior work^{24} and is based on three parameters the determine what neighboring ordered particles are clustered: the neighbor cutoff distance *r*_{c}, the orientational correlation cutoff *P*_{2,c}, and the connection threshold *ζ*_{c}. The neighbor cutoff *r*_{c} is defined as $1.4\xd7min(2,1+AR2)$, while the other parameters are optimized based on the mislabeling criteria.^{41} The specific parameter values used for various aspect ratios are given in the supplementary material.

### F. Umbrella sampling

Umbrella sampling (US) simulations were performed to determine free energy barriers for various ordering transitions. The recipe remains largely identical to that used in our previous publications;^{24} here, we specify key details of the procedure. The size of the largest nucleus (*n*) in the system was used as the order parameter along which the free energy calculation was performed. The transition path along *n* was divided into overlapping equal-sized windows. The process is started with a window where a range of nucleus sizes are readily sampled in an unbiased simulation, i.e., around the metastable basin. A configuration with one of the largest nuclei thus found is used to launch the next US window that is now set up to restrict the nucleus to remain in a range of larger *n* values (than in the previous window). The process is repeated until the free-energy (found as explained shortly) starts decreasing, signaling that the barrier top has been surpassed. The size of each window varied depending on the conditions and the size of the system (see the supplementary material). Unbiased isothermal-isobaric simulations were performed at each window to collect histograms of the relative frequency with which different *n* values are visited (rejecting any move that would take the system outside the window bounds); ** n** values were sampled every two MC cycles. Free energy curves extracted from the frequency histograms collected from individual sections are stitched together by matching values in the middle of the overlapping regions of neighboring windows. The values Δ

*G*and the final barrier Δ

*G** are calculated in reference to the free energy of the homogeneous disordered phase using unbiased simulations to estimate the probability to form a crystal nucleus of size

*n*, as described in earlier studies.

^{37,42}We also made sure that the nuclei did not become too large compared to the system size (preferably no greater than

*N*/10) and interact with their periodic images, the latter being particularly likely for highly anisotropic nuclei.

## III. RESULTS AND DISCUSSION

### A. Anisotropic nuclei for GBF

In this section, we investigate the origins of the anisotropic shape of the GBF crystal nucleus during the isotropic–crystal transition. Our earlier studies revealed that the aspect ratio of the crystal nucleus in regular GBF did not appreciably change with size, with the shape getting increasingly better defined for larger nuclei. This motivates us to investigate the existence of an intrinsic shape of the nucleus, defined by the properties of the crystal–isotropic interface. In this section, we predict a simple nucleus’ polyhedral geometry by using a Wulff construction based on calculated interfacial free energies of distinct crystals planes.

#### 1. Crystal planes of interest

For simplicity, we assume that the surface of the crystal nucleus is comprised of low index crystal planes that are fundamentally flat. A given crystal may have several closed packed crystal planes, close packing being defined by the absence of steps and sheet-like arrangement of particles along the plane.^{43} Since interfacial free energy calculations can be expensive, we narrow down the relevant crystal planes to those prominently present in the simulation-observed nucleus shape (Fig. 2). We hence choose the three crystal planes (100), (001), and $(1201)$. The last plane arises from the lattice geometry and has a flatter profile than the other two. This plane will henceforth be referred to as the *s* (slant) plane.

#### 2. Interfacial free energy calculations

We used different system sizes to simulate different crystal planes to avoid artifacts during the cleaving walls’ calculation. In particular, the (001) plane required the use of a longer crystal phase to minimize the warping that this phase tends to exhibit upon compression. All interfacial free energies were calculated at the coexistence pressure *p*_{ic,GBF} = 10.8. Final merged configurations and pressure vs wall position plots for the three crystal planes are shown in Figs. 3–5, and results are summarized in Table I. The pressure profiles can be interpreted in relation to the associated reversible work performed upon the system, which is equal to the area under the curves in each step. Thus, the difference in the area under the solid line and the dotted line correlates with the net reversible work done to get a phase from a bulk, uncleaved state to a state where it is interfacing with the other phase. This reversible work per unit interfacial area created is equal to *γ*. Note that since the pressure measured for these virtual cleaving planes can vary depending upon the arrangement of particles at the interface, one may not expect equal pressures on either side of the cleaving plane upon merging. This apparent discrepancy was shown to have no effect on the calculated interfacial free energy when compared to more realistic walls that interact with the particles’ facets and edges rather than with their centroids.^{34}

#### 3. Relation of interfacial roughness and interfacial free energy

In an earlier study,^{34} we conjectured that the interfacial free energy correlates with the roughness of the crystal plane. The results for the GBF crystal planes substantiate that conjecture as visual inspection [Fig. 6(a)] places the apparent roughness of facets in the order: (001) > (100) > *s*, which is the same as that of decreasing *γ*. This could be understood by considering how much excluded volume for the isotropic phase a given crystal interface creates due to its roughness. This effect can be captured by calculating the interfacial potential of mean force (PMF) experienced by a free particle at a given distance *x* from a perfect crystal interface. The interface location is assigned based on the limiting proximity that a particle can get to the surface (being forbidden from occupying any space closer than that). Such a limiting proximity would correspond to the point where *PMF* → ∞. However, here we use PMF = 7 *k*_{B}*T* as an effective exclusion threshold and as surrogate of an impractical infinite *PMF* value. More details of the method are provided in an earlier study.^{34} The interfacial PMF for all three facets is calculated as a function of distance from the interface, as shown in Fig. 6(b). We find that the *s* facet unequivocally has the shortest range of influence. (001) is has higher *PMF* at short distances (*x* < 0.5) than (100) but decays to zero quicker at longer distances, indicating a shorter range of disruption overall.

#### 4. Wulff construction

We can utilize the interfacial free energies calculated for the three crystal planes to perform a Wulff construction.^{43–45} The construction applies to the steady state geometry of a crystalline nucleus, and is based on the relation that the perpendicular distance of a crystal facet *i* from the center of mass (*h*_{i}) is proportional to its interfacial free energy *γ*_{i}, i.e., *h*_{i} ∝ *γ*_{i}. This relation can be used to determine the relative positions of crystal planes, with the enclosed volume representing the nucleus shape. The Wulff construction for a GBF crystal nucleus is shown in Fig. 7(a), where we take the perpendicular distance to the least energetically expensive facet *s* as *h*_{s} = 1, with subsequent planes placed at a distance proportional to their corresponding *γ*. We also leverage the symmetry of the crystal to place equivalent crystal planes, which allows us to find a closed polyhedron [Fig. 7(b)]. We find that the crystal plane (100) does not appear in the predicted nucleus shape, resulting in a shape we could describe as an isosceles square bipyramid with two vertices truncated along the long axis. The shape is remarkably similar to the largest nuclei observed in our simulations (see, e.g., Fig. 2). Furthermore, the aspect ratio of this predicted nucleus geometry is calculated to be 2.05, which is consistent with that from the nucleus size-pinning simulations, which ranged between 1.9 and 2.1.^{24} Based on relative areas of *s* and (001) facets in the Wulff construction, we calculated the average interfacial free energy of the nucleus at coexistence to be *γ*_{avg} = 0.79, which is also in agreement with the previous predictions based on classical nucleation theory.^{24} These calculations reveal that significant anisotropy in interfacial free energy along different crystal planes can result in anisotropic crystalline nuclei for hard polyhedral particles. Ultimately, the crystal grows adopting a shape that minimizes the free-energy cost associated with the crystal–liquid interface, which for GBF means a shape with a preponderance of the *s* facet.

### B. Effect of particle aspect ratio

#### 1. Phase behavior

The phase behavior for the GBF family of shapes is illustrated in Fig. 8 where we plot regions of various phases bounded by their coexistence volume fractions *ϕ* as a function of aspect ratio. For all aspect ratios considered, we observe the existence of an isotropic and an ABCD crystalline phase. For a broad range of aspect ratios, the phase behavior remains qualitatively identical to that of regular GBF with the crystal phase being the only ordered phase that forms. For extreme aspect ratios, we observe the stabilization of a nematic mesophase at intermediate packing fractions by both the compression of the isotropic phase and the expansion of the crystalline phase (albeit with some hysteresis around the coexistence pressure). This is expected, since at higher aspect ratios, particles are more anisotropic and, hence, become closer to the geometric category of oblate or prolate shapes that form lyotropic liquid crystalline mesophases.^{7}

Further details and equations of state for individual aspect ratios are provided in the supplementary material.

#### 2. Crystallization transition kinetics

We elucidate the effect of AR on nucleation kinetics for isotropic–crystal transition for three select aspect ratios, which exhibit no mesophase behavior: $AR\u2208{1,3,23\u22121}$. As in prior studies, we use the nucleation barrier Δ*G** for a set of degrees of supersaturation (Δ*μ*_{od}) as a measure of how difficult it is to nucleate the crystal phase, with larger Δ*G** correlating with the slower nucleation rate. This approximation is rooted in the classical nucleation theory wherein the rate constant depends exponentially on −Δ*G** (in *k*_{B}*T* units) and only linearly on a prefactor, whose values, while different for different AR systems, are expected to be within a similar order of magnitude. We have previously reported the nucleation barriers for $AR=3$,^{24} and hence, here we computed the values for the other two aspect ratios. Figure 9 shows the free-energy profiles as a function of *n* for AR = $23\u22121$ for different Δ*μ*_{od} values, and Fig. 10(a) summarizes the results for Δ*G** for the three AR of interest. Figure 10(a) reveals a clear trend wherein at a given supersaturation longer particles have lower nucleation barrier. In fact, the barrier for $AR=23\u22121$ is reduced to the point that the isotropic phase transitions spontaneously to the crystal phase in an unbiased simulation.

Interestingly, these trends in nucleation barriers correlate with the absolute chemical potential at the isotropic–crystal phase coexistence calculated for various aspect ratios shown in Fig. 10(b): Systems with a larger coexistence chemical potential have a higher crystallization nucleation barrier (for comparable degree of supersaturation). In previous work,^{47} it was argued that in considering a family of particle shapes capable of forming the same target ordered phase, the one that has the lowest free-energy at the disorder-to-order transition has the optimal shape to achieve such a target. In other words, if one is coming from the isotropic phase, the crystal phase is most stable at the point where it can be reached with its lowest free-energy, which in this case occurs for the more elongated GBF. It is then conjectured here that the most stable crystal phase would reflect some intrinsic affinity of the particle shape to arrange into their crystalline structure and should hence be easier to reach or nucleate from the isotropic phase, i.e., exhibit smaller Δ*G**. In addition, note that in this case, the kinetic optimality in AR coincides with the system with the lowest packing fraction of both the isotropic and crystal phases at the transition, a fact that may be associated with a lower interfacial tension and easier rearrangement of particles at the nucleus interface.

The anisotropic nucleus geometry observed for regular GBF ($AR=3$) also systematically varies upon changes in the particle AR. In general, the crystalline nucleus aspect ratio increases with particle AR (Fig. 11). The nucleus aspect ratio is defined via the principal moments of the nucleus inertial tensor as the ratio of the longest axis to the average of the other two. We performed this calculation for nuclei close to the critical sizes obtained through US simulations, but our observations here and in the previous study^{23} indicate that this aspect ratio remains relatively consistent during nucleus growth. Also, the nuclei become more fragmented with increasing particle AR (and inversely more compact for smaller AR). These results also align with the conditions for kinetic optimality alluded to in the previous paragraph: in the longer GBF, the crystal nucleus interface would be more dominated by the crystal facet with the lowest interfacial tension, thus contributing to a lower Δ*G** (noting that, according to the classical nucleation theory, $\Delta G*\u221d\gamma 3$ and a tendency for nucleus fragmentation. Interestingly, our observations indicate that the critical nucleus size always have at least four layers of particles along the long axis of the nucleus; this is consistent with the facts that (i) the crystal nucleus shape tends to have an elongated aspect ratio (as suggested by the Wulff construction) and (ii) to form a minimal seed of the ABCD lattice one would need at least four layers stacked along the particle’s main axis.

While simulating the transition kinetics or Δ*G**’s associated with the isotropic–nematic or the nematic–crystal phase transitions for our most oblate (AR = 1/3) and most prolate (AR = 3) GBF shapes lies beyond the scope of this study, it is interesting to note that those transitions occurred spontaneously upon compression from the isotropic phase. This should be contrasted with the trend observed for increasing oblateness, e.g., the AR = 1 and AR = $1/3$ cases for which no intermediate nematic phase occurs, and the crystal state cannot be reached spontaneously from the isotropic phase (regardless of the degree of supersaturation). This supports the conjecture^{18} that the emergence of a mesophase (whose structural order is intermediate between those of the isotropic and crystal phases) acts as a transitional state that effectively facilitates (or catalyzes) the attainment of crystalline order by splitting a single large barrier (for the isotropic-to-crystal transition) into two smaller, easier to overcome barriers (for the isotropic-to-mesophase and the mesophase-to-crystal transitions). While increasing prolateness does reduce the isotropic-to-crystal Δ*G** as shown in Fig. 10(a), the crystal phase obtained spontaneously upon compression becomes more prone to forming multiple grains; in contrast, once a nematic phase occurs (e.g., for AR = 3), it templates the bulk orientational order and, hence, promotes the formation of single-grain crystals upon compression. Altogether, it is seen that the appearance of a nematic mesophase is helpful in realizing the crystalline order from the isotropic state. A more detailed account of these observations and the corresponding equations of state for different systems are given in the supplementary material.

## IV. CONCLUSIONS AND OUTLOOK

In this study, we have studied both thermodynamics and kinetic properties of the disorder-to-order transition in a family of hard GBF-shaped particles, taken as representative of faceted particles that exhibit a low order of rotational symmetry and form a non-trivial crystalline lattice structure. For the regular GBF, we found that the unusual anisotropy of crystal nuclei in a metastable isotropic phase can be largely explained by the anisotropy in interfacial energy along different crystal planes. Indeed, by directly measuring interfacial free energies of various GBF crystal planes with the isotropic phase, we are able to predict a prolate nucleus shape and average nuclei interfacial free energy in remarkable agreement with independent estimates from umbrella sampling and nucleus size pinning simulations. While it is not fully resolved why certain crystal planes have higher interfacial free energy, our interfacial PMF calculations point to the significance of a correlation between surface roughness of the crystal plane and its interfacial free energy. If this correlation holds more generally, one could then engineer particle shapes such that certain crystal planes are preferred at the nuclei’s interface or that prominent facets are flatter to enhance the transition kinetics.

The phase behavior of the regular GBF (AR = $3$) qualitatively extrapolates unchanged to the GBF-like shapes with not too dissimilar aspect ratios, but it significantly changes at extreme aspect ratios (e.g., for AR = 1/3 and 3) with the stabilization of liquid crystalline phases in between the isotropic and crystal states. This is in line with the results from other anisotropic hard particle shapes, and in line with a prior conjecture, we find that the emergence of the nematic phase is associated with the facilitation of spontaneous formation of the crystalline phase through gradual compression simulations that start from the isotropic phase. For the range of aspect ratios where the phase behavior comprises of isotropic and crystalline phases only, we find that the height of nucleation barriers (for a given degree of supersaturation) inversely correlates with the particle aspect ratio. Thus, one could conjecture that a similar trend may hold for other hard particle shapes, i.e., that the crystallization will be accelerated upon elongation of particles.

There are several open questions regarding the disorder to order transitions of anisotropic hard particles. While we noted that both the crystal nucleation free-energy barrier decreases and the nucleus shape elongates as the GBF becomes more prolate, observations suggestive of a reduction of nucleus’ interfacial free energy, it would be instructive to support such trends by computing interfacial free-energies (say at coexistence conditions) for systems with varying aspect ratios. While we have focused on a specific particle geometry (GBF), we expect that some of the trends in phase and kinetic behavior to be general; however, further studies are needed on other particle shapes and their kinetics as they may reveal case-specific peculiarities during a phase transition. Importantly, while we used nucleation free-energy barriers as surrogate metrics for the transition kinetics, comparing the exponentials of −Δ*G** (in *k*_{B}*T* units) only provides a sense of the relative nucleation rates across the different conditions and systems studied; methods that explicitly probe nucleation rates and time scales^{16,37,42,48} should also be implemented to provide data more directly testable by experiments.

It would also be interesting to explore to what extent other instances of anisotropic nuclei in hard particle crystallization can similarly be explained by a simple Wulff construction with a few crystal planes considered. Note that in our study, we ignored the effects of any dislocations or defects that may lead to an irregular presentation of a crystal interface that may not correspond to a crystal plane of the perfect crystal. Such deviations, often borne out of screw dislocations and defects, could be interesting to consider in future studies. Furthermore, it would be illuminating to find and analyze cases where the particle AR is negatively correlated with the nucleus aspect ratio and elucidate the role played by trends in interfacial tension. Finally, having noticed that the nuclei become increasingly fragmented with increasing aspect ratio, we wonder how elongated a nucleus could become before becoming dendritic. In our system, reaching such a scenario was precluded by a change in phase behavior, but other systems may not have this limitation.

While we noted some qualitative features of the nematic phases and their transitions for very oblate and prolate GBFs, it would be of interest to study in more detail the liquid crystal to crystal phase transitions for faceted particles. For example, some similarities and differences may exist in how the liquid crystalline phase crystallizes for prolate vs oblate particles. Such systems would provide an opportunity to quantify the overall reduction in the ordering transition barrier due to the presence of the mesophase, as has been posited for other soft matter systems.^{18,21} Work along these lines is currently underway. Some preliminary analysis of the ordering occurring for the AR = 3 GBF system is included in the supplementary material.

Experimental studies with faceted particles whose aspect ratios are tunable would allow testing of some of our predictions. Our finding that longer particles were easier to crystallize could mean that the crystallization pathways could be engineered to promote or hinder phase transitions. For example, it has recently been reported that GBF related lattices could be of special importance for achieving the colloidal diamond lattice,^{22} with the specific AR = 1 being particularly useful for this purpose. However, since AR = 1 GBFs are one of the harder systems to crystallize according to our simulations, and it could hence be useful to consider approaching the sought-after lattice structure by first crystallizing a longer shape and then shrinking it, e.g., by exploiting the anisotropic swelling–deswelling of elastomeric particles.^{49}

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional details provided pertaining to the calculations of equations of state, chemical potentials at the disorder–order phase transition, and umbrella sampling of crystal nucleation barriers for GBFs of different aspect ratios.

## ACKNOWLEDGMENTS

Funding support from NSF Award Nos. CBET-1907369 and CHE-2101829 is gratefully acknowledged. We would like to thank Professor Marjolein Djikstra (U. Utrecht) for useful discussions.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Abhishek K. Sharma**: Conceptualization (equal); Investigation (equal); Methodology (lead); Software (equal); Visualization (lead); Writing – original draft (equal). **Fernando A. Escobedo**: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.