{'fraction_non_alphanumeric': 0.03172670629520877, 'fraction_numerical': 0.016418391463152618, 'mean_word_length': 4.948769837043349, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Liquid-Liquid Phase Transition in Supercooled Silicon

We present a review on the study of metastable silicon, primarily focusing mainly on the aspects of liquid-liquid transition, critical point and phase behaviour, structural and dynamic properties of liquid phase as well as crystal nucleation. We begin with an extensive survey of the investigations of liquid silicon pursued over three decades, with salient experimental, theoretical and simulation results. Following which we present various scenarios put forward to rationalize the density and related anomalies often observed in water and other network forming liquids. After which we present the more recent investigations (both simulation and experimental works) of the phase behavior of Silicon. Since a significant part of metastable silicon work is on a classical empirical potential an important question to address is the reliability of this potential in describing the behavior of silicon. To provide a critical assessment of the applicability of classical simulation results to real silicon we present a comparison of the structural, dynamical, and thermodynamic quantities obtained from the SW potential with those from ab initio simulations and with available experimental data. We also discuss the sensitivity of the thermodynamic properties to model parameters.


I. INTRODUCTION
Silicon, the second most abundant element in the earth's crust, is ubiquitous in the form of silica and silicates in the natural world. In the elemental form, it is an essential component of the semiconductor technology. It was first prepared in its amorphous form by J.J. Berzelius and later the crystalline form by H.E. Sainte-Claire Deville [1] in the 1800's.
The crystalline and amorphous solid are the two most familiar forms of silicon, which have been studied extensively. The crystalline form of silicon is a tetravalent semiconductor (as is the amorphous solid) and upon melting at 1687K at ambient pressure, transforms to a metallic liquid with higher coordination number, around 6. Liquid silicon is relatively less studied, given the elevated temperatures at which it exists. Nevertheless, it has been a subject of substantial experimental, theoretical and computational investigation, both at temperatures above the melting temperature, and in the supercooled and stretched (negative pressure) states. The investigations of the metastable liquid have been motivated, as this review seeks to demonstrate, by fundamental questions regarding (i) the eventual fate of metastable liquids upon deep undercooling and stretching, (ii) the interest in the possibility of a novel transition between two distinct liquid forms in a class of "tetrahedral" liquids to which silicon belongs, and (iii) the role of the thermodynamics of metastable liquid states on the kinetics of phase transformations, particularly to the crystalline state.
Based on the extrapolated Gibbs free energies of amorphous solid and liquid phases for germanium, and a scaling of temperatures for the case of silicon, Bagley and Chen [2] and independently Spaepen and Turnbull [3] suggested a first order phase change from the amorphous solid to the liquid near T = 1349K for silicon, below the freezing point of liquid T m = 1687K. Subsequent experimental as well as simulations studies [4][5][6][7][8][9][10][11][12][13] of non-crystalline silicon supported this notion though the precise nature of the transition was unclear owing to the limited availability of information. Using a two-state model, Aptekar [14,15] in 1979 proposed a phase diagram that described the liquid and amorphous states as two states of one noncrystalline (liquid) phase, and further predicted a negative pressure critical point. The perspective that the transition from amorphous silicon to the liquid should be viewed as a liquid-liquid, rather than an (amorphous)solid-liquid transition, attracted renewed interest in light of an independent proposal of a liquid-liquid transition in the case of water [16], and the growing appreciation that such a possibility was in principle also applicable to other substances such as, e. g. silica [17], which exhibited thermodynamic and structural features similar to water.
The analogy was pursued by Angell and co-workers [18] who proposed a first order liquidliquid transition line as a feature in the pressure-temperature phase diagram of silicon, based on simulation evidence using the Stillinger-Weber (SW ) potential of silicon [19]. The experimental work of Deb and co-workers [20] observed a pressure induced amorphous-amorphous transition and speculated a possibility of an underlying liquid-liquid transition (LLT ), to be found at ambient pressure at around (1400K). From extensive simulations of silicon using the SW model potential Sastry and Angell [21] found evidence of a liquid-liquid transition at zero pressure at around 1060K, which was also subsequently supported by ab initio simulations [22,23]. Considerable recent simulation and experimental evidence lends support to the idea of a liquid-liquid phase transition at ambient pressures in silicon [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. However, such a transition is also consistent with a "critical point free" scenario, proposed recently by Angell [37], and hence, the existence of a critical point needs independent verification. The simulation evidence of a critical point at negative pressures has been provided by Vasisht et al. [32,38] for SW silicon and in ab initio simulations by Ganesh and Widom [23]. The work of Vasisht et al. also maps out the detailed phase diagram of supercooled silicon, finding phase behaviour that is similar in qualitative features to that found in similar studies of water and silica. Changes in the electronic structure accompanying this phase transition have also been studied computationally [23,26,29], demonstrating a transition from a metallic high temperature liquid to a low temperature liquid with substantially reduced conductivity.
The introductory outline above spans investigations pursued over more than three decades and is inevitably not comprehensive. In the following sections, we elaborate on some of these themes. While we present salient experimental, theoretical, and ab initio simulation results, the main body of results discussed will be from classical computer simulations of SW silicon.
We do not make an attempt to be exhaustive, and regret the omission of any significant material concerning the topic of this review. We broadly divide the investigations described into those that address an "amorphous solid to liquid" transition, and those that address a "liquid-liquid" transition in supercooled silicon. Correspondingly, the next section describes early work that largely falls into the first category, although including the work of Aptekar [14] and Ponyatovsky and Barkolov [15] which address a liquid-liquid transition. The parallel developments in the case of water (and later by extension, other tetrahedral liquids) were aimed at developing an explanatory framework for its anomalous properties. Therefore in Section III, we describe the various scenarios put forward for fluids exhibiting density and related anomalies. In Section IV, we discuss the more recent investigations, starting with the work of Angell and co-workers [18] which explore the similarities in the phase behaviour of silicon with water and other analogous substances. We discuss in this section the details of the phase diagram as revealed by simulations of SW silicon including a negative pressure critical point [32] as well as structural and dynamical properties of supercooled silicon. An important feature of the behaviour of silicon, not necessarily shared with other substances that may exhibit a liquid-liquid transition is the change in electronic structure, from a semiconductor at low temperatures to a metal in the high temperature liquid state.
We therefore discuss the electronic structure in Section V separately, including electronic structure calculations performed for structures obtained in classical simulations [26], first principle molecular dynamics (F P MD) simulation studies [22,23] and experimental work of Beye and co-workers [30] which employ the changes in the electronic density of states as an experimental probe for the liquid-liquid transition in silicon. Since a significant part of the simulation results presented are based on a classical empirical potential (performed to a large extent by the authors of this review), an important question to address is the reliability of this potential in describing the behaviour of silicon. To this end, in Section VI we compare the structural, dynamical and thermodynamic quantities obtained from the SW potential with those from ab initio simulations and with available experimental data, to provide a critical assessment of the applicability of classical simulation results to real silicon.
We also discuss the sensitivity of the thermodynamic properties to model parameters. We end with a summary in Section VII.

II. EARLY WORK ON METASTABLE SILICON
The earliest experimental work on metastable silicon includes studies carried out by Bagley and Chen [2], Spaepan and Turbull [3], Donovan et al. [4] and Thompson et al. [5]. Bagley and Chen [2], and independently Spaepen and Turnbull [3], used the available thermodynamic data for germanium (such as the heat capacities measured by Chen and Turnbull [39]), and the kinetics of crystal growth, to estimate the temperature dependence of the excess Gibbs free energy of amorphous and liquid silicon. The significantly different entropies of the two limiting states leads to the prediction (by extrapolation) of a slope discontinuity in the Gibbs free energy, or a first order phase transition from a four coordinated amorphous to a metallic liquid state. The estimated transition temperature was T al at 1349K for silicon (by scaling the melting points of germanium and silicon). Later Donovan and co-workers [4,40] performed differential scanning calorimeter (DSC) measurements on amorphous silicon, produced by ion implantation and based on Gibbs free energies deduced, estimated T al to be 1420K. We shown in FIG. 1 the Gibbs free energies estimated by Donovan et al. [40]. This transition was also confirmed by another experimental group of Thompson and co-workers [5] who performing pulsed-laser melting of amorphous silicon and estimated the T al as 1480 ± 50K.
The theoretical analysis of the thermodynamics of supercooled silicon, presented by Aptekar [14] treats the liquid as a pseudo-binary regular solution of two components, along lines explored in related contexts by Rappaport [41], Ponyatovsky and co-workers [15,42].
The two components are characterised by different local bonding environments (covalent or metallic). Correspondingly, the Gibbs free energy of the liquid is written as The parameter ω describes the degree of metallisation and is determined by the equilibrium condition of the liquid. Writing the free energy difference between the two pure liquids it is a relatively simpler task to explore a wide range of temperatures and pressures in simulations and thus bracket the region of interest for further experimental verification. For more than three decades computer simulation studies of silicon have been carried out using various empirical interaction potentials (a comprehensive study of six different empirical potential is given in the reference [43]) and also using first principle simulations [23,29,44,45]. One of most widely used potential [6,7,18,21,26] for studying silicon in computer simulations is the Stillinger-Weber (SW ) potential [19]. Using the SW potential in molecular dynamics (MD) simulations, Broughton and Li [6] performed one of the earliest studies of the liquid, crystal and amorphous phase diagram of silicon. In this work Broughton and Li found that the crystal and liquid phases are well represented by the SW potential, but the thermodynamics of the amorphous phase is poorly described and that the supercooled liquid phase does do not undergo a first order transition to an amorphous state upon cooling. Luedtke and Landman [7] showed that this failure to obtain amorphous Si via direct cooling of the melt in simulations is related to the quench rates employed. These authors [7] noted that upon cooling, a sharp change in the energy and density of the system occurred at T ∼ 1060K followed by a slow variation in these properties as cooling continues to T = 300K. These authors compared their system obtained from quenching with the amorphous phase obtained from an alternate method (involving tuning of coefficient of three body part of the SW potential) and found that both the systems have comparable structural composition. Angell and co-workers [18] using the SW potential explored a relatively wide range of temperature and pressure to chart out the phase diagram of metastable silicon. These authors also reported a well defined transition, at T ∼ 1060K, from a highly diffusive liquid states to a non-diffusive "amorphous phase", with a coordination number of 4.1. The phase diagram as suggested by Angell and co-workers is shown in the FIG. 2 (note that at higher pressures, the transition occurs at lower temperatures, leading to a negatively sloped transition line).
As noted earlier, the work of Angell and co-workers [18] made contact with the possibility of a liquid-liquid transition in the case of water, which had been proposed as one of the possible scenarios within which to understand the anomalous properties of water. Since these scenarios are relevant for our further discussion, we review them briefly in the next section before returning to more recent results for silicon.

LIES
It is well known that ice floats on water, owing to the solid form having a lower density than the liquid. It is also well known that liquid water has the maximum density at ambient pressure at 4 • C, below which the density decreases as temperature decreases, contrary to the normal behaviour wherein liquids become denser as the temperature decreases. Hence, the decrease of density of water with temperature is described as anomalous. The temperature at which the density is a maximum depends on pressure, and thus once has a locus of temperatures of maximum density (T MD) which for water at positive pressures is negatively sloped in the (P, T ) plane. Liquid water also shows anomalous behaviour in thermodynamic quantities like compressibility (K T ) and heat capacity (C P ).  [47]. Water also shows anomalous behaviour in its dynamical properties. The diffusivity of liquid water increases with increase in pressure which is abnormal for liquids. The work of Errington and Debenedetti [48] identifies a region in the phase diagram where the structure of the liquid behaves anomalously. This work also found that the anomalies in density and diffusivity occur within the structurally anomalous region (in the ρ − T phase diagram of water).
In the case of silicon, even though simulations have predicted anomalies similar to that of water, experimentally none of the anomalies have been verified. The density maximum in silicon, at zero pressure, as predicted by classical SW simulations is at T = 1350K [18,32,49] and ab initio simulations predict it at around T = 1200K [45]. The lowest temperature at which density measurements have been carried out in experiment is T = 1370K [50]. An extrapolation of a polynomial fit to the experimental data from [50] would suggest a density maximum around 1200K. In SW silicon simulations density minima have also been identified [32]. A very recent experimental measurement of isobaric heat capacity [51] has been carried out down to T = 1548K which shows an increase in heat capacity with a decrease in temperature similar to that of water below T = 309K. The SW silicon simulations predict a very weak minimum in C P at around T = 3350K. There are no experimental reports of compressibility but in simulations of SW silicon, both compressibility maxima and minima have been identified [32]. Along with these thermodynamic anomalies, silicon is also found to show anomalous behaviour of self diffusion. The diffusivity of supercooled liquid silicon is found to increase with increase in compression [32,52], which is yet to be verified in experiments. Some of these anomalies have also been seen in the case of silica [17,53,54] Various models and scenarios (based on thermodynamic constraints) have been developed to explain the thermodynamic anomalies of water [55][56][57], (and by extension, other liquids with water-like anomalous behaviour, including silicon) which are briefly discussed in this section. The observation of negative melting curves in various systems including water and silicon, means (from the Clausius-Clapeyron relation dP/dT = ∆S m /∆V m , and assuming that the entropy of the crystal is lower than that of the liquid) that the liquid density will be greater than that of the solid phase. This is a feature that is typical of the substances to which the considerations in this section apply.
Below, we describe some of the scenarios that have been explored as a way of rationalising the thermodynamics of liquids displaying anomalies, such as water. These are: (i) the stability limit conjecture [46,55,58,59], (ii) the liquid-liquid critical point scenario [16], (iii) the singularity free scenario [60] and (iv) the critical point free scenario [37].
The stability limit conjecture: The anomalous increase in water's heat capacity and compressibility with decrease in temperature, with apparent power law divergences at T s = 228K [46] was explained by Speedy [58] to be due to the approach to a spinodal line originating from the liquid-gas critical point. This spinodal has a positive slope in the (P, T) plane near the critical point, but upon intersection with the negatively sloped line of density maxima, goes through a zero slope according to the thermodynamic condition (dP/dT ) spinodal = (∂P/∂T ) isochore and retraces to higher pressures with a negative slope.
Hence in this scenario the spinodal constitutes both the superheating and supercooling limit of the liquid. Debenedetti and D'Antonio [59,61] further proposed that thermodynamic consistency also requires that a density maxima locus must necessarily have an end point.
The density maxima locus should either intersect a density minima locus and hence the liquid shows a normal behaviour in its density or terminate by intersecting at a spinodal curve (FIG.   3 (a)). Although some theoretical works have shown that [62][63][64] a re-entrant spinodal is present in models with water like properties, no compelling experimental experimental verification exists of this scenario [65].
The liquid-liquid critical point scenario: Poole et al. [16] investigated the retracing spinodal scenario using molecular dynamics simulations of the ST2 model of water. In their simulation study, the spinodal was found to be a monotonic function of T. The locus of density maxima or the T MD line, although having a negative slope at high pressures, changes to positive slope at low pressures (FIG. 3 (b)). Hence no intersection between the spinodal and the TMD line occurs. Instead the authors found evidence for a second critical point, between two forms of the liquid. Considerable simulation and theoretical investigations since the original work of Poole et al. support the possibility of a second critical point. Indeed earlier theoretical analyses using a two state description [14,41] also generically lead to this possibility [66]. There has been a substantial amount of experimental work to verify the possibility of a liquid-liquid transition in water that has lead to much evidence in support of this possibility, including recent work on confined water as a way of circumventing crystallisation in bulk water experiments. Such evidence has been critically reviewed in [65,67,68].
The singularity free hypothesis: Sastry et al. [60] proposed that a minimal scenario that was consistent with the salient anomalies did not require recourse to any thermodynamic singularities, such as a critical point or a retracing spinodal. They analysed the interrelationship between the locus of density and compressibility extrema, and showed that the change of slope of the locus of density maxima (TMD) was associated with an intersection with the locus of compressibility extrema (T EC) (FIG. 3 (c)). The relation between the temperature dependence of isothermal compressibility at the TMD and the slope of the TMD is given by where K T is the isothermal compressibility. The subscript P and TMD represents the slope at constant pressure and at the TMD at a given pressure. Since ∂ 2 v/∂T 2 > 0 at the TMD, the above relation shows that for an anomalous liquid exhibiting a negatively sloped TMD, the isothermal compressibility at constant pressure increases upon decreasing temperature and hence such increases in compressibility are not a priori an indication of singular behaviour. Calculations with a lattice model displaying singularity free scenario [60,69,70] reveal a line of compressibility maxima at low temperatures. The metastable critical point scenario may be considered to be a special case where the compressibility along the line of compressibility maxima diverges (at the critical point). Alternatively, it has been argued that the singularity free scenario is a limiting case where the critical point moves to zero T [71].
Critical point free scenario: Recently Angell [37] has discussed a possibility, related to some of the early observation of Speedy and Angell [46], in which the high temperature liquid encounters a spinodal at positive pressure, but this is a spinodal associated with a first order transition between two liquid states. Such a first order transition however does not terminate in a critical point, but may terminate at the liquid-gas spinodal. A weaker version of this picture is that no critical point may exist at positive pressures. Analysis of a model calculation by Stokely et al. [72] indicates that such a scenario may indeed arise in the limit of extreme cooperativity of hydrogen bond formation.

IV. RECENT STUDIES OF METASTABLE SILICON
As described earlier, the early studies of metastable silicon [2][3][4][5][6][7]18] probed the possibility of a liquid-amorphous transition. More recent work has attempted to find evidence that the transition is one between two liquid phases. In this section we present a brief discussion of such recent work.

A. Experimental studies
Experimental studies of supercooled silicon are very challenging because of high crystallisation rates. To explore supercooled states by cooling from the high temperature liquid, one would need to quench the liquid at rates exceeding 10 9 Ks −1 [25] to avoid crystallisation and hence using simple quenching techniques exploring deeply undercooled metastable liquid is not possible. Alternate methods like chemical vapour deposition and pressure induced techniques [25] have been employed to study the amorphous phase. Other studies have been performed using methods like aerodynamic levitation [73][74][75] or electromagnetic levitation [76][77][78][79][80][81] to avoid crystallisation induced by the containers during the experiments.
In silicon, the phase change from a low density liquid (LDL) to a high density liquid (HDL) involves a change in electrical conductivity (from a low temperature semiconducting to a high temperature metallic state), which in turn presents a number of measurable properties that can be used to detect the phase transition. Optical micrograph methods have been used to measure the change in optical reflectivity upon a change in phase [28], and the luminescence of the material is also used to detect the phase transition [20]. X-ray diffraction spectra and Raman spectra have also been used to observe the phase transition [20,28]. Experimental measurement of densities is quite difficult but in-situ measurement of structural quantities and electronic properties have been reported by various groups [30,50,76,[78][79][80][82][83][84].
Evidence for a pressure induced amorphous-amorphous transition in silicon was first shown by Deb et al. [20]. These authors studied porous silicon -π-Si (silicon having nano-porous holes in its microstructure and a large surface to volume ratio) because of its luminescence property. At ambient pressure π-Si exhibits red luminescence upon irradiating with an Argon laser. With the application of pressure (using a Diamond anvil cell) the luminescence shifted to longer wavelengths and became opaque at around P = 10GP a.
X-ray diffraction measurements showed that at around P = 7 to 8GP a the crystal diffraction pattern disappears and a broad diffraction pattern, a characteristic of an amorphous material, is observed. At around P = 10 to 12GP a, the crystalline peak disappears entirely. The authors performed Raman scattering measurements both during compression and decompression and found that upon compression to P = 13GP a, the sharp crystalline feature at around wavenumber 520cm −1 disappears, and a broad peak appears between 200 to 400cm −1 , distinct from the Raman signature of tetrahedral low density amorphous (LDA) silicon (a broad peak around 400 to 500cm −1 ). This feature is interpreted to be due to a high density amorphous (HDA) phase. Upon decompression, this feature disappears giving way to a broad peak around 400 to 500cm −1 , corresponding to LDA at low pressure.
These observation led to the conclusion that π-Si undergoes a pressure induced amorphousamorphous phase transition. In turn, this amorphous-amorphous transition was suggested to be related to a liquid-liquid transition, employing a theoretical model [85]. Direct optical observation and electrical resistance measurements carried out on amorphous silicon by McMillan and co-workers [28] showed that the HDA is highly reflective and LDA is non-reflective (see FIG. 5). From the electrical resistance measurements the authors found that there is an abrupt decrease in resistivity across the LDA-HDA transition around P = 11GP a, indicating transformation to metallic HDA. The sample was verified to be in its amorphous state (using Raman spectroscopy), since pressure induced crystallisation to β-Sn phase could also lead to a drop in resistivity.
The above experiments suggest the possibility of a liquid-liquid transition but are performed under conditions at which the amorphous forms of silicon are solid. In an attempt obtain a more direct evidence that the transition is between two liquids, Hedler et al. [25] performed ion bombardment experiments on amorphous silicon. The plastic deformations they observe of the samples are similar to the deformation seen in conventional glasses undergoing the glass transition, and the authors deduce a glass transition of around 1000K.
As described below, in both classical and ab initio simulations a clear evidence of liquidliquid phase transition has been found. Computer simulations [21,23,24,26,32,33,74] predict that the first order phase transition is characterised by a change in coordination number from 4, in LDL to greater than 5, in HDL. It has been also observed that the tetrahedrally coordinated LDL is less diffusive compared to HDL. The electronic structure calculations [23,26,74] in these simulations have shown that the LDL is less metallic than of Sastry and Angell [21] addressed these uncertainties, by seeking evidence of (a) phase co-existence, and (b) finite mobility in the low temperature phase. To probe whether a first-order transition exists, the authors carried out constant enthalpy (NPH) simulations.
A non-monotonic dependence of the enthalpy on temperature was found (FIG. 7), which is an indication of a first order phase transition. The transition temperature was found to be around T = 1060K at zero pressure. Similar behaviour is also observed in first principles simulations by Jakse and Pasturel [29]. The equilibration times in the low temperature liquid phase (see below) range from tens to hundreds of nanoseconds. The T-dependence of the diffusivity in the high temperature liquid phase (till T = 1070K) was found to be highly non-Arrhenius, characteristic of a fragile liquid [49].
It has been argued [86] that silicon, along with other tetrahedral liquids such as water and silica, should undergo a transition from fragile liquid behaviour (non-Arrhenius temperature dependence of viscosity and other transport coefficients) to strong liquid behaviour (Arrhenius temperature dependence) as the liquid makes a transition from the HDL to the LDL (either discontinuously or continuously at pressures below the critical pressure). Since data over a sufficient range of temperatures in the low temperature phase was not available to judge this matter directly for silicon, Sastry and Angell [21] took recourse to an empirical observation that the intermediate scattering function (F(k,t)) shows oscillatory behaviour in strong liquids that becomes more pronounced in small systems. In HDL phase (T = 1055K) no oscillations were observed in the F(k,t) (see FIG. 9) and upon transition to LDL phase (T = 1070K) oscillatory behaviour appears and it becomes more significant at lower system sizes (inset of FIG. 9). Independent evidence for such a transition also is obtained by the fact that the heat capacity C P drops to a value of 3.6Nk B in the low temperature phase [87].
The pair correlation function g(r), fifth neighbour distribution g 5 (r), bond angle correlation function G 3 and local bond orientation order parameter Q 3 (see IV B(4) for definitions were calculated to study the structural properties of the two liquids. Coordination numbers calculated by integrating the g(r) till its first minimum. The average coordination number was found to change from 5.12 to 4.61 in the high-temperature liquid and around 4.2 in the low-temperature liquid. The fraction of four coordinated atoms increased from about 50% (at high T) to 80% in the low T phase, indicating a greater degree of local tetrahedral order.
This change is also reflected in the local bond orientation order Q 3 (see FIG. 10 (a)). The Q 3 values for the low temperature liquid peaked at the crystal's Q 3 value, suggesting a tetrahedral local ordering whereas the high temperature liquid showed a broader peak in Q 3 . The fifth neighbour distribution (which is distribution of distances of the 5 th nearest neighbour to a given atom) for the high temperature liquid was found to be unimodal indicating that the fifth neighbour resides inside the first coordination shell, whereas the low temperature liquid showed a bimodal g 5 (r) distribution (see FIG. 10 (b)), with the larger peak shifting to the second neighbour shell in the low temperature phase.

Liquid-liquid critical point
Recently, extending the work of Sastry and Angell [21], which provided evidence of a first order liquid-liquid phase transition in silicon at zero pressure, Vasisht et al. [32,38] reported evidence of the existence of a negative pressure liquid-liquid critical point in SW silicon, based on extensive simulations of the SW model of silicon. These authors also constructed the complete phase diagram of supercooled silicon, which clearly demonstrates the interconnection between various thermodynamic anomalies and the phase behaviour of the liquid as analysed in previous works [16, 46, 58-60, 69, 88].
It is well known that the liquid-gas coexistence line does not extend to arbitrarily large T, but terminates at the critical point, where second derivatives of the free energy are singular (e. g., the heat capacity and the isothermal compressibility). Density along isotherms vs.
pressure are continuous at temperatures above the critical temperature, whereas they display a density discontinuity at lower temperatures, at the coexistence pressure. were run up to 10 relaxation times when possible and in all LDL cases, simulations were carried out for times required for the MSD to reach 1nm 2 (5 σ 2 , where σ is the atomic diameter) or for 100ns (300 million MD steps), whichever was larger. Equilibration of the system was monitored by the MSD and from the relaxation of self overlap function Q(t) defined as  . 11 (a). Below the critical temperature, the system phase separates and hence these isotherms (below T = 1133K), in the isothermal-isobaric (NP T ) MD simulations, showed jumps in densities for a small change in pressure (FIG. 11 (b)) The isothermal-isochoric (NV T ) MD simulations performed at these temperatures (T = 1133K to 1070K) (and at densities spanning the range of the density jumps in (NP T ) MD simulations) showed non-monotonic isotherms (FIG. 11 (b)). First principles simulations: In the cases of carbon and silica, computer simulations using classical empirical potentials have shown a liquid-liquid transition [17,92], but first principle MD (F P MD) simulations [93,94] show results that are not consistent with classical simulations. In silicon, Jakse and Pasturel [22] and independently Ganesh and Widom [23] have reported first principle simulation results, both of which support the proposed liquiddensity) and restricted ensemble MC simulations for locating the spinodal at low temperatures.
liquid transition in silicon. In the work of Ganesh and Widom, the authors report the emergence of a van der Waals-like loop (shown in FIG. 13), as signature of a first order phase transition at temperatures below 1182K. The maximum time span of these simulations is around 40ps [22], which seems to be very small compared to the relaxation times of LDL (tens to hundreds of nanoseconds; see below) obtained from simulations of SW silicon [21].
But the FPMD calculations are computationally very expensive compared to classical MD simulations. Hence it would be of interest to compare the equilibration times of the system simulated in FPMD and classical MD and also do a systematic study of relaxation processes in these two different methods of simulation. A comparison of properties obtained in different simulations are discussed in a later section.

Phase diagram
In order to obtain the complete phase behaviour of supercooled silicon Vasisht et al. Temperature of Minimum Density: The TMinD line is locus of density minima points, crossing which, the system returns to normal behaviour (density increases with decrease in temperature). Finding the TMinD line in supercooled silicon is challenging since one must simulate the system deep inside the supercooled region of the phase diagram (where crystallisation, slow equilibration and cavitation pose hurdles to obtaining equilibrated data).
In order to obtain equilibrated data NPT parallel tempering MC simulation technique was employed [91], in which copies of the system at different T and P are swapped periodically according to a Metropolis acceptance criterion, thereby avoiding the possibility shows high temperature spinodal isotherms (T > 2200K). These isotherms were obtained from constant volume temperature (NVT) MD simulations.
For T < 2200 K, cavitation was observed in the NVT simulations before the minima along the isotherms are reached, due to which a drastic increase in the pressure was observed. In an attempt to circumvent this problem, the authors performed REMC simulations wherein arbitrary bounds were imposed on the magnitude of the allowed density fluctuations by dividing the simulation box into a number of equal sub-cells and constraining the number of particles in each of these sub-cells [95]. However, even in the REMC simulations, the system was found to cavitate occasionally, with the formation of voids across sub-cell boundaries (with each sub-cell satisfying the applied constraint on number of particles). Hence the estimation of the spinodal at these state points were done from a quadratic extrapolation of the isotherms. The data points obtained from REMC simulations was fitted with a quadratic As a further check on the spinodal estimate, the tensile limit of the liquid was obtained by increasing the tensile pressure on the simulation cell at different constant rates. The tensile limit line is defined as the the locus of maximum tensile stress (negative pressure) a system can withstand before it fails. At a given temperature the system was first equilibrated at a high pressure value (for T < 1510K at P = −2.26GP a and for T > 1510K at P = 0GP a, by performing NPT MD simulation) and then a tensile pressure was applied which increases at a specified rate. Simulations were performed at four different constant rates of change of tensile pressure (0.1MP a/ps, 1.0MP a/ps, 10.0MP a/ps, 50.0MP a/ps). When the system reaches its limit of tensile strength, the system's density decreases drastically towards zero.
FIG . 21 show the applied pressure against the measured density for a range of temperatures, from which tensile limit line was obtained. At faster stretching rate (10MP a/ps) the tensile limit was found to be consistent with the spinodal estimates. For slow stretching rate (0.1MP a/ps) it was found that the system cavitates at higher (less negative) pressure values.
At the intermediate rate, (1.0MP a/ps) the estimated tensile limit line lies between the estimates obtained from the faster (10MP a/ps) and slow (0.01MP a/ps) stretching rates.
At very high rate of change of tensile pressure (50.0MP a/ps), the tensile limit was found to be below the spinodal estimates, indicating perhaps that the stretching rates are faster than microscopic relaxation time scales of the system. temperature T , ie not "reentrant" as predicted to be the case [58] if it intersects the TMD.
The TMD, changes slope upon intersection with the TMinC, as analysed in [60]. Evaluating the relevant equation of state data as the TMD approaches the spinodal is particularly challenging. From available data, the authors conclude that the TMinC appears to be smoothly joining with the TMC (line of compressibility maxima) that emanates from the liquid-liquid critical point. Using the PT MC simulations below the critical temperature and pressure a line of density minima was also recognised (which was recently observed in the case of water in experiments and computer simulations [88,96]) which appears to smoothly join the TMD line, as required by thermodynamic consistency.

Structural and dynamical properties
Various structural and dynamical properties of supercooled liquid silicon calculated as function of pressure and temperature are now summarised. These properties include relaxation times, structural properties and dynamical properties. In simulations, the structural relaxation times are calculated from the coherent intermediate scattering function (F (q, t)), which is defined as where δρ( q, t) is the fluctuation in local density (at time t) in Fourier space given by The F (q, t) is calculated at a q value corresponding to the first peak of structure factor S(q).
The alpha relaxation time (τ α ) is obtained usually as the time at which the F (q, t) decays by a factor of e.
The structural properties are quantified by the S(q) and its Fourier transform, the radial distribution function g(r). The structure factor is defined as The g(r) can be directly calcuated in simulation using the relation In simulation one can also measure the structural order parameters which provides local structural arrangements. The orientational order among the bonds present in the system is characterised by local orientational order Q l which is given by where or from the velocity autocorrelation function using the Green-Kubo formula Dynamics: One of the biggest challenges in the study of supercooled silicon is that, with deep under- the phase diagram and found that approaching the liquid-liquid transition line (below the critical temperature) and the compressibility maxima line (above the critical temperature), the relaxation time increases in a non-Arrhenius manner. In the FIG. 25 (a), the relaxation time as a function of temperature at two different pressure values (above -P = 0GP a and below -P = −1.88GP a critical point) is shown. FIG. 25 (b) shows relaxation times as a function of pressure for three different temperatures values.
Vasisht et al. [32] report the calculation the diffusivity at temperatures below and above the critical temperature. As the liquid transforms from HDL to LDL the diffusivity is found to change by two orders of magnitude. The diffusivity is found decrease with decrease in pressure, which is an anomalous behaviour similar that of water and silica [48,54]. At higher pressure values diffusivity goes through a maximum and a return to normal behaviour (FIG.   26). This has been also observed in ab inito simulations [52].

Structural properties:
The work of Sastry and Angell [21] showed that the liquid-liquid transition marks change in structural and dynamic properties. Structurally the high temperature, HDL is less tetra-  . 27 and FIG. 28. The authors find that the P = 0GP a isobar is above the critical pressure and hence one sees a discontinuous change of g(r), whereas the P = −1.88GP a isobar is below the critical pressure and hence a continuous evolution of structural change is seen.
The coordination number C nn is the number of atoms in the first coordination shell, and is calculated by integrating the g(r) till its first minimum (r c ) using the equation C nn = rc 0 4πr 2 ρg(r)dr. The coordination number was found to be very sensitive to the location of the first minimum of g(r). A cutoff of r c ∼ 2.96Å was found for temperatures less than 1259K and for T > 1259K, the minimum was found to shift towards higher values of r (FIG.   29 (a)). For a given temperature, with varying pressure, r c was found not to change much (FIG. 29 (b)). In FIG. 30

Relationship between Structure and Dynamics
From the data reported above of the coordination number and the diffusivity, it is evident at a qualitative level that the diffusivity in silicon is correlated with coordination number, with higher coordination number corresponding to larger diffusivity. It has been shown in [32] that the diffusivity depends quite strongly on coordination number and has only a weak temperature dependence (see FIG. 31). Scaling the diffusivity to its value at a fixed C nn in the HDL phase for different temperatures, a remarkable data collapse is obtained that spans two distinct phases, a wide range of temperature and pressure, and four decades of change in diffusivity (see inset of 31). The resulting master curve was found to fit well to a Vogel-Fulcher-Tammann (V F T ) form, D(C nn ) = D 0 exp(− A Cnn−n 0 ) with n 0 = 3.86, and also to a power law D(n) = D 0 (n − n 0 ) 3 , with n 0 = 4.06. These results suggest that the mobility of atoms is strongly tied to the presence of coordination larger than four, and that regions of higher coordination number act as "defects" that promote faster rearrangements of atomic positions. This observation is consistent with previous analysis of the role of bifurcated bonds or the fifth neighbour in determining molecular mobility in water [97,98], though ‡ Note that in [21] and [32] the integration was performed up to the first minimum of ρ4πr 2 g(r), rather than the g (r) directly, which leads to a small underestimate in the coordination number at high temperatures and pressures. seen here for a remarkably wide temperature and pressure range.

V. ELECTRONIC STRUCTURE
Among liquids that may exhibit a liquid-liquid phase transition, a feature that is special to silicon (though not uniquely so; see earlier discussions) is the change in electronic properties that accompany the liquid-liquid phase transition. Indeed, this is a feature that has been exploited in studies from early on in experimental probing of the transition. The amorphousamorphous transition in silicon has been also found to be a transition from a semiconducting low density state to a metallic high density state. The liquid form of these phases have shown similar change in the conductivity. Given that the change in electronic properties has a strong influence on the effective interatomic interactions, a question has been raised about the extent to which a classical empirical potential can capture the behaviour seen in silicon. To address a part of this question, Ashwin et al. [26], performed electronic structure calculations, using an empirical pseudo-potential, for atomic configurations obtained from classical MD simulations using the SW potential. The electronic density of states DOS(E),

VI. CRITICAL ASSESSMENT OF CLASSICAL SIMULATION RESULTS
The supercooled phase of silicon has been extensively studied in the past three decades and a vast amount of knowledge about supercooled silicon comes from computer simulations using the classical SW model potential, although there have also been many ab initio simulations performed in recent years. Any model interaction potential is parameterised so as to reproduce certain experimentally known properties like phase transition temperature, crystalline structure etc. It is difficult for a single empirical interaction potential to repro- and hence it is difficult to access large system sizes and time scales using them. These shortcomings becomes a real bottleneck when one needs to study systems at low temperatures or near a critical point, as is the case with supercooled silicon. Relaxation times increase with the lowering of temperature and also with the approach to a critical point. The spatial correlation in the system increases with the approach to a critical point which makes it necessary to study larger system sizes (though it must be mentioned that this shortcoming also applies to much of the classical simulation results presented in this review). In the studies involving supercooled phases one has to also confront the issue of crystallisation. Since crystallisation is inherent to a supercooled phase, it is not a feature that can be easily eliminated without introducing artefacts in the sampling. Hence, it may be necessary to perform simulations of multiple independent samples to obtain reasonable information on the metastable states.
Given these considerations, in the case of silicon, one must of necessity perform some of the simulations using a classical empirical potential such as the SW potential. It is therefore interesting to see to what extent the liquid state properties evaluated using the SW potential agree with available first principles and experimental properties.
To this end, in this section we compare various properties obtained in the simulations using SW potential with the available experimental as well as ab initio and other model potential results. These comparisons allow us to understand to what extent the SW silicon results are applicable to real silicon. Note that in the comprehensive study of six different empirical potentials carried out by Balamane et al. [43], the authors looked at various properties at T = 0K, but very little focus was given to the liquid state of silicon.
The Stillinger-Weber potential, by far the most widely used interaction potential for silicon, comprises of a two-and a three-body interaction potential. The Crystalline phase of silicon at low pressures is in the diamond cubic structure and it melts into a high density liquid phase. Stillinger and Weber, after a search through a class of interaction potentials with two and three body interactions, defined their empirical potential as follows: where σ is the diameter of the particles, r i is the position of particle i, and r ij is the distance between particles i and j. The two-body potential is short-ranged and has the form where A = 7.049556277, B = 0.6022245584, and a = 1.8. The repulsive three-body potential is also short-ranged, and is given by where θ jik is the angle formed by the vectors r ij and r ik and Regarding the second point raised by these authors, it is clear that a change in model parameters will shift the phase diagram around. Importantly, a shift by 5% in the value of λ changes the liquid-liquid phase transition temperature to higher values [100], to about T = 1390K at which Beaucage and Mousseau see crystallisation. It has been observed by many previous studies that crystal nucleation rates becomes very high when the transition boundary to LDL is crossed [18,21,32]. Hence it is not surprising that crystallisation is 2. The SW potential considerably underestimates the densities. However it is comparatively better than other estimates from other classical empirical potentials [102,103].
The density maximum as predicted by SW potential at zero pressure is at T = 1350K.
The EDIP simulations estimate for the density maximum is T = 1300K and MEAM potential estimates at T = 2500K.
3. Estimates of density from ab initio simulations compare better with the experimentally measured values. The density maximum as predicted by FPMD simulation is at T = 1200K for P = 0GP a.
We next look at the comparison of structural properties.
Structure factor S(q) and Pair correlation function g(r): Experiments measure the structure factor S(q) and the Fourier transform of the S(q) yields the radial distribution function g(r). By specifying the density as input and integrating the g(r), one obtains the coordination number. From our literature search we find that agreement between experiments and simulations is as good (or bad) as the agreement between different experiments.
We also find that agreement between experiments and MD simulations using SW potential is as good as between experiments and first principles simulations. A detailed comparison of experimental results, first principles results and our results using SW potential is shown below.
Experimental investigations have measured S(q) down to 1382K [81] and not below. We have chosen a range of temperatures (From T 1100 K to T 1700 K) at normal pressure.
For T = 1382K we compare the SW simulation data with the experimental data extracted We note the following points regarding the structure factor data: 1. The differences in structure factor values among different experimental reports as well as among different FPMD simulations are significant.
2. Even though all the experimental data (except Ansell et al. [73] at T 1542 K) show similar trends in the structure factor, we find a noticeable difference in the amplitudes of the structure factor. The trend captured by ab initio simulations is similar to that of SW simulations.
3. The feature of split peaks in the structure factor is seen in both ab initio and SW simulation. But in comparison with experiments, it seems that ab initio simulation captures better the first peak and the SW potential captures better the second peak of structure factor.
4. The minima and maxima of S(q) for Q > 4Å −1 is well captured by the SW potential.
5. The ratio of the first peak to the second peak is better captured by SW potential than by the ab initio simulations.
We now focus on radial distribution function obtained from experiments and simulations (shown in FIG. 43 and 44) and note the following points: 1. The most noticeable feature is the difference in the first minimum of g(r) within different experiments. The coordination number is very sensitive to the location of the first minimum of g(r). The results from Jakse et al. [74] and Krishnan et al. [75] are consistent with each other and shows the minimum of g(r) at around 3.3Å. The results from Kim [81] shows the minimum is at around 3Å.
2. Both ab initio and SW simulations have the first minimum of g(r) at around 3Å (except for Wang et al. [105] which shows the first minimum of g(r) to be at 3.3Å).
3. The second important difference is in the amplitude of the first peak of g(r). Both ab initio and SW simulations have similar amplitudes of the first peak and they are bigger than the observed amplitudes in experiments.
4. The intermediate peak between the first and the second prominent peak of g(r) as reported by Kim [81] is not found in the experiments of Jakse [74] and Krishnan [75].
But it is a prominent feature in both ab initio and SW simulations. We also compare the SW data with experimental S(q) and g(r) [107] at three different high pressure values at T = 1737K (FIG. 45 and FIG. 46). We find that at P = 4GP a, the SW data compare reasonably well with experimental data. But at around P = 14GP a the SW potential fails to capture the experimental S(q) and g(r).
From the above comparisons we find that there are noticeable differences between the g(r)'s obtained in different experiments as well as from different simulations. We also find from the comparison of simulation results with experiments, that the SW potential does as good a job as the ab initio simulations.

Coordination number:
The coordination number C nn as discussed before is calculated by integrating the g(r) till its first minimum (r c ) using the equation C nn = rc 0 4πr 2 ρg(r)dr. Hence the most important inputs that goes in to the calculation of the C nn is the r c and the density of the liquid. Hence the differences in the g(r) (both from experiments and simulations) and the difference in the densities measured and calculated by experiments and simulations, will reflect in the differences we find in the calculated coordination numbers as  [109]. We find that only FPMD simulation of Stich el al.
[44] predict the diffusivity close to the experimental data. All other FPMD as well as SW estimates of D (at T m ) is smaller than the experimental value. We find that the SW data is comparable with the FPMD data only down to T = 1250K.
In summary we find that the SW potential, while displaying significant differences with experimental data for the quantities we have discussed here, does so within deviations that are comparable to the spread between different experimental results, and these differences are comparable to those displayed by first principles simulation results. For the different quantities compared the summary of observations are as follows: 1. There are noticeable differences in structure factors and radial distribution functions reported by different experimental groups.

The density values reported by different experimental groups have a variation of about
10%.
3. The coordination numbers calculated are very sensitive to the location of the first minima of g(r) and the density. There is a large spread in the coordination numbers reported by different experimental groups, between 5 and 6.5.
4. The first principles simulations do not show a better agreement with experiments than SW simulations. 5. The SW potential underestimates the density by 5%, which leads to lower estimates of the coordination number by 4%. Even then these numbers are within the error bars of reported experimental data. 6. The SW potential estimate of diffusivity (which is comparable to the first principle simulations) is approximately a factor of 4 less than the experimental value.
It is evident that the most significant shortcoming of the SW potential appears to be in the estimates of the density. Given the detailed understanding of the SW potential, it may be interesting to fine tune the potential by varying the strength of the three body interactions and the range of interactions, which is a useful future direction to pursue.  [110,111]. A careful analysis of the free energy landscape, and finite size analysis, along the lines prompted by the work of Limmer and Chandler [112] will add valuable insights into the behaviour of metastable liquid silicon.    Below these densities, the system cavitates before the isochore passes through a minimum.          The inset in (a) shows the fifth neighbour distribution.