Dynamics and Elastic Properties of Glassy Metastable States

Abstract

By a molecular dynamics (MD) simulation method which ensures the system will be under hydrostatic pressure, dynamic and elastic properties of glassy metatstable states are investigated. In the MD method, the simulation cell fluctuates not only in volume but also in shape under constant hydrostatic pressure and temperature. As observed in experiments for many glass forming materials, metastable states in our simulation show a sharp increase in mean-square-displacement at certain temperatures $T_D$. Dynamic heterogeneity is also observed at $T_D$. Elastic properties are calculated from stress and strain relations obtained from the spontaneous fluctuation of internal stress tensor and simulation cell parameters. Each investigated state shows distinctive dynamics while maintaining solid-like elastic properties. The elastic properties stay intact even above $T_D$. It has been shown that the rigidity and mobility of glassy metastable states are compatible under dynamic heterogeneity.

1. Introduction

Virtually any material can be formed as a glass under the proper experimental conditions [1]. Thus, there is no apparent reason to excluded monodispersed simple molecular models as an exception. A simple model has an advantage that the knowledge of thermodynamic equilibrium states and phase transitions among them are abundant, making them a preferred reference state for a theoretical description. However, molecular dynamics (MD) simulations of glassy states to date frequently use binary models [2,3,4]. These binary mixtures are made with size ratio which cannot form a crystal and size polydispersity larger than the value where the solid–fluid co-existence region disappears in the phase diagram [5]. The state of lowest energy of such binary mixtures is not known. In addition, most of these simulation are constant volume studies.

There are several reasons why simulations at constant pressure are interesting. First of all, most of the experiments are done under constant pressure at both ambient and high pressures [6]. Second, inherent structure excitation profiles are steeper at constant pressure than at constant volume [7]. This suggests that glassy states might be confined in narrow basins for long enough to permit a through analysis at constant pressure. It is known that constant pressure MD simulations at high pressures lead to sharper phase transitions and the metastable states possess distinct physical properties, not only for soft spheres [8,9], but also for molecules with anisotropic shape [10,11]. In these systems with anisotropic constituents, multiple metastable states of hexatic smectic B (HexB) liquid crystal with discrete values of order parameter appear. At temperatures close to the 1st order phase transition to smectic A phase, barrier crossings from one HexB metastable state to another is often observed [11]. At lower temperatures, these distinct metastable states only show oscillation around the local minima [10].

Amorphous solids of monodispersed Weeks-Chandler-Andersen (WCA) spheres [12] trapped in metastable states at high pressure was first reported in [13]. The MD simulation method used in the above study is a constant hydrostatic pressure method where not only the size but also the shape of the system change to avoid non-hydrostatic stress under constant pressure. This method was developed due to the fact that systems of anisotropic molecules give rise to non-hydrostatic stress not only in constant volume simulations, but also in constant pressure methods [14]. It has been shown that not only quantitative deviation but also structural transformation occurs by non-hydrostatic stress. A series of stress-controlled MD simulation methods have been developed to avoid such artifacts feigned by non-hydrostatic stress [15,16,17,18]. Furthermore, it has recently been noticed that residual shear stress exists near inherent structures in constant volume MD simulations as a direct consequence of the imposed boundary [19,20]. These studies demonstrate that glassy states are also prone to non-hydrostatic stress and more broadly to system size confinement effects.

Since glass is ubiquitous, it is difficult to distinguish the characteristics of the material itself and those common to the glassy states. By investigating several different metastable states obtained by the same monodispersed system, we aim here to make clear the relation of diffusion dynamics and elastic properties. All particles interact through the same pair-wise potential, thus the differences spontaneously emerge solely from the state of particle assemblage.

It has been experimentally revealed that a wide range of glass forming materials show a sharp increase in mean square displacement (MSD) near the glass transition temperature; from proteins and biologically relevant systems to simple liquids as othroterphenyl and selenium [21,22,23,24,25,26,27]. In this paper, we pick up several metastable states of the same model system which show a sharp increase in the MSD and investigate not only the dynamics but the elastic properties of those states in isothermal–isobaric ensemble. In Section 2, dynamics properties are presented. In Section 3, elastic properties of the same metastable states are presented. In Section 4, we summarize and discuss future directions. Data are presented in reduced units. For simulation details, see Appendix A.1, Appendix A.2 and Appendix A.3.

