Mapping Disorder in Entropically Ordered Crystals

TL;DR

This paper investigates how entropy is distributed in crystals of hard hexagons, revealing that translational degrees of freedom dominate entropy near key phase transitions, through a mode decomposition approach.

Contribution

It introduces a method to decompose thermal fluctuations into collective modes, extending Onsager's theory to quantify entropy contributions in hard shape crystals.

Findings

Translational entropy dominates near densest packing.

Rotational entropy is less significant in the studied regimes.

The approach generalizes Onsager's theory to complex crystal systems.

Abstract

Systems of hard shapes crystallize due to entropy. How is entropy distributed among translational and rotational microscopic contributions? We answer this question by decomposing thermal fluctuation of crystals of hard hexagons into collective modes, a generalization and quantification of the Onsager picture of hard rod liquid crystals. We show that at densities both near densest packing and near the solid-hexatic melting transition, solids of hard regular hexagons hold most of their entropy in translational degrees of freedom.

Figures

Equations

$$\frac{{\partial }^2{\mathcal{U}}_{1,2}}{\partial r_{1,2}^2}=A,\frac{{\partial }^2{\mathcal{U}}_{1,2}}{\partial ({\theta }_1+{\theta }_2{)}^2}=B_{+},\frac{{\partial }^2{\mathcal{U}}_{1,2}}{\partial ({\theta }_1-{\theta }_2{)}^2}=B_{-},$$

$$u_{\mu }(\vec{k})=\frac{1}{\sqrt{N}}\sum _{\alpha }(r_{\mu }^{\alpha }-⟨r_{\mu }^{\alpha }⟩)e^{-i\vec{k}\cdot ⟨{\vec{r}}^{\alpha }⟩},$$

$$\theta (\vec{k})=\frac{1}{\sqrt{N}}\sum _{\alpha }\left[({\theta }^{\alpha }-⟨{\theta }^{\alpha }⟩)\mathrm{mod}\frac{2\pi }{6}\right]e^{-i\vec{k}\cdot ⟨{\vec{r}}^{\alpha }⟩}.$$

$$S=-\frac{1}{2}\sum _{i,\vec{k}}\ln ({\omega }_i^2(\vec{k})).$$

Full text

Plato’s order/chaos dichotomy Archer-Hind et al. (1888) has a long intellectual history that achieves its most precise mathematical formulation in the theory of phase transitions Landau (1937). However, phase transitions can be driven by entropy Onsager (1949); Alder and Wainwright (1957); Wood and Jacobson (1957), which is historically associated with disorder Clausius (1864). Entropy-driven phase transitions confront Plato’s dichotomy with an apparent paradox: a macroscopically ordered system can exhibit more microscopic disorder than a macroscopically disordered system Frenkel (1999, 2015). Entropy-driven liquid crystalline Frenkel (1999) and crystalline order van Anders et al. (2014) can be rationalized through intuitive arguments for simple enough model systems, but resolving the paradox for more complex types of order Haji-Akbari et al. (2009); Damasceno et al. (2012); Agarwal and Escobedo (2011); Henzie et al. (2012) requires a detailed description of the entropy in an ordered system. Quantifying entropy contributions is complicated by the fact that entropy arises from statistical considerations, and is inherently macroscopic unlike energy, which trivially decomposes into a sum of microscopic contributions.

Here, we quantify contributions to the entropy of hard particle crystals by decomposing thermal fluctuations in two-dimensional systems of hard regular polygons into collective modes. Collective vibrational (phonon) modes have been studied previously for hard spheres Stillinger (1967); Honda (1976), and propagating rotational (libron) modes have been studied in detail in molecular crystals Schomaker and Trueblood (1968); Venkataraman and Sahni (1970); Raich (1971). We compute the dispersion relations for phonon and libron modes for hard hexagons Anderson et al. (2017), and show that the macroscopic system entropy is partitioned into phonon and libron contributions. From these contributions, we construct a “map,” by wavevector, of the contributions of microscopic order that lead to macroscopic order. For macroscopic order in systems at densities just above crystallization and close to maximal packing, we find a predominance of translational disorder. We find that the reason for the predominance of translational order, however, is different at low and high densities. Moreover, we find that between density extremes systems receive primary contributions from both translational and orientational disorder at different wavevectors. Our results give a concrete picture of “disorder within order,” and offer an intuitive explanation of recent, counterintuitive findings on the behavior of densely packed anisotropic colloids.Cersonsky et al. (2018)

