# Three-state structural heterogeneity in a model two dimensional fluid

## Abstract

Three structural populations with distinct average mobility are identified within an equilibrium two-dimensional Lennard-Jones fluid simulated via molecular dynamics at a constant temperature and varying density. Quantifying the structure of the immediate neighborhood of particles by a tessellation allows us to identify three distinct structural states by the shapes of the tessellation cells. Irrespective of dynamic particle exchange among these populations, each is observed to maintain their own thermodynamic and average dynamic properties across the liquid-solid transition. We expect these findings to be valuable for better understanding the structural basis of dynamical heterogeneity in complex fluids defined in terms of local mobility fluctuations.

It is now widely appreciated that localized fluctuations in mobility are characteristic features of strongly interacting dense fluids. (1); (2) Commonly referred as dynamical heterogeneity (DH), (3); (4) this phenomena is observed in virtually all condensed matter systems, (5) ranging from microscopic to macroscopic scales, such as colloidal dispersions, molecular and polymeric liquids, metallic alloys, plastics (thermal systems) and granular materials-like sands, powders, foams and pastes (athermal systems). While dynamically heterogeneity is usually described in terms of the local mobility fluctuations of the fluid molecules or particles, (6) it is natural to expect some sort of complementary structural fluctuations that somehow facilitate such fluctuations, but the precise nature of these fluctuations is still unresolved. (7) The present work seeks to better quantify structural heterogeneity in a model strongly interacting dense fluids in order to better understand the origin of mobility fluctuations in this type of fluid.

The key question is how correlations in spatial structures influence the local dynamics of condensed materials.(8); (9); (10); (11); (12) Over the years, various local fluid properties have been proposed for this purpose: free volume defined by local cavity in fluids, (13) local Debye-Waller factor as a measure of local stiffness, (14) local elastic moduli, locally favored packing structures, (15); (16) local thermal energy, (17) to name a few proposed forms of structural heterogeneity. Some proposed structural measures are derived from the knowledge of particle position alone, while others require additional information about interatomic interaction and/or vibrational properties. While each of these local structural properties has had some success in understanding the dynamics of fluids, they often are system specific and lack consistency when compared with other structural indicators. Sometimes, the physical foundation of such indicators remains unspecified, as in recent application of machine learning approaches to glass-forming liquids. (18); (19)

The notion of DH was initially introduceed (20) to help understand the dynamical response of glass-forming liquids. In such materials, the existence of particles belonging to clusters of finite size and lifetime implies the existence of spatio-temporal correlations of some sort within the system. Certain long-lived fluctuations are then required to disintegrate such structures to allow their flow through structural relaxation. This suggests the coexistence of at least two dynamic populations of fast and slow particles (21) in this class of dynamically heterogeneous fluids. Similar broad distributions of particle mobility also arise in many other condensed materials at equilibrium, especially, near a liquid-solid phase transition. Given the stark difference in the long-time dynamics of liquids and crystalline media, a critical system undergoing such transition can easily be conceived as an admixture of dynamic populations with different mobility which shows average slow dynamics and increasing correlation length due to structural ordering, characteristic of DH. The obvious structural change across the liquid to solid transition also makes such structures suitable for understanding the structure-dynamics relations in a controlled setting other than the glass-forming liquids.

To this end, we employ a radial angular distance (RAD) (22) based space tessellation method to characterize the structure of a two-dimensional (2D) model single component equilibrium fluid approaching a liquid-solid phase transition by increasing density at a fixed temperature. This newly introduced technique is argued to better capture the local coordination of particles within a inhomogeneous polydispers system than the Voronoi tessellation. A 2D case is considered for the sake of simplicity and easy visualization. The local structures, quantified in terms of the shape of first coordination shells, then reveal at least three dynamic population contributing to DH which are structurally, thermodynamically and dynamically distinct from each other.

