Relationship between structural order and water-like anomalies in metastable liquid silicon: Ab initio molecular dynamics

The relationship between structural order and water-like anomalies in tetrahedral liquids is still open. Here, first-principle molecular dynamics are performed to study it in metastable liquid Si. It is found that in T-P phase diagram, there indeed exists a structural anomaly region, which encloses density anomaly but not diffusivity anomaly. This is consistent with that of SW Si and BKS SiO2 but different from that of SPC/E water. Two-body excess entropy anomaly can neither capture the diffusivity, structural, and density anomalies, as it can in a two-scale potential fluid. In structural anomaly region, tetrahedrality order qtetra (measuring the extent to which an atom and its four nearest neighbours adopt tetrahedral arrangement) and translational order ttrans (measuring the tendency of two atoms to adopt preferential separation) are not perfectly correlated, which is different from that in SW Si and renders it impossible to use the isotaxis line to quantify the degree of structural order needed for water-like anomalies to occur. Along the isotherm of critical temperature Tc, ttrans/qtetra is approximately linear with pressure. With decreasing pressure along the isotherm below Tc, ttrans/qtetra departs downward from the line, while it is the opposite case above Tc.

can not anticipate the occurrence of the diffusivity anomaly and hence can not provide a satisfactory microscopic picture of the anomalous behaviour in liquid silica. The later results of liquid BeF 2 with the transferable rigid ion model (TRIM) by Agarwal et al. 22 shows the similar result with BKS SiO 2 , but that of the Oeffner-Elliot (OE) model of GeO 2 by Jabes et al. 23 is similar to SPC/E water. The systematic comparison for different model potentials of water (mTIP3P, TIP4P, TIP4P/2005, TIP5P, SPC/E), ionic liquids (SiO 2 and BeF 2 ), and liquids characterised by the Stillinger-Weber (SW) potential (including Si) [24][25][26][27] suggests that a strong correlation between q tetra and t trans exists only in rigid-body model potentials for water and SW liquids (in a limited range of tetrahedrality strength), but not in ionic melts. In addition, to obtain a more general picture of the origins of these anomalies, the excess entropy has also been attempted to understand these anomalies and has been tested to be able to predict the regions of anomalies in the phase diagram [22][23][24][25][26][27][28][29] .
It should be noted that the water-like cascading of anomalous regions in the phase diagram has been observed not only in systems having directional interactions but also in the simulation studies of spherically symmetric potentials [29][30][31][32][33] . Especially, the study of Jagla model by Yan et al. 30 suggests the water-like relationship between structural order and anomalies is related to the presence of two different length scales in the potential. Moreover, the study of a core-softened model by Fomin et al. 33 shows that the order of the cascading regions of anomaly can be changed by increasing the depth of the attractive part of the potential.
Although so many computer simulation studies of tetrahedral liquids have focused on the relationship between structural order and water-like anomalies, they all depended on different model potentials or empirical potentials. To our knowledge, first-principle computer simulations have not been performed to study this question so far. It is well-known that the first-principle force neither makes assumptions such as empirical model nor includes fitting parameter to experimental data, so it is expected to give more realistic results. In the past few decades, first-principle simulations had achieved great success in investigating structural and dynamic properties of liquids. Based on first-principle calculations, we have studied the inherent structures of high-density liquid (HDL) and low-density liquid (LDL) Si and found that the first-order LLPT in metastable liquid Si is a transition between a sp 3 -hybridization LDL and a white-tin-like HDL 34 . Further, we found that the liquid-liquid critical point (LLCP) is hidden, blended in with the continuous reentrant spinodal of HDL, suggesting that the phase behaviour of metastable liquid Si tends to be a critical-point-free scenario rather than a second-critical-point one based on SW potential 35 . In this work, based on the prior research, we turn our attention to the relationship between structural order and water-like anomalies (including the well-known first-order LLPT) in metastable liquid silicon.