Dispersion relations relate the frequency and wavelength of vibrational and librational modes in crystals. We computed dispersion relations for phonons and librons in a crystal of hard hexagons from NVT simulation using the hard-particle Monte Carlo (HPMC) plugin Anderson et al. (2016) for HOOMD-blue Anderson et al. (2008); Glaser et al. (2015). We use systems of $N=10,000$ hexagons near the hexatic-solid transition Anderson et al. (2017), and $2,500$ hexagons for higher densities. In addition, we validated our analytical and numerical approaches with systems of hard disks. In all cases, we initialized systems into hexagonal lattices, thermalized and equilibrated at a selected density, and sampled over $1\times 10^{7}$ Monte Carlo (MC) sweeps at high density, and $5\times 10^{7}$ MC sweeps at low density. Four independent samples were run at each density. From these MC data we computed dynamical matrix elements at all wavevectors $\vec{k}$ commensurate with the periodic boundary conditions using ensemble averages of statistically independent samples.

To calculate dynamical matrix elements, we constructed a linear elastic model (see SI for details) of the system that assumes that anisotropic hard core interactions can be modeled via effective, entropic pair interactions van Anders et al. (2014), $\mathcal{U}_{1,2}(\left|{\vec{r}_{1}-\vec{r}_{2}}\right|,\theta_{1},\theta_{2})$, which yield the parameters

[TABLE]

where $\vec{r}_{i}$ and $\theta_{i}$ describe particle positions and orientations. $A$ and $B_{\pm}$ determine the eigenvalues of the dynamical matrix (see Fig. 1 for illustration and SI for mathematical form) and thus function as effective, entropic stiffness parameters Bates and Frenkel (2000). We extracted them from dispersion relations computed with MC data using least squares fitting.

At high densities, when dislocations are absent, we use displacement covariance analysis (DCA) Hess (2000) to compute dispersion relations. DCA uses equipartition to relate small, thermal deviations $x_{i}$ to a dynamical matrix $K_{ij}$ according to $\langle x_{i}x_{j}\rangle=k_{B}TK_{ij}^{-1},$ where the $x_{i}$ of interest are phonon modes $u_{\mu}(\vec{k})$,

[TABLE]

and libron modes $\theta(\vec{k})$,

[TABLE]

Here, $r^{\alpha}_{\mu}$ is the $\mu$ spatial component of the position of the center of particle $\alpha$, $\theta^{\alpha}$ is the orientation of the body of the particle with respect to the positive $x$-axis, and $N$ is the total number of particles. Particles’ rotation symmetry is removed by the modulo $2\pi/6$ operation. We fix units by taking the particle circumdiameter and the thermal energy $k_{B}T$ to be unity. At lower densities, where dislocations become important and particles are able to diffuse from their lattice positions, we adapt a method that decomposes density fluctuations into contributions from transverse and longitudinal hydrodynamic wave modes Walz and Fuchs (2010). This method is tolerant of dislocations and does not assume particles vibrate closely around well-defined lattice positions, but is more computationally intensive than simple DCA. See SI for details.

We determined the density-dependent, effective elastic constants for both phonon and libron modes in hard hexagon crystals. We validated our results for $A$ by extracting the melting transition density, at which the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory Nelson and Halperin (1979) predicts $2A/\sqrt{3}$ should change discontinuously from the universal value of $16\pi$ to [math], and find good agreement with melting densities reported in Anderson et al. (2017) (see Fig. 2). Results for hard disks are also shown for comparison. For both hexagons and disks Fig. 2 data are given at densities $\phi$ expressed with respect to dense packing density $\phi_{d}$ (hexagons: 1, disks: $\pi/(2\sqrt{3})$) and melting density $\phi_{m}$ (hexagons: $0.710$, disks: $0.719$ Anderson et al. (2017)) in the form $(\phi-\phi_{m})/(\phi_{d}-\phi_{m})$. The near collapse of the data observed for the phonon elastic constant between hexagons and disks reveals the interesting result that using ordinary, translational sound, it is difficult at moderate densities and above to distinguish the difference between hexagon- and disk-shaped particles if density is scaled appropriately.

We present in Fig. 3 phonon and libron dispersion relations for a range of densities for high symmetry points in the Brillouin zone. These frequencies are computed from eigenvalues of the dynamical matrix multiplied by the inverse mass matrix, which has diagonal elements $1/m$ for phonon modes and $1/I=48/(5m)$ in reduced units for libron modes. We find several features of the dispersion relations that are common at all densities. First, we find deviations between our simulation data and the linear elastic model are only appreciable near the edge of the Brillouin zone. The edges correspond to short distances, and is precisely where we would expect a linear elastic model of hard interactions to break down. Although deviations increase with density as expected, we find that the linear elastic model provides a good approximation for most of the Brillouin zone, even for relatively high densities of 0.9. Second, we find that phonon modes are always ungapped, which is guaranteed by Goldstone’s theorem. Third, we find that the libron modes are gapped at all densities, but that their group velocity increases rapidly with system density.