2. Dynamic Properties

Several metastable states which show a sharp increase in MSD are chosen for simulation of a larger system size $N=5376$ to investigate the stability, dynamics, and elastic properties. The system size investigated in this paper is four times larger than the original work [8,9]. All the states have distinctive values of densities and radial distribution functions. In these larger systems, it is confirmed that the metastable states are “stable enough” to warrant further investigation. We use the same label of each state as in the original report [8,9] for clarity. In the first subsection, we observe the temperature dependence of dynamic properties and show that MSD drastically increase beyond an inflection temperature $T_D$. Next, we observe the MSD of each state one by one, since characteristics of each state are drastically different. In the final subsection, the dynamical heterogeneity is discussed. For comparison, properties of crystals, and thermodynamic equilibrium liquids as well as supercooled liquids, are calculated in addition to the glassy metastable states.

2.1. Temperature Dependence

In Figure 1, the temperature dependence of the diffusion coefficient D (time average of MSD) for different states is shown. The melting temperature is $T_m=153$ in the investigated system. For non diffusive states, the MSD is often stepwise/oscillatory as shown in the following subsections. The determination of D will be difficult in such situation. In principle, D should be determined in a time scale (or length scale) where MSD becomes linear against time. Nevertheless, we make a linear fit to the MSD data against time shown in the MSD figures in the following Section 2.2, Section 2.3, Section 2.4 and Section 2.5 and calculate D (except otherwise stated). Because of this strong non-linearity, especially at low temperatures, the variance of D are temperature dependent and larger for lower temperatures in general. When a clear change in the nature of MSD occurs during a simulation at a certain fixed temperature, the values of MSD after that clear change are only used for calculating D. Such a situation occurs in state J at $T=102$, 106, and 118 (see Section 2.3). See Appendix A.4 for calculation details of MSD.

For the metastable states, an inflection temperature is observed where D increases drastically. We define the temperature where D starts to increase drastically as $T_D$ and use it as a reference temperature to denote the dynamical change. In experimental findings, $T_D$ is correlated to the glass transition temperature [21,22,23,24,25,26,27]. For crystals at temperatures $T\le 130$, slight diffusion $D<4.0\times {10}^{-7}$ is observed. However, in this temperature range, no particles escape from the region beyond the 1st nearest neighbor in a certain time duration (see Figure 2). For temperatures $T>135$, hopping occurs in the crystal. However, it occurs in a spatially isotropic manner in contrast to the metastable glassy states which often show anisotropic hopping. Note that the diffusion resulting from the hopping dynamics in crystal become notably large near the melting transition at $T=153$.

In Figure 2, the temperature dependence of the fraction of mobile particles for each state is shown. The number of particles which move more than the average excluded volume length $d_c$ (effective core) in time duration $\Delta t=1000$ are defined as mobile. See Appendix A for the definition of $d_c$. Note that the measurement time for Figure 2 is limited compared to measurement of MSD and of D in Figure 1.

In state O (symbol blue +), there is a slow oscillation in MSDs at low temperatures where collective movements of particles occur beyond the measuring time. The fluctuation in the values of state O in Figure 2 is due to the choice of measuring time because of the slow oscillating nature of MSD of state O. See Section 2.5 for details.

2.2. State G

State G has the lowest diffusion constant D among the metastable states for a wide range of temperature.

In Figure 3, MSD of state G versus time at different temperatures is shown.

At temperatures $T=114$ and $T=122$, MSD changes only slightly. A stepwise increase in MSD is observed at $T\ge 126$. This is because the values of MSD are extremely small and the hopping nature of individual particles appears in Figure 3. Most of the glassy states of monodispersed softcore systems show hopping diffusion [8,9].

2.3. State J

State J shows intermittent structural rearrangements.