Model System and Simulation Details – We use large-scale molecular dynamics simulation (23) to prepare a two-dimensional equilibrium system of about particles in a canonical constant number-area-temperature (NAT) ensemble. Particles of size and unit mass interact via a Lennard-Jones potential where is unit energy and is the distance between a pair of particles. Temperature (in unit) is fixed by a DPD (24) thermostat as implemented within LAMMPS (25) where is the Boltzmann constant set to unity. Measuring time in units of , numerical integration of the equations of motion is done in steps of with a velocity Verlet algorithm. Systems with different , densities are considered ranging from to within a square box of suitable size with periodic boundary conditions in all directions. After equilibrating each system over , particle trajectories of length are stored for analysis. Equilibration is confirmed by the total energy fluctuation which is of the order of over the observation period. Further details of the simulation and thermodynamic characterization of the system is presented elsewhere in detail.(26) For the choice of densities, the model system shows (27) liquid, crystalline and liquid-solid coexistence behavior at . Below, we present and discuss the main findings of our study.

Equilibrium Dynamic Heterogeneity – Following the simulated trajectories of particles, it is easy to quantify their relative mobility in terms of their displacement over certain time period. We define individual particle mobility as where is the position vector of -th particle at time . Although the choice of time period is somewhat arbitrary, it is long enough to observe overall diffusive behavior in our model system. Color-coding the instantaneous particle configuration for immediately shows (Fig.1) the presence of fast (red) and slow (blue) particles with mobilities differing by the orders of magnitude for all densities considered. We see that the number of highly mobile particles decreases as particle crowding occurs with the increase in density and more particles become immobile. Importantly, the spatial organization of particles with similar mobility become apparent from such visualization which reveal the characteristic string-like arrangement of the most mobile particles,(28); (29) a hallmark of DH. Now that we have confirmed the existance of DH in our equilibrium model system, we focus on quantifying the structural heterogeneity of the same fluid.

Geometrical Neighborhood as a Measure of Structural Dynamical Heterogeneity – Identification of the immediate neighborhood of individual particles is a natural first step to gather structural information about a disordered system. As such systems enjoy the full translational and rotational symmetry, obtaining complete knowledge of both radial and angular arrangement of particles around a central one is a difficult task and no unique method to do so has been prescribed so far. Two widely used methods in this regard are the radial distribution function (RDF) based method and Voronoi tessellation (VT) technique. (30); (31) RDF based methods suffer from the arbitrariness of predefined cutoff to identify the closest neighbors. Also, this method which involves ensemble averaging over angles and configurations, is only capable of providing average properties relating to fluctuations in the particle arrangement so that this method does not provide insight into the dynamics of a particular particle arrangement. The VT provides more detailed information of particles’ neighborhood by defining space-filling non-overlapping polygons around each of the particles. The area of such polygons in two dimensions and volumes in three dimensions are known to follow a well-defined distributions, (32) universal to a wide variety of systems such as glass-forming liquids,(33) granular materials, (34); (35) etc. However, it is now well appreciated that free volume defined by this local measure of local packing is inadequate for predicting local dynamics. (36); (37) It is possible that other measures of local structures based on the Voronoi construction, such as, number of faces or edges of the Voronoi cells, average shape of the cells, etc. might be more informative about local dynamics. There is the additional problem that the Voronoi based methods are also computationally expensive and often require further modifications for highly inhomogeneous materials such as aggregate forming systems, polymer liquids etc . The matter of computational efficiency motivated us to consider another tessellation method, radial angular displacement (RAD) (22) to quantify the geometrical neighborhood of finite-size particles in any arrangement. This tessellation is particularly good at quantifying changes in local coordination number which we intutively expect to correlate to changes in local fluid dynamics since this physical property should alter the local cohesive potential energy and local barrier for molecular transport.

The RAD method identifies a particle as nearest neighbor to particle if it is not blocked by any other particle and not further away from any other blocked particle. The blocking condition is ensured by the following geometric inequality: where is the three-body angle made by and at and is the solid angle subtended by particle at the center of particle . Making a small angle approximation for , the criterion for equal size particles reads , is the distance between particles and and similarly, . Starting from a list of possible neighbors arranged in an order of ascending distance from the central particle, this parameter free method can quickly identify the first neighbor shells of all particles in a considerably large system as ours.

