# Tuning Pairwise Potential Can Control the Fragility of Glass-Forming Liquids: From Tetrahedral Network to Isotropic Soft Sphere Models

## Abstract

We perform molecular dynamics simulations for a glass former model proposed by Coslovich and Pastore (CP) over a wide range of density. The density variation can be mapped onto the change of the potential depth between Si and O interactions of the CP model. By reducing the potential depth (or increasing the density), the anisotropic tetrahedral network structure observed in the original CP model transforms into the isotropic structure with the purely repulsive soft-sphere potential. Correspondingly, the temperature dependence of the relaxation time exhibits the crossover from the Arrhenius to the super-Arrhenius behavior. Being able to control the fragility over a wide range by tuning the potential of a single model system helps us to bridge the gap between the network and isotropic glass formers and to obtain the insight into the underlying mechanism of the fragility. We study the relationship between the fragility and dynamical properties such as the magnitude of the Stokes-Einstein violation and the stretch exponent in the density correlation function. We also demonstrate that the peak of the specific heat systematically shifts as the density increases, hinting that the fragility is correlated with the hidden thermodynamic anomalies of the system.

## 1 Introduction

The current understanding of the mechanism behind the dramatic slowing down of supercooled liquids near the glass transition temperatures still remains far from complete. Although underlying mechanism of the glass transition is believed to be universal, many dynamical properties of the glass formers are diverse and system-dependent. Among them, the concept of the fragility, or the temperature dependence of the relaxation time and transport coefficients such as the viscosity and diffusion constant, is one of the most important but the least understood problems. For some glass formers, which are called the strong liquids according to the Angell’s classification, the relaxation times obey the Arrhenius law, while others, called the fragile liquids, show a strong departure from the Arrhenius law in their temperature dependence of the relaxation times [1, 2]. Representative strong glass formers are and , whose local molecular configurations are characterized by anisotropic network structures. On the other hand, the fragile liquids such as -terphenyl and toluene tend to have more isotropic and compact local structures. The concept of the fragility has been playing a key role in the study of the glass transition. Many experimental studies have reported that there exist correlations between the fragility and various thermodynamic, structural, mechanical, and dynamical properties of the glass formers [2]. For example, it has been observed that the entropy and specific heat [1, 3], the elastic constants [4, 5], and the spatially heterogeneous dynamics [6, 7] sensitively depend on the fragility of the systems. Theoretical understanding of the fragility, on the other hand, is left at the phenomenological level. Virtually all theoretical scenarios of the glass transition proposed so far, including the Adam-Gibbs theory[8], the random first order transition theory [9], energy landscape picture [10], geometric frustration scenarios [11, 12, 13, 14], elastic models [15, 16], soft modes [17], and the dynamic facilitation scenario [18], have their own explanations of the fragility but the first-principles and microscopic theory to quantitatively describe the fragility is still lacking [10, 19, 20, 21, 22].

Difficulty to gain an unified picture of the fragility out of the experimental observations is obviously due to the diversity and complexity of real molecular or polymeric glass formers. The computer simulation is an ideal tool to avoid this difficulty because one can gain all microscopic information of dynamics for relatively simple model systems. Popular glass former models in the simulation studies, such as the soft-sphere (SS) [23] and the Lennard-Jones (LJ) potential liquids [12, 24, 25], are regarded as typical fragile glass formers. On the other hand, strong glass former models with anisotropic networks, where the particles are connected by covalent bonds, such as and , have also been studied extensively by simulations [26, 27, 28, 29].

There have been several attempts to understand the origin of the fragility systematically by simulations. In many studies, the isotropic potential systems such as the SS and LJ potentials have been employed. The fragility has been controlled using various methods, for example, by varying the density or pressure of the system [30, 31, 32, 33, 34], by changing the polydispersity [35, 36, 37] and the size ratio of the multi-component mixtures [12, 38], by tuning the interaction potential [31, 39, 40], or by truncating the attractive part of the LJ potential [41, 42]. The fragility was also found to be sensitive to many-body interactions which influence the local geometrical frustration [43, 44]. An idea to sort out thermodynamic and purely kinetic contributions to the fragility has been discussed recently [45]. The external parameters such as the impurities [46, 47] and the curvature of the non-Euclidean space [48] have been proposed as a new method to control the fragility. The findings of these studies, however, are still not sufficient to unravel the entire picture of the fragility. Aside from the obvious drawback that the time window which the current simulation can cover is narrow compared with experiments, the one of the obstacles of computational studies is that the variation of the fragility which the simple glass former models can encompass is small [30, 32, 35, 36, 12, 38, 37, 31, 39, 40, 41, 42]. Furthermore, the variation of the fragility is often masked by an apparent density-temperature scaling law, which makes it difficult to extract the generic mechanism controlling the fragility [31, 32]. There also remains the nagging question of how to decompose the kinetic and purely thermodynamic contributions of the fragility [45]. Even more important is to bridge the gap between the fragile isotropic systems and the strong network glass formers, which has been studied almost in different arenas in the past. Given these circumstances, it is beneficial to consider a simple model system which can cover a wide range of fragility, spanning from the strongest network glass former down to the very fragile isotropic system, simply by tuning a system parameter.

In this paper, we numerically study a simple model glass former originally introduced by Coslovich and Pastore (CP) as a model strong glass former mimicking [29]. The CP model is a binary mixture of spherical atoms interacting by the isotropic Lennard-Jones potential with a very strong attraction between Si and O atoms, which allows them to form anisotropic tetrahedral network structures. In this study, the fragility of the CP model is examined over an extremely wide range of density from the order of unity up to virtually infinity. Varying the density is equivalent with changing the potential depth between Si and O atoms of the CP model. By reducing this potential depth (or increasing the density), the anisotropic tetrahedral network structure observed in the CP model is eventually transformed into the isotropic structure with the purely repulsive soft-sphere potential. Thus, it enables one to control the fragility from strong to fragile systematically.

