Microscopic structure and dynamics of glass like domains in a low density binary mixture

We present results from extensive equilibrium molecular dynamic simulations of a binary Lennard Jones 80:20 mixture. Our study focuses on the changes in structural and dynamic properties of this system during its transition from a stable phase to an unstable phase as temperature reduces. This system, phase separate into a mixture of glass like domains and gas at lowest temperatures. We compute coarse-grained density(CGD) to identify density distribution in the system, which shows formation of very dense domains. Studies of bond order parameter shows a mixture of crystal structures in these domains at low temperatures. These results show that the structure of the dense domains formed at low density is akin to glassy domains that is found in studies of glass transition at high density, at low temperatures. Our studies on mean square displacement of particles at various state points shows pronounced sub-diffusive regime and plateau that marks the existence of dynamic arrest due to cage formation known in the studies of glassy and supercooled liquids.


Introduction
In simple model of liquids, such as Lennard-Jones (LJ) potential, show universal phases of gas, liquid and solid. It is a well studied fact that a binary mixture of LJ system that consists of particles differ in their radii and interaction strength, at state points of lower temperatures and high density, form a glassy or amorphous state [1,2,3]. A computationally well studied example of this system is 80:20 binary mixture used in the studies of Kob and Andersen which is extensively used for computational studies of glass transition [4,5,6,7,8,9,10]. One of the outstanding problem of the glass forming liquids is to understand nature and relaxation properties of metastable states formed at state points that is close to the glass transition temperature T g [11,12,7,13]. There are many computational studies of glass transition in the binary mixture which are done at high densities [14,15,16,17,18,19,20,21]; where, it is difficult for the structural rearrangement to occur at low temperatures. Then one of the interesting question to ask is, what is the role of density in the formation of glassy domains or how the cooperative rearrangements that results in glassy dynamics is effected by density of the system. To investigate the role played by the density of the system in the formation of glassy or crystalline phases, let us first look into phase changes of a single component LJ system at low densities.
In many regions of the phase diagram of LJ liquid, see point (d) in Fig.1, at low densities and low temperatures where system cannot exist as a single phase: the system exists as a mixture of domains of vapor and liquid phases. These state points are considered as unstable states in the phase diagram of the system. In between (d) and (h) in Fig.1, where the system phase separate into random shaped gas and liquid domains [22,23,24]. In state points where temperature is further lowered, see point (h) in Fig.1, the liquid domains freeze into crystallites and exist as a mixture of crystallites and dilute gas. The crystallites formed can be considered as local potential energy minimized structures. Unlike a single component LJ system which at high densities and low temperatures, freezes into a crystal state, the glass forming binary mixture at these state points forms droplets of glassy domains [25,24]. In non-equilibrium simulation studies of this binary mixture by Testard et al. [25,24], they showed the temperature dependency of phase separation. Their study in binary mixture showed gel-like structures exist at low densities and temperatures. It is now interesting to study the structure and dynamics of these dense domains For BLJ, d shows the gas+liquid regime, e and f represents binodal and spinodal respectively, g represents liquid, h is a mixture of gas and glass, i represents glass and j represents the glass transition line. Schematic phase diagram is based on Ref. [26,27,25].
in the equilibrium state of this binary mixture at low temperatures and low densities. As these systems have excess free volume [28], therefore, the relaxation of these dense domains may differ from the glassy phase formed at high densities. In addition, it is also very interesting to see what are the local structures formed in a binary mixture at low densities. The structural changes of this binary fluid, in the state points of phase separation, while temperature is lowered, result in the structural arrest. We have studied these structural changes using the radial distribution function g(r). A useful way of distinguishing dense domains is by using coarse grain density (CGD) defined in the studies of Testard et al. [24]. We have used this method to characterize the domains obtained in our simulations. Our work shows coexistence of dense amorphous domain as well as gas in the low density and low temperature state points of the binary mixture we have studied.
Another feature that characterizes clusters formed is their local structure. In an important work that characterize local structures of LJ liquid by Steinhardt et al. [29], they used bond order parameters to identify various possible ordered structures. Bond order parameter has been used in many studies of glass forming liquids to identify the local structures [30,31,32]. Studies of bond order parameter of 80:20 binary mixture at state points near glass transition at high densities shows coexisting collection of structures like hcp, fcc, bcc, etc. [33]. To identify the existence of local structures in the system we looked into the bond order parameter of nearest neighbors of each molecule at various state points we have studied which also show coexistence of the mixture of structures found at high density. To further identify the cage relaxation process in the low density binary mixture we have computed mean squared displacement which clearly shows the sub-diffusive regime that get prolonged as temperature reduces. At lowest temperatures we find a plateau region in the mean squared displacements computed, which is a characteristic property of a glass forming system that is found in many studies [4,34,35,36,37].
This paper is organized as follows: In next Sec. model and computation details are presented. In section III we discuss the results from the simulations. Summary and conclusions are presented in section IV.