State J has been first obtained as a result of spontaneous structural transformation of state H at $T=100$ in the original system size [9]. The initial configuration of state H is the amorphous solid state reported in Figure 12 of [13]. During the calculation of state H at $T=100$, the volume and the total potential energy decreased and transformed to J state [9]. The physical properties of H and J states are are completely distinct and the pair distribution function $g\left(r\right)$ shows a clear structural transformation. State H is stable for temperature regions $65\le T\le 95$ and $105\le T\le 115$. See Figure 3 of [9] for other structural transformations observed in WCA spheres. These facts show that the connectivity of basins and their depth for these metastable states are temperature dependent.

In Figure 4, the MSD of state J versus time at different temperatures is shown. For system size $N=5376$ reported here, state J seems to be relatively unstable at temperatures $102\le T\le 118$.

At $T=118$ (cyan with larger MSD values), MSD rapidly increases at the beginning; however, it slows down to a constant rate. The diffusion constant D at $T=118$ in Figure 1 is calculated for $4080<t\le$ 41,000. At $T=116$ (larger red values), a stepwise nature due to hopping appears in the MSD. At $T=106$ (magenta, also shown in the inset of Figure 4), string-like hopping motion occurs; however, the rate decreases at t > 24,000. In Figure 1, D is calculated at 24,000 < $t\le$ 41,000 for $T=106$. At $T=102$ (green), MSD drastically increases upto $t\simeq 3.235\times {10}^4$ but both the diffusion and the fluctuation become much smaller at larger times. A landslide in the stacking layer occurs before that time and results in a structural transformation. At $T=102$, the diffusion constant D shown in Figure 1 is calculated at 32,500 < $t\le$ 51,000. For other temperatures not mentioned above, the diffusion constant D in Figure 1 is obtained by averaging through the whole measuring time. At temperature $T=120$, MSD increases drastically compared to lower temperatures.

In Figure 5, displacement plots of state J at $T=102$ and at times (a) $t=1.11\times {10}^4$, and (b) $t=3.25\times {10}^4$ are shown. The MSD shows large fluctuation and increase before $t=3.25\times {10}^4$ (green line in Figure 4). Displacement plots show positions relative from the initial position of all particles projected on a plane. Thus, displacement plots reveal spatially cooperative movements. A straight string-like hopping motion can be clearly discerned in Figure 5. In addition, the center crowded position where particles are fluctuating around the original position, splits into two as shown in Figure 5b. A landslide occurred between the time shown in Figure 5a,b. The length scale of the landslide is about six times larger than the distance of the neighboring particles (see Supplementary Material).

The phenomenon where the displacement center splits into two (along with string-like motion) is also observed at temperature $T=106$ of state J for system size $N=2688$. Thus, state J in this temperature range seems to be susceptible to landslides, due to the small elastic modulus (see Section 3 for further discussion).

2.4. State M

State M is the most diffusive state among the metastable states studied in this paper. A clear temperature dependence of MSD is observed for state M as shown in Figure 6.

The inflection point in the increment of D observed in logscale (Figure 1) is not so clear compared to other states.

2.5. State O

State O is a state with a long trajectory reversal time and distance.

In Figure 7, MSD of state O is shown. The overall temperature dependency is shown in Figure 7a. At $T=128$ and 132, the system is quite diffusive. For temperatures $T\le 120$, MSD of state O shows a long time fluctuation with a large return distance as clearly seen in Figure 7b.

A large oscillation in the MSD shows that there is a collective movement among the particles. This collective movement clearly appears in the self-van Hove correlation function,

$$F(d,t)=\frac{1}{N}\langle \sum _{i=1}^N\delta [d-|{\mathbf{r}}_i\left(t\right)-{\mathbf{r}}_i\left(0\right)\left|\right]\rangle ,$$

(1)

which is the probability of displacement d in duration of time $\Delta t$, where ${\mathbf{r}}_i\left(t\right)$ and ${\mathbf{r}}_i\left(0\right)$ are the position of particle i at time t and at $t=0$, the origin of the measurement. In Figure 8, the self-van Hove correlation functions measured at several time duration against the distance scaled by the effective core length $d_c$ are shown.