Here, we summarize the results of this paper. By extensive simulations, dependence of the structure on the density is monitored using the radial distribution functions, the static structure factors, and the coordination numbers. By measuring the temperature dependence of the relaxation time obtained from the density-density time correlation function, we observe that the fragility systematically changes from strong to fragile with the change of the structures. We carefully examine the density-temperature scaling for the temperature dependence of the relaxation time [49, 50, 51, 52, 53, 54, 55]. It is confirmed that the relaxation time does not collapse onto a master curve. This implies that the observed fragility variation is not superficial but it is due to a qualitative change of the underlying mechanism controlling the slow dynamics. We investigate correlation between the fragility and several dynamical observables, i.e., the magnitude of the Stokes-Einstein (SE) violation and the exponent of the stretched exponential relaxation of the density correlation function. It has been argued that these observables are manifestations of spatially heterogeneous dynamics, or the dynamical heterogeneities, and they are intimately correlated to the fragility [56]. We confirm that more fragile systems show stronger SE violation, whereas a clear-cut correlation between the stretch exponent and the fragility is not observed. Finally, we also discuss the possibility that the fragility is related to the peak position of the specific heat observed in its temperature dependence.

## 2 Simulation Methods

In this section, we describe the simulation methods.

### 2.1 The Coslovich-Pastore model

Coslovich and Pastore (CP) have introduced a simple (strong glass former) model [29]. Its dynamical and structural properties have been intensively investigated for fixed densities [57, 58, 59, 60]. The advantage of this model is its simplicity compared to other established models [26, 27, 28]. Although its potential is given by a combination of the soft-sphere (SS) and the Lennard-Jones (LJ) potentials, the atoms in the CP model form an anisotropic tetrahedral network structure similar to those of the realistic models [29]. According to References [61, 62], there are three conditions for a binary mixture with spherical potential to generate the tetrahedral network structure; (1) The composition ratio is , where and are the number of particles for species 1 and 2, respectively. (2) The potential is non-additive, i.e, the range of the interaction between the species 1 and 2 is not a simple sum of their diameters. (3) The attractive interaction between different species is strong. The CP model satisfies these conditions.

The CP model is a binary mixture in three dimensions whose composition ratio is and the mass ratio is . The species 1 and 2 correspond to and atoms of , respectively. The interaction potentials between two particles are given by

(1) |

where is a constant which was set to unity in the original CP model [29]. We specify the particles by the Roman indices and the species by the Greek indices . The parameters are set to , , , and . The last condition warrants that the potential is non-additive, i.e., . is truncated at . In order to ensure the continuity of the potential, we add a switching function, , up to . The switching function connects two points, and , continuously and is defined as

(2) |

We used and .

### 2.2 Tuning of the potential depth

In this study, we introduce a parameter in the second term (the attraction part) of in Equation (1). is the strength of the attraction between different species. The depth of the potential well is written in the unit of as . Hereafter, we use instead of . In Figure 1, we show for various ’s. () corresponds to the original CP model and () corresponds to a simple soft-sphere (SS) potential. Note that the potential for is slightly different from the conventional SS potential model which has been extensively studied in the past [23], since the range of the interaction is still non-additive. Our model can seamlessly connect the original CP model to the purely repulsive SS potential model. As mentioned above, the atoms with large , i.e., the strong attractive interaction, tend to form the tetrahedral network structure. This structure is broken as reduces and eventually the local structure becomes isotropic and compact in the limit of . Mathematically, changing at a constant density is equivalent to changing at a constant . We explain the relation between and using a simple scaling argument, following the similar argument for the SS potential model [63]. The Hamiltonian of the system is written as

(3) |

where , , , and are the Boltzmann constant, momentum, mass, and the potential given by Equation (1), respectively. Introducing the unit of the length , this Hamiltonian can be rewritten as

(4) |

where , , , and are reduced parameters scaled by , , , and , respectively. The time unit is defined by . is a dimensionless parameter commonly used for the SS potential [63]. Note that there exist two dimensionless parameters which control the system: and . This means that, for a fixed (or the temperature), changing the density with a fixed is equivalent with changing for a fixed density. Therefore, if one chooses the purely Lennard-Jones system () with a density as a reference system, one can infer thermodynamic and dynamical properties at arbitrary densities by tuning for the fixed . The mapping from to is given by , or equivalently in terms of by

(5) |

From this mapping, one can cover the whole range of density up to by tuning (or ) from a finite value down to 0.

### 2.3 Details of molecular dynamics simulations

We perform the molecular dynamics (MD) simulations for various ’s. We use , , , , , , , and at a constant density . Correspondence of various ’s for a fixed density to various densities for a fixed is summarized in Table 1. By tuning , we explore a wide range of density from up to infinity of the original CP model. The number of the particles is with . The temperature range which we perform the simulation are listed in Table 1. Hereafter, we use , , and , as the units of length, time, and temperature, respectively. The micro-canonical MD method with the periodic boundary condition is used to produce the trajectories. A time step is chosen throughout this study. For the calculations of the dynamic and thermodynamic quantities, we average over four independent simulation runs.

## 3 Results

### 3.1 Structural properties