Results
The density anomaly and the first-order LLPT have been shown in our previous work 35 . In the studied pressure range from − 27.25 to 6.50 kbar, an obvious maximum can be observed along each ρ-T isobar line 35 and in T-P phase diagram, the line of density maxima tends to pass through the LLCP, located at T c ~ 1420 ± 10 K and P c ~ − 30.5 ± 1.0 kbar 35 , which is different from the result in ref. 17 (T c ~ 1120 K and P c ~ 0.60 GPa). We speculate that such a difference maybe origins from the different potentials adopted. Although the LLPT is proved to exist for SW liquid Si, it has been reported that it strongly depends on the parameters of the potential and an arbitrary variation, albeit small, of the parameters may even lead to its disappearance 36 . Here, we begin by presenting our results related to the diffusivity anomaly. Diffusivity anomaly. In Fig. 1, we show our calculated diffusion coefficient D as a function of pressure along eight different isotherms ranging from 1032 to 1800 K. Along each isotherm, the pressure is gradually decreased to the stability limit of HDL and an obvious maximum can be found in the whole temperature range studied here, which marks the upper bound of diffusivity anomaly region. Along each isotherm below T c , the locus of diffusivity minimum, which marks the lower bound of diffusivity anomaly region, can not be found and is expected to be located in LDL, where system's dynamics is slow and there is considerable error in computing D and locating the diffusivity minima. So we have not attempted to find this locus in our study. Above T c , the diffusivity minimum is expected to be located at the HDL-V spinodal because of the larger D in vapor than that in liquid. Figure 2 presents the evolution of structural order (t trans and q tetra ) with decreasing pressure along the isotherm. From Fig. 2(a), it can be found that q tetra increases with the decrease of pressure in the simulated temperature range from 1032 to 1800 K, which is in contrast to what is observed for simple fluids and so is termed as structural anomaly. Below T c , we do not find a maximum of q tetra in our simulated pressure range, which marks the lower bound of structural anomaly region and is expected to be found in LDL. At temperatures above T c , the maximum of q tetra is expected to be located at the HDL-V spinodal, such as for T = 1532 and 1800 K. In Fig. 2(b), in the temperature range from 1032 to 1700 K, t trans decreases with decreasing pressure and goes through a minimum, which marks the upper bound of structural anomaly region as suggested by Errington and Debenedetti 6 . For temperature T = 1800 K, the minimum cannot be found until the HDL-V spinodal is encountered.

Structural anomaly.
Excess entropy anomaly. The excess entropy S ex is defined as the entropy of the liquid relative to that of the ideal gas at the same temperature and pressure and can be expressed in terms of an expansion of n-body correlation functions. The first term, namely the two-body contribution S 2 , accounts for a large proportion, such as 80-90% in Lennard-Jones fluids 37 . Because S ex characterizes the reduction in the number of states (relative to an ideal gas) accessible to a system due to translational interparticle correlations, it can also be thought of as a metric for translational order. Errington and co-workers 29 argued that the region in the phase diagram where the excess entropy anomalously increases with decrease in density is also the region where translational order anomaly is observed. Indeed, many studies mentioned above have showed that there also exists a S ex or S 2 anomaly region in the T-ρ phase diagram, where S ex or S 2 anomalously increases with increasing pressure. In this work, we calculated S 2 as a function of pressure along eight different isotherms, shown in Fig. 3. In the temperature range from 1032 to 1432 K, along each isotherm a maximum occurs but the minimum can not be found until the HDL spinodal is encountered, which is expected to be found in LDL. From 1532 to 1800 K, the maximum can not be found until the HDL-V spinodal is encountered.   Figure 2. The average tetrahedrality order q tetra (a) and translational order t trans (b) as a function of pressure for different temperatures from 1032 to 1800 K obtained from ab initio molecular dynamics simulations of liquid Si. q tetra increases with decreasing pressure. From 1032 to 1700 K, t trans decreases with decreasing pressure and goes through a minimum, which marks an onset of structural anomaly.
It is known that there is a strong connection between transport properties of fluids and excess entropy, such as the Rosenfeld scaling relation 38,39 and the Dzugutov scaling relation D z * = a 0 exp(S 2 ) 40 . In the former, the reduced diffusion coefficient D R * = Dρ 1/3 /(k B T/m) 1/2 , m is the mass of the particle, and ρ is the number density. In the later, D z * = DΓ −1 σ −2 , Γ is the collision frequency and equals 4σ 2 g(σ)ρ(πk B T/m) 1/2 , σ is the hard-sphere diameter, and g(σ) is the value of radial distribution function at the constant distance. In Fig. 3 (inset), we also test the scaling relation between D R * and S 2 . The relation between D z * and S 2 is similar to that presented in Fig. 3 (inset) and not shown here. However, it can be seen that the relation between ln(D R *) and S 2 is not linear, showing a different behaviour from that of Rosenfeld and Dzugutov. Based on the results in Figs 1 and 3, it is due to the fact that along an isotherm, the change of diffusivity with pressure is out of step with that of S 2 versus pressure, which is different from the results in ref. 41 based on a core-softened potential (i.e., a Lennard-Jones potential plus a Gaussian repulsion). Thus we can comment that, based on first-principle calculations, the Rosenfeld and Dzugutov scaling relations are invalid along an isotherm in the anomalous region. This is consistent with that observed in supercooled liquid Si based on SW potential 27 and the systems based on some two-scale potentials 42,43 .

