Dynamics and Elastic Properties of Glassy Metastable States

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 TD. Dynamic heterogeneity is also observed at TD. 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 TD. It has been shown that the rigidity and mobility of glassy metastable states are compatible under dynamic heterogeneity.


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 Appendixes A.1-A.3.

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.

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 Sections 2.2-2.5 and calculate D (except otherwise stated). Because of this strong non-linearity, especially at low temperatures, the variance of 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 ≤ 130, slight diffusion D < 4.0 × 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 ∆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.

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 ≥ 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].

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(r) shows a clear structural transformation. State H is stable for temperature regions 65 ≤ T ≤ 95 and 105 ≤ T ≤ 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 ≤ T ≤ 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 ≤ 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 ≤ 41,000 for T = 106. At T = 102 (green), MSD drastically increases upto t ≃ 3.235 × 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 ≤ 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 × 10 4 , and (b) t = 3.25 × 10 4 are shown. The MSD shows large fluctuation and increase before t = 3.25 × 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).

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.

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 ≤ 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, which is the probability of displacement d in duration of time ∆t, where r i (t) and r i (0) 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 ∆t, but also on the absolute time where the measurement is taken. To show the effect of time duration ∆t alone, the finishing time of the measurement is common in Figure 8 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 ≤ 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.

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 where l is the dimension, ∆r = r(t) − r(0) is the displacement of the particles in duration of time t and ⟨∆r 2 ⟩ 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 × 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.

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

Temperature Dependence of Modulus of Spontaneous Elastic Tension/Compression
In Figure 10, the temperature dependence of the longitudinal E ∥ and transverse E ⊥ elastic modulus is shown. All states, except state J, elastic modulus E ∥ and E ⊥ 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 ∥ and E ⊥ 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 ∥ and E ⊥ 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 ∥ and E ⊥ is not affected even at temperatures higher than T D where the particles get mobile. For supercooled liquids, E ∥ and E ⊥ are at least two orders of magnitude smaller than other metastable states, thus not shown in Figure 10.

Temperature Dependence Spontaneous Strain Ratio
In Figure 11, the ratio of spontaneous strains µ = ⊥ are plotted against temperature. For comparison, some data for the supercooled and thermodynamic equilibrium liquids are also shown. The strain ratios are µ ≃ 0.5 for the liquids, like the Poisson's ratio [33][34][35]. The relation between elastic modulus for isotropic medium is This equation shows that as µ → 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 → 0 at melting.

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 η to pressure. The Birch-Murnaghan equation of state [37] in the second degree is, where η = 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 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].

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.

Appendix A. Computer Simulation Methods
The principle of corresponding states [46] holds for systems with simple pair-wise interaction, i.e., the inter-particle potential having the functional form φ ij = f (r σ). By introducing the mass m, length σ, and energy , as fundamental units, other physical quantities may be defined by them; such as pressure P * = Pσ 3 , temperature T * = k B T , and molecular volume v * = V (Nσ 3 ). Under these reduced units, a thermodynamic state point determines a set of corresponding states which applies to thermodynamics, structural, and dynamic properties [47].
Furthermore, if particles interact by repulsive force alone, a single reduced variable defines the excess properties. For inverse power potential φ ij = r ij σ −n , the reduced density ρ * = ρσ 3 ( k B T) 3 n determine the state. A single isotherm, isochore, or isobar is sufficient to determine the entire phase diagram of thermodynamic equilibrium [48]. This can be understood as the average excluded distance d c of a pair of particles (the effective core) equals the average kinetic energy [14] defined by the systems temperature T, It reveals that the temperature dependence of the effective excluded volume d 3 c = σ 3 ( k B T) 3 n is the crucial factor. Similarly, for the WCA potential, the relation will determine the effective core d c at temperature T. Again, the temperature dependence of excluded volume d 3 c is the crucial factor. Even when there is anisotropy in the shape of the particle, similar comprehension leads to the scaling law of the system [49]. The equation of state under the reduced units for these simple models will be given by ; expressed by the inverse function of the pair potential φ −1 ij .

Appendix A.1. Symplectic Integrator for Soft Matter
To properly investigate thermodynamic and dynamic properties of condensed matter, the system should not be under non-hydrostatic stress. To achieve simulations under hydrostatic pressure, not only the volume but also the shape of the simulation cell must change according to the stress tensor caused by the constituent particles. In the method used in this work, the shape and size of the simulation cell is described by an anisotropic factor α and an isotropic volume factor Q where the volume of the simulation cell is V = αQ [13]. The anisotropic factor is the ratio of the simulation cell lengths in z and xydirections (length in x and y-directions are the same in this study), i.e., α = L z L xy . The isotropic volume factor Q expands and contracts isotropically, while α expresses the change in shape. These variables α and Q change spontaneously, evolving the system towards hydrostatic pressure with the scalar pressure balancing with the given value P ex . Periodic boundary conditions are applied to x, y, and z directions. In this method, the Hamiltonian is solved by a symplectic integrator which preserves the structure of the Hamilton's equations of motion. Thus, high precision and stable trajectory is guarantied. System size effects are small in this method and notable only near the phase transition. See Figure 6 of [50].