Model and Computational Details
We simulated 80:20 binary mixture widely used in studies of glass transition. Interaction potential between particles of binary Lennard Jones system is given by [4,5,6].
where, α, β ∈ {A, B}, AA = 1.0, σ AA = 1.0, AB = 1.5, σ AB = 0.8, BB = 0.5, and σ BB = 0.88 [4,5,6]. The interaction potential was truncated at a cutoff distance of 2.5σ αβ . All the quantities presented in reduced units of AA interaction, the Lennard Jones well depth AA is the unit of energy, while the Lennard Jones diameter σ AA is the unit of length, time in units of t * = mσ 2 AA / AA and Boltzmann constant is set as k B = 1. The equation of motion are determined by integrating Verlet algorithm with time step of 0.01 [24].
To study the properties of this system, we prepare equilibrated homogeneous system at high temperature T = 3.0 and at desired density. Then the system is quenched to state points of desired temperature and equilibrated before the production runs. Equilibration time used in our simulation is time that required for domain growth to saturate this is determined from the average time dependent chord length [24] reaches a plateau; it is nearly 10 3 t * at lowest temperatures. State points that we investigated are at T = 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 3.0 and at density ρ = 0.6. Our system has 10 4 particles (8000 A and 2000 B ) and the total length of the simulation is 10 6 number of steps.