Discussion
Based on the above results, now we turn to discuss the relationship of these anomalies in T-P phase diagram, shown in Fig. 4. The spinodals, LLCP, HDL-LDL coexistence line, and density maxima are obtained from our  previous work 35 . The LDL-V and the HDL-V spinodals form a continuous liquid-vapor spinodal, which meets the liquid-vapor critical point at higher temperature and higher pressure (not shown here), and the LLCP seems to just be hidden, blended in with the liquid-vapour spinodal. This is consistent with the critical-point-free scenario suggested by Angell in refs 44 and 45, but is different with the two-critical-point scenario, in which the LLCP occurs above the liquid-vapor spinodal, similar to that presented in Fig. 2 in ref. 46. From Fig. 4, it can be seen that the region of diffusivity anomaly encloses the region of structural anomaly, which in turn encloses the region of density anomaly. This result is consistent with that of SW Si 27 and is similar to that of BKS SiO 2 21 but different from that of SPC/E model of water 6 . In addition, it should be noticed that the two-body excess entropy anomaly region does not enclose the diffusivity anomaly as it does in a fluid interacting through a two-scale potential 29 . On the contrary, it is also included in the structural anomaly region determined by the translational order minima, with the S 2 maxima line going across the density maxima line. So S 2 anomaly can not be used as a criteria of structural anomaly to capture the region of diffusivity and density anomalies. This may origin in the fact that the relation between ln(D R *) and S 2 shows an anomalous behaviour, presented in Fig. 3 (inset).
Recently, it has been shown for the two-scale potential that different criteria of structural anomaly can lead to completely different results 47 . It was found that the structural anomaly region determined by the S ex anomaly lies much deeper than that obtained from the approximate S 2 anomaly. This means that the structural anomaly calculated from the S ex criterion may be inside the density anomaly shown in Fig. 4, which is in contradiction with the thermodynamically consistent relation that the region of anomalous density is always inside the region of anomalous S ex 47 . We think that there are two main reasons accountable for this. The first is that there is a larger error in determining the temperature of density maximum, the higher the pressure, the greater the error (see Fig. 3 in ref. 35). The second is that, although the structural anomaly region determined by S ex lies much deeper than that obtained from S 2 for the two-scale potential 47 , it is not necessarily the case for other systems. For example, in SW Si the onset pressure of S ex anomaly is the same with that of S 2 anomaly and at lower temperatures the former is even slightly larger than the latter (see Fig. 12(a) in ref. 27).
As mentioned above, in the structural anomaly region, a strong correlation between q tetra and t trans exists in rigid-body model potentials for water and SW liquids (in a limited range of tetrahedrality strength), but not in ionic melts. Here, we also present the correlation between q tetra and t trans , shown in Fig. 5. Similar to water, our Si order parameter map also suggests an inaccessible region in the high-q tetra , low-t trans quadrant. In other words, for any given value of t trans , there is a maximum value of q tetra that the system can attain. Likewise, for any given value of q tetra , there is a minimum attainable t trans . However, in the structural anomaly region, q tetra and t trans do not collapse onto a single line, as is the case for water, suggesting that q tetra and t trans are not perfectly correlated. With the decrease of temperature, the correlation becomes weaker and weaker. These results are different from those observed in SW Si [25][26][27] , in which q tetra and t trans are perfectly correlated in the structural anomaly region. This degree of independence renders it impossible to use the isotaxis line to quantify the degree of structural order needed for the water-like anomalies to occur.
To obtain the relationship between structural order and the LLPT, in Fig. 5 (inset), we also present the value of t trans /q tetra as a function of pressure. It can be found that at higher pressure above 20 kbar, t trans /q tetra is independent of temperature and is linear with pressure. Along the isotherm of T c , t trans /q tetra can be approximately expressed as a linear relationship with pressure, t trans /q tetra = 0.0015P + 0.494. With decreasing pressure along the isotherm below T c , t trans /q tetra departs downward from the straight line until the LLPT is encountered, while it is the opposite case above T c until the liquid-vapor transition occurs. Figure 5. Parametric plot of translational order t trans against tetrahedrality order q tetra for different temperatures from 1032 to 1800 K obtained from ab initio molecular dynamics simulations of liquid Si. q tetra and t trans are not perfectly correlated in the structural anomaly region and with decreasing temperature, the correlation becomes weaker and weaker. Inset: The value of t trans /q tetra against pressure for different isotherms. Along the isotherm of T c , t trans /q tetra is approximately linear with pressure, t trans /q tetra = 0.0015P + 0.494. With decreasing pressure along the isotherm below T c , t trans /q tetra departs downward from the line until the LLPT occurs, while it is the opposite case above T c .
Scientific RepoRts | 7:39952 | DOI: 10.1038/srep39952 In summary, first-principle molecular dynamics are performed for the first time to study the relationship between structural order and water-like anomalies in metastable liquid Si. Our results show that in T-P phase diagram, there indeed exists a structural anomaly region, which encloses the density anomaly but not the diffusivity anomaly. This is consistent with that of SW Si and BKS SiO 2 but different from that of SPC/E model of water. Two-body excess entropy anomaly can neither capture the diffusivity, structural, and density anomalies, as it can in a fluid interacting through a two-scale potential. Moreover, in the structural anomaly region, q tetra and t trans are not perfectly correlated, which is different from that in SW Si and renders it impossible to use the isotaxis line to quantify the degree of structural order needed for water-like anomalies to occur. Along the isotherm of critical temperature T c , t trans /q tetra is approximately linear with pressure, t trans /q tetra = 0.0015P + 0.494. With decreasing pressure along the isotherm below T c , t trans /q tetra departs downward from the line until the liquid-liquid phase transition occurs, while it is the opposite case above T c .