We first study the structural properties of the CP model for various potential depth . We start with the radial distribution function . Figure 2 shows for each at the lowest temperatures in our simulation which are listed in Table 1. For (see Figure 2 (d)), exhibits a very sharp peak at followed by the very broad and small first minimum which persists up to beyond which the small and broad second peak appears. In contrast, the first peaks of and are observed at and , respectively. This indicates that the particles of species 1 and 2 are strongly bonded [29]. As decreases, the first peak of gradually decreases and broadens. Concomitantly the height of the first minimum rises to a finite value at around . Finally, at (see Figure 2 (a)), the profiles of , , and become similar to those of typical fragile glass formers such as the SS and LJ binary mixtures [25, 23]. This structural change reflects the fact that the bonds between species 1 and 2 are weakened by decreasing .

To understand the effect of on the structure in more detail, we calculate the coordination number . is the averaged number of the particles in the first neighbor shell of the particle and defined by

(6) |

where is the number density of species and is the position of the first minimum of . In particular, and are useful to characterize the tetrahedral network structure. In , a single Si atom (species 1) is surrounded by four oxygen atoms (species 2), thus, , , and . In Figure 2 (f), we show and at the lowest temperatures for several ’s. We use instead of as the horizontal axis for the sake of visual clarity. At (), we obtain and , corresponding to the perfect tetrahedral network structure [29]. As decreases, both and increase. Eventually, we observe and at . This indicates that the tetrahedral network structure is broken and that the local configuration of the system becomes more isotropic and compact. We find that, in a narrow range of , (or ), the system crystallizes at low temperatures. Some of the samples in our simulation runs also crystallized at the lowest at and , when we carried out simulation for very long time (). We show results for the crystalline state by open symbols in Figure 2 (f). We confirm by carefully inspecting the real space snapshots that the systems crystallize completely over the entire simulation box. The coordination numbers for crystallized samples are and . This suggests that the observed crystalline structure is Stishovite type. Indeed, it has been reported that forms the Stishovite crystal under very high pressures [64]. We remark that an exceptionally large value of at () is due to the ambiguity to define the first coordination shell. As shown in Figure 2 (e), the profile of is broadened around , which makes it difficult to define clearly. Thus, can be regarded as the crossover regime from the network dominated region to the isotropic/compact structure region.

Next, we calculate the static structure factor . We show at the lowest temperatures in Figure 3 for several ’s. At , the maximal peaks are observed around in all components of . These peaks originate from the tetrahedral structure whose bond length is estimated to be . We also find weaker peaks at smaller wave numbers around . This is the so-called first sharp diffraction peak (FSDP) [65], which is known as an indication of the medium-range order of the length scale larger than the neighbor shell. In particular, the FSDP has been attributed to the formation of the hierarchical clusters of the tetrahedra with larger lengths [66]. When the potential depth is reduced, these two peaks merge to a single peak. This result implies that the hierarchical tetrahedral structures disappear gradually as is decreased. At , the shape of the static structure factors are analogous to those observed in typical fragile glass formers such as the SS and LJ binary mixtures [67, 68].

0 | 0 | 0.07-0.34 | 5.8 | 0.18 | 0.0548 | 0.0547 | 0.085 | 0.372 | 0.475 | |
---|---|---|---|---|---|---|---|---|---|---|

0.042 | 0.204 | 5.733 | 0.06-0.34 | 5.9 | 0.16 | 0.0467 | 0.0465 | 0.065 | 0.357 | 0.459 |

0.167 | 0.408 | 4.054 | 0.057-0.34 | 5.9 | 0.12 | 0.0472 | 0.0461 | Non | 0.652 | 0.759 |

1.042 | 1.021 | 2.564 | 0.065-0.34 | 5.6 | 0.14 | 0.046 | 0.0456 | 0.075 | 0.241 | 0.286 |

1.5 | 1.225 | 2.341 | 0.085-0.34 | 5.5 | 0.18 | 0.0599 | 0.059 | 0.11 | 0.2 | 0.229 |

2.667 | 1.633 | 2.027 | 0.14-0.6 | 5.2 | 0.26 | 0.0776 | 0.0793 | 0.2 | 0.1 | 0.126 |

6 | 2.449 | 1.655 | 0.3-1.25 | 5.0 | 0.48 | 0.161 | 0.173 | 0.4 | 0.087 | 0.123 |

### 3.2 Dynamical properties

As demonstrated above, the tetrahedral network structures are broken when the potential depth is reduced to zero. In this section, we analyze the dynamical properties of the CP model. Specifically, the dependence of the fragility is quantified.

First, we calculate the self part of the intermediate scattering function for species 1 defined by

(7) |

where denotes the position of the -th particle at time . In Figure 4, we show the temperature dependence of for several ’s. The wave number is chosen at the peak position of the static structure factor, (see Figure 3). The values of are listed in Table 1. At high temperatures, shows the exponential decay with the short relaxation time but, at low temperatures, dynamics dramatically slows down and exhibits a two-step relaxation. This is the sign of the onset of the glassy dynamics.

We define the relaxation time by . We fit the observed by the Vogel-Fulcher-Tammann (VFT) equation;

(8) |

where and are fitting parameters. The parameter is referred to as the fragility index which is regularly used to quantify the degree of the super-Arrhenius temperature dependence [30]. The larger (smaller) values of correspond to the fragile (strong) glass formers. In this study, the VFT fitting is applied below , where the two-step relaxation for sets in. for each is presented in Table 1. We show the temperature dependence of for several ’s in Figure 6 (a). The temperature is scaled by . It is clearly seen that the results for and follow the Arrhenius law at low temperatures [29]. For smaller ’s, the temperature dependence of deviates from the Arrhenius behavior. In other words, the system changes from strong to fragile glass formers with decreasing .

