Thermodynamic anomalies in silicon and the relationship to the phase diagram

The evolution of thermodynamic anomalies are investigated in the pressure–temperature (pT) plane for silicon using the well-established Stillinger–Weber potential. Anomalies are observed in the density, compressibility and heat capacity. The relationships between them and with the liquid stability limit are investigated and related to the known thermodynamic constraints. The investigations are extended into the deeply supercooled regime using replica exchange techniques. Thermodynamic arguments are presented to justify the extension to low temperature, although a region of phase space is found to remain inaccessible due to unsuppressible crystallisation. The locus corresponding to the temperature of minimum compressibility is shown to display a characteristic ‘S’-shape in the pT projection which appears correlated with the underlying crystalline phase diagram. The progression of the anomalies is compared to the known underlying phase diagrams for both the crystal/liquid and amorphous/liquid states. The locations of the anomalies are also compared to those obtained from previous simulation work and (limited) experimental observations.


Introduction
Silicon remains a system of huge technological importance. Even so, an understanding of the fundamental properties of this key system, and how they inter-relate, remains elusive. For example, silicon shows a complex crystalline phase diagram, displaying changes in both coordination environment and conductivity. In the disordered state silicon shows lowand high-density amorphous forms (termed LDA and HDA respectively) characterised by different densities and entropies and which are separated by a first-order transformation [1] and which display a change in conductivity [2]. Furthermore, there is increasing evidence from both experiment and simulation which highlights how silicon may be more 'waterlike' in displaying additional, very-high density, amorphous forms [3,4]. In addition, and again mirroring the properties of water, silicon displays thermodynamic anomalies (most commonly known for the density, but also in related properties such as the compressibility and heat capacity). The structural origin of both the distinct amorphous forms and the thermodynamic anomalies share key similarities, arising from a balance between relatively close-packed and more open local coordination environments.
Effective experimental investigation of the thermodynamic anomalies remains problematic. Limited information is available, focussed on the behaviour of the density maximum, and mostly restricted to ambient pressure conditions. Furthermore, significant anomalous thermodynamic behaviour occurs in the supercooled regime (i.e. below the thermodynamic melting temperature) which makes experimental investigation, requiring highly effective suppression of crystallisation, difficult.
Levitation experiments offer the possibility of removing nucleation centres associated with the presence of surfaces and hence may allow deeper supercooling, but preclude the study of materials at anything other than ambient pressure.
A number of theories aim to explain the presence of the thermodynamic anomalies (and related stability limits). For example, Speedy [5] proposed a scenario for water in which the density anomaly locus intercepted the stability limit resulting in a tensile instability, the nature of which has subsequently been thoroughly studied by Debenedetti and co-workers [6][7][8][9]. An alternative explanation [10] invoked the existence of a second (liquid-liquid) critical point with the density anomaly locus avoiding collision with the liquid-vapour spinodal line. Subsequently, the interactions between the density and compressibility anomaly loci were rationalised using a singularityfree interpretation [11].
Simulation models offer potentially unique insight into the factors controlling the anomalous properties. Rapid cooling (by experimental standards) coupled with the ability to greatly suppress crystallisation, means that the supercooled regime may be accessed to a certain extent. Furthermore, methods such as replica exchange may effectively extend the accessible degree of supercooling. However, the rapid supercooling, coupled with the relatively short accessible simulation timescales, renders results from such simulations open to interpretation. For example, in reference [12] a liquid-liquid critical point is detected, while in reference [13] the authors identify it as an artefact of the replica-exchange method employed, specifically the use of relatively short equilibration times. However, more recent studies have identified potential problems with the determination of the free energy surfaces [14]. Even more recent studies claim the observation of an LLCP results from slow crystallisation and loss of equilibrium [15]. In the present work we use a definition of what is effectively equilibratable in order to identify the regions of the phase space in which crystallisation becomes prohibitive.
Further potential problems may be associated with the sensitivity of the anomalies to the details of the models themselves. The models applied are often relatively simple (which usually equates to a relatively small number of parameters) which allows for an efficient exploration of the phase space with the large systems required. However, potentials may lack sufficient transferability away from the state points on which they were parameterised. Electronic structure methods, in particular those based on density functional theory (DFT) offer potentially more accurate results at higher computational cost, although issues of model details (for example, the functionals applied) remain. Zhao et al, for example, have applied DFT to investigate the anomalies in silicon and have uncovered a hidden critical point [16]. The implication is that the details of the model may not only affect the presence and locations of any anomalies in phase space, but may even control the scenario by which they may be classified. More recent work suggests an even more complex landscape, with DFT calculations indicating the presence of three liquid-liquid phase transitions [17]. Machine-learning algorithms may help to bridge this gap by offering the potential for accessing DFT-level accuracy at a fraction of the computational cost [3,[18][19][20].
Simulation investigations generally utilise one of two approaches. In the first, models are developed which focus on specific (experimentally-accessible) systems (for silicon, for example, see references [12,[21][22][23][24][25][26]). In the second, generic potentials are employed in which key potential parameters may be varied systematically in order to follow the evolution of key properties (see, for example, references [23,24,[27][28][29]). Recent theoretical and computational work has focussed on two key aspects regarding these anomalies. In one case work has focussed on mapping how the anomalies shift in phase space as a potential model is systematically modified [30]. In the second case the mathematical constraints, which govern how these anomalies may intersect and interact with each other and stability limits, have been studied [29]. The latter builds on the previous work which linked the compressibility and density anomaly loci [11] and the heat capacity and density analogues [29,31,32].
In this paper we employ a Stillinger-Weber (SW) potential as an example of a well-studied (and computationallytractable) model [33]. We employ a combination of 'standard' molecular dynamics (MD) and replica exchange MD to allow the deeply supercooled regime to be effectively accessed and sampled. We map the thermodynamic anomalies in the density, compressibility and heat capacity as well as mapping the stability limit. This allows conclusions to be drawn regarding the nature of the underlying scenarios which are driving the behaviour of the anomalies. We also highlight the relationships to the underlying phase diagram of silicon (which has been detailed using the same potential model) considering both crystalline and amorphous forms. In addition, we highlight (using simulations and by reference to previous results) how relatively small changes in potential model significantly alter the underlying anomalies.

Extending the accessible phase space
In the present work we attempt to determine the thermodynamics of the disordered states beyond the limits of their stability with respect to phases with inherent symmetry. The primary motivation is that, although crystallisation cannot be suppressed in all regions of the thermodynamic phase space, often enough thermodynamic information can be extracted prior to crystallisation to allow meaningful conclusions to be drawn. The origin and nature of the anomalous properties are strictly tied to the disordered states, regardless of whether they are thermodynamically stable at a given state point. For example, many studies depict the density anomaly (and other related anomalies, such as those in the heat capacity and isothermal compressibility) as terminating at the homogeneous nucleation limit with uncontrolled gradients. However, such scenarios are at variance with basic thermodynamic relations (see, for example, reference [29]) which dictate the specific properties and shapes required for anomaly loci to terminate. As a result we look to extend the work from reference [12] in terms of all possible thermodynamic anomalies and consider the high pressure regions to look for consistencies between the thermodynamic anomalies in order to help verify our approach. The main purpose of the replica-exchange simulations (see below), therefore, is to attempt to estimate the thermodynamics of the disordered states beyond their stability limit with respect to the crystalline phases.
The main issue with sampling thermodynamic quantities from simulations in unstable regimes is that the existence of such states is in perfect accordance with the mathematical apparatus of thermodynamics but the laws of statistical mechanics are strictly violated [34].
In a physical sense such thermodynamic states are not accessible because the system is unstable (neither stable nor meta-stable) with respect to a form with a different chemical potential. At those state points the system 'flows' in the direction of the composition variable. Formally, the free energy, F, of a system with a single component (in the absence of external fields) is given, in the canonical ensemble, by where the symbols have their usual thermodynamic meanings. It is physically impossible to constrain the system to consider only phases with different chemical potentials that have no symmetry. However, the formal mathematical apparatus does not prohibit such states and so one can write, Such states are (in reality) unobtainable but can still be studied, for example purely in terms of the equation of state (EOS), a well-known example being the section of the EOS between spinodals separating liquid and gas often referred to as a van der Waals loop. Inside this region the compressibility becomes negative, signalling an unphysical behaviour which is often 'fixed' by using a Maxwell construction. A further step is to write the strict single component, single phase free energy surface: A van der Waals loop shows negative compressibility which is unphysical and yet are in accordance with the equations written above in a purely mathematical sense. Furthermore, simulating systems inside such loops is potentially problematic. Binder et al [35], for example, highlight how (in a simulation) the loops are mean-field artefacts resulting from finite system sizes.
As the laws of thermodynamics require strict equilibrium (i.e. no meta-stable states) well-defined boundaries can be drawn in the pressure-temperature (pT) projection between phases which define the phase diagram (see for example figure 1(b)-discussed below). The same holds for statistical mechanics, where ergodicity has to be guaranteed in order for it to be strictly applied. Both cases are in agreement with equation (1). However, strict equilibrium can be violated in real systems (leading to, for example, supercooled liquid and glassy states) and metastability can be achieved over significant time-scales for which 'thermodynamic properties' can be determined reproducibly (the 'traditional thermodynamics picture' [34]). In those cases the phase diagram can be extended into these meta-stable regimes and meta-stable equilibrium can be achieved (see for example figure 1(c)discussed below) until the limit of stability is reached. One might also be interested in the regions where the phases are unstable (see, for example figure 1(a)-discussed below). Considering only disordered phases, the extension of the unstable disordered phase diagram should be viable up to the point at which the heat capacity or compressibility start diverging and become negative, highlighting the limit of stability of a disordered phase with respect to another disordered phase (seen, for example, in a liquid-gas transition). Although strictly nonstable states do not conform to the equilibrium axiom of thermodynamics, the mathematical apparatus can, in practice, be extended to them without loss of generality. Furthermore, the thermodynamic functions obtained can be usefully assessed, for example in terms of their continuity between the different regions of stability and in terms of their required interactions (see, for example, reference [29]).
Given that the simulation (and experimental investigation) of metastable states are well accepted, the natural first step in studying unstable thermodynamic states is to simply neglect the fact that the postulates of statistical mechanics are violated (as they are while simulating metastable states) and simply perform simulations and attempt to extract any meaningful information. It should be noted that an important difference here is that, for metastable states, the ergodicity can usually be satisfied locally, meaning that within the part of the phase space in which system is trapped ergodicity can usually be satisfied. For unstable states this cannot be guaranteed-one could increase the phase space sampling by employing different enhanced sampling methods to get as much sampling in the unstable regime as possible before finally collapsing to a (meta)stable phase at best.

Potential models and methods
We employ a SW potential, originally parameterised for silicon [33]. The SW potential is chosen as an example of a relatively simple model with the system energy decomposed into a sum of two-and three-body terms, which favour local close-packed and tetrahedral environments respectively. The relative magnitude of these two energetic terms is controlled through a single parameter, λ. Varying this parameter allows this potential to effectively model a number of related systems, for example, silicon (λ = 21 [33]), phosphorus (λ 16.5 [36,37]), germanium (λ = 20 [38]), carbon (λ = 26.2 [39]), and even water (λ = 23.15 [40]).
We apply MD in the NVT ensemble using LAMMPS [41] with temperatures maintained using Nosé-Hoover thermostats [42,43] and with a system size of 512 atoms. A range of isochores are investigated using temperature increments of 50 K. The volume is increased by changing the simulation cell lengths in increments of 0.04 Å. System equilibration is achieved using simulations of t ∼ 2 ns with production runs of To emphasize this the liquid/crystal coexistence curves are shown for the diamond, SC16 and clathrate structures (grey lines from reference [50] and the brown line from reference [49]). The liquid/gas coexistence curve (orange line) is at p ∼ 0 for the temperature range shown. The figure is shown on an expanded pressure scale in order to highlight the 'S' shape of the compressibility locus. The TminC locus is shown extended below the liquid/gas coexistence curve (dashed green line) in order to highlight the 'S' shape. Panel (c) shows the full set of anomalies and related properties along with the region of phase space for which crystallisation is unavoidable. The locus T eq highlights the temperature below which equilibration is impossible on the simulation time-scale. The panel also highlights the temperature of maximum crystallisation (T X ) and homogeneous nucleation temperature (T hom -red triangle) scaled from references [24,56] as described in the text. The panel also highlights the region of the phase space for which crystallisation appears unavoidable (region labelled 'crys.'). In panels (a) and (c) the dashed lines show extrapolations of the observed loci to regions of phase space in which they are not resolved in our simulations though are consistent with thermodynamic analysis [29]. Panel (d) shows the set of anomalies from panel (a) but in the VT projection. All panels also show the pressure and temperature in reduced units. t ∼ 2.5 ns. We sweep the thermodynamic state space adjacent to the liquid-vapour stability limit using significantly longer production runs of t ∼ 150 ns and with increased temperature and volume sampling [12]. At state points at which dynamical arrest is observed on the simulation time-scale replica exchange molecular dynamics (REMD) [44,45] is employed in which the system is coupled to higher temperature states, enabling a more complete exploration of the disordered state configurational space. Temperatures are exchanged every 500 time steps (Δt ∼ 0.2 ps) with total simulation times of t ∼ 4 ns. The starting configurations are taken from the 'standard' MD simulations. The REMD is tested using simulations of t ∼ 2 ns and t ∼ 8 ns, observing a slight increase in the sampling quality with increased simulation time to t ∼4 ns but very little further improvement to t ∼ 8 ns. The 'crossover' temperature from MD to REMD is selected by identifying the temperature at which the system effectively diffuses on the simulation time-scale at this temperature across the whole volume range considered.
In this paper the temperature of maximum and minimum densities are denoted TMD and TminD respectively. The analogous anomalies in the compressibility and heat capacity are denoted TMC, TminC, TMH and TminH respectively. Thermodynamic data is analysed using NumPy and SciPy [46]. Points on an isochore are filtered using a second order Butterworth filter with a 4th order univariate spline used to interpolate the filtered data. The first derivative is then calculated and identified as maxima or minima and hence corresponding to points on the TMD or TminD loci. An analogous procedure is employed to locate the compressibility anomalies with isotherms interpolated using a 4th order spline fit and the derivative taken to obtain the isothermal compressibility. The compressibility data is then recast on a pressure/temperature grid and extrema are identified with respect to temperature at constant pressure using the same procedure as that used to identify the density anomalies. Filtering is avoided in identifying compressibility maxima in order to avoid potential underestimation close to the critical points. An analogous procedure is employed to determine the heat capacity anomalies using The liquid-vapour stability limits are calculated by determining the lowest pressure point on an isotherm for each temperature. This approach generates an upper bound estimate of the liquid-vapour spinodal line in the pT plane. As a result we refer to this line throughout as the 'stability limit' (labelled 'SL') rather than the 'spinodal line'. Technically we could refer to this as a 'cavitation line' as the instability manifests as a cavitation event (and as we have in previous work [29] and which is in agreement with reference [12]). The LLCPs are estimated from the liquid-liquid spinodals which are calculated by identifying pairs of minima and maxima on the isotherms. The existence of LLCPs in SW models remains controversial. The most recent high precision free energy studies confirm that the region of phase space where LLCPs were previously reported is unstable with respect to the crystal [15].
The results presented here are in agreement with this viewpoint and we where also unable to obtain a stable equilibrium in that region of the phase space. The unequilibrated data which we used to represent the thermodynamically-unstable continuation of the disordered manifold (as discussed previously) suggests that two LLCPs exist in this unstable region. These data would be in agreement with the concept of 'virtual criticality' suggested by Holten et al for water [47]. In that work the location of LLCP was extrapolated from a two state model beyond the limit of the stability for the supercooled water with respect to ice, with the result that the LLCP could sit in the unstable region. Finally, it is also possible that the onset of this virtual criticality is related to the peculiar shape of the instability region of the model (see the cut black region in figure 1(c)). The existence of two LLCPs does not affect the general discussions presented here regarding the behaviour of the thermodynamic anomalies. As discussed above, criteria are defined which allows us to extract thermodynamic information from underneath the stability limit. By considering the time evolution of the local tetrahedrality (q 4 ) the temperature locus, T eq , can be identified, below which full thermodynamic equilibrium cannot be established over the simulation time-scale. It is, of course, possible some states are could be equilibrated with longer simulation times. For each simulation the relaxation time of q 4 was determined from the corresponding autocorrelation function. Simulations are considered equilibrated when the total simulation time exceeds 100 of these relaxation times. The temperatures at which this equilibration is attained defines T eq . All state points at higher temperature (at a given pressure) can be fully equilibrated. To obtained more accurate values T eq was estimated from simulations of around 30 ns. In the majority of the phase space at T < T eq the system attains a quasi-equilibrium (constant q 4 ) from which thermodynamic properties can be derived. However, the region of phase space at temperatures below the identified critical points show crystallisation and, as a result, thermodynamic properties cannot be effectively extracted.

Results and discussion
We begin by considering the locations and behaviour of the key thermodynamic anomaly loci along with the liquid stability limit. We then continue to explore the relationships between these features and the underlying crystalline and amorphous phase diagrams. Figure 1 shows the complete set of anomalies studied along with the stability limit for the potential parameter λ = 21 (most appropriate to model liquid silicon) and acts to highlight the known relationships between them [29,31,32,48] and which we shall elaborate on below. Figure 1(a) shows the complete set of anomalies (and related properties) obtained here. Figure 1(b) shows only the properties obtained at thermodynamic equilibrium and highlights the liquid/crystal coexistence curves for the clathrate, diamond and SC16 phases (from references [49,50] respectively). This figure is shown on an expanded pressure scale to highlight the behaviour of the compressibility anomaly locus at high pressure (see below). Figure 1(c) highlights the region of the phase space for which crystallisation cannot be suppressed. In figure 1(a) the TMD line (by far the most studied of these anomalies) extends from the positive to negative pressure regimes. As the TMD line becomes re-entrant it collides with the stability limit as shown, re-emerging as the TminD line at lower temperature which then extends from negative to positive pressure. The collision of both the TMD and TminD loci with the stability limit occur at zero gradient for the latter as required [29], with the collisions at {p, T} = {−4.2 GPa, 1780 K} for the TMD and {p, T} = {−3.5 GPa, 1130 K} for the TminD. The shape of the stability limit shows a local maximum and minimum as a result of the intersection with the set of density anomaly loci. The TMD and TminD lines do not quite form a complete loop at this value of λ, although may do so at higher values for which the loop avoids suppression by the stability limit [29]. At negative pressures the possible transformation between the two is intercepted by the stability limit while at high pressures the features are not resolved. We note, however, that the existence of both the TminH ↔ TMH and TminD ↔ TMD crossovers is supported by thermodynamic analysis [29]. The compressibility locus intersects the density locus at the highest temperature accessed by the TMD and the lowest temperature accessed by the TminD and corresponding to ∂T ∂P (TminD or TMD) = 0 as required [11,29]. The TMD is intersected by the TminC at {p, T} = {−3.95 GPa, 1820 K} and the TminD by the TMC at {p, T} = {1. 70 GPa, 740 K}. The TminC ↔ TMC crossover occurs inside the region enclosed by the density anomalies. The singularity free interpretation [11] suggests the possibility of three scenarios which differ in terms of the minimum compressibility anomaly locus as it intersects the maximum density anomaly locus. The case observed here corresponds to TEC1 in reference [11], at which point the TminC intercepts the TMD at a positive TminC gradient. The actual crossover appears obscured by the stability limit but can be estimated as {p, T} = {−4.30 GPa, 1530 K}. The TMC collides with the stability limit at {p, T} = {3.87 GPa, 1430 K} and re-emerges at higher pressures and lower temperatures as a TminC {p, T} = {4.08 GPa, 1720 K}. The TminC and TMC loci persist to high pressures and appear to 'repel' at extreme pressures with a concomitant change in the sign of the gradient. This gives a characteristic 'S' shape to both compressibility anomalies as can be seen in figure 1(b). At lower values of λ = 18, the TMC and TminC loci may merge at high pressure [29], while high values of λ = 27 the anomalies weaken and the 'S' shape of the TminC locus is replaced by a line with only positive gradient in the pT projection [30,51]. A similar progression has been found in models for BeF 2 in which the anion polariability is used as a free parameter (in the analogue of the use of the SW λ parameter) [52]. The 'S' shape in the TminC locus in itself is a highly unusual feature of these systems. Until recently [51,52] it has been largely overlooked due to the focus on lower pressure behaviour, in particular in the region of phase space occupied by the thermodynamic anomalies and their associated critical points. However, simple analytic models, such as the two-state model, fail to predict an 'S' shaped compressibility anomaly locus for any of the main 'scenarios' commonly used to interpret anomalies in water and water-like systems [24,[53][54][55]. Figure 1(b) hints at a correlation between the liquid/crystal coexistence curves and the 'S' shape in the TminC locus (as been discussed previously in, for example, reference [30]). This will be discussed further in section 3.2. A recent DFT study identifies the liquid-liquid phase transition between high-and low-density silicon at T = 1200 K and at densities between 2.27 and 2.57 g cm −3 while the LLCPs in this study suggest the phase transition would occur between the densities of 2.356 and 2.426 g cm −3 [17]. The TMH locus emanates from the LLCP as required. The heat capacity locus anomalies intersect the density anomalies when the latter has zero gradient [29]. At both the high and low pressure limits these correspond to the TMD ↔ TminD crossovers (with the lower pressure crossover suppressed beneath the stability limit as discussed above). In addition, the TMH transforms to the TminH near the TMD  [29]. Figure 1(c) shows the anomalies and related properties generated across the whole pT plane studied, with the region of the phase space for which no equilibratable disordered states could be stabilised highlighted. This region terminates at high temperature (at around the identified LLCPs) indicating that the inability to equilibrate is naturally correlated with the expected uncontrollable fluctuations present near the critical points. Figure 1(c) also shows the temperature locus, T eq , below which equilibration is impossible over the simulation time-scales. The figure also shows estimates of the temperature of maximum crystallisation rate (T X , from reference [56]) and the homogeneous nucleation point (T hom , from reference [24]). The former is taken from reference [56] in which a T X locus is determined for the mW water model (effectively a modified SW potential for which λ = 23.15). The T X locus appears to mirror the TMD locus (with T X /T TMD ∼ 0.8-0.9) and so here T X is estimated from the TMD locus as T X ∼ 0.8T TMD . The homogeneous nucleation temperature is determined in reference [24] for λ = 20.75, 22.75 and 23.15. These values are used here to extrapolate a value for λ = 21 of T ∼ 1060 K. As would be expected, these extrapolated temperatures are broadly at the high temperature limit of the instability region identified here and correlate near-perfectly with T eq . Finally, figure 1(d) shows the complete set of anomalies highlighted in figure 1(a), shown in the alternative VT projection and acting to further highlight the key relationships between anomalies discussed above.

Comparison to previous work
The locations of the anomalies in the thermodynamic phase space has been the subject of a number of significant theoretical and computational investigations. Many investigations have focussed on particular systems (for silicon see, for example, references [12,15,[21][22][23][24][25][26]). Further studies have systematically varied key potential parameters to follow the evolution of key properties (see, for example, references [15,23,24,[27][28][29]). While the majority of work has focussed on the location and behaviour of the TMD locus (as the most well-known anomaly) some investigations have included the anomalies in the compressibility and heat capacity. Many of these anomalies are inter-linked in the sense that their intersection must follow well-defined thermodynamic rules (see, for example, reference [29] and references therein). However, experimental data is relatively limited, again reflecting the inherent difficulty of working in the supercooled regime (in particular, successfully suppressing crystallisation). Figure 2(a) shows the comparison of the density anomaly loci (and stability limits) obtained from the present work with that obtained from previous simulation [12,28,[57][58][59] and experiment [60,61]. For the TMD the present results are in agreement with those obtained by Vasisht et al [12], Beaucage and Mousseau [57] and Singh et al [28] (who used the same SW potential as applied here). The TMD at ambient pressure Figure 2. The anomaly loci and stability limit obtained from the present work compared to previous computational and experimental data from references [12,16,28,[57][58][59] and [60,61] respectively. In all panels the black and red lines show the TMD and TminD loci from the present work. Panel (a) presents data for the TMD and TminD loci. TMD and TminD loci from previous simulations are shown by open and closed circles respectively from the sources indicated on the legend, from Zhao et al [16], Vasisht et al [12], Beaucage and Mousseau [57], Wang et al [58], Angell et al [59] and Singh et al [28]. The spinodal line from Vasisht et al [12] is shown as a purple line. The black crosses from experimental TMD data from references [60,61]. Panel (b) shows the effect of changing the SW parameter λ. The orange → green → black → magenta → cyan lines show the effect of changing λ from 23 → 22 → 21 → 20 → 19. The black crosses from experimental TMD data from references [60,61]. Panel (c) presents data for the TMC, TminC, TMH, and TminH. The anomaly loci from the present work are shown as indicated in the legend (using the same colour scheme as for figure 1). TMC, TminC and TMH data from previous simulations (from Singh et al [28] and Vasisht et al [12]) are shown for comparison. Panel (d) presents the data compared to previous estimations of the low and high density coexistence curves. The green solid and dashed lines show the LDA/HDA coexistence curve, and related spinodals from Daisenberger et al [49]. The cyan lines show the corresponding data from Zhao et al [16] while the magenta star shows the LLCP from Vasisht et al [12]. The purple and light green lines show the HDL-V and LDL-V spinodals from Zhao et al [16]. All panels also show the pressure and temperature in reduced units. obtained by Wang et al (calculated using a modified SW potential) is suppressed by ΔT ∼ 400 K [58]. In that case the original SW potential was modified by using a smaller value for the parameter, σ, taking σ = 2.0618 • A compared with the original value of σ = 2.0951 • A, σ being related to the atom radius. The difference in the location of the respective TMDs in the pT phase space highlights the sensitivity of these loci to the details of the underlying potential model (here a change in effective radius of ∼1.6% resulting in a temperature change of around 400 K). The TMD locus obtained by Zhao et al [16] using a DFT approach, is suppressed compared to that obtained from the original SW potential by ΔT ∼ 400 K.
To help emphasize the sensitivity of location of the TMD locus to the details of the potential model, figure 2(b) shows the TMD locus determined for the original SW potential except varying the λ-parameter (from reference [33] and which effectively controls the balance between the two-and three-body interaction strengths). TMD loci are shown for λ = 23 to 19 inclusive (taken from reference [30]). Again, a relatively small change in a potential parameter, here from λ = 21 to λ = 20 suppresses the TMD loci by ΔT ∼ 400 K [62].
The observed range of temperatures for the TMD obtained from experiment at ambient pressure is T = 996 − 1221 K obtained using laser-heated samples suspended by application of an electrostatic levitator [60,61]. The large range in the temperatures reflects the required extrapolation of the experimental data into the deeply supercooled regime. Comparisons for the observed TminD are more sparse with no experimental data available. The locus determined here agrees well with that from Beaucage and Mousseau [57] while the data from Vasisht et al [12] and Singh et al [28] are ΔT ∼ 100 K higher and lower respectively. The greater spread of temperatures observed for the TminD with respect to the TMD reflects the greater difficulty in working in the even more deeply supercooled regime. Figure 2(a) also shows the stability limit determined here compared with the spinodal line from reference [12] estimated from a quadratic extrapolation of isotherms obtained from restricted ensemble Monte Carlo simulations. As discussed above, here the stability limit is estimated as an upper pressure bound of liquid stability from the isotherms and so is not strictly equivalent to the spinodal line. It is worth noting that the evaluation of spinodal lines outside of the mean-field limit is impossible [63] and, as such, all approaches strictly generate estimates. In addition Vasisht et al [12] calculate a tensile limit line as the most negative pressure obtainable from simulation in which the simulation cell is slowly stretched. Both the spinodal and tensile limits point towards the density anomaly loci avoiding collision with the true spinodal, although neither of the two results guarantee such behaviour. We can, however, show that the observed collision between the stability limit and the density anomalies (and for that matter any collisions at lower pressures) are not governed by speedy's stability limit conjecture scenario [5] by following the methodology set out in our previous work [29] which includes analysis of the signs of the second derivatives of the isochores, isotherms and density anomaly loci in different projections. Finally, the fact that stability limit (cavitation line) behaves according to thermodynamic rules as if it was a spinodal line [29] is quite odd. Such behaviour is not directly obvious from a purely theoretical viewpoint, but remains interesting nevertheless. The application of the REMD methods allow us to determine the stability limit to the deeply supercooled regime and hence allows us to uncover a more complex shape (associated with the intersection with both the TMD and TminD loci as described above). Figure 2(c) shows the comparison of the anomalies observed in the compressibility and heat capacity. The former are in excellent agreement with those obtained by Vasisht et al [12] using the same SW potential, while the latter is in good agreement with the single point obtained by Singh et al [28]. Note again that the TminC locus has a characteristic 'S' shape at high T and p as shown on the expanded pressure scale in figure 1(b) [29]. Considering the heat capacity anomalies, experimental data is again limited, with that from Rhim et al showing a weak maximum and minimum in C p at T ∼ 1380 K and 1850 K respectively (see figure 5 in reference [61]), at the very extremes of the experimental range. A later study by Zhou et al shows no low temperature maximum but does show a weak minimum at T ∼ 1836 K (see figure 6 in reference [60] and an extrapolation equation (11) in the same work, although the latter is just outside of the stated range of applicable temperatures). Neither the observed weak minimum or maxima appear to correspond to any of the values predicted from the models. Figure 2(d) shows the comparison with the critical points and associated amorphous (liquid) coexistence curves. The higher temperature LLCP located here is in good agreement with that identified by Vasisht et al [12] and close to that identified prior to that [2,49]. These critical points are at significantly less negative pressure and lower temperature compared to that identified by Zhao et al [16]. In that case the LLCP appear at {p, T} {−3.1 GPa, 1420 K} compared to {p, T} {−0.5 GPa, 1050 K} for those obtained from the SW potential. Differences are also seen in the high and low density state spinodals with those obtained by Zhao et al more widely spaced than those obtained previously [2,16]. The figure also shows the low and high density state/vapour spinodals from reference [16]. Finally, figure 2(d) also shows the range of experimental melting points observed for thin films of amorphous silicon; T al = 1480 ± 50 K [64] and T al = 1295 − 1420 K [65].

The relationship to the phase diagram
A value of λ = 21 was originally determined as ideal to model liquid silicon [33] (with 'enhanced' values applied to mimic the formation of amorphous states [66,67]). In order to help place the anomaly loci shown in figure 1 into the context of the broader phase diagram, figures 1(b) and 3 shows the same loci along with a number of previously determined thermodynamic features. Figure 1(b) is shown on a more expanded pressure scale in order to highlight the characteristic 'S' shape of the TminC locus (see section 3.1) whilst figure 3 is divided into (a) crystalline phase and (b) amorphous form diagrams [49,68]. Figure 1(b) highlights how the 'S'-shape of the TminC anomaly locus effectively traces along the liquid/crystal coexistence curves. In fact, the 'S'-shape of this locus becomes less pronounced as λ increases which is consistent with the diamond phase becoming the only thermodynamically stable crystalline phase [30]. At λ = 21 (as used here to model Si) the thermodynamically stable crystalline phases are clathrate, diamond and SC16. However, at lower values of λ structures such as the β -Sn and recently characterised 'X'phase [69,70] will become significant. The negative Clapeyron slope (dp/dT) associated with the diamond/liquid coexistence curve, coupled with the positive slope associated with the liquid/crystal coexistence curves for higher pressure crystal phases, leads to the characteristic 'V' shape (when displayed with the pressure on the abscissa) [71]. Figure 1(b) shows that this 'V"' shape is clearly correlated with the 'S' shape of the TminC anomaly locus. Figure 3(a) shows the diamond/liquid and clathrate/liquid coexistence curves. As discussed in references [49,68] the diamond crystal is thermodynamically favoured for p −1.36 GPa while, for more negative pressures, the more open (less dense) clathrate structure becomes stable. The figure shows the coexistence curve for the Si 136 clathrate/liquid interface, although we note that the analogous curve for the Si 46 clathrate follows approximately the same path [49]. Interestingly, the TMD line determined in the The anomaly loci (and related properties) shown in the pT projection for λ = 21 (corresponding to the SW model parametrisation for silicon) and shown alongside coexistence data from references [49,68]. Panel (a) shows the crystalline coexistence curves for the diamond/liquid and Si 136 clathrate/liquid interfaces and the previously determined liquid/gas spinodal. Panel (b) shows the LDA/HDA coexistence curve from reference [49]. Both panels also show the pressure and temperature in reduced units.
present work approximately follows the clathrate/liquid coexistence curve, but supercooled by ∼150 K. The TMD remains in the supercooled regime with respect to the clathrate/liquid coexistence curve across the whole pressure range. The TMD line does, however, intercept the diamond/liquid coexistence curve at {P, T} = {−3.5 GPa, 1800 K}. This observation further highlights the potential problems in observing the anomalous behaviour from direct experiments. Crystallisation to the diamond crystal phase may become thermodynamically unfavourable (and 'uncover' the TMD locus) but crystallisation to an alternative (here clathrate) crystal structure nearsimultaneously becomes favourable. It is interesting to note that the TMD locus appears to follow a similar pathway to the clathrate/liquid coexistence curve, albeit suppressed by ΔT ∼ 100 K. At ambient pressure the TMD and TminD are supercooled (with respect to the diamond/liquid phase transition temperature) by ΔT = 325 K and 830 K respectively. Figure 3(a) also shows the liquid/gas spinodal from reference [49] compared to the present, more accurate calculation which highlights, for example, the curvature of the stability limit locus as discussed above. The remaining anomalies are present in the supercooled regime with the exception of the TminC which is ΔT ∼ 300 K above the melting temperature at ambient pressure. In the positive pressure regime the TminC locus approximately follows the path of the diamond/liquid coexistence curve. In the negative pressure regime the TminC becomes reentrant, crossing the clathrate/liquid phase boundary at {p, T} {−3.5 GPa, 1900 K} as it approaches the TMD. The high pressure section of the TminH locus crosses the diamond/liquid coexistence curve (at {p, T} {6.0 GPa, 1250 K}) and appears to correspond to the limit at which this boundary could be detected from simulation.
For completeness figure 3(b) shows the amorphous phase diagram including highlighting the location of the LLCP from reference [49] (in good agreement with that determined here more precisely) and the associated LDA/HDA coexistence curves and associated spinodals which emanate from this point as determined in reference [2].

Discussion and conclusions
In this paper the evolution of thermodynamic anomalies in the density, compressibility and heat capacity have been investigated for silicon using the SW potential. We have highlighted how the anomalies inter-relate and relate to the stability limit. Replica exchange methods have allowed the deeply supercooled phase space regime to be accessed and, although arguments can be made as to the validity (or not) of this methodology, the results obtained are thermodynamically consistent. We have related the locations of the anomaly loci to the underlying crystalline and amorphous phase diagrams. The temperature of maximum density anomaly appears to follow the liquid/clathrate coexistence curve, while the temperature of minimum compressibility follows the analogous liquid/diamond curve at lower pressure and the liquid/SC16 coexistence curve at higher pressure, giving this locus a characteristic 'S' shape. A comparison between the anomaly loci determined here and those obtained previously highlights how the locations of the anomalies in phase space are highly sensitive to the model details. More surprisingly, comparison to recent DFT studies indicates that even the scenario (which describes the broad behaviour of the anomalies in the phase space) may also be highly model dependent. Large scale simulations, using machine-learning protocols to bridge between the DFT and potential model simulations, are in progress to further evaluate these differences.
RCUK data management requirements. DF is grateful to St. Edmund Hall and the Clarendon Fund (University of Oxford) for financial support.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.