The progression of thermodynamic anomalies in MX2 networks with local tetrahedral geometries

Key thermodynamic anomalies in density and compressibility, as well as the related stability limits, are determined using an ionic model for BeF2 which includes many-body polarization terms. BeF2 is chosen as an example of an archetypal network-forming system whose structure can be rationalised in terms of connected local tetrahedral coordination polyhedra. The anion dipole polarizability (which effectively controls the bond angles linking neighbouring tetrahedra) is used as a single free parameter in order to help rationalise the changes in the anomaly locations in phase space, whilst all other potential parameters remain fixed. The anomalies and stability limits systematically shift to lower temperature and higher pressure as the anion polarizability is increased. At high dipole polarizabilities the temperature of maximum density anomaly locus becomes suppressed into the supercooled regime of the phase space. The movements of the anomaly loci are analysed in terms of the network structure and the correlation with the inter-tetrahedral bond angles is considered. The high sensitivity of the anomalies to the details of the potential models applied is discussed with reference to previous works on related systems. The relationship to analogous studies on Stillinger–Weber liquids is discussed.


Introduction
Materials which are generally considered as network-forming (for example, atomistic systems such as C, Si, or Ge, and molecular systems such as H 2 O, SiO 2 , GeO 2 or BeF 2 ) may display so-called anomalous properties which may extend to both structure and dynamics. H 2 O is perhaps the most wellknown example as it displays an anomaly in the density at 1 Present address: Institute of Industrial Science, University of Tokyo, 4 Chome-6-1 Komaba, Meguro City, Tokyo 153-8508, Japan. 2 The author to whom any correspondence should be addressed.
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
T ∼ 4 • C at ambient pressure (see, for example, [1] and references therein). Determining the location of the anomaly as a function of pressure allows for the construction of a locus of the turning points (often displayed in the pT or ρT planes) and termed a temperature of maximum density (TMD) line. Whilst the density anomaly is perhaps the most widely-studied, further anomalies (both maxima and minima) are observed in other related thermodynamic properties such as the heat capacity and isothermal compressibilities as well as in the diffusivities [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. The origins and relationships between these anomalies are complex. At one level fundamental thermodynamic 'rules' may be derived which govern how specific anomalies interact with each other and with related properties such as stability limits or critical points [3,4,[20][21][22][23][24][25][26][27][28][29]. However, these rules say nothing about the structural origin of the anomalies themselves. In the systems listed above these appear to arise from the subtle disruption of a basically tetrahedral network as a function of temperature and/or pressure. To try to understand the structural origin of the anomalies two basic modelling strategies have been employed. In the first accurate models are employed which model specific systems of interest, for example, SiO 2 [6,7,15,18], BeF 2 [16][17][18], GeO 2 [18] and H 2 O [1, 3-5, 18, 30-32]. In the second generic models are employed, for example using a Stillinger-Weber potential [33,34], or a ramp [15,[35][36][37] or core-softened [38] potential. In the Stillinger-Weber (SW) potential, for example, the energy of interaction is deconvoluted into a sum of twoand three-body interactions whose relative magnitudes may be systematically varied. Recent work has also highlighted the sensitivity of the location of critical points to the potential model, focusing on both colloidal systems [39] and modified water potentials [40].
Understanding the structure of systems dominated by electrostatic (or Coulombic) interactions has a rich history [41][42][43][44]. At short length-scales the ionic ordering imposed by the electrostatic interactions promotes a relatively simple structure in which ions with opposite charges sit in nearest-neighbour 'shells'. As a result, the coordination number is effectively determined by the relative ion sizes (commonly referred to as radius ratio rules). The three systems listed above, for example, are dominated by tetrahedral local coordination polyhedra. However, the manner in which these coordination polyhedra connect (via M-X-M anion 'bridges') may differ with many-body (ion polarization) interactions controlling the key M-X-M bond angles, θ MXM [45,46]. The primary mechanism for this control may be summarised as follows. Dipoles induced on anions connecting a pair of cations introduces negative charge bisecting those cations and effectively screening the cation-cation repulsive electrostatic interaction [47]. The magnitude of the induced dipoles depends on the local electric field and the dipole polarizability. If the anion polarizability is increased the mean M-X-M bond angle becomes more acute until edge-sharing tetrahedra become stabilised [45,46,48]. This control of the bond angles then effectively controls the system topology. Furthermore, in relatively simple ionic models this means that these bond angles are controlled by a single model parameter, namely the (anion) dipole polarizability, α (the response of the ion electron density to an applied electric field). As a result, generic models can be interrogated in which α is systematically varied allowing the relationships between (chemically-related) systems to be more widely appreciated [48][49][50][51]. The control of the bridging bond angle affects the system topology (for example, as seen in the distribution of ring sizes [48]). BeF 2 and SiO 2 may reasonably be considered near-isomorphous in the sense that these key bridging bond angles are very similar.
The control of the overall topology means that these systems may display a wide range of dynamic, mechanical and structural properties. For example, BeF 2 , GeO 2 and SiO 2 are strong liquids (Angell's classification [52], with fragility indices of m = 20-28 [53]) as their viscosities change on cooling following an Arrhenius behaviour and show relatively obtuse M-X-M bond angle distributions dominated by corner-sharing unitsθ MXM = 155-120 • [54,55]. Systems such as ZnCl 2 , ZnBr 2 [56], GeS 2 [57] and GeSe 2 [58] are more fragile liquids (m = 30-60) with more acute bond angles, θ MXM = 120-85 • [54,[59][60][61][62]. The inclusion of the manybody effects may also control the system dynamics, acting to significantly soften (in particular) bending modes centred on the anion bridges (see, for example, [63] which considers the effects on these modes in the context of the infrared spectrum of SiO 2 ).
In this paper we make use of an electrostatic model for BeF 2 whilst systematically varying the anion polarizability in order to follow the progression of any observed thermodynamic anomalies. The use of a generic model allows us to relate any changes in structure to any shift on the phase space of the anomalies and will help connect with previous work which considered separate (more accurate) models for BeF 2 , GeO 2 and SiO 2 [18].

A summary of anomalous behaviour
A mathematical analysis confirms that there exists a number of constraints which govern the manner in which anomaly loci (and related properties such as critical points and stability limits) may interact or intersect (see, for example, [20] and references therein). In the present work the density and compressibility anomalies are considered along with the stability limits, and so we shall summarise their key interactions here. Maxima in temperature are referred to as TMD and maximum compressibility (TMC) respectively, whilst the corresponding minima are termed temperature of minimum density and compressibility (TminD and TminC) respectively.
In summary, • The density anomalies may or may not collide with the stability limit. If there is a collision then the gradient (dP/dT) of the stability limit in the pT plane at the collision state point is zero [20]. • A compressibility anomaly intersects a density anomaly when the latter has an infinite gradient in the pT projection ((dP/dT) TMD = ∞) [4].
In addition, we note that the singularity free interpretation suggests three alternative scenarios [4] which differ in terms of the relationship of the compressibility anomaly locus (specifically the point at which the TMC and TminC loci merge) as it intersects the density anomaly locus. The case in which compressibility anomaly locus intercepts the TMD at a positive TminC gradient is termed TEC-I, whilst intersection with a negative gradient is termed TEC-III. If the TMD is intersected at a point of inflection where the TMC and TminC merge, then this is termed TEC-II.

Potential models
The potential model used is fully described in [64]. The short-range interactions are accounted for using a modified Born-Mayer [65] potential (see [66] and references therein) of the form, where B ij and a ij are parameters controlling the ion radii and the rate of decay of the repulsive wall respectively, C i j n and f n (r ij ) are the dispersion parameters and respective damping functions (the latter mimicking the loss of asymptotic behaviour in the dispersive interactions at short range [67]). Damping functions of the form suggested by Tang and Toennies [68] are used throughout.
The pair potential shown above may be augmented by a description of ion polarization which introduces a manybody character. Historically, this has been incorporated using a shell model [69] which is a mechanical representation of the dipolar charge displacement comprising a charged shell connected by a harmonic spring to a charged core. In the polarizable-ion model (PIM) used here the induced moments are incorporated directly as simple point entities [47], requiring two sets of parameters: dipole polarizabilities, α i , and short-range damping parameters (SRDPs). The dipole polarizability describes the response of the ion electron density to the internal electric field arising from the presence of both formal ionic charges and other induced dipoles. The SRDP controls the effect of the overlap of nearest-neighbour electron density on the induced dipole moments [70][71][72]. The additional degrees of freedom representing the induced moments can be integrated in parallel with the ion equations of motion via a Car-Parrinello [73] method or by continuous energy minimisation [74]. Dipole polarizabilities can be obtained from ab initio electronic structure calculations (in which ions are studied in a condensed phase-an 'in-crystal' calculation) [74,75] or by reference to experimentally-determined refractive indices [76,77]. SRDPs can also be obtained from ab initio calculations by applying specific distortions to the nearestneighbour shell of counter-ions [70][71][72]. Alternatively, more extensive parameter sets may be obtained by 'force-fitting' methods. In these systems properties, such as atomic forces, cell stresses and induced multipoles, can be obtained from empirical (usually density-functional-based) electronic structure calculations. The potential parameters can then be systematically varied in order to obtain the best fit for the model to these properties.
For the MX 2 stoichiometry PIMs have been developed for GeSe 2 [78], GeO 2 [79], ZnCl 2 [78,80] and SiO 2 [81]. For the present work we use the model for BeF 2 [64,82] which was parameterised using a force-fitting method using typical crystal and liquid configurations with forces and multipoles obtained from density-functional calculations. The anion polarizability, α F − , in the full PIM is α F − = 7.09au [64]. Anion polarizabilities are dependent upon the local coordination environment. Ions such as F − are stable as free ions (with α FREE F − ∼16.8 au [83]) and are compressed in the condensed phase environments by both the electrostatic potential and overlap with nearest-neighbour ion electron density [70]. The value of α F − = 7.09au represents a mean average value for key fluorides at their respective equilibrium lattice parameters in the ground state crystal structure [71,84]. The polarizability of the Be 2+ cation is vanishingly small (α Be 2+ = 0.052au [85]). Furthermore, the cation sits at the centre of a tetrahedron of anions which precludes the formation of electric fields by symmetry. As a result, cation polarization can be safely neglected and a single (anion) polarizability applied (hereafter simply referred to as α). Dispersion terms are parameterised from the ion polarizabilites with the small cation value allowing the anion-cation and cation-cation terms to be set to zero.

Methods
Molecular dynamics simulation are performed in the NVT ensemble using a system of 576 MX 2 molecules. Initial configurations are generated by melting an α-quartz crystal configuration constructed from 4 × 2 × 4 unit cells. The equations of motion were integrated with a time-step of 25 au (∼0.6 fs) with constant temperature maintained using the Nosé-Hoover thermostats [86,87]. The starting crystal configuration was melted at a temperature of T = 5000 K, well above the estimated melting point, and then reduced to the required temperature and equilibrated for 50 000 time steps (≈30 ps). The production runs of 1700 000 time steps (≈1 ns) for the rigid ion model (RIM) and 600 000 time steps (≈360 ps) for the PIM were started after equilibration. The PIM simulations are roughly an order of magnitude slower than the RIM counterpart as a result of evaluating the many-body interactions associated with the induced dipole moments. In the PIM the potential energy must remain minimised with respect to the additional degrees of freedom which describe the induced dipoles. Here the potential energy was minimised using a steepest descent algorithm every five time steps in order to return the dipoles to the adiabatic surface. A frequency of 500 time steps was used for sampling thermodynamic data.
The thermodynamic phase space is explored by changing the system volume and temperature at different values of the anion polarizability, α. Simulations are performed using the 'full' ab initio 'in-crystal' polarizability of α = 7.09 au [71] as well as at lower values of α = 0, 1, 2, 3, 4, 5 and 6 au. Recall that the anion polarizability can be employed as a free parameter which effectively controls the M-X-M bond angles and hence potentially controls the network topology [48][49][50][51]. Setting α = 0 corresponds to the RIM in which the interactions are strictly pair-wise additive.
The thermodynamic data is analysed using the NumPy and SciPy packages [88]. Points on each isochore are filtered using a second order Butterworth filter and interpolated with a fourth order univariate spline allowing the first derivative to be calculated and identified as a maxima or minima, and which then corresponds to points on the TMD or TminD lines. An analogous procedure is employed to locate the compressibility anomalies. The isotherms are interpolated using a fourth order spline and the derivative determined to obtain the isothermal compressibility. The compressibility data is then re-cast on a pressure and temperature grid and extrema with respect to   figure 1 to highlight the shift in the anomaly loci relative to the stability limit as the anion polarizability, α, is increased. The black line shows Δp calculated from the minimum in the compressibility anomaly locus (at which point the TMC and TminC inter-change) to the stability limit at the same temperature. The red line shows the analogous pressure difference calculated with respect to the maximum temperature along the TMD locus. temperature at constant pressure are identified using an analogous procedure to that used to identify the density anomalies. No filtering was employed in order to avoid a potential underestimation of the compressibility peaks close to any critical points. The liquid-vapour stability limits were obtained by determining the lowest pressure point on each isotherm. This gives an upper bound of the liquid-vapour spinodal line in the pT plane. Figure 1 shows the anomaly loci and stability limits observed over the range of polarizabilities studied, in the pT projection. To allow for a more direct comparison, figures 2(a)-(c) shows the collective evolution of the TMD and TminC/TMC anomaly loci, as well as the stability limits, as a function of anion polarizability. Both the density and compressibility anomaly loci progress towards lower temperatures and slightly higher pressures as the polarizability is increased, with the increase in pressure almost negligible for the compressibility anomalies, but somewhat larger for the density anomaly. The density anomaly and TMC lines vanish for α = 7.09 au but are present at all the lower polarizabilities studied (see below). The stability limits show significant shifts to both lower temperatures and higher pressures as the polarizability increases. The general shape of the TMD loci remains approximately unchanged as the polarizability is increased showing the standard negative gradient in the pT projection at high pressure which turns at lower pressures towards a positive slope. The maximum temperature along the TMD loci (at which point ∂ p ∂T TMD = ∞) decreases with α (see below).

Results
The TMD line avoids collision with the stability limit for all values of the polarizability studied. Figure 2(d) shows the pressure difference, Δp, between specific anomaly loci and the stability limit. Δp is determined both with respect to the maximum temperature along the TMD loci and the lowest pressure observed for the compressibility anomaly (at which point ∂ p ∂T TMC/TminC = 0) at which point the TMC and TminC loci interchange. In both cases the pressure difference increases with polarizability showing the anomalies to be shifting away from the stability limit despite the shift to higher pressure (and lower temperature) seen with the stability limit over the same parameter space.
The TminC line is present for all polarizabilities studied, slowly becoming completely positively sloped as the polarizability increases. At low polarizability the TminC locus shows a characteristic 'S' shape observed in SW liquids [89]. The collision between the density anomaly and compressibility anomaly occurs at the infinitely sloped part of the TMD in the pT projection (figure 1) as required [4,20]. Sastry et al [4] defined the three possible scenarios which follow from the intersection of the density and compressibility anomalies and the crossover TMC ↔ TminC. For the RIM (α = 0) the crossover at a temperature just below the TMD locus (corresponding to TEC-I in [4]) is observed. As the polarizability increases the TMC ↔ TminC crossover shifts to lower temperature (see discussion). Figures 1 and 2 indicate that the locations of the anomaly loci (and related stability limits) in the pT plane are highly sensitive to the magnitude of the anion polarizability. One possibility is that the inclusion of the anion polarizability is significantly changing the network structure for the reasons discussed in the introduction. Figure 3 shows the evolution of the pair distribution functions, g αβ (r) (αβ = FF, BeF, BeBe), as a function of the anion polarizability. To try to form a meaningful comparison each distribution function is evaluated at the respective state point corresponding to the maximum temperature in the TMD anomaly locus. The inclusion of the anion polarizability has a relatively small effect on the static structure (as highlighted previously for the approximately isomorphous SiO 2 [63]). The most notable changes observed here is in the position of the cation-cation pair distribution function, g BeBe (r), which shows a small shift to a shorter length-scale as the anion polarizability is increased. The inset to figure 3 shows the position of the first peak in g BeBe (r) as a function of the anion polarizability. The peak shifts from r BeBe ∼ 3.10Å at the rigid-ion limit to r BeBe ∼ 3.01Å at α = 7.09 au, corresponding to a change in the Be-F-Be bond angle from θ BeFBe ∼ 154 • to θ BeFBe ∼ 143 • , assuming a fixed Be-F nearest-neighbour length-scale of r BeF ∼ 1.59Å. BeF 2 is effectively a charge-ordered system in that the basic structures adopted in the condensed phases can be rationalised in terms of a simple electrostatic ionic model in which the ions take their full formal valence charges of Be 2+ and F − . Both the crystalline and typical liquid (and amorphous) structures are dominated by nearest-neighbour anion-cation interactions (for which the radius-ratio rules translated to local tetrahedral order at ambient pressure). Furthermore, the cation-cation nearest-neighbour length-scale is longer than the anion-anion analogue, simply reflecting the larger formal cation charge and hence stronger repulsive forces between the pair of like ions. As noted earlier in the PIM the anion polarizability effectively controls the cation-cation nearest-neighbour length-scale. Furthermore, even for relatively small anion polarizabilities, in which the static structure may be relatively unchanged when compared to that generated with an RIM, the vibrational modes (and associated ion dynamics) may be significantly affected (see, for example, [63] for a discussion of the effects in the infrared spectrum of SiO 2 ).
Previous work has used the tetrahedral order parameter, Q 4 , as a useful metric [18], which is given by with the sum over the four nearest-neighbour cations about a given central cation. A local tetrahedral configuration gives Q 4 = 1. Figure 4 shows the mean average Q 4 over a range of isotherms as a function of the reduced density and the anion polarizability, α. The curves at each value of α are qualitatively similar showing maxima at ρ * ∼ 0.6-0.8 and which shift to higher density as the temperature is increased. The densities are shown in a reduced form, with ρ * = n 0 σ 3 , where n 0 is the number density. The magnitude of the shift becomes greater with increasing α. Figure 5(a) shows the values of Q 4 along isotherms corresponding to the highest temperature at which a TMD is observed at each polarizability, that is, corresponding to the temperature at which ∂T ∂ p TMD = 0 in figure 1 and ranging from T = 2700 K for the RIM to T = 1830 K for α = 6 au. As the anion polarizability increases the maximum in Q 4 shifts to higher density and increases in value, consistent with the simple reduction in temperature at which the TMD loci are observed. Figure 5(b) shows the density at the maximum  themselves both shown as a function of the reduced temperature as described in the text. In both panels the solid lines (no symbols) are for the present work with each line corresponding to a single value of the anion polarizability. Key: α = 0 (RIM), 1, 2, 3, 4, 5, 6, 7.09 au-black, red, green, blue, yellow, cyan, magenta and orange respectively. Both panels also show data from [18] for BeF 2 (black line, ×), SiO 2 (green line, ) and GeO 2 (red line-•).
value of Q 4 (< Q max 4 >) as a function of temperature reduced by the experimental melting point (T m = 813 K-see discussion). At low α there is a minimum in ρ * (< Q max 4 >) as a function of T which shifts to lower T as α increases before becoming suppressed into the supercooled regime at α ∼ 4 au. The density at which the maximum is observed also shifts to higher values as α increases. Figure 5(c) shows the value of < Q max 4 > as a function of the (reduced) temperature for each polarizability. The value of < Q max 4 > shows a weak maximum at low α which vanishes at α ∼ 3 au. The values fall as the temperature increases (as would be expected given the greater degree of local disorder).

Discussion
The results outlined in the previous section highlight how sensitive the locations of the thermodynamic anomalies are to the fine details of the potential models. We now consider these results in the context of previous work, both focused on the MX 2 stoichiometry, and more broadly to monatomic systems  [18,93] respectively. The lower panel shows the results from [18] for the related network-forming systems SiO 2 , GeO 2 and H 2 O, as well as for the BeF 2 TRIM from both [18,93] (black and cyan respectively). modelled with a Stillinger-Weber potential, the latter having been more extensively studied.

Comparison to previous MX 2 results
The results presented in the previous section indicate that the inclusion of the anion polarization has the potential to greatly influence the behaviour of the anomalous loci in ionic melts. Here, for example, the maximum temperature along the TMD locus falls from T ∼ 2750 K for the RIM to T ∼ 1800 K for α = 6 au. Above this polarizability no further TMD loci can be detected. It should be noted that a possibility of finding the TMD locus at α = 7.09 au still exists as we were unable to avoid the onset of amorphisation using standard MD simulations. It is possible that the TMD locus is suppressed inside this amorphous region, but to prove this we would have to employ enhanced sampling methods such as replica exchange [90,91] and avoid crystallisation at the same time [20]. Furthermore, as shown in figure 1, even at α = 7.09 au the compressibility anomaly locus is still observed (as this appears at higher temperature than the TMD). Given the known relationships between the compressibility and density anomalies then the observation of the former reflects an 'echo' of the latter, albeit present at too low a temperature to be easily observed. However, we note that it is perfectly possible to observe compressibility anomaly loci without an associated density anomaly [20]. Indeed, at the highest α studied in figure 1 it appears that the TminC locus may not be heading to zero gradient which would indicate that the density anomalies have ceased at this α.
In previous work on BeF 2 a transferable rigid-ion model (TRIM) was employed [92]. The TRIM and RIM employed here differ only in the details of the pair potential parameters (the former, for example, omits any dispersion terms, although these are relatively small given the small polarizabilities of the two ions) but are otherwise equivalent, for example, employing full formal ion charges. As a result, therefore, the TRIM and RIM employed here differ only in the (seemingly minor) details of the effective pair potentials. Even so, these subtle differences lead to a relatively large shift in the location of the TMD lines. Figure 6 shows the present TMD anomaly loci (presented in the ρT plane to make contact both with previous work and other chemically-related systems). The shape of the density anomaly generated by the TRIM matches that reported here. The TRIM shows a maximum TMD temperature of T * = 2.79 compared with T * = 3.34 for the RIM, where T * is the temperature reduced by the experimental ambient pressure melting point. The temperature range reported for the TRIM (from [18,93]) is comparable to our results, with best agreement shown at α = 3 au figure 6 also shows TMD data for other chemically-related systems, SiO 2 , GeO 2 and H 2 O (again from [18]). For each system the temperature is shown reduced by the respective experimental ambient pressure melting points (T m = 813 K, 1996 K, 1389 K and 273 K for BeF 2 , SiO 2 , GeO 2 and H 2 O respectively-see the references given in [18]). The densities are again shown in a reduced form, with ρ * = n 0 σ 3 , where n 0 is the number density. The maximum reduced temperatures attained along each TMD falls in the order BeF 2 → SiO 2 → GeO 2 → H 2 O from T * ∼ 2.8 for BeF 2 to T * ∼ 0.9 for H 2 O. Figure 7 shows the maximum TMD temperature (again normalised by the respective melting points) as a function of the M-X-M bond angles, θ MXM , for the systems studied here as well as data for SiO 2 , GeO 2 and BeF 2 from [18]. Bond angle data is taken from [54]. The error bars shown reflect the difficulties in obtaining accurate and reliable values for the bond angles. For BeF 2 , for example, θ MXM has been determined as θ MXM ∼ 155.6 • for the amorphous system at T = 25 • C [94] to θ MXM ∼ 144.6 • in the α−quartz crystal structure [95]. H 2 O is included in figure 6 for completeness although the network structure in that case is controlled by more complex interactions including hydrogen bonding in which quantum effects may become significant (for example, the temperature of the density maximum at ambient pressure varies from 4 • C to 11.4 • C to 13.4 • for H 2 O, D 2 O and T 2 O respectively).
There appears a clear and strong correlation between this temperature and the bond angle. As the bond angle is reduced below θ MXM ∼ 130 • the TMD locus is predicted to become suppressed firmly into the supercooled (metastable) regime. These results correlate with the lack of an observed TMD in, for example, ZnCl 2 (θ MXM ∼ 110 • [54]) although we again note the possibility of detecting echoes of the TMD via, for example, detection of a compressibility anomaly at higher temperature. Furthermore, these three systems are strong liquids (Angell's clssification) showing fragility indices in the range m = 20-28 [53]. As a result, the system viscosities would be expected to reduce only slowly on heating above the melting point. The above analysis indicates that the topology of the network (as probed, for example, by the bridging anion bond angles) is key in determining the location of the anomalies in phase space. Increasing the anion polarizability in the PIM is equivalent (in a coarse-grained sense) to moving from BeF 2 → SiO 2 → GeO 2 , with BeF 2 and SiO 2 near-isomorphous (the latter showing a marginally more acute bridging bond angle). Figure 5 shows the densities at the respective maximum values of the mean tetrahedral order parameter, Q 4 , (panel(b)) and the values at that maximum, < Q max 4 > (panel(c)), as a function of the anion polarizability. The figure also shows the results from [18] for BeF 2 (obtained using a TRIM), SiO 2 and GeO 2 . The BeF 2 TRIM density matches that of the current PIM with α ∼ 3 au (corresponding to the near-matching of the respective TMDs-see figure 7). The SiO 2 data is best matched by the PIM at a slightly higher α (α ∼ 4-5 au) consistent with the slightly more acute bridging bond angles when compared with BeF 2 . GeO 2 behaves like a PIM with a significantly higher α albeit shifted to lower T. In terms of the values for < Q max 4 > the TRIM appears to favour more ordered structures (higher Q 4 ) than the RIM, again highlighting how subtle differences in the potential models can significantly affect the location of the anomalies in phase space. The reduction in < Q max 4 > for BeF 2 → SiO 2 → GeO 2 is matched by the PIM on increasing α. Overall, therefore, whilst the detail of the local coordination polyhedra (as controlled in models of the type employed here by the radius-ratio rules) must have some effect on the TMD (and broader anomaly) loci locations (as suggested in [18]) it is the connectivity between the local coordination polyhedra which is critical.

Comparison to Stillinger-Weber liquids
Another potentially useful comparison is to contrast the systems studied here with the (typically monatomic) systems studied using the SW potential. In the present work the anion polarizability is used as a free parameter which is shown to control the location in phase space of the anomaly loci. In the SW potential a single parameter may be varied in order to control the balance between the two-and three-body interactions and, as a result, control the location in phase space of the anomaly loci (see, for example, [89] and references therein). The most common parameter to vary is λ (see [20] and references therein), although varying the ideal local polyhedra angle, θ 0 , has an analogous influence on the underlying system properties [96]. Figure 8, for example, shows the evolution of the TMD loci as a function of (a) λ and (b) θ 0 (data from [89,96]). In the case of the MX 2 system studied here, an increase in the anion polarizability acts to suppress the TMD loci (shift to low temperature) whilst, for the SW liquids, a similar increase in both λ or θ 0 acts to shift the TMD loci to higher temperature and pressure. At high θ 0 the TMD becomes 're-entrant' in the sense that the highest TMD temperature achieved shifts back to lower T. Figure 9(a) shows the evolution of the maximum TMD temperature as a function of both α and λ, highlighting the similar behaviour observed for reducing α or increasing λ. Note that T/T m is typically lower for the SW liquid indicative of the greater fragility of these systems (compared with the relatively strong MX 2 liquids).
The use of a single melting point to normalise the temperature scales is, of course, an approximation as we might reasonably expect the melting point for the MX 2 crystals to be a Figure 9. (a) Maximum TMD temperature, reduced by melting point, for (left) the MX 2 system studied here, and (right) the SW liquid. For the left-hand panel the melting point for BeF 2 (T m = 813 K) is used, whist for the right-hand panel the corresponding temperature for silicon (T m = 1687 K) is used (red circles) and the melting curve obtained from [33] (blue triangles). The inset to the right-hand panel shows the liquid/crystal melting curve from [33]. (b) Temperature difference, ΔT, between the temperature at which the TMC and TminC interchange and the maximum TMD temperature shown for (left-hand panel) the MX 2 system considered here as a function of the anion polarizability, and (right-hand panel) for the SW liquid as a function of λ. The blue line shows the limiting value at which TEC-I scenario changes to TEC-III case via TEC-II case [4]. function (albeit as yet unknown) of the anion polarizability. Determining the full melting curve for the MX 2 systems is beyond the scope of the present work. We can, however, highlight the analogous change when considering the SW liquids for which the melting curve is known [33]. Figure 9(a) shows the maximum TMD temperature normalised by the melting point as a function of λ (from [33]-the melting curve determined in this reference is shown in the inset to the figure). The melting point rises strongly as a function of λ resulting in a very different behaviour to that assuming a fixed value. A potential additional difficulty is the potentially complex underlying crystalline phase diagram [33,96]. At low λ, for example, close-packed crystal structures become thermodynamically stable over the diamond [33]. Figure 9(b) shows the comparative behaviour of the temperature difference between the temperature at which the TMC and TminC interchange and the maximum TMD temperature. In the present work the interchange occurs exclusively at lower temperature (negative ΔT) commensurate with the TEC-I scenario [4]. As the anion polarizability increases ΔT becomes more negative, again analogous to the behaviour observed for decreasing λ for the SW liquids. In the latter case, as shown in figure 9(b), ΔT can change sign and so all three TEC scenarios are observed.
Further points of contact are as follows. The gradient of the high pressure region of the density anomaly in the pT projection progresses towards a more negative value as the polarizability is increased analogous to the progression observed in the SW system as a function of λ. In addition, the density anomaly weakens in both the PIM and SW liquids concomitant with the gradient of the TminC anomaly locus losing its characteristic 'S' shape and developing an infinite gradient in the pT projection (figure 2). The weakest density anomalies detected (observed for α = 6.0 au and λ = 27 respectively) coincide with the gradient of the TminC locus crossing infinity in the pT projection. Finally, at the lowest polarizability studied (α = 0, corresponding to an RIM) a second TminC branch emerges at high pressure (not shown) which replicates the behaviour observed in the SW liquid at low λ [20]. This observation might suggest that the PIM systems also show the 'loop-like' anomaly behaviour similar to the one seen in SW liquids. The low fragility of these systems may open the way to future replica exchange studies which may help to clarify the observed changes.

Conclusions
In this paper we have studied thermodynamic anomalies in both the density and compressibility, as well as the related stability limits, using an ionic model for BeF 2 which includes many-body polarization terms (confined to the anions only). BeF 2 represents an archetypal network-forming system whose structure can be rationalised in terms of connected local tetrahedral coordination polyhedra. By varying the anion dipole polarizability we have mapped the changes in the anomaly locations in phase space. As the anion polarizability is increased the anomalies and stability limits systematically shift to lower temperature and higher pressure whilst, at high dipole polarizabilities the density anomaly locus becomes suppressed into the supercooled regime. The locations of the anomalies appears most highly correlated with the (small) changes in the bond angles between neighbouring (linked) tetrahedra. These angles are effectively controlled by the anion polarizability in models of the sort applied here. Furthermore, the systematic evolution of a single parameter has allowed us to isolate the effect of the changes in system topology as opposed to applying more detailed models to a range of different (although physically-related systems) whose models may differ in a number of different ways making isolating the root cause of any changes problematic.
Comparison to previous work, however, makes it clear that relatively small differences in the potential models may lead to relatively large shifts in the anomaly locations. Finally, we note that the anion polarizability behaves somewhat analogously to the Stillinger-Weber parameter, λ, which controls the balance of the two-and three-body interactions in that model. Similar comments apply to other models. For example, systematic modification of the angular flexibility inherent in the ST2 model for water strongly affects the location of the liquid-liquid critical point, with the necessary changes in the locations of the related anomalous properties [40]. Analogous comments apply to the study of 'patchy' colloids with the angular terms again critical [39]. Overall, therefore, the ability to control both the network 'rigidity' (here effectively meaning the bond lengths) and flexibility (the angle between well-defined structural units) is critical in controlling the locations in the thermodynamic phases space of any anomalies and related properties.