We also quantify the self diffusion constant () from the long-time behavior of the mean squared displacement (MSD). The MSD for species 1 is presented in Figure 5. Figure 6 (b) shows as a function of .

The dependence of the fragility index is demonstrated in Figure 7. At , the value of is small (), which is consistent with the results of the previous study [29]. As shown in Figure 7, increases with decreasing . At (), exceeds the value of the Kob-Andersen (KA) mixture, which is a typical fragile glass former () [40]. We address that the observed fragility index covers a wide range from to . Especially, the maximal value of at is comparable with those of the most fragile glass formers studied in simulations [69, 14]. is also obtained by using the VFT equation for . The results are plotted with filled circles in Figure 7. Although obtained from are larger than those obtained from for all ’s, the overall behavior is qualitatively the same. We remark that the increase of the fragility index with increasing the density has been reported for the BKS model [28], although the investigated densities were limited [70].

The observation that the system becomes the most fragile ( becomes the largest) just next to the crystalline state may be related to the frustration scenario which claims that the fragility is controlled by the frustration against crystallization [71]. In fact, it has been reported that the system with a weaker frustration against crystallization tends to have larger fragility indices [36, 43]. Further investigation of the dynamics around this Stishovite-like crystallization regime would be valuable to verify the scenario.

In order to characterize the difference of dynamics between the strong (large ) and fragile (small ) regimes, we evaluate the ratio of the diffusion constants between species 1 and 2. In the experiments for , it has been reported that diffusion of silicon and oxygen atoms decouples at low temperatures and the ratio of the diffusion constants for Si and O, , increases with decreasing the temperature [72]. Similar behavior has been also demonstrated in simulations [73, 74]. As mentioned in Reference [74], this decoupling is attributed to the rotational motion of the oxygen atom in the tetrahedral structure. In Figure 8, we show as a function of for several ’s. At , increases with increasing (with decreasing the temperature). This is consistent with the results for the original CP model [29]. As decreases, variation of becomes milder and eventually becomes almost a flat line at . As decreases further below , the trend is reversed and becomes an increasing function of again and the slope keeps increasing until reaches , where the slope becomes maximum. The degree of the decoupling at , which corresponds to the SS model, is even larger than that of the network glass former at . Similar strong decoupling between and has been reported for the conventional SS model [75]. This non-monotonic dependence of is another signature that the dynamics of the CP model is changed qualitatively by tuning the potential depth . It is natural to expect that the underlying mechanisms behind the strong decoupling of and at the two extreme ends of and are completely different. It would be worthwhile to seek for the origin of this strong decoupling for .

### 3.3 Absence of the density-temperature scaling

The density-temperature (DT) scaling is an useful way to single out the parameter which controls the thermodynamic and dynamical properties of liquids. The simplest example is the inverse power law (IPL) potential, , where the DT scaling exactly holds and the system is characterized by a single parameter, , where is the spatial dimension [63]. It has been known that a broad class of liquids can be scaled by a single parameter , where is a constant [50, 52, 53, 54, 55]. It has been argued that the DT scaling holds if there is a strong correlation between the virial and the potential energy . The liquids for which the correlation between and is more than , are called the strongly correlating or Roskilde-simple liquids [52, 53, 54, 55]. Since the CP model with is the purely IPL system, the DT scaling exactly holds. On the other hand, the previous study has demonstrated that at (non-IPL system), the correlation coefficient between and is less than (non-strongly correlating liquid), thus the DT scaling does not hold [76]. Consequently, the validity of the DT scaling for the CP model varies depending on .

There is another type of the DT scaling which is convenient to characterize dynamical quantities such as the relaxation time . It has been argued that for many liquids is scaled by the characteristic time and energy at high temperatures, and . and are determined by fitting at high temperatures using the Arrhenius law;

(9) |

This scaling is demonstrated to work well also in many supercooled liquids and polymer systems [49, 50, 51].

As mentioned in Section 1, the density or pressure dependence of the fragility has been studied in various systems [30, 31, 32, 33]. However, some of those results can be collapsed onto a master curve using the DT scaling, which means that the observed variation of the fragility is only superficial with no physical significance [31, 32].

Here, we examine the DT scaling using Equation (9), for the relaxation time of the CP model. In Figure 9, scaled by is plotted as a function of . In this figure, the potential depth is used instead of the density . At the small regime, the relaxation times are collapsed onto a master curve, whereas they systematically deviate from the master curve as increases. This result eloquently demonstrates that the dynamical properties of the CP model do not follow the DT scaling. From these results, it is concluded that the observed variation of the fragility is not superficial but it is a genuine manifestation of the changeover of the physical mechanism.

### 3.4 Stokes-Einstein violation and stretch exponent

The advantage of the CP model is that the correlation of physical quantities with the fragility can be studied systematically in a single system over a wide range of fragility. In particular, it enables one to examine the much debated issues on the relation of the fragility with dynamical properties, such as the SE violation and the stretched exponential relaxation. In normal liquids, it is expected that the SE relation, const. holds, where and are the diffusion constant and viscosity, respectively. However, the SE relation is violated in most supercooled liquids near the glass transition temperature. This SE violation is often regarded as a manifestation of the spatially heterogeneous dynamics, or the dynamical heterogeneities [56]. Previously, the SE violation of the original CP model () has been investigated [59, 60]. These studies highlight the qualitative difference between the CP model () and other fragile glass formers in their dynamical behavior [77]. The aim of this section is to elucidate the dependence of the SE violation on the fragility in a systematic way.