The geometry of each tessellation cell is then quantified by its shape factor, (38) where are the circumference and area of the shell, respectively. We note that is the inverse square root of a well-known shape descriptor, circularity. (39) Variations of basically the same shape descriptor have already been used in previous morphometric studies (40); (41) under different names such as reduced area, (42) isoperimetric ratio, (43) etc. While is designed to be unity for a circle, the analytic values for a regular -gon can be obtained by . For example, perfect hexagon, pentagon and square has values , and , respectively. The distribution shows a trimodal distribution with fixed minima at and (Fig.2(a)) for all studied. We note that this trimodal feature resulting from the RAD algorithm used is in contrast with the VT method. computed using Voronoi neighborhoods is reported to show a bimodal to unimodal transition across the phase transition in model hard-disk fluid (38) granular media (45) and colloidal systems. (44) The existence of three structural populations, proposed in these works based on physical intuition, can now be clearly identified using our tessellation method. The significance of these peaks for thermodynamic and dynamical properties is considered below.

Definitions of Three Structural Populations – The prominent peak close to at high densities assures crystalline packing with expected hexagonal symmetry. With decreasing density, this peak shifts and broadens as distorted hexagons start to appear allowing more disorder in the system. Perfect five- and four-fold symmetries are not allowed in two dimensions which is consistent with the position of fixed minima in . Particles with and can then be statistically identified as distorted pentagons and squares with increasing number as decreases. This well-behaved trimodal distribution allows us to unambiguously identify three distinct geometric populations: (i) , (ii) and (iii) . Similar structural characterization is also possible based on the area of local neighborhood as discussed elsewhere. (26) The shape factor criterion employed in this work is more rigorous as the differentiation of populations is based on well-defined analytical values rather than relying on the numerical fitting of the local area distribution adapted in the other method. Nonetheless, the structural populations identified by either of the criteria are essentially the same.

Thermodynamics of Three Interconverting Structural States – The number density of each geometric population computed as a ratio of the number of particles within a population over the total number of particles in the system is shown in Fig.2(b). We view for population as being a kind of “order parameter” that exhibits a sigmoidal variation from a very low value at low and to unity at high . A complementary variation is observed for which has a finite value at low and decreases in a sigmoidal way to a very small positive value with increasing . Interestingly, for both and near which marks the onset of liquid-solid coexistence phase, as signaled by a prominent inflection point in , the specific heat at constant volume computed from the potential energy fluctuation as a function of . (Inset Fig.2(a)) We see that for is negligible at high , but we find a linear increase of with decreasing starting from , where we observe another sharp inflection point in . A significant relation between local geometry and thermodynamics of the system is indicated as we observe that all three populations show well-separated mean energy values E at all . (Fig.2(c)) We mention that the number of particles involved in computation of is not restricted to just the nearest neighbors but includes all particles within , a cut-off set by the potential. Similarly, the hydrostatic pressure computed from the trace of virial stress tensor of each population is different from each other for all (see Fig.2(d)). The pressure for different populations plotted together with the bulk pressure of the system (shown by dashed line in Fig.2(d)) indicates that the overall thermodynamics for high systems is dominated by population, particles with hexagon-like neighborhood while for low systems with mostly pentagonal neighborhood i.e. population, consistent with previous observations of . Particles within the population, having square-like neighborhoods represent locally high pressure regions as their average pressure is always larger than the bulk. The well-defined differences in thermodynamic properties among these populations of particles naturally lead to a consideration of the time evolution of these classes of particles.

Time Evolution of Structural State Populations –

