A possible four-phase coexistence in a single-component system

For different phases to coexist in equilibrium at constant temperature T and pressure P, the condition of equal chemical potential μ must be satisfied. This condition dictates that, for a single-component system, the maximum number of phases that can coexist is three. Historically this is known as the Gibbs phase rule, and is one of the oldest and venerable rules of thermodynamics. Here we make use of the fact that, by varying model parameters, the Gibbs phase rule can be generalized so that four phases can coexist even in single-component systems. To systematically search for the quadruple point, we use a monoatomic system interacting with a Stillinger–Weber potential with variable tetrahedrality. Our study indicates that the quadruple point provides flexibility in controlling multiple equilibrium phases and may be realized in systems with tunable interactions, which are nowadays feasible in several soft matter systems such as patchy colloids.

W hen different phases are in thermodynamic equilibrium with each other at constant temperature T and pressure P, the chemical potentials of the phases must be equal. The number of equality relationships determines the number of degrees of freedom F. This leads to the famous Gibbs phase rule 1 : F ¼ C À N þ 2, where C is the number of chemically independent constituents of the system, and N is the number of phases. The rule should be valid, provided that the equilibrium between phases is not influenced by external fields and there is no spatial constraint on the phases. The latter condition is known to be violated for coherent solids 2 . This rule tells us that for a pure substance, it is only possible that three phases can exist together in equilibrium (N ¼ 3). For a one-component system, there are no degrees of freedom (F ¼ 0) when there are three phases (A, B and C), and the three-phase mixture can only exist at a single temperature and pressure, which is known as a triple point. The two equations m A (T, P) ¼ m B (T, P) ¼ m C (T, P) are sufficient to uniquely determine the two thermodynamic variables, T and P. Four-phase coexistence should then be absent, as three chemical potential equations admit no solutions when there are only two independent variables T and P. Mathematically, however, this does not necessarily rule out the possibility that the set of equations may be solved in a special case. Here we seek such a possibility in a systematic manner by tuning the interaction potential, or the Hamiltonian of the system. Extending the dimensionality of the system will allow us to investigate what are the conditions for the existence of a quadruple point.
The presence of a point where different phases coexist provides an interesting possibility of switching materials properties, including electric, magnetic, optical and mechanical properties, by a weak thermodynamic perturbation such as stressing or heating/cooling. The technological importance of a triple point has recently been shown for a popular candidate material for ultrafast optical and electrical switching applications, vanadium oxide (VO 2 ) (ref. 3): it has been revealed that the well-known metal-insulator transition in this material actually takes place exactly at the triple point. Large piezoelectricity near a morphotropic phase boundary is another important example of the importance of multi-phase coexistence [4][5][6] . In these systems, structural transformations in lattice order are coupled with other orders such as dipole, spin, charge and orbital, which can be used for applications such as electromechanical or magnetoelectronic devices. Although the role of multi-phase coexistence in the ease of the transition is not so clear, the minimization of the volume change associated with a phase transition may be realized by combined nucleation of two phases with different signs of the volume change on the transition, which has been reported for a transition near a ferroelectric-anitiferroelectric-paraelectric triple point 5 . So, the presence of a multiple point may provide a novel kinetic pathway of phase transition, for which the barrier for phase transformation is much lower than an ordinary phase transition between two phases. Thus, the fundamental understanding of multiple-phase coexistence is not only of scientific interest but also of technological importance.
To study the basics of multi-phase coexistence, we need a model system that shows rich polymorphism. In this context, it is well known that water exhibits a rich variety of crystal polymorphs (at least, 16 types of crystals 7 ). Motivated by this, here we study systems interacting with tetrahedral interactions (for example, covalent bonding and hydrogen bonding). Tetrahedral interactions are the most important category of directional interactions found in nature, both in terms of abundance, and in terms of unique physical properties. They are ubiquitous in terrestrial and biological environments, and fundamental for technological applications. The disordered (liquid) phases of tetrahedral materials show unique thermodynamic properties, the most important being water's anomalies 8,9 , like the density maximum, the isothermal compressibility and specific heat anomaly and so on. Ordered phases of tetrahedral materials are of fundamental importance in industrial application, as they include open crystalline structures, like the diamond cubic (dc) crystal, or the quartz crystal, with unique mechanical, optical and electronic properties. For example, in Si and Ge, the diamond cubic (dc) crystal is a semiconductor, whereas the liquid and body-centred cubic (BCC) crystal are metals. Furthermore, dc crystals of mesoscopic particles (like colloids) are also a promising candidate for photonic crystals 10 . It is thus not surprising that tetrahedral interactions are one of the focus of nanotechnology, with the aim of producing new generation of materials with properties that can be finely controlled by design [10][11][12] .
To understand the bulk behaviour of materials with tetrahedral interactions, several coarse-grained models have been introduced. Among them, probably the most important and successful model is the Stillinger-Weber (SW) potential 13 , in which tetrahedrality is enforced with the use of three-body force terms. Originally devised as a potential for silicon, the model has found widespread applicability for several materials, especially group XIV elements, like germanium and carbon. The key parameter controlling the tetrahedrality of the model is the ratio between the strength of three-body interactions over two-body interactions, often referred to as l. As tetrahedrality becomes less strong with increasing atomic number, the basic idea is that group XIV elements, apart from energy and length-scale differences, can be modelled by simply varying l (ref. 14). Even more importantly, the modified SW potential has found general application as a coarse-grained model for molecular and supramolecular systems. The most important example is water, whose structural properties have been accurately reproduced with a parametrization of the SW potential (called mW water) with a precision that is competitive (if not superior) to the best classical molecular models available to date [15][16][17] .
Despite the importance and widespread applicability of the SW model, our knowledge of its phase diagram is still lacking. Determining the phase diagram is challenging because of the three-dimensional parameter space (temperature, pressure and l), and the fact that all calculations are multiplied by a large number of crystalline structures with local tetrahedral symmetry that have to be tested for thermodynamic stability. Previous attempts have considered stable crystalline structures taken from the elements that the SW potential is ought to describe, as for example silicon. The first study of the model as a function of l was introduced in ref. 14, where three crystalline structures were identified at zero pressure P: BCC, b-tin, and dc, respectively for low (lt18), intermediate (18tlt19) and high values (l\19) of l. Interestingly, the intermediate region showed increased glass-forming ability 14 . The b-tin phase was believed to be the high-pressure phase for SW silicon (l ¼ 21) (ref. 18). So, according to the current view, starting from a perfectly tetrahedral diamond (dc) phase, the SW system would transform into b-tin by reducing the amount of global tetrahedrality, either by applying pressure or decreasing the value of l. But this view was recently proven wrong when a new crystal of SW silicon was found, sc16, which replaces b-tin as the stable phase 19 at high pressures. sc16 is a new crystal with a simple cubic unit cell and 16 atoms in the unit cell. This calls for a new understanding of the phase behaviour of the SW model, with new behaviour that should emerge in between the low P-l region (where b-tin is stable) and the high P-l region (where sc16 is stable).
In this Article, we show that indeed the SW model exhibits a behaviour, which takes the form of a 'quadruple point' (QP), where the fluid and three different crystalline structures (dc, b-tin and sc16) have the same chemical potential. We also show that the newly found quadruple point is stable against liquid-gas phase separation, by computing the liquid-gas coexistence line and critical point. At the quadruple point, l takes the value l QP B20.08 and the model describes a one-component system with four-phase coexistence.