In Figure 10 (a), we show the Stokes-Einstein ratio, , normalized by the values at , as a function of . For all ’s, the SE relation is violated, i.e., the SE ratio increases with increasing (or decreasing the temperature). We remark that the deviation of the SE ratio at high temperatures (small ) is an artifact caused by the use of instead of [78, 79, 59]. Dependence of the SE violation on the fragility is shown in Figure 10 (b). Here we plot the SE ratio at low temperatures (at ) against the fragility index . A clear correlation between the SE violation and the fragility is observed, that is, the more fragile systems tend to exhibit the stronger SE violation.

Next, we study the stretch exponent of the nonexponential relaxation observed in the density correlation function. Experimental studies have shown that the fragile systems tend to have smaller stretch exponents than the strong systems [6, 7]. Some theories explain this observation [80], while others argue that there is no direct correlation between the fragility and [81, 82].

Here, we show the fragility dependence of the stretch exponents of the CP model. We determine from using the following fitting function,

(10) |

where , and are fitting parameters [78]. Figure 11 (a) shows as a function of . For all ’s, decreases gradually with increasing (decreasing the temperature). In Figure 11 (b), we plot at as a function of . Surprisingly, we observe no clear correlation between the fragility and , although becomes somewhat larger at strong-liquid regime of and . This observation is not inconsistent with the argument in Reference [82]. A similar trend has been reported for the harmonic sphere model [33].

### 3.5 Relationship with specific heat peak

Finally, we discuss the possible relationship between the thermodynamic properties and the fragility variation observed in the CP model. Several simulation studies using the BKS model [28] for have demonstrated that there exists the so-called fragile-to-strong (FS) crossover; the temperature dependence of the relaxation time changes from the super-Arrhenius (fragile) to the Arrhenius (strong) behavior with decreasing the temperature [73, 83, 84]. It has been argued that this FS crossover is related with a thermodynamic anomaly signaled by appearance of the peak of the specific heat [83, 84, 85] which tends to shift to lower temperatures if the density is increased [83, 84, 86]. This anomaly might be connected to a hidden thermodynamic singularity such as the liquid-liquid transition. But the overall picture is still unclear. Under these circumstances, it is of great interest to examine whether the FS crossover and the shift of the specific heat are also observed in the CP model by tuning the potential depth .

The specific heat is calculated by . Figure 12 shows the temperature dependence of for various potential depth . The temperature is normalized by the onset temperature . As observed, there exists the broad but clear peak of at a finite temperature for all ’s except for . We cannot access the lower temperatures for because the system crystallizes. The temperature for each is listed in Table 1. We find that shifts to lower temperatures with decreasing (increasing the density).

Contrary to the distinct peaks of the specific heat, the FS crossover is harder to detect in the simulation data. Due to the limited range of the temperature in the Arrhenius plot of Figure 6, one hardly observe any sort of distinct crossover. In this figure, we marked the temperatures at which the peak of the specific heat are observed with the empty symbols for the guide of eyes. For the data of and 6, it is not impossible to fit the data of with the Arrhenius and super-Arrhenius laws on the low and high temperature sides across , but the fitting range is too narrow to call it convincing. Similar fitting was also possible for smaller values of ’s, but it is less trustworthy. Therefore, it would be fair to say that the existence of the correlation between the FS crossover and the peak of the specific heat is not conclusive in current numerical simulations [83, 84, 86]. More extensive simulations are required to clarify the relation between the FS crossover and the specific heat peak.

## 4 Summary and Discussion

We have numerically investigated the CP model by tuning the depth of the potential . Changing corresponds to changing the density of the system. This approach has enabled us to perform simulations over a wide range of density from the order of unity up to infinity. From the radial distribution functions, the coordination numbers, and the static structure factors, we found that the anisotropic tetrahedral network structures are broken and then the isotropic structures are formed by decreasing (or increasing the density). We have also calculated various dynamical quantities such as the relaxation time and diffusion constant. These calculations have revealed that their temperature dependence seamlessly changes from the Arrhenius to the super-Arrhenius behavior. We also confirmed that the temperature dependence of the dynamics is not collapsed by the DT scaling, assuring that the observed fragility change is not a consequence of the trivial density scaling but due to the generic change of the mechanism of the glassy dynamics. We have studied the relationship between the fragility and two dynamical quantities, the magnitude of the SE violation and the stretch exponent of the density correlation function. The magnitude of the SE violation, which is believed to be a manifestation of the dynamical heterogeneities, correlates well with the fragility of the CP model. On the contrary, the clear correlation between the fragility and the stretch exponent has not been observed. Finally, the peak of the specific heat and the possible correlation with the fragile-to-strong (FS) crossover have been argued. The peak has been observed in the temperature dependence of the specific heat for most ’s in the CP model. In addition, its peak position shifts to lower temperatures with decreasing the potential depth (or increasing the density). These observations are related with the recent numerical results obtained in the BKS model [83, 84, 86]. However, it was difficult to identify the FS crossover temperature from the relaxation times because the simulation time windows are too narrow and we conclude that the correlation between the FS crossover and thermodynamic anomaly is still elusive. Further numerical investigation is required.

We address that the CP model exhibits a very broad variation of the fragility, the highest and lowest values of which are comparable to those of the most fragile and strongest liquids studied numerically in the past, and that it can serve as an ideal bench for numerical study of the fragility. Note that the harmonic sphere model also shows a very wide variation of the fragility as the density is varied [33, 87]. This model also does not satisfy the DT scaling. However, it should be emphasized that this model is qualitatively different from the CP model. In the harmonic sphere model, the potential is truncated at a certain cutoff distance. At low densities, the model becomes effectively a hard sphere system. Therefore, the temperature and energy are not relevant but, instead, the pressure or density become the controlling thermodynamic parameters. This means that at low densities the temperature dependence of the relaxation time is trivially Arrhenius, whereas at higher densities, the temperature plays a pivotal role, leading to more diverging (or fragile) behavior. For the CP model, on the contrary, the temperature plays a role as a relevant control parameter for all densities and the fragility is continuously controlled.