Methods
Ab initio molecular dynamics. The simulations are performed by the Vienna ab initio simulation package (VASP) 48,49 together with projector augmented wave (PAW) potential 50,51 in the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) 52 . A system of 216 Si atoms in a cubic supercell with periodical boundary conditions is used. The electronic wave functions are expanded in the plane wave basis set with an energy cutoff of 245 eV and a Pulay stress 53 of 1.65 kbar is added to offset the incompleteness of the plane wave basis set. Only Γ -point is used to sample the supercell Brillouin zone. The canonical (NVT) ensemble simulations are performed with a Nosé thermostat for temperature control 54 . Newton's equations of motion are integrated using Verlet's algorithm in the velocity form with a time step of 2 fs.
The initial configuration, a diamond cubic crystal structure with the volume of 17.2 Å 3 /atom, is heated up to 2000 K. After a run of 20 ps with constant temperature 2000 K, the system arrives at an equilibrium liquid state. The temperature is then gradually reduced to the desired temperature with a cooling rate of 0.5 × 10 15 K/s and a subsequent run of 15-20 ps is performed at this temperature. Then, the simulations are performed with a series of successively larger volumes along the isotherm up to the stability limit of HDL. At each simulation point, we performed a simulation about 15-30 ps. It is found that for HDL 5 ps is enough to equilibrate and the statistical averages are collected for the rest of the simulation. Along each isotherm, the end point of the previous simulation is scaled to obtain the initial conditions for the next.
Calculation of self-diffusion coefficient D. The self-diffusion coefficient D in liquids was obtained by fitting the long time mean square displacement (MSD). MSD is defined as where ξ = rρ 1/3 is the distance r between two Si atoms divided by the mean separation between a pair of atoms at the given number density ρ; g(ξ) is the pair correlation function, and ξ c is a cutoff distance beyond which the system's pair correlation function cannot be distinguished from its asymptotic value of 1. In this work, ξ c is chosen to be 9.0 × ρ 1/3 . Scaled coordinates are used so that the above integral sums over an equivalent number of coordinate shells at each density. For an ideal gas, since g(r) = 1.0, t trans is zero and for a crystal, it has a finite value depending on the cutoff distance. In the liquid phase, t trans will have a value in between that of the ideal gas and a crystal. The tetrahedrality order parameter q tetra is defined as ∑ ∑ where ψ jk is the angle formed by the lines joining a reference atom i and its nearest neighbours j and k. The average q tetra varies between 0 (in the case of an ideal gas) and 1 (in the case of a cubic diamond crystal).
Excess entropy. Excess entropy S E is defined as S E = S − S id , where S is the total entropy of the system and S id is the entropy of an ideal gas system. In present work, we calculate the two body approximation (S 2 ) to S E from the pair correlation function, and is given by where g(r) is the pair correlation function, and ρ is the number density.