Results
Modified Stillinger-Weber model. To compute the phase diagram, we run Monte Carlo simulations in the isothermal-isobaric NPT ensemble 20 . The SW potential can be written as the sum of a pairwise term U 2 and three-body interaction term U 3 : Here U 2 models a steep repulsion at short distances and a shortrange attraction, while U 3 is a directional repulsive interaction that promotes tetrahedral angles between triplets of particles (for the analytic expressions of these terms, see Methods). l is a dimensionless parameter controlling the relative strength between pairwise and three-body term. Free energy calculations of all relevant crystalline structures are conducted with the Einstein crystal method 20 , and both Gibbs-Duhem integration 21 and Hamiltonian integration 22 are employed to compute coexistence planes and triple lines. Critical points are estimated with grand canonical simulations and histogram reweighting techniques 23 .
A description of all methods can be found in Methods and from here we use internal units as explained there.
Phase behaviour of the modified SW model. We start from a liquid phase and four crystalline phases which are known to be stable for the SW model. The crystalline phases are body-centred cubic (BCC), b-tin, diamond cubic (dc) and sc16. BCC is known as a stable crystalline phase of the SW model at lower l (refs 14,15). b-tin crystal has a body-centred-tetragonal structure with two atoms per cell and is known as a stable crystalline phase for silicon at intermediate pressure 24 . dc is known as a stable crystalline phase for group XIV elements. sc16 is a crystal which has recently been found to be stable at intermediate and high pressure 19 Fig. 1a. To aid the visualization, we also plot in Fig. 1b a projection of the coexistence surfaces on the (P, l) plane. Each surface represents a coexistence surface between the liquid and the corresponding crystal. Thick lines are triple lines, where two crystalline phases and the liquid phase coexist. The order of the different crystals is as follows: BCC at low l; b-tin at intermediate l; dc and sc16, for low and high pressures, respectively, at high l. The dot in Fig. 1 highlights a quadruple point at the intersection of three triple lines. At the quadruple point dc, b-tin, sc16 and the liquid phase all coexist at the same T QP and P QP . The coordinates are approximately: l QP ¼ 20.08, T QP ¼ 0.042 and P QP ¼ 0.120. Incidentally, we note that l QP is very close to the value of l for the SW model of Germanium (l ¼ 20). We have checked our results with direct-coexistence simulations 25 , in which each crystalline phase is placed in contact with the fluid phase and shown to be at coexistence.
Our results confirm that the sc16 is indeed the stable crystalline phase at high (l, P) and show that it shares a triple line with the previously known b-tin phase down to the quadruple point, where the b-tin transforms directly into dc. We have further confirmed the stability and relevance of the sc16 crystalline phase by direct nucleation events, and showed that the fluid phase directly crystallizes in the sc16 phase at l ¼ 21, P ¼ 0. Unlike T and P, l is not a thermodynamic variable but a parameter of the Hamiltonian (equation (1)). By choosing l ¼ l QP we thus have a system with a stable quadruple point in its phase diagram. To show this we compute the phase diagrams for l ¼ l QP in the P-T and r-T planes, plotted respectively in Fig. 2a,b. The dc phase is stable at lower P and the b-tin phase at intermediate P. The stable region of the sc16 is instead split into two regions, at lower and higher P. Later, we discuss how the quadruple point emerges when l-l QP .
In Fig. 2b, we show the densities of the different crystalline states. Diagonal lines represent coexistence regions between two different phases, while horizontal lines are plotted at the temperatures of the triple points and the quadruple point. dc is the phase with lowest density, lower than the fluid's density, as already could be inferred by the slope of the P, T coexistence line in Fig. 2a. The sc16 crystal can coexist with the fluid phase at two different densities: at the quadruple point with a density lower than the liquid's density, and at higher P with a density higher than the liquid's density. To summarize our study of the liquid-solid phase diagram, we report all thermodynamic values at triple and quadruple points in Table 1. Figure 2 shows also the presence of a new phase (denoted as X) that we found while computing the coexisting line between b-tin and sc16 at the l QP . The phase spontaneously forms from b-tin, with which it shares similar densities. As this is a low T phase that does not coexist with the liquid, and that lies well below the quadruple point, we have not focused on identifying it. In Supplementary Fig. 3 and Supplementary Note 2, we report our preliminary studies on this new phase, and leave to a future work the task of determining which crystal it represents.
Location of the gas-liquid phase transition. The phase diagrams obtained so far only include solid and liquid phases, so there is the possibility that the quadruple point is metastable with respect to liquid-gas phase separation. To exclude this possibility, we have computed the liquid-gas critical point and coexistence lines for l ¼ l QP . To obtain the critical point, we conducted grand canonical simulations to get the distribution function of the mixing order parameter M (M ¼ r þ su; r is density and u is internal energy per particle and s is mixing parameter), and use histogram reweighting to find the state point where this distribution matches the one from the Ising universality class 23,26 . The results for the critical point is P CP ¼ 0.004, T CP ¼ 0.321. The gas-liquid phase diagram (the critical point and the coexistence line) is shown in Fig. 3, where it is clear that the liquid-gas critical point is located at pressures two orders of magnitude lower than P QP . Therefore, the quadruple point is indeed a stable thermodynamic point of the model.
Emergence of the quadruple point. It is of particular interest to observe how the quadruple point emerges as a function of the parameter l in two dimensions, where the Hamiltonian is given by equation (1), and P and T are the only intensive thermodynamic variables. To do this, we consider the Gibbs surface, which expresses all the thermodynamic information on the system as an energy surface u ¼ u(s, r) in the entropy s and density r plane. In this representation each thermodynamic phase is represented as a surface, whose tangents are the temperature, T ¼ (qu/qs) r , and pressure, P ¼ r 2 (qu/qr) s /N. Triple points are represented as triangles whose vertices lie on the pure phases surfaces. Slightly below l QP , we have two triple points (fluid/dc/b-tin triple point and dc/sc16/b-tin triple point). Slightly above l QP , on the other hand, we have two different triple points (fluid/dc/sc16 triple point and fluid/b-tin/sc16 triple point). By increasing l continuously from below to above l QP , the Gibbs surface shows the following change. Below l QP (Fig. 4a), two triangle areas corresponding to the fluid/dc/b-tin and dc/sc16/b-tin triple points sandwich the dc/b-tin coexistence surface. With an increase in l, this surface continuously becomes narrower, eventually becomes a line, until the two triangles merge forming a quadrangle surface. This quadrangle is made of coplanar points and is the Gibbs representation of a quadruple point. Figure 4b shows the quadruple point as obtained from thermodynamic calculations. A further increase in l, leads to the splitting of the quadrangle to two triangles corresponding to the fluid/dc/sc16 triple and fluid/b-tin/sc16 triple points (Fig. 4c). The splitting of the quadrangle into two triangles now occurs along the fluid/sc16 coexistence surface. The entire process is schematically depicted in Fig. 4d. In two dimensions, the quadruple point is thus formed by the merging of two pairs of triple points located on the same plane in the Gibbs surface. The degeneracy in the number of degrees of freedom (four phases, but only three equations to determine the volume of each phase) means that there is no lever rule, and the volume of each phase is not determined by bulk properties alone. The degeneracy is removed if we consider the system in three dimensions, where l is   treated like a thermodynamic parameter (equivalent to an external field). In this case the Gibbs surface is defined in the three dimensional space of r, s and u 3 (which is the sum of all three body contributions in equation (1)). The quadruple point will be represented as a triangular pyramid, and the volume of each phase at four-phase coexistence is not only determined by the total volume and entropy, but also by the total three-body energy, if we constrain it.
Physical properties associated with the quadruple point. We next examine some of the unique physical properties associated with the quadruple point. The simulation shown in Fig. 5a was started by interfacing the cubic face of each crystal with a fluid slab, and equilibrated in the isobaric ensemble where the pressure perpendicular to the interfaces is kept fixed at P QP . It is important to ensure that the size of the interfacial plane is commensurate with an integer number of unit cells for all three crystals at their equilibrium unit cell size. In our simulations, we use 11, 14 and 8 unit cells in each direction of the interfacial plane for b-tin, dc and sc16, respectively. With this choice we ensure that all crystals have the correct unit cell size within B2%. The simulations show that the system indeed feels the underlying proximity of a quadruple point, in a way that the free-energy differences between all four phases become negligible, and thus the four bulk phases coexist over practical time scales with free diffusion of all crystal/liquid interfaces (Fig. 5a).
Next, we show that the properties of a quadruple point can be exploited to gain a very fine control over the stability and number of crystalline phases. Figure 5e shows the liquid/solid lines for all crystalline phases, both stable (continuous lines) and metastable (dashed lines) ones. The quadruple point is the point where all these lines cross. This means that, by choosing thermodynamic conditions arbitrarily close to the quadruple point, we can stabilize systems with one, two or three crystalline phases. This is demonstrated in Fig. 5b-d, where a small change in temperature allows us to equilibrate systems with varying number of crystalline phases. The ability to tune the stability of a varying number of crystalline phases with just small changes in thermodynamic parameters is one of the most interesting properties of a quadruple point. Our simulations also point to the importance of the liquid layer during the solid-solid transition. By placing any couple of crystalline phases directly in contact we always observe the development of several interfacial liquid layers, despite the fact that the liquid is not a stable phase in these state points (Fig. 5c-e). Solid-solid transitions in our systems always proceed in two-steps, where the first step is the melting of the interfacial particles of one solid, and the second step is the nucleation of the second solid phase from this liquid layer. This scenario is similar to what recently observed in solid-solid phase transitions of colloidal systems 27 , and is due to the high interfacial cost of forming a solid-solid interface.