There are other physical quantities that are expected to be related with the fragility [2]; the excess entropy [3], the elastic constants [4, 5], the non-ergodicity parameter [88], the boson peak frequency [5], and the temperature dependence of microscopic structure [89]. In addition, the relationship between the fragility and the dynamical heterogeneities is a topic which attracts much attention recently [90, 91, 58]. The comprehensive information on the correlation of these properties with the fragility will be obtained from the systematic study of the CP model by tuning a single parameter . Further investigations along this line are currently underway.

## Acknowledgments

We thank D. Coslovich, A. Ikeda, H. Ikeda, and T. Kawasaki for helpful discussions. M. O. acknowledges the financial support by Grant-in-Aid for Japan Society for the Promotion of Science Fellows (26.1878). This work is partially supported by KAKENHI Grants No. 24340098 (K.M.), 25103005 “Fluctuation & Structure” (K.M.), 25000002 (K.M.), 26400428 (K.K.), 16H00829 “Soft Molecular Systems” (K.K.), and 26400428 (K.K.). K. M. also thanks the JSPS Core-to-Core program. The numerical calculations have been performed at Research Center for Computational Science, Okazaki, Japan.

## References

### References

- Angell C A Formation of glasses from liquids and biopolymers 1995 Science 267, 1924
- A Lindsay G, Kenneth F K and Srikanth S 2014 Fragility of Glass-Forming Liquids (Hindustan Book Agency)
- Martinez L M and Angell C A thermodynamic connection to the fragility of glass-forming liquids 2001 Nature 410, 663
- Novikov V and Sokolov A Poisson’s ratio and the fragility of glass-forming liquids 2004 Nature 431, 961
- Novikov V, Ding Y and Sokolov A Correlation of fragility of supercooled liquids with elastic properties of glasses 2005 Phys. Rev. E 71, 061501
- Böhmer R, Ngai K, Angell C and Plazek D Nonexponential relaxations in strong and fragile glass formers 1993 J. Chem. Phys. 99, 4201
- Niss K, Dalle-Ferrier C, Tarjus G and Alba-Simionesco C On the correlation between fragility and stretching in glass-forming liquids 2007 J. Phys.: Condens. Matter 19, 076102
- Adam G and Gibbs J H On the temperature dependence of cooperative relaxation properties in glass-forming liquids 1965 J. Chem. Phys. 43, 139
- Xia X and Wolynes P G Fragilities of liquids predicted from the random first order transition theory of glasses 2000 Proc. Natl. Acad. Sci. U.S.A. 97, 2990
- Debenedetti P G and Stillinger F H Supercooled liquids and the glass transition 2001 Nature 410, 259
- Tarjus G, Kivelson D and Viot P The viscous slowing down of supercooled liquids as a temperature-controlled super-Arrhenius activated process: a description in terms of frustration-limited domains 2000 J. Phys.: Condens. Matter 12, 6497
- Coslovich D and Pastore G Understanding fragility in supercooled Lennard-Jones mixtures. I. Locally preferred structures 2007 J. Chem. Phys. 127, 124504
- Tanaka H Bond orientational order in liquids: Towards a unified description of water-like anomalies, liquid-liquid transition, glass transition, and crystallization 2012 Eur. Phys. J. E 35, 1
- Royall C P, Malins A, Dunleavy A J and Pinney R Strong geometric frustration in model glassformers 2015 J. Non-Cryst. Solids 407, 34
- Dyre J C Colloquium: The glass transition and elastic models of glass-forming liquids 2006 Rev. Mod. Phys. 78, 953
- Krausser J, Samwer K H and Zaccone A Interatomic repulsion softness directly controls the fragility of supercooled metallic melts 2015 Proc. Natl. Acad. Sci. U.S.A. 112, 13762
- Yan L, Düring G and Wyart M Why glass elasticity affects the thermodynamics and fragility of supercooled liquids 2013 Proc. Natl. Acad. Sci. U.S.A. 110, 6307
- Garrahan J P and Chandler D Coarse-grained microscopic model of glass formers 2003 Proc. Natl. Acad. Sci. U.S.A. 100, 9710
- Cavagna A Supercooled liquids for pedestrians 2009 Phys. Rep. 476, 51
- Berthier L and Biroli G Theoretical perspective on the glass transition and amorphous materials 2011 Rev. Mod. Phys. 83, 587
- Binder K and Kob W 2011 Glassy Materials and Disordered Solids: An Introduction to their Statistical Mechanics (World Scientific)
- Royall C P and Williams S R The role of local structure in dynamical arrest 2015 Phys. Rep. 560, 1
- Bernu B, Hansen J, Hiwatari Y and Pastore G Soft-sphere model for the glass transition in binary alloys: Pair structure and self-diffusion 1987 Phys. Rev. A 36, 4891
- Wahnström G Molecular-dynamics study of a supercooled two-component Lennard-Jones system 1991 Phys. Rev. A 44, 3752
- Kob W and Andersen H C Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function 1995 Phys. Rev. E 51, 4626
- Woodcock L, Angell C and Cheeseman P Molecular dynamics studies of the vitreous state: Simple ionic systems and silica 1976 J. Chem. Phys. 65, 1565
- Tsuneyuki S, Tsukada M, Aoki H and Matsui Y First-principles interatomic potential of silica applied to molecular dynamics 1988 Phys. Rev. Lett. 61, 869
- van Beest B, Kramer G J and Van Santen R Force fields for silicas and aluminophosphates based on ab initio calculations 1990 Phys. Rev. Lett. 64, 1955
- Coslovich D and Pastore G Dynamics and energy landscape in a tetrahedral network glass-former: direct comparison with models of fragile liquids 2009 J. Phys.: Condens. Matter 21, 285107
- Sastry S The relationship between fragility, configurational entropy and the potential energy landscape of glass-forming liquids 2001 Nature 409, 164
- De Michele C, Sciortino F and Coniglio A Scaling in soft spheres: fragility invariance on the repulsive potential softness 2004 J. Phys.: Condens. Matter 16, L489
- Sengupta S, Schrøder T B and Sastry S Density-temperature scaling of the fragility in a model glass-former 2013 Eur. Phys. J.. E 36, 1
- Berthier L and Witten T A Compressing nearly hard sphere fluids increases glass fragility 2009 Europhys. Lett. 86, 10001
- Wang L, Duan Y and Xu N Non-monotonic pressure dependence of the dynamics of soft glass-formers at high compressions 2012 Soft Matter 8, 11831
- Kawasaki T, Araki T and Tanaka H Correlation between dynamic heterogeneity and medium-range order in two-dimensional glass-forming liquids 2007 Phys. Rev. Lett. 99, 215701
- Kawasaki T and Tanaka H Structural origin of dynamic heterogeneity in three-dimensional colloidal glass formers and its link to crystal nucleation 2010 J. Phys.: Condens. Matter 22, 232102
- Abraham S E, Bhattacharrya S M and Bagchi B Energy landscape, antiplasticization, and polydispersity induced crossover of heterogeneity in supercooled polydisperse liquids 2008 Phys. Rev. Lett. 100, 167801
- Kurita R and Weeks E R Glass transition of two-dimensional binary soft-disk mixtures with large size ratios 2010 Phys. Rev. E 82, 041402
- Bordat P, Affouard F, Descamps M and Ngai K Does the interaction potential determine both the fragility of a liquid and the vibrational properties of its glassy state? 2004 Phys. Rev. Lett. 93, 105502
- Sengupta S, Vasconcelos F, Affouard F and Sastry S Dependence of the fragility of a glass former on the softness of interparticle interactions 2011 J. Chem. Phys. 135, 194503
- Berthier L and Tarjus G Nonperturbative effect of attractive forces in viscous liquids 2009 Phys. Rev. Lett. 103, 170601
- Coslovich D Locally preferred structures and many-body static correlations in viscous liquids 2011 Phys. Rev. E 83, 051505
- Shintani H and Tanaka H Frustration on the way to crystallization in glass 2006 Nat. Phys. 2, 200
- Molinero V, Sastry S and Angell C A Tuning of tetrahedrality in a silicon potential yields a series of monatomic (metal-like) glass formers of very high fragility 2006 Phys. Rev. Lett. 97, 075701
- Parmar A D and Sastry S Kinetic and Thermodynamic Fragilities of Square Well Fluids with Tunable Barriers to Bond Breaking 2015 J. Phys. Chem. B 119, 11243
- Kim K, Miyazaki K and Saito S Slow dynamics, dynamic heterogeneities, and fragility of supercooled liquids confined in random media 2011 J. Phys.: Condens. Matter 23, 234123
- Chakrabarty S, Karmakar S and Dasgupta C Dynamics of glass forming liquids with randomly pinned particles 2015 Scientific reports 5, 12577
- Sausset F, Tarjus G and Viot P Tuning the fragility of a glass-forming liquid by curving space 2008 Phys. Rev. Lett. 101, 155701
- Alba-Simionesco C, Kivelson D and Tarjus G Temperature, density, and pressure dependence of relaxation times in supercooled liquids 2002 J. Chem. Phys. 116, 5033
- Casalini R and Roland C Thermodynamical scaling of the glass transition dynamics 2004 Phys. Rev. E 69, 062501
- Alba-Simionesco C and Tarjus G Temperature versus density effects in glassforming liquids and polymers: A scaling hypothesis and its consequences 2006 J. Non-Cryst. Solids 352, 4888
- Pedersen U R, Bailey N P, Schrøder T B and Dyre J C Strong pressure-energy correlations in van der Waals liquids 2008 Phys. Rev. Lett. 100, 015701
- Schrøder T B, Bailey N P, Pedersen U R, Gnan N and Dyre J C Pressure-energy correlations in liquids. III. Statistical mechanics and thermodynamics of liquids with hidden scale invariance 2009 J. Chem. Phys. 131, 234503
- Gnan N, Schrøder T B, Pedersen U R, Bailey N P and Dyre J C Pressure-energy correlations in liquids. IV.“Isomorphs” in liquid phase diagrams 2009 J. Chem. Phys. 131, 234504
- Schrøder T B, Gnan N, Pedersen U R, Bailey N P and Dyre J C Pressure-energy correlations in liquids. V. Isomorphs in generalized Lennard-Jones systems 2011 J. Chem. Phys. 134, 164505
- Ediger M D Spatially heterogeneous dynamics in supercooled liquids 2000 Annu. Rev. Phys. Chem. 51, 99
- Berthier L, Biroli G, Coslovich D, Kob W and Toninelli C Finite-size effects in the dynamics of glass-forming liquids 2012 Phys. Rev. E 86, 031502
- Kim K and Saito S Multiple length and time scales of dynamic heterogeneities in model glass-forming liquids: A systematic analysis of multi-point and multi-time correlations 2013 J. Chem. Phys. 138, 12A506
- Kawasaki T, Kim K and Onuki A Dynamics in a tetrahedral network glassformer: Vibrations, network rearrangements, and diffusion 2014 J. Chem. Phys. 140, 184502
- Staley H, Flenner E and Szamel G Reduced strength and extent of dynamic heterogeneity in a strong glass former as compared to fragile glass formers 2015 J. Chem. Phys. 143, 244501
- Zaccarelli E, Sciortino F and Tartaglia P A spherical model with directional interactions. I. Static properties 2007 J. Chem. Phys. 127, 174501
- Mayer C, Sciortino F, Tartaglia P and Zaccarelli E A spherical model with directional interactions: II. Dynamics and landscape properties 2010 J. Phys.: Condens. Matter 22, 104110
- Hiwatari Y, Matsuda H, Ogawa T, Ogita N and Ueda A Molecular dynamics studies on the soft-core model 1974 Progr. Theor. Phys. 52, 1105
- Keskar N R and Chelikowsky J R Structural properties of nine silica polymorphs 1992 Phys. Rev. B 46, 1
- Elliott S Origin of the first sharp diffraction peak in the structure factor of covalent glasses 1991 Phys. Rev. Lett. 67, 711
- Nakamura T, Hiraoka Y, Hirata A, Escolar G E, Matsue K and Nishiura Y Description of Medium-Range Order in Amorphous Structures by Persistent Homology 2015 arXiv:1501.03611
- Kob W and Andersen H C Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture. II. Intermediate scattering function and dynamic susceptibility 1995 Phys. Rev. E 52, 4134
- Kob W, Roldan-Vargas S and Berthier L Spatial correlations in glass-forming liquids across the mode-coupling crossover 2012 Phys. Procedia 34, 70
- Malins A, Eggers J, Royall C P, Williams S R and Tanaka H Identification of long-lived clusters and their link to slow dynamics in a model glass former 2013 J. Chem. Phys. 138, 12A535
- Barrat J L, Badro J and Gillet P A strong to fragile transition in a model of liquid silica 1997 Mol. Simul. 20, 17
- Tanaka H Roles of bond orientational ordering in glass transition and crystallization 2011 J. Phys.: Condens. Matter 23, 284115
- Mikkelsen Jr J Self-diffusivity of network oxygen in vitreous SiO2 1984 Appl. Phys. Lett. 45, 1187
- Horbach J and Kob W Static and dynamic properties of a viscous silica melt 1999 Phys. Rev. B 60, 3169
- Saksaengwijit A and Heuer A Origin of the decoupling of oxygen and silicon dynamics in liquid silica as expressed by its potential energy landscape 2006 Phys. Rev. E 74, 051502
- Kawasaki T and Onuki A Slow relaxations and stringlike jump motions in fragile glass-forming liquids: Breakdown of the Stokes-Einstein relation 2013 Phys. Rev. E 87, 012312
- Coslovich D and Roland C Heterogeneous slow dynamics and the interaction potential of glass-forming liquids 2011 J. Non-Cryst. Solids 357, 397
- Flenner E, Staley H and Szamel G Universal features of dynamic heterogeneity in supercooled liquids 2014 Phys. Rev. Lett. 112, 097801
- Sengupta S, Karmakar S, Dasgupta C and Sastry S Breakdown of the Stokes-Einstein relation in two, three, and four dimensions 2013 J. Chem. Phys. 138, 12A548
- Shi Z, Debenedetti P G and Stillinger F H Relaxation processes in liquids: Variations on a theme by Stokes and Einstein 2013 J. Chem. Phys. 138, 12A526
- Xia X and Wolynes P G Microscopic theory of heterogeneity and nonexponential relaxations in supercooled liquids 2001 Phys. Rev. Lett. 86, 5526
- Dyre J C Ten themes of viscous liquid dynamics 2007 J. Phys.: Condens. Matter 19, 205105
- Heuer A Exploring the potential energy landscape of glass-forming systems: from inherent structures via metabasins to macroscopic transport 2008 J. Phys.: Condens. Matter 20, 373101
- Saika-Voivod I, Poole P H and Sciortino F Fragile-to-strong transition and polyamorphism in the energy landscape of liquid silica 2001 Nature 412, 514
- Saika-Voivod I, Sciortino F and Poole P H Free energy and configurational entropy of liquid silica: fragile-to-strong crossover and polyamorphism 2004 Phys. Rev. E 69, 041503
- Speck T, Royall C P and Williams S R Liquid-liquid phase transition in an atomistic model glass former 2014 arXiv:1409.0751
- Saika-Voivod I, Sciortino F, Grande T and Poole P H Phase diagram of silica from computer simulation 2004 Phys. Rev. E 70, 061507
- Berthier L and Tarjus G The role of attractive forces in viscous liquids 2011 J. Chem. Phys. 134, 214503
- Scopigno T, Ruocco G, Sette F and Monaco G Is the fragility of a liquid embedded in the properties of its glass? 2003 Science 302, 849
- Mauro N, Blodgett M, Johnson M, Vogt A and Kelton K A structural signature of liquid fragility 2014 Nat. Commun. 5, 4616
- Qiu X and Ediger M Length scale of dynamic heterogeneity in supercooled D-sorbitol: comparison to model predictions 2003 J. Phys. Chem. B 107, 459
- Berthier L, Biroli G, Bouchaud J P, Cipelletti L, El Masri D, L’Hôte D, Ladieu F and Pierno M Direct experimental evidence of a growing length scale accompanying the glass transition 2005 Science 310, 1797