We first attempt to visualize the local geometry as it evolves in time. First, we consider a set of particles chosen along a line cross-sectional cut of Fig.1 where the particles on this line are colored at each time step according to their population identity: blue (), green () and red (). A line segment of same color then represents the time period over which the shape of the neighborhood of an individual particle is maintained. These kymographs [REF.] are presented for three different : Fig.3(a) , (b) and (c) . At low density, the neighborhood changes its shape quite frequently. However, the population identity of adjacent particles seem to maintain some degree of correlation at least for and populations. The persistence of the blue lines increases with increasing pointing, that the particles maintains their neighborhood shape much longer than other populations. This is expected since the hexagonal shaped cells, characteristic of , emerge in the system with increasing . From this exercise, we can now quantify the average life time of the local geometry around a particle. We find that each population generally has different at a fixed density.(Fig.3(d)) For high (crystalline phase), characteristic for each population is well separated by values of different order of magnitudes. as a function follows the same behavior of which is expected. We note that small is indicative of a fast changing local environment than a particle with larger . This brings us back to the question of DH in the mobility sense. A calculation of the mobility computed separately for each population, correspondingly shows distinctly different values at each .(Fig.3(e)) Evidently, and correspond to the immobile and highly mobile particle populations across all , while maintains an intermediate mobility. Our geometric characterization of structure is then successful in identifying the key structural population across the phase transition.

Concluding Remarks – In summary, a scale- and parameter-free analysis for simulated particle configurations of a model two-dimensional equilibrium liquid is presented. We emphasize that both the identification of neighborhood and its shape analysis are crucial to recognize the key structural motifs with distinctly different thermodynamics and dynamical properties. Since both of the methods have already been applied to the three dimensional (3D) systems separately, their combination as in our case should easily be translated to 3D systems. In fact, as a purely geometric approach which uses only the positional information of the constituents of a system, this method should be readily applicable to practically any inhomogeneous arrangement of constituent particles, even those that have significant polydispersity in 3D. Being relatively computationally efficient, our method can be used for on-the-fly structural characterization of equilibrium systems and non-equilibrium fluids like granular media, colloidal suspensions, nanoparticle composites etc. We note that competition in three population systems have also been intensively studied in the context of population dynamics of living systems to explain the spatial organization and movement of interacting species.(46); (47) Although the structural populations revealed by our study are dynamical interchangeable and thus, different from those studied in population dynamics, constructing a model dynamical systems theory by taking inputs from our study can be instructive. We believe that our approach to defining structural heterogeneity in strongly interacting fluids provides a promising first step forward in development of a rigorous structure-dynamics relation that should be applicable to a broad class of condensed matter systems. However, further analysis of the model will be required to demonstrate the generality and practical utility of this relationship.

Acknowledgement TD acknowledges support under the Cooperative Research Agreement between the University of Maryland and the National Institute of Standards and Technology Center for Nanoscale Science and Technology, Award 70NANB10H193, through the University of Maryland.

Official work of the US Government- Not subject to copyright in United States.

### References