Results and discussion
We have used radial distribution function g(r) to identify the formation of secondary structures while temperature is reduced. The partial radial distribution function for the binary mixture is defined as [38] where, α ∈ {A, B}, r ij is position vectors of the particles. The radial distribution function of the particles of three distinct pairs, AA, BB, and AB are given by g AA (r), g BB (r), and g AB (r). These radial distribution functions gives how respective species of the particles are distributed in the system. In the simulation of the binary mixture at high densities, these radial distribution functions show splitting of the second peak due to the formation of multiple structures at state points of low temperature [4,39,40,31].
Figs. 2 shows radial distribution functions g AA (r), g AB (r) and g BB (r) at ρ = 0.6 and at different temperatures. For better representation, curves have been displaced vertically in Fig. 2. As state points of lower temperatures peaks become sharper in all three radial distribution functions presented. In addition, g AA (r), g AB (r) and g BB (r) show multiple splitting of second peaks at state points of low temperatures which is the signature of different local structures present in the system; this splitting of the secondary peaks similar to what found in the simulation studies of binary mixture at high densities [40,4]. In contrast results of simulation studies at high density system, in our simulations, splitting of second peak at low density shows multiple peaks, which indicates the existence of the various structures that differ from each other in this system. Another characteristic of g AA (r), g AB (r) and g BB (r) is , at low temperatures their mean go above value 1 which shows clustering of particles the system, which is the expected behavior of the system in unstable region given in Fig.1. It is now interesting to see the characteristics of these dense domains formed and the corresponding structures.
To identify the nature of high density domains we have used method of CGD proposed by Testard et al. [24]. This method can directly measure the average density of a domain and it also this method helps to identify various domains and their interfaces. The methodology can be summarized as follows. First discretize the simulation box into multiple cells of linear size ξ b , so a total of V /ξ 3 b sites [24]. Then for each cell we compute discrete density defined as is the Heaviside function, ξ s is the radius of sphere from center of each cell. The final CGD computed by considering the neighboring cells too and it is given bȳ (2) Then the probability distribution of CGD is computed.      The CGD distribution we get is similar to that obtained by Testard et al. [24] at state points T = 0.1 and ρ = 0.6. Fig. 3 shows CGD at ρ = 0.6 for various temperatures. At low temperatures CGD shows two peaks: first belong to regions of gas phase; second belongs to dense phase. Between these phases there are regions of interface which also appears in the CGD. At high temperatures, T = 3.0 and density ρ = 0.6, CGD shows a homogeneous distribution of densities in the cells that appears as a single peak. We see that at low temperatures, the system prefer to exist in a dense domains whose CGD ranges from 0.8 to 1.4, even though simulation are performed at a reference density of 0.6; in a homogeneous system CGD peaks around reference density. It is an indication that there exist a preferred range of densities in in these dense domains. Many simulations of glass transition in this binary mixture is performed at ρ = 1.2, CGD computed at this density near glass transition temperature in our model is given by the solid line in Fig. 3, which show same CGD as that found at low density. In order to look into various structures formed in these dense domains, we compute the bond order parameter which indicate the existence of crystalline as well as amorphous structures [29,41]. In the liquid where there is only a transient orientational order, bond order parameters are small. When local structures stabilize as temperature is reduced there is considerable increase in the value of bond-order parameters. The local bond order parameter is defined as [29,42] q l (i) = 4π 2l+1 l m=−l |q lm (i)| 2 where, q lm (i) is the local orientational order vectors. The averaged form of local bond order is given byq l (i) = 4π 2l+1 l m=−l |q lm (i)| 2 where,q lm (i) is the local orientational order vectors averaged over particle i and its surroundings. We computed the average local bond order, see Fig. 4 for the system and compared with standard values q 4 =0.08, q 4 =0.03 for bcc, q 4 =0.01,q 4 =0.08 for hcp, See Tab.I of Ref. [42]. A similar result has been obtained in the calculation of bond order parameters at high densities in the same system [33]. This shows structural similarities of dense domains obtained in this simulations with density distribution of a glass forming system near its glass transition temperature.
Many studies support the view that glass transition occurs due to the dynamical arrest or cage formation in glass systems. One way to characterize the single particle cage is by studying the mean square displacement r 2 (t) = 1 is the particle displacement of i th particle over a time t. The cage formation in the supercooled and glassy system is indicated by the long sub-diffusive regime that appears at lower temperatures [4,34,35,36,37] near T g . Fig. 5 shows plot of mean square displacement for ρ = 0.6 for various temperatures. At high temperatures r 2 (t) of the particles changes from ballistic region to diffusive region. Below temperature T =0.6, r 2 (t) shows a sub-diffusive region between ballistic and diffusive regions as shown in the Fig. 5. A sub-diffusive region which extends when temperature reduces called plateau region, suggest cage formation in this system [4,16], that results in slow relaxation, which is found in the studies of many glass forming liquids at high density ρ = 1.2 [4,34,35,36,37] and in colloids [43].

Summary and Conclusions
We have done long time molecular dynamic simulation study on the formation of glass like domains at the unstable region of the 80:20 binary mixture which is widely used in the studies of glass transition at high densities [4,5,6]. The radial distribution functions g AA (r), g BB (r) and g AB (r) computed have shown the existence of multiple local structures at low temperatures. In comparison to studies of same system at high densities [40,4], the second peak of the radial distribution functions split into multiple peaks. The local structure of these dense domains are studied using method of CGD given by Testard et. al. [24]. We found that uniform phase at high temperature splits into gas and liquid phases at intermediate temperature. At very low temperatures the system show coexistence of dense domains whose density distribution is comparable to that near the state points used in simulation studies of glass transition at high density [4,5,6]. Multiple structures that appear in these dense domains are characterized by the bond order parameters [29]. We find the existence of mixture of bcc, hcp structures in this system which is also found at high density [33]. The cage relaxation properties of these dense systems are studied using mean square displacement of particles of the system. The mean squared displacement shows clearly the growth of the sub-diffusive regime which indicate the formation and caging that is found in the computational studies of models of glassy and supercooled liquids [4,34,35,36,37]. Thus our studies indicate presence high density domains in phase separating low density regime of 80:20 binary mixture, whose structure is comparable to what found high density studies of same system near glass transition. The nature of diffusive motion of particle at these state points is also comparable to that found in high density studies of glass transition in the same system. PPJ acknowledge DRDO for partial financial support and HPC facility at IIT Mandi for computational support.