What do these features imply the contributions are to the macroscopic entropy of the system? It can be shown (see SI) that the fluctuation spectrum yields a mode-by-mode decomposition of the configurational entropy (up to overall additive constants) according to

[TABLE]

Here, $\omega_{i}(\vec{k})$ are the dispersion frequencies, and the sum extends over all mode types $i$ corresponding to longitudinal and transverse phonon modes and libron modes at all wavevectors $\vec{k}$ in the first Brillouin zone. Brillouin zones are shown in Fig. 3, as “entropy maps,” in which all points are colored according to the mode that provides the dominant contribution to the system entropy at that wavevector, using Eq. 4.

This analysis produces a similar entropy map at densities just above the hexatic-solid transition (0.722) and close to dense packing (0.9) in which entropy at all wavevectors comes predominantly from translational disorder. The predominance of translational disorder just above the hexatic-solid transition arises because although the libron group velocity is small, the gap in the libron spectrum is sufficiently large to put the libron dispersion above the phonon dispersion throughout the Brillouin zone (save for wavevectors at the very edge of the Brillouin zone where we expect the harmonic approximation to breakdown). This is in contrast to what occurs at high densities. At high densities the phonon group velocity itself has increased by an order of magnitude (cf. vertical scales in Fig. 3a vs. d). In addition, we find there is an increase in the libron gap, but that this increase is much less dramatic than the increase in the libron group velocity, which becomes similar in scale to the phonon group velocity.

Our observation of ungapped phonon modes is a consequence of Goldstone’s theorem, but our observation of a gapped libron spectrum is the result of more subtle field-theoretic effects. For systems of isotropic particles, it has been shown that although crystals break both translational and rotational symmetry, the absence of Goldstone modes associated with rotational symmetry breaking is due to linear relations among translation and rotation generators that render the Goldstone modes associated with rotational symmetry breaking redundant Watanabe and Murayama (2013). For systems of anisotropic particles, rotation generators contain additional terms that correspond to particle body rotations which introduce additional degrees of freedom. These degrees of freedom can be shown (see SI) to always couple to the global rotation generator as a mass term in a linear elastic model. The existence of coupling between the generator of global rotations and the generator of body rotations implies that, in generic situations, the libron modes for anisotropic colloids or molecules should be gapped, whereas Goldstone’s theorem requires that phonon modes remain ungapped Brauner and Watanabe (2014). Thus, librational gaps should appear in any crystal or hexatic phase of anisotropic particles with global orientational order. Because of this, we expect that systems of anisotropic particles that melt from crystals with quasi-long range translational and long range orientational order into hexatic or other phases with quasi-long range orientational order but without translational order will exhibit entropy maps similar to those found here. We note that, practically, the identification of the signatures of solid-hexatic transition found here were far less computationally intensive to produce than approaches that were used to confirm the transition elsewhere Anderson et al. (2017). In contrast, if a system were to “melt” from a crystal with well-defined orientational order to a plastic crystal, we would expect a dramatically different entropy map.

At high densities, the predominance of translational disorder at all wavevectors arose through a combination of the gap in the libron spectrum and the comparable group velocity for phonon and libron modes. In our model systems we found that translational disorder predominated at all wavevectors at all densities at or above 0.9. Because the libron gap is a consequence of global orientational order, and because dense packing of highly anisotropic shapes generically leads to strong constraints on both particle positions and orientations, we would expect that similar computations in other systems would lead to similar phonon and libron group velocities. We therefore expect the predominance of translational disorder to be a generic feature of densely packed colloids. If this is so, it would help to explain recently reported, unexpected results from inverse-design van Anders et al. (2015) indicating that non-space-filling shapes are thermodynamically preferred in model systems of hard colloids even at densities above 0.99.Cersonsky et al. (2018) In the context of the present results, this finding can be rationalized because systems stand to gain more entropy through rotational modes – which the entropy map shows are restricted at all wavevectors at high densities – than through translational modes. This is most easily achieved by slight rounding of the particles into non-space-filling shapes, as observed in Ref.  Cersonsky et al. (2018).

We decomposed the total entropy of an entropically stabilized solid into microscopic contributions from collective modes of particle motion. We showed for a simple hard hexagon solid that translational and rotational degrees of freedom contribute different amounts of entropy at different wavevectors, extending the Onsager picture of hard rod liquid crystals Onsager (1949); Frenkel (1987, 1999) to hard particle crystals. Our method is applicable to any entropically stabilized solid in which the particles exhibit harmonic motion about equilibrium positions, and it is generalizable to three dimensions.