Appendix A.2. Calculation of Metastable States
All the metastable states investigated in this paper are morphologically changed from the pristine forms of the amorphous solid or supercooled liquid reported in [13] to settle in a basin under isothermal and hydrostatic pressure. That is the basin of the free energy profile of the total system where the crystal state is in the global minimum and metastable states are in local minima. The value of the given pressure P ex in these simulation has an effect on the sharpness of the transition. When P ex is high, the melting transition temperature T m becomes high, thus leading to a smaller d c near T m . Since the potential becomes steeper when d c becomes smaller, the particles are effectively harder at higher pressures. The potential landscape becomes steeper as well.
To seek the stable temperature range of a certain metastable state, a series of calculations at different temperatures are done independently from a common initial configuration. During these calculations, we occasionally observe the system transform to a different state. Some examples of such structural transformation are given in Figure 3 of [9]. The transformation occurs in a short time duration, compared to the calculated time. Such transformations are caused by a rare event knocking the system over the energy barrier to a new state. Information of the free energy profile of WCA phase space, such as the approximate depth of the basin, is obtained through such calculations. Occurrence of such transformation not only appears in the static values, but also in the dynamics such as the MSD.

Appendix A.3. Model and Initial Configurations
To investigate to what extent the excluded volume effect alone can describe the properties of glass, soft isotropic monodispersed spheres are employed as a model system. Pairwise potential of WCA [12] is used. Reduced units are employed throughout our work. In addition to the metastable states explained below, the dynamic and elastic properties of crystalline solid are investigated for comparison. The given pressure is set to P ex = 1.0 × 10 4 where the melting temperature from equilibrium crystalline solid to liquid is T m = 153. The effective particle size at the melting temperature is d c = 0.653.
Metastable states for monodispersed WCA spheres are initially reported in [13] as amorphous solids. A total of seventeen states, including the supercooled liquid, are reported in [8,9]. Among them, four states in [9] are chosen, i.e., states G, J, M, and O; the labels of each state will remain the same as the original article to avoid confusion. The original calculations were done for number of particles N = 1344. In this work, the configuration of these states were stacked four times and used as the initial configuration for a new system size N = 5376. The initial configurations to construct the N = 5376 system, are those of N = 1344; states G, J, M, and O at temperatures T = 96, 100, 104, and 106, respectively, at which these states were first identified [13].
Note that the supercooled liquids reported here are metastable states obtained independently at each temperature where the thermal fluctuation cannot overcome the energy barrier at each temperature in the calculated time. The supercooled liquid exists as a thermodynamic metastable state at T ≥ 109 for the original system size N = 1344. At lower temperatures, solidification occur and marked by a clear change in MSD and the volume/density [8]. The supercooled liquid is calculated for T ≥ 110 for system size N = 5376 where the typical total calculation time is t = 4000 to 6000. The initial configuration was the original system (of N = 1344 at T = 110) stacked four times. For crystalline solid, data at T = 65 of N = 1344 was used to construct the initial configuration. The crystalline state of N = 5376 was calculated to t = 3000 at T = 152 and did not melt.
The mass of the barostat is M = 1.0 × 10 −3 and the mass of the thermostat is K = 1000 throughout this work. All calculations were done with time step dt = 1.0 × 10 −4 .

Appendix A.4. Calculation of Transport Coefficient
A linear fit to the mean square displacement against time gives the diffusion constant D.
Since the velocity v i of each particle i is calculated at every time step for the molecular dynamics, the displacements ∆x i (t) = ∑ t v x dt calculated from the velocities v β (in each direction β = x, y, z) are tracked along the time evolution and used for calculating the transport properties. This allows the calculation of transport properties unaffected by the periodic boundary conditions.

Appendix A.5. Calculation of Elastic Properties
The spontaneous elastic modulus for each state are calculated from the stress versus strain plot as in [51]. Although we do not have any uniaxial load in our simulation, the geometry of the spontaneous fluctuations is those of Poisson's ratio measurements [33,34]. Instead of a uniaxial load, we have a spontaneous uniaxial shape fluctuation (in addition to the volume fluctuation) of the simulation cell where the strains can be directly measured. Simultaneously, the internal stress tensor is tracked by the simulation method. The longitudinal and transverse strains/stresses can be measured regardless of the state of matter in our method. The isotropic volume factor Q is related to the cell length in x and y-direction by L xy = Q 1 3 , and will give a modulus of spontaneous elastic tension/compression. This modulus E ⊥ is calculated by plotting the transverse stress P ⊥ = (P xx + P yy ) 2 against the transverse strain ⊥ = L xy ⟨L xy ⟩ = Q 1 3 ⟨Q 1 3 ⟩. Similarly, the uniaxial elastic modulus E ∥ is obtained from longitudinal stress P ∥ = P zz and longitudinal strain ∥ = L z ⟨L z ⟩. Values of the transverse strain ⊥ against longitudinal strain ∥ give the spontaneous stress ratio. The instantaneous stress and strain are measured at each ∆t = 1.0 × 10 −2 (100 time steps) for a time duration of ∆t = 1000 (1.0 × 10 7 time steps) at the end of each simulation. The modulus of volume expansion/contraction (the spontaneous bulk modulus) is calculated from collecting the instantaneous values of mean stress P = (P xx + P yy + P zz ) 3 and strain v = V ⟨V⟩ = αQ ⟨αQ⟩ of the simulation cell. The elastic modulus and the strain ratio are calculated when the cell shape is relatively stable for time duration ∆t = 1.0 × 10 3 , otherwise we only calculate the bulk modulus.
In Figure A1, the spontaneous stress and strain are plotted for state G at temperature T = 132. Each dot represents an instantaneous value measured at 20,000 < t ≤ 21,000. Red lines are the linear fit to all points. In the strain-stress plot, the data points form a ellipse showing that the state is viscoelastic to some extent. Nevertheless, we make a linear fit to all data points. For Figure A1, the red line which is the linear fit to all points does not go through the middle of area of dots. This is because data points are more distributed on the high-strain/low-stress area which shows that it is more difficult to compress than to expand in this system.