Since the large fluctuation of MSD seems quite random but is confined in a certain width for state O (Figure 7b), the self-van Hove correlation function depends not only on the time duration $\Delta t$, but also on the absolute time where the measurement is taken. To show the effect of time duration $\Delta t$ alone, the finishing time of the measurement is common in Figure 8 except for the red dot-dashed line. The fraction of mobile particles (Figure 2) and elastic properties (Section 3) for state O are taken at the same time duration t = 20,000–21,000 shown by the red dot-dashed line. It is clear from Figure 8 why the fraction of mobile particles are high at $T=88$ and 116 in Figure 2. At the time duration t = 20,000–21,000 where the measurements of Figure 2 are done, particles are scattered to relatively large distances, which tend to return close to their original positions at longer time duration as shown by the solid line in Figure 8a,b. Figure 8 clearly shows that the distance for trajectory reversal (rattling distance) extends beyond twice the distance $d_c$. The time scale differences between the trajectory reversal motion and diffusive motion clearly appear in the MSD at low temperatures $T\le T_D$ (Figure 7b). State O has a large return distance regardless of system size and the same characteristics are persistently observed in other system sizes $N=1344$ and 2688.

2.6. Dynamical Heterogeneity at $T_D$

The dynamical heterogeneity is a phenomenon where dynamically active and inactive regions appear in the same system. To measure the degree of dynamical heterogeneity, the non-Gaussian parameter(NGP) is calculated. NGP is defined as

$$\alpha =\left(\frac{l}{l+2}\right)\frac{\langle \Delta {\mathbf{r}}^4\rangle }{{\langle \Delta {\mathbf{r}}^2\rangle }^2}-1$$

(2)

where l is the dimension, $\Delta \mathbf{r}=\mathbf{r}\left(t\right)-\mathbf{r}\left(0\right)$ is the displacement of the particles in duration of time t and $\langle \Delta {\mathbf{r}}^2\rangle$ is the MSD. The NGP plotted against time will show the degree and time scale of the dynamical heterogeneity.

In Figure 9, the three dimensional ($l=3$) NGP near $T_D$ is shown. In state G (Figure 9a), NGP become notably positive after $t>{10}^3$, where a sharp asymmetric peak appear near $t=6\times {10}^3$. The calculated time is not enough to observe whether the decay is linearly falling off. For state J (Figure 9b), the decay of NGP is very slow and show large fluctuation. The NGP curve of state M (Figure 9c) resembles that of colloidal glasses in experimental observations [28,29,30,31,32]. Reflecting the large particle reversal movements in state O, the NGP curve (Figure 9d) shows two separated large peaks, consisting of smaller peaks. All the states show distinctively different characteristics in the NGP curves, which reveals that the nature of dynamical heterogeneity is different for each state.

3. Elastic Properties

To calculate elastic properties, spontaneous stress and strain in longitudinal and transverse directions are measured. For calculation details, see Appendix A.5.

3.1. Temperature Dependence of Modulus of Spontaneous Elastic Tension/Compression

In Figure 10, the temperature dependence of the longitudinal $E_{\Vert }$ and transverse $E_{\perp }$ elastic modulus is shown. All states, except state J, elastic modulus $E_{\Vert }$ and $E_{\perp }$ in log-scale show approximately a linear dependence on T. Only state O shows a tendency of softening at $T=132$. Note that the sharp increase in D occurs around $T_D=115$ for state O (Figure 1), revealing that the elastic properties stay intact even at higher temperatures where the particles get mobile.

For state J, the values are scattered in Figure 10, and difficult to discern a relation against T. Note that at the calculated lowest and highest temperatures ($T=99$ and 120) of state J, the values of $E_{\Vert }$ and $E_{\perp }$ are higher than other temperatures and overlap with the value of state M at $T=99$. As discussed in Section 2.3, MSD at $T=102$, 106, and 118 showed peculiar characteristics and D was calculated for a limited time at the end of the simulation. At $T=102$ which mark the 2nd smallest values of $E_{\Vert }$ and $E_{\perp }$ of state J, a landslide was observed (Figure 5). The fraction of mobile particles of state J was high only for $T=120$. For other states, the temperature dependence of $E_{\Vert }$ and $E_{\perp }$ is not affected even at temperatures higher than $T_D$ where the particles get mobile. For supercooled liquids, $E_{\Vert }$ and $E_{\perp }$ are at least two orders of magnitude smaller than other metastable states, thus not shown in Figure 10.

3.2. Temperature Dependence Spontaneous Strain Ratio