Discussion
Quadruple points have not been found before in singlecomponent systems, and the reason behind their improbability is embodied in the Gibbs phase rule. Each phase corresponds to a different thermodynamic function (in the isothermal-isobaric ensemble these functions are the chemical potentials of each phase), and a quadruple point requires the equality of these functions at a single (T, P) state point. Since the equality of four different chemical potentials is equivalent to a system of three equations in two variables (T and P), the system should have no solution, and genuine quadruple points should not occur in onecomponent systems. The way this limitation is circumvented in our system is by promoting the parameter l to an independent thermodynamic variable, so extending the dimensionality of space to three dimensions. This allowed us to explore the full (T, P, l) space (Fig. 1), where we were able to find the condition on l for the existence of the quadruple point, where the four phases have the same chemical potential. But, besides the extension of the dimensionality of state space, a quadruple point also requires three different stable crystalline phases in a close region of state space, a condition which is rarely met in one-component systems. We note that the SW model itself did not meet this requirement until the recent discovery of sc16 (ref. 19).
Our work opens for the possibility of realizing a onecomponent system with a quadruple point in a practical sense, by finding the conditions at which the strength of the tetrahedral interaction matches the conditions we highlighted in this work. The most likely candidate for such a system are patchy particles, which are colloidal systems with functionalized patches on their surface [28][29][30][31][32] . Tuning the angular width of the directional interactions plays a role similar to the parameter l in the SW potential 33,34 , and allows for the tuning necessary to unveil four-phase coexistence.
Here we have determined the full three-dimensional phase diagram for the SW model. The relevant crystalline phases are dc, b-tin, BCC, sc16 (and the yet unidentified crystalline phase X). Apart from sc16 and X, the remaining phases have been confirmed experimentally for group XIV elements, for which the SW potential is a good coarse-grained model. The model also displays a quadruple point in the sense that the four phases have the same chemical potential there. For l ¼ l QP we have fully determined the (P, T) and (T, r) phase diagrams, and also computed the liquid-gas critical point, showing that the quadruple point is a stable feature of the phase diagram. Thanks to the development of technology, a systematic control of the Hamiltonian of a system is now realistic even in experiments: for examples, ordinary and patchy colloids with more than two types of interactions (see, for example, refs 35,36), proteins (see, for example, ref. 37), and application of optical and magnetic fields in quantum systems (see, for example, refs 38,39). In such a case, the Gibbs phase rule, which has been considered for two independent thermodynamic variables, should be extended by including an additional variable linked to the interaction potential: F ¼ C À N þ M, where M is a number of independent thermodynamic and Hamiltonian-related variables 40 . Although we have studied a case of a pure substance (C ¼ 1 and M ¼ 3), the above extension of the Gibbs phase rule is not limited to single-component systems.
Furthermore, by computing the phase diagram as a function of the tetrahedral parameter l our results can lead to better modelling of atomic systems or enable the design of novel coarsegrained potential for new generation materials with directional interactions (for example, patchy particles) [29][30][31][32][33] . Our results can also lead to a better understanding of tetrahedral materials which are arguably the most important class of materials in nature and technology. For example, our results reinforce the parallelism between pressure and frustration effects, which is an important principle in the understanding of water mixtures and their glass-forming ability (see, for example, ref. 41). In the SW model, l controls the degree of deviation from tetrahedrality, and the phase diagram in l has many points in common with V-shaped phase diagram of real elements in pressure 42 . For example, the l-T phase diagram at ambient pressure, resembles the the P-T phase diagram of tin (Sn), with its succession of grey-tin (dc), white-tin (b-tin) and BCC phases.
As recently shown in VO 2 (ref. 3) for the ultrafast insulatormetal transition 43 , high-order points can be used to design phasechange materials, where properties change rapidly by applying mechanical stress, heating/cooling, or even tuning the interaction potential (as shown here) by modifying internal degrees of freedom such as spin and electronic states with electromagnetic excitation (see, for example, ref. shown that we can indeed produce any of four different phases by arbitrarily small variations in T or P. Furthermore, the quadruple point in a single component system should provide great flexibility in controlling multiple phases. Here we study phase transitions where the density is the relevant order parameter, but they can be any types of order including dipole, spin, charge, and orbital order (see, for example, refs 45,46), which are important in functional materials.

Methods
Details of the model. The details of the Stillinger-Weber model are as follows.
The pairwise and three-body interaction terms are given by, i exp s r À as ð2Þ U 3 r ij ; r ik À Á ¼ E cosy ijk À cosy 0 Â Ã 2 exp gs r ij À as exp gs r ik À as ð3Þ where A ¼ 7.049556277, B ¼ 0.6022245584, p ¼ 4, q ¼ 0, cosy 0 ¼ À 1/3, g ¼ 1.2, a ¼ 1.8. The parameter E sets the energy scale and s the length scale. They correspond to the depth of the two-body interaction potential and the particle diameter, respectively, and determined by materials for which the model is used. We use internal units where E and s are the units of energy and length, respectively. Therefore, l is only parameter which differentiates the models.
Simulation methods. We compute internal energies and densities at each equilibrium state by performing Monte Carlo simulations. The size and shape of the simulation box can fluctuate so as to allow crystalline phases to change their structures 47,48 . A volume-change attempt occurs on every N translation attempts. The number of particles in the box is 1,024, which is large enough so that finite-size effects are negligible. We obtained the phase diagram by extending coexisting lines from phase diagrams at zero pressure and at l ¼ 23. 15 (ref. 19). We extend two-phase coexisting lines in two directions, along pressure axis and along l axis. We perform Gibbs-Duhem integration and Hamiltonian Gibbs-Duhem integration 21,22 to obtain coexisting lines along pressure axis and along l axis, respectively.
When we compute triple lines, we use the following relationships. For three phases 1, 2, 3 at coexistence, u 1 dp À s 1 dT þ @g 1 @l dl ¼ u 2 dp À s 2 dT þ @g 2 @l dl & ð4Þ u 1 dp À s 1 dT þ @g 1 @l dl ¼ u 3 dp À s 3 dT þ @g 3 @l where u, s and g are the volume, entropy and Gibbs free energy per particle. By solving this set of equations, we get Here, we write o ¼ @g @l : When the internal energy per particle u can be written as u ¼ u a þ lu b , we can use 20,22 o ¼ @g @l Here hyi N,p,T,l means ensemble average with constant N, p, T, l, which can be determined within an NPT simulation. By integrating these equations, we extend triple points.
To obtain critical points, we compute distribution functions of densities and energies with using histogram reweighting method 23 and fit them according to the universal Ising universality class 26 .
Data availability. The data that support the findings of this study are available from the corresponding authors upon request.