- L. Larini, A. Ottochian, C. de Michele and D. Leporini, Nat. Phys. 4, 42 (2008).
- G. M. Hocky, D. Coslovich, A. Ikeda and D. R. Reichman, Phys. Rev. Lett. bf 113, 157801 (2014).
- M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000).
- L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).
- L. Berthier, G. Biroli, J.-P. Bouchod, L. Cipellitti and W. van Saarloos, Dynamical Heterogenities in Glasses, Colloids and Granular Media (Oxford University Press, New York, 2011).
- B. A. Pazmiño Betancourt, P. A. Hanakata, F. W. Starr and J. F. Douglas, Proc. Nat. Acad. Sc. 112, 1418654 (2015).
- C. P. Royall and S. R. Williams, Phys. Rep. 560, 1-75 (2015).
- G. Adam and J. H. Gibbs, J. Chem. Phys. 43, 139 (1965).
- R. W. Hall and P.G. Wolynes, J. Chem. Phys. 86, 2943 (1987).
- P. G. Debenedetti and F. H. Stillinger, Nature 410 259 (2001).
- F. W. Starr, J. F. Douglas and S. Sastry, J. Chem. Phys. 138, 12A541 (2013).
- K. F. Freed, J. Chem. Phys. 141, 141102 (2014).
- F. Spaepen, Scr Mater 54, 363 (2006).
- A. Widmer-Cooper and P. Harowell, Phys. Rev. Lett. 96, 185701 (2006).
- D. Coslovich and G. Pastore, J. Chem. Phys. 127, 124504 (2007).
- C. P. Royall, S. R. Williams, T. Ohtsuka and H. Tanaka, Nat. Mater. 7, 556 (2008).
- J. Zylberg, E. Lerner, Y. Bar-Sinai and E. Bouchbinder, Proc. Nat. Acad. Sc. 114, 7289 (2017).
- E. D. Cubuk, S. S. Schoenholz, J. M. Rieser, B. D. Malone, J. Rotler, D. J. Durian, E. Kaxiras, A. J. Liu, Phys. Rev. Lett. 114 108001 (2015).
- S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, A. J. Liu, Nat. Phys 12 469 (2016).
- J. H. Hildebrand, Proc. Nat. Acad. Sc. 72 1970, (1975).
- H. Sillescu, J. Non-Cryst. Solids 243, 81 (1999).
- J. Higham and R. H. Henchman, J. Chem. Phys. 145 084108 (2016).
- D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, 1996).
- R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997).
- Freely available at www.lammps.sandia.gov
- T. Das and J. F. Douglas, (to be published).
- J. A. Barker, D. Henderson and F. F. Abraham, Physica 106A, 226 (1981).
- H. Zhang, P. Kalvapalle and J. F. Douglas, Soft Matter 6, 5944 (2010).
- B. A. Pazmiño Betancourt, J. F. Douglas and F. W. Starr, Soft Matter 9, 241 (2012).
- G. Voronoi and J. Reine, Angew. Math. 134, 198 (1908).
- Marina L. Gavrilova (Ed.), Generalized Voronoi Diagram: A Geometry-Based Approach to Computational Intelligence (Springer-Verlag Berlin Heidelberg, 2008).
- V. S. Kumar and V. Kumaran, J. Chem. Phys. 123, 114501 (2005).
- F. W. Starr, S. Sastry, J. F. Douglas and S. C. Glotzer, Phys. Rev. Lett. 89, 125501 (2002).
- T. Aste, T. D. Matteo, M. Saadatfar, T. J. Senden, M. Schröter and H. L. Swinney, Europhys. Lett. 79, 24003 (2007).
- S. Slotterback, M. Toiya, L. Goff, J. F. Douglas and W. Losert, Phys. Rev. Lett. 101, 258001 (2008).
- F. W. Starr, S. Sastry, J. F. Douglas, and S. C. Glotzer, Phys. Rev. Lett. 89, 125501(2002).
- A. Widmer-Cooper and P. Harrowell, J. Non-cryst. Solids 352 5098 (2006).
- F. Moučka and I. Nezbeda, Phys. Rev. Lett. 94, 040601 (2005).
- A. M. Bouwman, J. C. Bosma, P. Vonk, J. A. Wesselingh, H. W. Frijlink, Powder Technol. 146, 66 (2004).
- J. Wasen, R. Warren, Mater. Sci. Tech. 5, 222 (1989).
- R. McAfee, I. Nettleship, Acta Mater. 51, 4603 (2003).
- A. Hočevar, S. El Shawish and P. Ziherl, Eur. Phys. J. E 33, 369 (2010).
- V. Lucarini, Symmetry 1, 21 (2009).
- P. M. Reis, R. A. Ingale and M. D. Shattuck, Phys. Rev. Lett. 96, 258001 (2006).
- Z. Wang, A. M. Alsayed, A. G. Yodh and Y. Han, J. Chem. Phys. 132, 154501 (2010).
- P. P. Avelino, D. Bazeia, J. Menezes and B. F. Oliveira, Physics Letters A 378, 393 (2014).
- P. P. Avelino, D. Bazeia, L. Losano, J. Menezes and B. F. Oliveira, Physics Letters A 381, 1014 (2017).