In Figure 11, the ratio of spontaneous strains $\mu ={ϵ}_{\perp }/{ϵ}_{\left|\right|}$ are plotted against temperature.

For comparison, some data for the supercooled and thermodynamic equilibrium liquids are also shown. The strain ratios are $\mu \simeq 0.5$ for the liquids, like the Poisson’s ratio [33,34,35]. The relation between elastic modulus for isotropic medium is

$$\frac{S}{B}=\frac{3(1-2\mu )}{2(1+\mu )}.$$

(3)

This equation shows that as $\mu \to 1/2$, the shear modulus S drastically decreases if the change in bulk modulus is small [36]. We show, in the next section, that the bulk moduli in supercooled and equilibrium liquids have a smaller temperature dependence compared to the metastable states, thus indicate that shear modulus $S\to 0$ at melting.

3.3. Temperature Dependence of Spontaneous Bulk Modulus

In Figure 12, temperature dependence of bulk modulus for the glassy states and liquids is shown.

Finite-strain theory has been applied extensively to problems in geophysics. The theory relates compression/expansion $\eta$ to pressure. The Birch–Murnaghan equation of state [37] in the second degree is,

$$P=\frac{3}{2}{B}_0\left({\eta }^{-7/3}-{\eta }^{-5/3}\right),$$

(4)

where $\eta =V/V_0$, V is the volume at P, and the bulk modulus $B_0$ is that of volume $V_0$ at a reference pressure $P_0$. This semi-empirical Equation (4) requires a reference value of the bulk modulus $B_0$ for volume $V_0$ to be determined from experiment. In geophysics, the reference values are usually those at the ambient pressure. Using the above Relation (4), the temperature dependence of the bulk modulus is

$$B_T=-V{\left(\frac{\partial P}{\partial V}\right)}_T=\frac{1}{2}{B}_0\left(7{\eta }^{-7/3}-5{\eta }^{-5/3}\right).$$

(5)

The temperature dependence of bulk modulus B of glassy states are fit to Equation (5). The reference values $B_0$ and $V_0$ are taken from the values at $T=80$ for crystal (red line) and state M (green line).

Most of the glassy states conform to Equation (5) obtained from the Birch–Murnaghan equation of state and fit between the region inside the red and green lines. However, when the temperatures get closer to the melting transition, the crystal state (•) becomes drastically softer than Equation (5) predict. In contrast, the temperature dependence of the liquids (blue line) is small and do not decrease drastically for temperatures higher than the melting temperature $T_m=153$. This conforms to the fact that when the Poisson’s ratio is close to $1/2$, the change in bulk modulus is small [36].

4. Concluding Remarks

Molecular dynamic simulations which allow anisotropic fluctuations of the cell shape to preserve hydrostatic pressure have been conducted. Using a simple monodispersed model, we have shown that below the melting temperature $T_m$, there exists multiple metastable states with distinctive dynamics. They are glassy metastable states where diffusion constants increase sharply at a certain temperature $T_D$. Around $T_D$, dynamic heterogeneity is observed. The drastic increase in diffusion constants at $T_D$ does not induce change in the temperature dependence of the elastic modulus. The elastic properties show that these metastable states are basically solids even at temperatures higher than $T_D$. The rigidity and mobility of glassy metastable states are compatible even at higher temperatures than $T_D$. There exists a rigid supporting framework which allows other parts to flow. Because of this framework, dynamical heterogeneity is sustained. How to quantitatively assess these structures, especially the supporting framework, will be addressed in future work. The bulk moduli for the supercooled/equilibrium liquids only show a slight temperature dependence compared to glassy metastable states.

Once we admit that there exist multiple metastable solid states with glassy dynamics, the scenario behind polyamorphous and liquid–liquid transition [38,39,40,41,42,43,44,45] will attain many possibilities. The temperature range of each metastable state might be different; thus, beyond those bounded temperatures, the state is destined to transform to another state. These transformations may also occur among two diffusive states. In the original system, we observed a clear change in MSD at a transformation among two different states (see Figure 11 of [8]). To address these important subjects, not only simulations of larger systems near $T_D$ as performed here, but also in the whole temperature range, especially towards $T_m$, is needed and remain for future work.