Polymorphism and Perfection in Crystallization of Hard Sphere Polymers

We present results on polymorphism and perfection, as observed in the spontaneous crystallization of freely jointed polymers of hard spheres, obtained in an unprecedentedly long Monte Carlo (MC) simulation on a system of 54 chains of 1000 monomers. Starting from a purely amorphous configuration, after an initial dominance of the hexagonal closed packed (HCP) polymorph and a transitory random hexagonal close packed (rHCP) morphology, the system crystallizes in a final, stable, face centered cubic (FCC) crystal of very high perfection. An analysis of chain conformational characteristics, of the spatial distribution of monomers and of the volume accessible to them shows that the phase transition is caused by an increase in translational entropy that is larger than the loss of conformational entropy of the chains in the crystal, compared to the amorphous state. In spite of the significant local re-arrangements, as reflected in the bending and torsion angle distributions, the average chain size remains unaltered during crystallization. Polymers in the crystal adopt ideal random walk statistics as their great length renders local conformational details, imposed by the geometry of the FCC crystal, irrelevant.


Introduction
The importance of crystallization in physics, materials, life science, and technology cannot be sufficiently emphasized. It plays a major role in the pharmaceutical industry [1], the storage of clean fuels such as hydrogen [2], and catalytic processes [3], among other industrial processes [4]. Furthermore, ordered structures of polymer-based physical and chemical systems are of paramount importance in the development of novel solar cells [5], semiconductors [6], biological materials [7], or conventional plastics [8]. Due to its high complexity, crystallization remains a topic surrounded frequently by scientific controversy and debate, especially since the classical view is considered to be too simple to cover a wide range of physical systems and conditions [9][10][11]. The basics of crystallization have been studied, among other systems, on simple spherical hard colloids, due to their suitability as macroscopic models, as they can be easily probed, and they bear properties and characteristics that can be tailored within wide ranges [12,13].
Extensive and careful experimental work has unmistakably demonstrated the difficulty of obtaining quite perfect crystals [14][15][16]. The identification of competing ordered structures of colloidal crystals has been carried out using methods such as fluorescence confocal scanning laser microscopy [17], small-angle synchrotron X-ray diffraction [18], laser scanning confocal microscopy [15], or light scattering [19]. In many of these experimental studies, the random hexagonal close packed (rHCP) structure prevails as the final ordered structure. It has also been shown that gravity and microgravity play a significant role in the crystallization of colloids [14,[20][21][22][23]. In many cases, the aging of the structures of colloidal spheres leads to a slow transition from the rHCP morphology to the theoretically expected, and thermodynamically more stable, face centered cubic (FCC) structure of varying degrees of perfection [24][25][26].
The available estimates of the entropy (directly proportional to free energy for hard spheres) demonstate an advantage of FCC against HCP that ranges between 9 × 10 −4 and 50 × 10 −4 per particle (expressed in terms of Boltzmann's constant, k) [38,[40][41][42]46]. The variation is attributable primarily to the methods used, and to a lesser extent, to the conditions under which it is calculated: ∆S seems to vary by about 25-30% between the melting transition and the maximum density. It is the smallness of this value that is responsible for the difficulty of obtaining neat crystals of the stable FCC polymorph, both experimentally and computationally. As a consequence, it is not surprising that the vast majority of isochoric simulations, starting from predominately amorphous monomeric HS packings, result in highly defective ordered structures of rHCP character [28,30,49,50]. The body of simulation works where FCC-like crystals are obtained is very limited [51,52]. The investigation of the FCC-HCP competition has been extended in simulations to the study of the effect of gravity [53], crystallization from seeds [36], sedimentation [54], templateassisted crystallization [55][56][57], or crystallization on surfaces and interfaces [58][59][60] .
Research studies addressing the crystallization of hard spheres forming linear sequences of chains are very sparse compared to the monomeric (single) HS analogs. Experiments on linear polymers of hard spheres are very challenging. In spite of the granular, colloidal, or droplet polymers being significantly less well explored than the "traditional" ones, over the years, there have been significant advances in their synthesis and characterization [61][62][63][64][65]. From the perspective of theory and simulation [66,67], emphasis is placed on the phase transition of semi-flexible [68][69][70] or flexible [71,72] chains of hard spheres, while recently, it was demonstrated that the use of block copolymers leads to HCP stable phases [73][74][75]. The stabilization of HCP colloidal structures has been investigated through the insertion of polymers [76,77]. In spite of these advances, crystal perfection and the relative stability of crystals of long, entangled HS polymers are still considered to be uncharted territory.
The very long relaxation times are responsible for the difficulty of preparing crystals of polymers of reasonable quality in the laboratory. It is no exaggeration to say that high-quality bulk polymer crystals with a well defined habit remain a laboratory curiosity, although a very valuable one, for they provide fundamental information on the unique characteristics of polymer crystallization. Most crystalline and semicrystalline polymers consist of huge assemblies of imperfect crystals, often and unavoidably combined with amorphous regions. Focusing on the computational work, the large size and concomitant sluggish dynamics of polymers present researchers with unusual challenges. Major conformational rearrangements involving the slowest modes, which play a key role in the formation of polymer crystals, can be well beyond the reach of deterministic methods.
Monte Carlo (MC) methods do not suffer from the slow dynamics associated with large molecular size, as in Molecular Dynamics (MD). The most advanced MC methods are precisely based on highly non-physical moves that allow for rapid equilibration and robust sampling in the configuration space [78,79]. This advantage is obviously offset by the loss of dynamic information, although methods such as Kinetic Monte Carlo (KMC) [80,81] do offer a reasonable compromise when the rates of individual events are known in advance.
Developing MC methods that correctly and efficiently sample polymer conformational space is not a trivial matter, especially for very long chains at very high concentrations. In spite of these difficulties, a wide variety of increasingly more efficient MC methods have been developed over the last few decades [78,79,[82][83][84][85][86][87][88]. The work to be reported in this manuscript is based on a powerful suite of advanced MC moves [89], which in the past has enabled us to observe the entropy-driven athermal polymer crystallization for the first time [71], and to identify and analyze the factors that affect the phenomenon, including chain length and its distribution [71,90], the presence of bond gaps or tangency [91], and confinement in one [92] or all three [93] dimensions.
In the present contribution, unprecedentedly long Monte Carlo simulations allow for the study of polymorphism and perfection in very long chains of hard spheres. These isochoric simulations start from an isotropic amorphous packing, and after a transient dominance of the HCP polymorph and the successive establishment of rHCP morphologies of various levels of ordering, they eventually reach a final FCC crystalline state of remarkable perfection. We gauge the established crystal structures and the fivefold local symmetry, and analyze the bond geometry and global sizes of the hard-sphere chains in each region of the phase transition. The entropic origins of crystallization, strongly related to the structural rearrangements of the local environment around each site, are also described and quantified.
In a companion paper [94] we support the computationally observed stability of the FCC polymorph by means of quantitative analytic calculations.

Methodology
We adopt the freely jointed model of linear chains comprising hard sphere monomers with uniform diameter σ, which is taken as the unit length. The pair-wise energy, u HS (r ij ), is given by the equation: where r ij is the distance between the centers of monomers i and j. This is the only type of interaction; neither bending nor torsional potentials are applied to successive monomers along the chain backbone. For numerical convenience, bond lengths, b len , are allowed to vary uniformly in the interval b len ∈ [σ, σ + db], where db is the maximum bond tolerance (gap) between two successive monomers, and it is set equal to 0.13σ. b len corresponds to the average bond length, where denote the average over all bonds for a given set of system configurations (frames).
The simulations are conducted using the Simu-D simulator-descriptor suite [89] in the isochoric, semigrand ensemble [VTN sites µ * ]. In the [VTN sites µ * ] ensemble [95] underlying the present calculation, chain lengths obey a given distribution, which is enforced via means of the chemical potentials µ * , and are allowed to fluctuate within a predetermined range. The practical implementation for the uniform and Flory chain length distributions is explained in detail in the Appendix of Ref. [96]. For athermal systems, the formal conjugate pair of the variables T, U is inactive (see [95] for details). The system under study comprises N = 54 chains of average size l av = 1000, leading to a total of N sites = 54,000 monomers or sites. Here, a flat (uniform) chain length distribution is chosen within the interval [l min , l max ] with l min = 600 and l max = 1400, as a requirement for the application of specific algorithms (see below). As will be shown in the results section, this chain length range lies deep in the polymeric regime. The uniform distribution is selected over the Flory one, as it allows for the robust sampling of the long-range chain characteristics in the whole interval [l min , l max ].
The initial configuration is generated through the progressive shrinkage of a very dilute configuration that meets the constraints imposed by chain connectivity until a desired packing density (volume fraction) of ϕ = 0.56 is reached [89], or ϕ * = ϕ ϕ FCC = 0.756 relative to the maximum compacity of the FCC or HCP crystals (ϕ FCC = ϕ HCP = 0.7405). The cubic simulation cell is subjected to periodic boundary conditions in all directions. Due to the high density, a configurational MC bias pattern is used throughout, with n dis = 50 candidate configurations being attempted per local move, as explained in detail in [89,99]. The isochoric simulation is carried out over 1.4 × 10 12 MC steps, and system configurations are recorded every 10 8 MC steps. Due to the very large system size, the time required to establish crystal perfection was close to 4 years of continuous wall-clock CPU time on a single Intel i5 processor with 32 Gb of memory.
The core post-processing tool of the work is the Characteristic Crystallographic Element (CCE norm), as implemented in the descriptor part of the Simu-D software [89]. This descriptor of structural order gauges both the radial and orientational similarity of the local environment of a monomer with respect to the reference crystals in two or three dimensions [101,102]. To this end, a tessellation into Voronoi cells (VC) around all monomers of the system is carried out (by means of the voro++ software [103]), and as a first step, the closest neighbors are identified for each monomer. The CCE norm is then computed for each monomer by carrying out the geometric symmetry operations of the point group of several candidate reference crystals (e.g., FCC and HCP) and evaluating the deviation from geometric invariance. The CCE norm, ε X j is assigned to each monomer, j, which is a measure of the deviation of its environment with respect to each reference crystal, X. The closer the value of the norm to zero, the greater the similarity of the monomer's environment to the reference crystal. A type is assigned to each monomer, depending on its structural similarity to a given characteristic point group (crystalline or not), by checking its CCE-norm against a threshold value, set here at ε thres = 0.3. This value was found to ensure selective discrimination among competing polymorphs for the system at hand. We calculate this CCE norm for each monomer at each MC frame, independently of which chain they belong to, in order to quantitatively analyze the evolution of crystallinity and the competition among polymorphs [101,102].
System configurations are tested against all reference 3D crystals implemented in Simu-D: hexagonal close packed (HCP), face centered cubic (FCC), body centered cubic (BCC), simple hexagonal (HEX), but also with respect to the non-crystallographic fivefold (FIV) local symmetry. The salient differences in the crystallographic operations and point groups, as well as in the shape and size of the corresponding Voronoi polyhedra for each reference crystal, are analyzed in detail in Ref. [102]. Only four different types of monomers are detected along the entire simulation: • Monomers with first neighbors whose arrangement conforms to the point symmetry group m3m of the face centered cubic crystal (FCC sites), • Monomers with first neighbors whose positions conform to the point group6m2 of the hexagonal close packed crystal (HCP sites), • Monomers whose first neighbors conform to the point group 5 of a non-crystallographic fivefold axis (FIV sites), • Monomers with first neighbors that are not arranged according to any kind of symmetry (apart from the identity E), either crystallographic or non-crystallographic (amorphous, denoted AMO).
In all of the analyzed sites and over all of the recorded configurations (frames), no instances of sites with BCC or HEX symmetry are found. Based on the CCE-norm, an order parameter, S X ∈ [0, 1] for each reference point group type X ∈ [HCP, FCC, FIV] is defined as: where P(ε X ) is the probability distribution function for the CCE norm of the X reference point group symmetry. In addition, given that only FCC and HCP crystalline sites exist in the system configurations S AMO = 1 − S HCP − S FCC − S FIV . As neither AMO nor FIV sites contribute, the total crystallinity is given simply by:

Results
As stated earlier, the simulation reaches 1.4 × 10 12 MC steps, with one configuration for every 10 8 MC steps being stored, for a total of 14,000 configurations (or frames or snapshots). This separation between successive frames was found to be sufficient to ensure frame decorrelation. The number of 10 8 MC steps between successive frames corresponds to an average of 1.85 × 10 6 MC moves per chain between successive frames.

Evolution of Crystallinity
Starting from a purely amorphous, statistically homogeneous configuration, the system evolves through intermediate states until a stable polymorph of remarkable perfection is formed. The evolution of the individual order parameters and total crystalinity as a function of the MC steps is shown in Figures 1 and 2. The intermediate states, which appear before the MC calculation settles in an equilibrated FCC polymorph, are characterized by a spatial distribution of regions of varied crystallographic nature, which evolve both in size and perfection along the simulation. Figure 1 covers the whole trajectory while Figure 2 focuses on the early part. The fraction of amorphous sites, S AMO , quantifies the degree of disorder in the system, to be contrasted against the degree of crystallinity, τ c . The system evolution can be conveniently split in four qualitatively different regions (Figure 1), numbered I through IV. Region I: MC steps 1 × 10 8 -1.5 × 10 11 (or equivalently frames 1-1500), II: MC steps 1.5 × 10 11 -5.3 × 10 11 (frames 1500-5300), III: 5.3 × 10 11 -9.1 × 10 11 (frames 5300-9100), IV: 9.1 × 10 11 -1.4 × 10 12 (frames 9100-14,000). From now on, except if otherwise stated, the frame numbers are quoted as rounded to the nearest 100, and are understood to be approximate, just as the boundaries between regions are not strictly sharp. Throughout the manuscript we will use kinetic terms such as "rate", "fast", etc., for brevity and simplicity, in order to describe the evolution of the system measured in terms of MC steps, but without claiming any truly kinetic or dynamic meaning.
In Region I, starting from the initial amorphous configuration, a rapidly growing number of monomers spontaneously develop crystalline characters, as evidenced by the steep drop in the AMO curve and the simultaneous growth of the number of FCC and HCP sites. As a matter of fact, 8.4% of all sites already have definite HCP or FCC character, and 11% show FIV character in the very first frame (10 8 MC steps). This is the reason for the curves in Figures 1 and 2 apparently not starting at 0 (HCP, FCC, and FIV) or 1 (AMO), and is a consequence of the great length of the calculation: alone in the first frame, 1.85 × 10 3 MC moves have been carried out, on average, for each of the N sites = 54,000 monomers. These observations are also qualitatively consistent with past simulations of dense random packings of monomers [49,50] and polymers [91,104,105], where the population of fivefold sites can exceed the sum of sites with HCP or FCC character.
After growing to approximately ≈10% within the first frame, the noncrystallographic FIV population remains largest in the first few frames, and then it drops steadily in parallel with the population of AMO sites, in favor of the FCC and HCP sites. The initial nucleation sites are transformed into compact assemblies of four sites of very approximate tetrahedral shape. The appearance of FIV sites is favored by the ease with which these four-monomer tetrahedra can be arranged in a variety of almost compact clusters that fill space efficiently at a small scale, like the pentagonal bipyramid shown on the right panel of Figure 3. The reason for the metastability of such arrangements is primarily geometric: at the smallest scale, five four-site tetrahedra can be arranged around a common axis so that they share an edge and build a seven-monomer bipyramid of almost, but not quite, exact fivefold symmetry, the fivefold axis being the common edge. By the ordered accretion of further elemental tetrahedra, structures of pentagonal symmetry and varying complexity [106][107][108][109] appear as metastable, long-lived morphologies in simulations. The left panel of Figure 3 shows the first neighbor shell of a representative site of FIV character (ε FIV 14717 = 0.05), as is spontaneously formed in Region I. Two of the pentagonal bipyramids described above share one site (marked in red in Figure 3) and also have a common, single fivefold axis, running vertically on the paper in this figure. This 13-site structure and other similar arrangements [110] are favored with respect to local disorder, because they lead to a localized increase in available volume and thus in translational entropy.   Initially, FIV sites appear abundantly and are randomly distributed in the simulation cell, but as their population subsequently decreases, they develop a remarkably nonhomogeneous distribution in space. The left panel of Figure 4 shows this spatial distribution of sites with FIV symmetry at a configuration within Region I (frame 125). At this stage, most of the numerous FIV sites (S FIV = 0.129) are still isolated, although the first traces of the linear assemblies of FIV sites are already visible. At a later stage of the simulation (frame 1622, early in Region II), as can be seen in the right panel of Figure 4, the FIV population is significantly lower (S FIV = 0.014) and the majority is organized in specific geometric patterns that are entirely analogous to those observed in the past [108,109]. Pentatwin formation is especially favorable due to the similarity between the pentagonal 72 • angle and the value of the dihedral angle 70.53 • between two faces of the regular tetrahedron. A gap of 360 • − 5 × 70.53 • = 7.35 • is left when five equally sized tetrahedra share an edge. Such tetrahedral nuclei appear easily via site accretion around the initial, minimal, four-monomer tetrahedra of Region I, and five such clusters form an almost perfect pentagonal dipyramid. Interestingly, the small angular mismatch of 7.35 • is not taken up by non-crystalline monomers as five equal, wedge-like gaps ≈ 1 5 7.35 • placed among the five twins, but the whole mismatch remains confined between two adjacent tetrahedra that are less perfect than the remaining three (see the rightmost two sectors in the middle panel of Figure 5). Because of the impossibility of compactly tesselating the 3D space using regular tetrahedra [112], noncrystallographic FIV sites are a transient phenomenon and tend to disappear in favor of compact structures of crystallographic symmetry. The simulation does indeed evolve the system into configurations of higher crystallographic FCC or HCP character. Consequently, a marked decrease in the number of FIV sites is observed in Region I, while total crystallinity increases rapidly. A video showing the FIV sites' evolution in Regions I and II (first 5000 frames) can be found in the Supplementary Material.  [34,[38][39][40][41][42], i.e., there is no initial overwhelming preference for the formation of the one over the other. As a consequence, and as long as a sufficient pool of amorphous sites is available, both polymorphs have roughly similar probabilities of being formed. The very large size of the simulation cell offers abundant volume for the independent growth of the nuclei of both types of crystals, i.e., there is plenty of available space for both the FCC and HCP crystalline regions to grow, as described above, without direct mutual competition.
In Region II, the rapid decrease in the number of AMO and FIV, and the acompanying growth in the crystalline FCC and HCP sites slows down appreciably, especially when compared to the trends observed in Region I. In this region, the numbers of FCC, FIV, and especially HCP sites remain much more constant, while the total crystallinity already exceeds ≈ 60%. Region II is characterized by the competition between the HCP and FCC crystalline regions. Region II can be seen as an induction period, and is, in a sense, a crucial phase and a stringent test of the efficiency of the MC scheme to properly and efficiently evolve the system towards the stable FCC polymorph, which is separated from HCP and the various rHCP structures by tiny free energy differences. The small fraction (O(10 −2 )) of FIV sites remaining in Region II have a different character from those in Region I, and they are not homogeneously distributed in space any more. After the rapid disappearance of the initial isolated FIV sites of Region I, the remaining FIV sites associate in almost perfectly linear clusters of ≈O(10) sites. These linear assemblies of FIV sites are the axes of cyclically twinned crystalline domains (see the schematic in Figure 6). In the ideal case, the pentatwins consist of five tetrahedral sectors, each made up of between a few dozens and a hundred monomers. Five such tetrahedra assemble in cyclic twin structures, the axis of which consists of sites with FIV symmetry, while the sectors can have mixed FCC and/or HCP character. Figure 6. In a cyclic pentatwin (pentagonal bipyramid), chains that span the boundary between twin sectors (A) or that pass through the twin axis (B) lose configurational entropy with respect to those in a bulk, untwinned crystal, as explained in detail in [113]. The white (empty) horizontal wedge is the angular gap of 7.35 • left when five equally sized tetrahedra are arranged cyclically around the twin axis (perpendicular to the paper in this figure).
In Region II, the number of HCP sites reaches a plateau, while the number of FCC monomers grows slowly and preferentially at the expense of both AMO and FIV sites. The total crystallinity in Region II is already appreciable, and most HCP and FCC clusters are contiguous, so that a conversion of crystalline HCP sites into FCC ones would be geometrically possible. In spite of this, it turns out to be easier to generate new FCC sites from the pool of available AMO sites than to convert neighboring HCP sites into FCC ones. This is consistent with the entropy difference between the amorphous state and both FCC and HCP crystals being larger than the entropy difference between FCC and HCP [94].
The difficulty of the FCC↔HCP interconversion is also a consequence of the appreciable metastability of the pentatwins: so long as they exist, the composition (HCP and FCC) of their sectors remains basically unchanged, so that HCP sites residing in pentatwin sectors are as immune to transformation into FCC as the twin itself [72]. This is the reason for why the proportion of HCP sites does not vary appreciably, as long as pentatwins exist.
There is also an extra factor, absent in monomeric systems and specific to chain molecules, that renders most mechanisms of FCC↔HCP interconversion for chains more difficult than for monomeric systems. At first sight, the most straightforward mechanism for HCP→FCC conversion is the rearrangement of the stacking . . . A-B-A-B-A-B-. . . →. . . A-B-C-A-B-C-. . . , which would allow FCC crystalline domains to grow at the expense of contiguous HCP domains. This rearrangement is, however, totally suppressed in polymer systems, because it entails a lateral slide of the (0001) HCP planes along the 12 · 0 directions. This slide, which is energetically expensive but allowed in monomeric systems [39], is not possible, even in principle in a polymer system, without breaking the chain backbones.
Easy pentatwin formation, metastability and slide suppression are the main reasons for the relative long life of the HCP polymorph in Region II. In spite of the efficiency of the MC scheme, it takes a considerable number of steps in the induction regime (plateau Region II in Figure 1), which correspond to months in wall-clock time, to initiate the final HCP→FCC conversion. The axes of these pentatwins appear as linear assemblies of sites of FIV symmetry, such as those in the right panel of Figure 4, which show the spatial distribution of sites with FIV symmetry at frame 1622 in Region II. These compact, non-crystallographic aggregates of tetrahedra also show up conspicuously as long-lived metastable states, even for monomeric systems. They represent intermediate states of appreciable stability. Their stability stems from the fact that the tetrahedra that compose the pentatwin can all grow simultaneously and independently of each other by accreting further sites from the periphery of the pentatwin. For every new surface layer accreted on the pentatwin, its axis, i.e., the linear assembly of FIV sites, grows by one site, and the number of sites in the pentatwin by 5n FIV (n FIV − 1) + 1 sites, where n FIV is the number of FIV sites along the twin axis, without necessitating any change in the existing structure. As a consequence of their easy formation and growth, ordered, noncrystallographic, metastable aggregates of FIV sites are found not to disappear, even in very long MD simulations of monomeric or polymeric packings [49,50,104,105].
The length of the induction period (Region II) is primarily determined by the slow disappearance of metastable structures having linear assemblies of FIV sites as their axes. The disappearance of metastable pentatwins is dictated by the impossibility to indefinitely grow into a periodic crystalline structure. For completeness' sake: the exotic possibility of the appearance of quasicrystalline structures has never been observed in our simulations of HS packings. In the present work, the difficulty of evolving the system out of this metastable situation, i.e., of spontaneously resolving the HCP-FCC competition, is marked by the slowdown in the evolution of all types of sites, i.e., the relative flatness of all curves in Region II, whose duration, measured in wallclock time, was about two years. The disappearance of FIV sites signals the end of Region II. It is only when the last linear FIV assemblies disappear that crystals resume their growth. The HCP-FCC competition is by then resolved, and the HCP regions start transforming into the stable FCC polymorph.
The mechanism of the disappearance of linear FIV assemblies in crystals of polymers is not just due to fluctuations brought about by the localized rearrangements of sites, as is the case in systems of monomeric hard spheres. In the case of chain molecules, the cyclic pentatwin incurs an extra entropic penalty with respect to the untwinned HCP or FCC crystals, because its sites are not individual spheres, but monomers of a polymeric chain. There exists an entropic penalty for chains in a cyclic pentatwin caused by the reduced number of available chain conformations with respect to those in a untwinned crystal. In the pentatwin, some chains are forced to span the (111) boundaries between sectors (marked with "A" in Figure 6), or to occupy sites in the [110] twinning axis (the linear FIV assembly, perpendicular to the plane of the paper marked with "B" in Figure 6). These two geometric situations force some of the torsional angles along the chain to adopt highly improbable values, which entails a decrease in entropy.
In [113], we evaluated the entropic difference between a twinned and a crystalline, but untwinned (bulk) system of hard-sphere chains, and found it to quantitatively explain the almost complete absence of FIV sites in properly equilibrated crystals of hard-sphere polymeric chains. The magnitudes of these entropic penalties in a pentatwin have been estimated to be −0.002k per monomer for chains spanning (111) sector boundaries, and −0.074k per monomer for chains that pass through the [110] twin axis. The former value is comparable to the entropy difference between the HCP and FCC polymorphs for equal monomeric spheres [34,[38][39][40][41][42], while the latter is ≈35 times larger. These two entropic penalties are thus sufficient to destabilize cyclic twins for chains, in agreement with previous simulation results on hard-sphere chains [71,90,91,114].
Once the induction phase is completed, the numbers of FIV and HCP sites in Region III both drop to zero, while the number of FCC monomers grows monotonically, up to a maximum of ≈94%. At the end of Region III, the populations of all site types reach an equilibrium state, which signals the onset of the production phase (Region IV), in which all of the properties fluctuate about well defined average values. Each of the frames in Region IV is a microstate of the macroscopic, crystalline, stable polymorph, generated with the correct probability from the [VTN sites µ * ] ensemble. The difference in entropy (and thus, the stability) between the stable polymorph and all other macrostates comes exclusively from the number of microstates that correspond to each macrostate, and not from differences in energy between microstates. Although also present in thermal systems, this entropic aspect is especially clear in the simulation of athermal ones: the entire evolution and selection of a given polymorph is not related in any way to differences in internal energy U (which can only be 0 or ∞, for feasible or infeasible configurations, respectively), but solely to the number of feasible configurations. The stable polymorph turns out to be the one with the (overwhelmingly) larger number of microstates.

Chain Statistics and Conformations
An inspection of system frames shows that as the simulation progresses, monomers tend to occupy (on average) the sites of crystalline clusters. This requirement has an impact on the bending and torsional angles along the backbone (see a sketch of Figure A1 for the definition of bending and torsion angles), for chains with unrestricted torsional and bending angles cannot efficiently occupy the sites of a crystal.
The distribution of bending angles in the crystal displays a good deal more structure than in the initial, amorphous state. Individual peaks in the distribution in Region IV can easily be assigned to angles formed by two successive bonds joining three sites in the ideal FCC crystal. In particular, the peaks at 60 • and 120 • are due to parts of the chains laying on sites of the 2D compact hexagonal layers from which the FCC crystal is built up. The broad peaks in Figure 7 that correspond to specially simple arrangements of three consecutive monomers which occur with high probability, are shown as sketches in the same figure. The arrows in Figure 7 mark the values of the bending angles for these simple arrangements. The bending angle values ≈ 30 • , 75 • , and 105 • are strongly suppressed in the crystal because they do not match the geometry of three consecutive sites along a chain on the FCC crystal. The angles of θ = 33.5 and 70.5 • have been found to play an important role in the crystallization of bend-core trimers interacting with the Lennard-Jones potential [115]. These angles are also unique for the HCP crystal, compared to the angles of θ = 0, 60, 90, and 120 • , which are common to both the HCP and FCC crystals. In our simulations at high densities, these bending conformations are improbable also in the initial amorphous phase, as they correspond to obtuse angles that do not minimize the local volume for triplets of monomers along the chain backbone. As crystallization proceeds, their population diminishes as a combined effect of the hard sphere model for the monomer interactions and the freely jointed model for chain connectivity. Values of θ > 120 • are forbidden, due to an overlap of monomers 1-3. Three-site arrangements that correspond to specific bending angles are sketched as indicated by the arrows. The following color and style format is used throughout the manuscript: Region I (dashed black); Region II (dotted green); Region III (dashed-dotted blue); Region IV (solid red).
As Figure 8 shows, the distribution of torsional angles (φ defined by monomers 1-2-3-4 in Figure A1) in the crystal also strongly deviates from the flat distribution of individual, isolated chains, and also from the distribution of the initial, amorphous configuration. Wide, bimodal peaks result from the overlap of closely spaced, favored torsion angles, like 54 • and 60 • . The peaks in the torsion angle distribution also correspond to specific conformations of four consecutive monomers (three bonds) whose positions correspond to neighboring sites in an FCC crystal. The torsion angles compatible with an ideal FCC crystal correspond to the values of φ = 0, 54.7, 70.5, 90, 109.5, 125.3, and 180 • . Some of these four-monomer sequences on the sites of a perfect FCC crystal are sketched along the torsion angle distribution in Figure 8. The virtually perfect symmetry about φ = 0 • of the distribution in Figure 8 is further proof of the robust configurational sampling due to the efficiency of the MC protocol and the great length of the simulation. It is further evident that for both the bending and torsion angle distributions, the minima and maxima become progressively more pronounced as crystal perfection advances, which is consistent with the observation that monomers tend to occupy, within fluctuations, the sites of an FCC crystal. A key question that arises is whether chains with torsional and bending angles distributions, that so markedly deviate from those in the amorphous state, also have unusual overall conformational properties. Ideal chain behavior in the crystal may at first sight seem counterintuitive and incompatible with the fact that monomers occupy the sites of a crystal and with the far from uniform distributions of bending and torsion angles. On the one hand, the great length of the chains renders local conformational details imposed by the geometry of the FCC crystal irrelevant. Still, there is another reason for the chains to display ideal behavior. As will be qualitatively discussed in the next Section (and quantitatively in a companion paper [94]), the stable macrostate for the polymer system is a crystal because it maximizes the sum of monomer translational entropy, because they occupy the sites of a crystal, and chain conformational entropy. This is because chains adopt ideal conformations compatible with monomers occupying the sites of the crystal. Thus, the crystal results from a combination of maximal positional freedom for the individual monomers about the sites of an ideal crystal, and maximal conformational variability for the entire chains. This result is independent of chain length, even for chains of moderate length (l ≈ 10) [71]. The chains in the present study are far longer. Confirmation that all chains in the system are deep in the polymeric regime is also provided by a plot of Kuhn length b 0 in Figure A2. Figure 9 shows the probability distribution function of the end-to-end vector modulus |R ee | for chains of lengths l = 800 ± 50 and l = 1200 ± 50 (scatter symbols), together with the expected distributions (dashed lines) for the same chain lengths, given by 4π|R ee | 2 f (R ee ) with: being the distribution of the end-to-end vector, and σ = 1 has already been taken into account in (2). In addition, the dashed lines are curves of the same functional form (2), fitted to the distributions obtained from the simulation. Gaussian behavior is clearly observed for both chain lengths in the FCC crystal in accordance with the ideal random-walk behavior.

Entropic Origins of Crystallization
The space available to monomers, its homogeneity and its isotropy, can conveniently be analyzed by means of the Dirichlet tesselation of space in Voronoi polyhedra, based on the coordinates of the monomers, as described in the context of the CCE-norm.
The increase in the regularity of the volume available to the monomers as the crystal is formed can be observed in Figure 10. This figure shows the probability distribution of the local density (measured as the ratio of the volume of a monomer (V sph = π 6 ) to the volume of its Voronoi cell (V VC ) for Regions I through IV. Since the [VTN sites µ * ] ensemble is isochoric, the curves in Figure 10 have the same mean value, but their shapes and variances do change in the course of the simulation. The wider initial distribution reflects the greater abundance of monomers with either low or high local volume fraction (i.e., high or low volume of their Voronoi cell), while the narrower distribution in Region IV highlights the greater spatial homogeneity (in the sense of the distribution having smaller variance) of the volume of the Voronoi cells in the FCC state. The distribution in Region IV (crystal) follows a Gaussian distribution almost perfectly (dashed line in Figure 10), which is consistent with the fluctuations of the monomers about the sites of an FCC crystal [116,117]. For reference purposes, the local volume fractions of the liquid and solid phases of the monomeric HS fluid at equilibrium are 0.667 and 0.736 times the maximum close packed density (0.7405), i.e., ϕ HS l = 0.494 (ϕ * l HS = 0.667) and ϕ HS s = 0.545 (ϕ * s HS = 0.736), respectively [42,118] (marked with vertical lines in Figure 10).
A clearer insight into the mechanism of the phase transition is given by the evolution of the volume accessible (V ac ) to a monomer within its own Voronoi cell. A simple yet convenient definition of accessible volume is that of the region of the Voronoi cell within which the center of the spherical monomer can be placed while keeping its distance to the nearest face of the Voronoi polyhedron < 0.5. The accessible volume in each Voronoi cell is computed here via straightforward MC integration. Although much more refined methods exist, like the analytical treatment in [119], the present approach is sufficiently precise for the calculations. The shape of the accessible volume is a smaller version of the Voronoi polyhedron, but with slightly rounded vertices, and whose faces are parallel to the faces of the Voronoi polyhedron. The logarithm of its size gives, to first order [120][121][122], a quantitative measure of the translational entropy of the monomer, under the neglection of multibody correlation and conformational entropy; that is, ignoring the fact that it belongs to a polymeric chain. Figure 11 shows a typical Voronoi cell as a wireframe representation, and the accessible volume as the red solid shape within the Voronoi cell. Unlike the total system volume V, the total accessible volume V ac is not conserved and is very sensitively related to changes in translational entropy. Figure 12 shows the distribution of the accessible volume for all N sites = 54,000 monomers, in the initial amorphous state at the beginning of the MC calculation, early in Region I and in the last frame of Region IV, i.e., when the FCC polymorph has fully developed. The distribution of accessible volume becomes clearly narrower, and its mean value shifts to the right. The mean relative accessible volume V ac /V VC almost triples from 9.0 × 10 −5 to 2.5 × 10 −4 as the system crystallizes. The whole distribution shows this increase vividly (Figure 12), and is direct evidence that, by adopting a spatially more homogeneous and geometrically more isotropic arrangement in the crystal, monomers increase their translational freedom, and hence, the system entropy. This is the main driving force for the crystallization of the chains, just as in a system of monomeric hard spheres.  The entropic driving force for crystallization can also be gauged by the correlation between local density and crystal quality, as measured using the CCE norm. Low values of this norm ε FCC , i.e., high local similarity to the perfect FCC crystal, are correlated with high local density, as Figure 13 demonstrates. In Figure 13, data points with ε FCC values below 0.30 have quite a perfect FCC crystalline environment. This numerical value separates the monomers into two non-overlapping clusters, one for non-FCC sites, in which no correlation between density and crystal quality is apparent, and a much more numerous second one, in which the correlation between density and crystal quality is very marked. As expected, monomers with an almost perfect FCC environment have significantly higher local density. The bottom panel in Figure 13 documents the wide range of accessible volumes for FCC sites: the accessible volume V ac for sites of very similar and quite perfect FCC character (with an average ε FCC ≈ 0.2) spans a range of two and a half orders of magnitude. This distribution is of course not static: the accessible volume of a given monomer fluctuates along the simulation within the approximate range [2 × 10 −5 , 6 × 10 −3 ] while keeping an almost constant value of ε FCC , i.e., a quite perfect FCC character. This is in complete agreement with the entropic origin of the transition: even though local density (inversely related to Voronoi cell volume) increases upon crystallization, monomer translational entropy increases because the accessible volume increases on average. The broad distribution for Region I in Figure 10 (black line) implies that an appreciable number of sites have very small Voronoi cells and hence very little accessible volume (black line in Figure 12).
During crystallization, Voronoi polyhedra not only become more homogenous in size, (as evidenced by the narrower distribution in Figure 10 and the shifted distribution in Figure 12 for Region IV), but also more isotropic. We have monitored the change in Voronoi cell isotropy in two complementary ways. A first straightforward measure is given by the modulus |R| of the vector joining the position of a given monomer and the centroid of its Voronoi cell. As crystallization takes place, the monomer and the centroid of its Voronoi cell tend to be closer together on average (shift to the left in the mean of the distribution of |R| in Figure 14), but also the distribution of |R| becomes narrower.
A resolved view in the change of the Voronoi cells is further provided by the evolution of the descriptors quantifying the shape of the Voronoi polyhedra, as derived from the second-order gyration tensor (see definition and corresponding results in Appendix A.3). In a rather literal sense, the increase in the accessible volume within Voronoi polyhedra, accompanied by greater uniformity and isotropy, is the force driving the transition from the amorphous to the perfect FCC crystal state, in spite of the opposition due to a loss of chain conformational entropy.

Conclusions
We have investigated the stability of crystal polymorphs for freely jointed chains of hard spheres deep in the polymeric regime through unprecedentedly long simulations, based on Monte Carlo algorithms. These simulations shed light on the competition between crystal polymorphs and on crystal perfection. Structural characterization of the local environment around each site allowed for a precise identification of crystal morphologies, which range from the early, defect-ridden, cyclic twin structures of mixed HCP/FCC type to the final, highly perfect, stable FCC crystal. We find the stable polymorph macrostate to be highly degenerate: all realizations (microstates) of the stable FCC macrostate consist of chains whose monomers tend to occupy, within spatial fluctuations, the sites of the ideal FCC crystal, while maintaining the constraints of chain connectivity and bond length. Polymer chains nevertheless display ideal behaviors, as their great lengths render local conformational details imposed by the geometry of the FCC crystal irrelevant.
The present simulations allow us to identify the entropic origin of the phase transition: the loss of chain conformational entropy is more than compensated for by the increase in translational entropy as the accessible volume within the Voronoi polyhedra increases. We have also carried out an analysis of accessible volume that provides further insights into the changes in homogeneity and isotropy of Voronoi cells responsible for the increase in monomer translational entropy. This increase in translational entropy is still large enough to drive the phase transition, in spite of chain conformational loss.
In addition, the present MC results also provide a basis for the calculation [94] of the entropic advantage of the FCC with respect to the HCP polymorph. The proposed methodology is currently being extended to tackle semi-flexible athermal polymers in the bulk and under confinement, as well as in composites with colloidal nanoparticles in the form of spheres and cylinders of varied size. . The authors acknowledge support through projects "RTI2018-097338-B-I00" and "PID2021-127533NB-I00" of MINECO/FEDER (Ministerio de Economia y Competitividad, Fondo Europeo de Desarrollo Regional). M.H. acknowledges financial support through "Programa Propio UPM Santander" of UPM and Santander Bank. The authors gratefully acknowledge the Universidad Politécnica de Madrid for providing computing resources on the Magerit Supercomputer through projects "r553","r727", "s341", and "u242".

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:  Figure A1 schematically shows the definition of bending and torsion angles, as adopted in the present work. This equivalent phantom chain has the same contour length and the same squared end-toend distance as the original chain (with excluded volume through Equation (1)). In order to fulfill these two conditions, the equivalent chain has two degrees of freedom: the number of equivalent segments l and its bond length, also known as the Kuhn length, b 0 .
From these two conditions: where |R ee | 2 is obtained for all chain lengths. Figure A2 shows a plot of Kuhn length b 0 as a function of chain length in Regions I and IV. The independence of Kuhn length b 0 on chain length indicates that even the shortest (l = l min = 600) chains in the system behave as ideal and fully flexible chains, and that all chains in the system can be characterized by the same Kuhn length b 0 = 1.52 ± 0.05 (taking the bond length as a unit).
The very small value of the Kuhn length (≈1.5 times the bond length) indicates an extremely high degree of chain coiling. Interestingly, no particular deviation in b 0 is detected between the initial, amorphous state (early part of Region I) and the perfected FCC chain crystal (Region IV), even if the bonded geometry undergoes significant changes during the transition, as demonstrated in Figures 7 and 8. The only appreciable difference is the slightly larger variance in the distribution of b 0 (larger fluctuations of the dashed curve in Figure A2) for different chain lengths, as expected from the greater spatial heterogeneity of the system in the amorphous state. i.e., the number-averaged analog of the moment of inertia tensor of a set of points of unit mass placed at the vertices of the Voronoi polyhedron.
Several scalar shape measures can be extracted from G, such as the asphericity b, acylindricity c, and the relative shape anisotropy κ 2 coefficient, which are defined, respectively by: where λ 2 1 , λ 2 2 , and λ 2 3 are the three eigenvalues of the positive definite G, ordered such that λ 2 1 > λ 2 2 > λ 2 3 . Notice that while κ 2 is dimensionless, b and c have dimensions of length squared; they are usually made dimensionless through division by 1 2 G. The lower the values of b, c, and κ 2 , the more isotropic the shape of the Voronoi polyhedron. An axisymmetric point distribution has c = 0, while a spherically symmetric one has b = c = κ 2 = 0. The Voronoi polyhedron for the FCC crystal, the rhombic dodecahedron, and for the HCP crystal, the trapezo-rhombic dodecahedron, both also have b = c = κ 2 = 0.
In Figure A3, we show the evolution of the (dimensionless) b, c, and κ 2 along the entire MC simulation. All shape descriptors (A2) in Figure A3 decrease monotonically from the initial values as the simulation evolves (the relative drops are b ≈ 33%, c ≈ 17%, and κ 2 ≈ 43%), and remain constant (within fluctuations) in Region IV of the stable polymorph. It is remarkable how faithfully all shape descriptors track the evolution of the systems through the four regions and the final onset of crystallinity. As with accessible volume, the distributions of the shape descriptors ( Figure A4) also become narrower, indicating a greater homogeneity of Voronoi cells. Figure A3. Shape descriptors asphericity b, acylindricity c, and relative shape anisotropy κ 2 . The first two are made dimensionless by dividing by half the trace of G. Figure A4. Probability distributions of shape descriptors b, c, and κ 2 for Regions I and IV. The first two are made dimensionless by dividing by half the trace of G.