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

Abstract

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.

ARTICLE

Received 26 Apr 2016 | Accepted 15 Jul 2016 | Published 25 Aug 2016

For different phases to coexist in equilibrium at constant temperature T and pressure P, the condition of equal chemical potential m must be satised. 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 StillingerWeber potential with variable tetrahedrality. Our study indicates that the quadruple point provides exibility 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.

DOI: 10.1038/ncomms12599 OPEN

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

Kenji Akahane1, John Russo1,2 & Hajime Tanaka1

1 Institute of Industrial Science, University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan. 2 School of Mathematics, University of Bristol, Bristol BS8 1TW, United Kingdom. Correspondence and requests for materials should be addressed to J.R. (email: mailto:[email protected]

Web End [email protected] ) or to H.T.(email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599

When 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 rule1: 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 inuenced by external elds and there is no spatial constraint on the phases. The latter condition is known to be violated for coherent solids2. 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 mA(T, P) mB(T, P) mC(T, P) are sufcient 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 (VO2) (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 coexistence46. 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 ferroelectricanitiferroelectricparaelectric triple point5. 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 scientic 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 crystals7). 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 waters anomalies8,9, like the density maximum, the isothermal compressibility and specic 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 crystals10. 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 nely controlled by design1012.

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 StillingerWeber (SW) potential13, 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 modied

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 date1517.

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 rst study of the model as a function of l was introduced in ref. 14, where three crystalline structures were identied 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 ability14. 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 phase19 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 Pl region (where b-tin is stable) and the high Pl 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 uid and three different crystalline structures (dc, b-tin

2 NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599 ARTICLE

a b

Fluid

sc16

1

0.120.100.080.060.04

0.020.00

22 20

T

P

0.5

[afii9826]tin

dc

BCC

18

16

14

12

1.0

[afii9838]

0

0.5

14 16 18 20 22 [afii9838]

P

0.0

Figure 1 | Phase diagrams of a system interacting with the SW potential. (a) The lPT phase diagram. The green, pink, turquoise and orange surfaces are liquid-BCC, liquid-b-tin, liquid-dc and liquid-sc16 coexistence surfaces. The purple, red, yellow and blue lines are liquid-BCCb-tin, liquid-b-tindc, liquid-dcsc16 and liquid-b-tinsc16 coexisting lines. The brown point is a four-phase coexistence point for liquid, b-tin, dc and sc16. (b) The projection of the coexisting regions into lP plane. The green, pink, turquoise and orange regions are the projection of BCC-liquid, b-tin-liquid, dc-liquid, sc16-liquid surfaces into lP plane, respectively. The brown dot is the quadruple point.

and sc16) have the same chemical potential. We also show that the newly found quadruple point is stable against liquidgas phase separation, by computing the liquidgas coexistence line and critical point. At the quadruple point, l takes the value lQPB20.08 and the model describes a one-component system with four-phase coexistence.

ResultsModied StillingerWeber model. To compute the phase diagram, we run Monte Carlo simulations in the isothermalisobaric NPT ensemble20. The SW potential can be written as the sum of a pairwise term U2 and three-body interaction term U3:

U X

i

Xj4iU2 rij

l

Xi

:

1

Here U2 models a steep repulsion at short distances and a short-range attraction, while U3 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 method20, and both GibbsDuhem integration21 and Hamiltonian integration22 are employed to compute coexistence planes and triple lines. Critical points are estimated with grand canonical simulations and histogram reweighting techniques23. A description of all methods can be found in Methods and from here we use internal units as explained there.

Phase behaviour of the modied 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 pressure24. 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 pressure19. The sc16 crystal has a simple cubic unit cell with 16 atoms per cell. The space group of the sc16 crystal is Pa 3.

First, we show the three-dimensional phase diagram of the SW model, for lA[12.2:23.15], TA[0:1.12] and PA[0:1], in 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 TQP and PQP. The coordinates are approximately: lQP 20.08, TQP 0.042 and PQP 0.120.

Incidentally, we note that lQP is very close to the value of l for the SW model of Germanium (l 20). We have checked our

results with direct-coexistence simulations25, in which each crystalline phase is placed in contact with the uid phase and shown to be at coexistence.

Our results conrm 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 conrmed the stability and relevance of the sc16 crystalline phase by direct nucleation events, and showed that the uid phase directly crystallizes in the sc16 phase at l 21, P 0.5, and

T 0.0395. These results are shown in (Supplementary Figs 1

and 2) and discussed in (Supplementary Note 1).

Unlike T and P, l is not a thermodynamic variable but a parameter of the Hamiltonian (equation (1)). By choosing l lQP we thus have a system with a stable quadruple point in

its phase diagram. To show this we compute the phase diagrams for l lQP in the PT and rT 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-lQP.

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 uids 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 uid phase at two different densities: at the quadruple point with a density lower than the liquids density, and at higher P with a density

Xj 6 i

Xk4jU3 rij; rjk

NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599

a b

0.08

0.08

0.06

0.06

sc16

T

T

0.04

0.04

0.02

0.02

0 0 0.2 0.4 0.6 0.8 1 1.2P

0 0.4 0.45 0.5 0.55 0.6 0.65 [afii9845]

Figure 2 | PT and qT phase diagrams of the SW potential at kCkQP The dc, b-tin, sc16 and X phases are stable in blue, pink, orange and grey regions, respectively. (a) PT phase diagram. Here circle points are triple points, and a square point is a quadruple point. (b) rT phase diagram. Regions with diagonal lines are the coexisting regions of two phases. The horizontal lines indicate temperatures of four triple points and one quadruple point.

Table 1 | Quadruple and triple points.

0.40

0.35

T P qdc qbtin qsc16 qliquid qX

0.0093 0.171 0.486 0.536 0.5600.0211 0.199 0.553 0.535 0.5580.0244 0.831 0.583 0.609 0.5800.0420 0.120 0.475 0.544 0.518 0.5250.0680 1.063 0.585 0.610 0.581

T, P and r at triple and quadruple points for the system with l 20.08, as shown also in Fig. 2.

0.30

Liquid

T

0.25

0.20

0.15

higher than the liquids density. To summarize our study of the liquidsolid 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 lQP. 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 gasliquid 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 liquidgas phase separation. To exclude this possibility, we have computed the liquidgas critical point and coexistence lines for l lQP. 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 nd the state point where this distribution matches the one from the Ising universality class23,26. The results for the critical point is PCP 0.004, TCP 0.321. The

gasliquid phase diagram (the critical point and the coexistence line) is shown in Fig. 3, where it is clear that the liquidgas critical point is located at pressures two orders of magnitude lower than PQP. 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 r2(qu/qr)s/N. Triple points

are represented as triangles whose vertices lie on the pure phases surfaces. Slightly below lQP, we have two triple points (uid/dc/b-tin triple point and dc/sc16/b-tin triple point). Slightly above lQP, on the other hand, we have two different triple points (uid/dc/sc16 triple point and uid/b-tin/sc16 triple point). By increasing l continuously from below to above lQP, the

Gibbs surface shows the following change. Below lQP (Fig. 4a), two triangle areas corresponding to the uid/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 uid/dc/sc16 triple and uid/b-tin/sc16 triple points (Fig. 4c).

The splitting of the quadrangle into two triangles now occurs along the uid/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

0.10

108 107 106 105 104 103 102 101

P

Figure 3 | The liquidgas phase diagram of the SW potential at kCkQP A

circle point is a critical point. A broken brown line indicates the pressure of the quadruple point.

4 NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599 ARTICLE

a c

[afii9838]=19.95

[afii9838]=[afii9838]QP=20.08

[afii9845] [afii9845]

[afii9845]

[afii9845]

[afii9838]=20.75

Fluid

u

u

u

Fluid

1.7 1.75

1.8 1.85

1.9 1.95

2

1.7 1.75

1.8 1.85

1.9 1.95

2

1.75

[afii9826]-tin

1.8 1.85

1.9 1.95

2

[afii9826]-tin

sc16

sc16

2

3

4

5

6

dc

dc

2 34

5

S S

6

0.54

0.58

0.54

0.58

7

8

0.46

0.5

0.46

0.5

b d

S

S

S

[afii9838]<[afii9838]QP

[afii9838]=[afii9838]QP

[afii9838]>[afii9838]QP

Fluid -tin

-tin

[afii9826]-tin

sc16

2

3

4

dc

5

S 6

7

0.54

0.58

[afii9845]

0.46

0.5

Figure 4 | The emergence of the quadruple point in a single-component system. (ac) The Gibbs surface below (a), at (b) and above (c) the quadruple point. Skew surfaces correspond to pure phases, ruled surfaces to phase coexistence regions, and triangles to triple points. In b, we show that a quadrangle represents four-phase coexistence on the Gibbs surface. (d) Schematic gure showing how two triple points merge into a quadruplepoint and how it splits to form other two triple points as a function of l.

treated like a thermodynamic parameter (equivalent to an external eld). In this case the Gibbs surface is dened in the three dimensional space of r, s and u3 (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 uid slab, and equilibrated in the isobaric ensemble where the pressure perpendicular to the interfaces is kept xed at PQP. 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 ne 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. 5bd, 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. 5ce). Solidsolid transitions in our systems always proceed in two-steps, where the rst 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 solidsolid phase transitions of colloidal systems27, and is due to the high interfacial cost of forming a solidsolid interface.

NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599

a b

e

0.06

c

0.05

T

0.04

d

0.03

0 0.05 0.1 0.15 0.2 0.25 0.3P

Figure 5 | Direct coexistence near the quadruple point. (a) Simulation snapshot of four phases coexisting at the quadruple point: liquid (white), dc (cyan), b-tin (magenta) and sc16 (orange). (bd) Snapshots of congurations with one (dc), two (dc and sc16) and three (dc, sc16, b-tin) crystalline phases, respectively, obtained at the following state points near the quadruple point, P 0.1 and T 0.0445 (b), T 0.041 (c) and T 0.040 (d). (e) Phase

diagram at l lQP, where continuous lines are the stable thermodynamic lines of Fig. 2a, while dashed lines represent their metastable extensions. The

black dot is the quadruple point. The triangle symbols indicate the state points at which snapshots (bd) were taken from top to bottom. Note that the liquid layers remain between different types of crystals in c,d despite the liquid is not a stable phase there.

DiscussionQuadruple points have not been found before in single-component systems, and the reason behind their improbability is embodied in the Gibbs phase rule. Each phase corresponds to a different thermodynamic function (in the isothermalisobaric 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 one-component 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 nd 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 one-component system with a quadruple point in a practical sense, by nding 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 surface2832. Tuning the angular width of the directional interactions plays a role similar to the parameter l in the SW potential33,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 unidentied crystalline phase X).

Apart from sc16 and X, the remaining phases have been conrmed 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 lQP we have fully

determined the (P, T) and (T, r) phase diagrams, and also computed the liquidgas 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 elds 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 variables40. 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 coarse-grained potential for new generation materials with directional interactions (for example, patchy particles)2933. 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 pressure42. For example, the lT phase diagram at ambient pressure, resembles the the PT phase diagram of tin (Sn), with its succession of grey-tin (dc), white-tin (b-tin) and BCC phases.

As recently shown in VO2 (ref. 3) for the ultrafast insulator-metal transition43, high-order points can be used to design phase-change 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. 44). Thus high-order points open new directions in the control of materials properties. We have

6 NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599 ARTICLE

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 exibility 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 StillingerWeber model are as follows. The pairwise and three-body interaction terms are given by,

U2 r AE B

h i

exp s

r as

s r

p

s r

q

2

E cosyijk cosy0

2exp

3

where A 7.049556277, B 0.6022245584, p 4, q 0, cosy0 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 uctuate so as to allow crystalline phases to change their structures47,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 nite-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 GibbsDuhem integration and Hamiltonian GibbsDuhem integration21,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,

u1dp s1dT

@g1

@l

U3 rij; rik

exp

gs rik as

gsrij as

dl

4

u1dp s1dT

@g1

@l

dl u2dp s2dT

@g2

@l

dl

5

where u, s and g are the volume, entropy and Gibbs free energy per particle. By solving this set of equations, we get

dT dl

u1

dl u3dp s3dT

@g3

@l

6

dp dl

s1 s3

o1 o2

s1 s2

o1 o3

s1 s2

u1 u3

s1 s3

u1 u2

u3

o1 o2

u1 u2

o1 o3

u1

u3

s1 s2

u1 u2

s1 s3

7 Here, we write

o

@g @l

:

When the internal energy per particle u can be written as u ua lub,

we can use20,22

o

@g @l

ub l

h iN;p;T;l: 8

Here hyiN,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 method23 and t them according to the universal Ising universality class26.

Data availability. The data that support the ndings of this study are available from the corresponding authors upon request.

References

1. Gibbs, J. W. Scientic Papers: Thermodynamics Vol. 1 (Dover Publications, 1961).

2. Johnson, W. C. On the inapplicability of Gibbs phase rule to coherent solids. Metal. Trans. A 18, 10931097 (1991).

3. Park, J. H. et al. Measurement of a solid-state triple point at the metal-insulator transition in VO2. Nature 500, 431434 (2013).

4. Cox, D. E. et al. Universal phase diagram for high-piezoelectric perovskite systems. Appl. Phys. Lett. 79, 400402 (2001).

5. Ishchuk, V. M. & Sobolev, V. L. Investigation of two-phase nucleation in paraelectric phase of ferroelectrics with ferroelectric-antiferroelectricparaelectric triple point. J. Appl. Phys. 92, 20862093 (2002).

6. Ahart, M. et al. Origin of morphotropic phase boundaries in ferroelectrics. Nature 451, 545548 (2008).

7. Russo, J., Romano, F. & Tanaka, H. New metastable form of ice and its role in the homogeneous crystallization of water. Nat. Mater. 13, 733739 (2014).

8. Debenedetti, P. G. Supercooled and glassy water. J. Phys. Condens. Matter 15, R1669R1726 (2003).

9. Russo, J. & Tanaka, H. Understanding waters anomalies with locally favoured structures. Nat. Commun. 5, 3556 (2014).

10. Maldovan, M. & Thomas, E. L. Diamond-structured photonic crystals. Nat. Mater. 3, 593600 (2004).

11. Zhang, Z., Keys, A. S., Chen, T. & Glotzer, S. C. Self-assembly of patchy particles into diamond structures through molecular mimicry. Langmuir 21, 1154711551 (2005).

12. Bif, S. et al. Phase behavior and critical activated dynamics of limited-valence DNA nanostars. Proc. Natl Acad. Sci. USA 110, 1563315637 (2013).

13. Stillinger, F. H. & Weber, T. A. Computer simulation of local order in condensed phases of silicon. Phys. Rev. B 31, 52625271 (1985).

14. Molinero, V., Sastry, S. & Angell, C. A. Tuning of tetrahedrality in a silicon potential yields a series of monatomic (metal-like) glass formers of very high fragility. Phys. Rev. Lett. 97, 075701 (2006).

15. Molinero, V. & Moore, E. B. Water modeled as an intermediate element between carbon and silicon. J. Phys. Chem. B 113, 40084016 (2008).16. Moore, E. & Molinero, V. Structural transformation in supercooled water controls the crystallization rate of ice. Nature 479, 506508 (2011).

17. Lu, J., Qiu, Y., Baron, R. & Molinero, V. Coarse-graining of TIP4P/2005, TIP4P-Ew, SPC/E, and TIP3P to monatomic anisotropic water modelsusing relative entropy minimization. J. Chem. Theory Comput. 10, 41044120 (2014).

18. Kaczmarski, M., Bedoya-Martinez, O. N. & Hernndez, E. R. Phase diagram of silicon from atomistic simulations. Phys. Rev. Lett. 94, 095701 (2005).19. Romano, F., Russo, J. & Tanaka, H. Novel stable crystalline phase for the Stillinger-Weber potential. Phys. Rev. B 90, 014204 (2014).

20. Frenkel, D. & Smit, B. Understanding Molecular Simulation: from Algorithms to Applications Vol. 1 (Academic Press, 2001).

21. Kofke, D. A. Direct evaluation of phase coexistence by molecular simulation via integration along the saturation line. J. Chem. Phys. 98, 41494162 (1993).

22. Vega, C., Sanz, E., Abascal, J. & Noya, E. Determination of phase diagrams via computer simulation: methodology and applications to water, electrolytes and proteins. J. Phys. Condens. Matter 20, 153101 (2008).

23. Panagiotopoulos, A. Z. Monte Carlo methods for phase equilibria of uids.J. Phys. Condens. Matter 12, R25R52 (2000).24. Hu, J. Z., Merkle, L. D., Menoni, C. S. & Spain, I. L. Crystal data for high-pressure phases of silicon. Phys. Rev. B 34, 46794684 (1986).

25. Vega, C., Sanz, E., Abascal, J. L. F. & Noya, E. G. Determination of phase diagrams via computer simulation: methodology and applications to water, electrolytes and proteins. J. Phys. Condens. Matter 20, 153101 (2008).

26. Tsypin, M. & Blte, H. Probability distribution of the order parameter for the three-dimensional Ising-model universality class: a high-precision Monte Carlo study. Phys. Rev. E 62, 7376 (2000).

27. Peng, Y. et al. Two-step nucleation mechanism in solid-solid phase transitions. Nat. Mater. 14, 101108 (2015).

28. Wang, Y. et al. Colloids with valence and specic directional bonding. Nature 491, 5155 (2012).

29. Romano, F. & Sciortino, F. Patterning symmetry in the rational design of colloidal crystals. Nat. Commun. 3, 975 (2012).

30. Rovigatti, L., Bomboi, F. & Sciortino, F. Accurate phase diagram of tetravalent DNA nanostars. J. Chem. Phys. 140, 154903 (2014).

31. Starr, F. W. & Sciortino, F. Crystal-clear liquid-liquid transition in a tetrahedral uid. Soft Matter 10, 94139422 (2014).

32. Smallenburg, F. & Sciortino, F. Tuning the liquid-liquid transition by modulating the hydrogen-bond angular exibility in a model for water. Phys. Rev. Lett. 115, 015701 (2015).

33. Romano, F., Sanz, E. & Sciortino, F. Phase diagram of a tetrahedral patchy particle model for different interaction ranges. J. Chem. Phys. 132, 184501 (2010).

34. Saika-Voivod, I., Smallenburg, F. & Sciortino, F. Understanding tetrahedral liquids through patchy colloids. J. Chem. Phys. 139, 234901 (2013).

35. Yethiraj, A. Tunable colloids: control of colloidal phase transitions with tunable interactions. Soft Matter. 3, 10991115 (2007).

NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12599

36. Mnsson, L. K., Immink, J. N., Mihut, A. M., Schurtenberger, P. & Crassous, J. J. A new route towards colloidal molecules with externally tunable interaction sites. Faraday Discuss. 181, 4969 (2015).

37. Bhardwaj, N., Clarke, D. & Gerstein, M. Systematic control of protein interactions for systems biology. Proc. Natl Acad. Sci. USA 108, 2027920280 (2011).

38. Lukin, M. D. & Hemmer, P. R. Quantum entanglement via optical control of atom-atom interactions. Phys. Rev. Lett. 84, 28182821 (2000).

39. Zaccanti, M. et al. Control of the interaction in a Fermi-Bose mixture. Phys. Rev. A 74, 041605(R) (2006).

40. Vega, C. & Monson, P. A. Plastic crystal phases of hard dumbbells and hard spherocylinders. J. Chem. Phys. 107, 26962697 (1997).

41. Kobayashi, M. & Tanaka, H. Possible link of the V-shaped phase diagram to the glass-forming ability and fragility in a water-salt mixture. Phys. Rev. Lett. 106, 125703 (2011).

42. Tanaka, H. Simple view of waterlike anomalies of atomic liquids with directional bonding. Phys Rev. B 66, 064202 (2002).

43. Cavalleri, A. et al. Femtosecond structural dynamics in VO 2 during an ultrafast solid-solid phase transition. Phys. Rev. Lett. 87, 237401 (2001).

44. Kirilyuk, A., Kimel, A. V. & Rasing, T. Ultrafast optical manipulation of magnetic order. Rev. Mod. Phys. 82, 27312784 (2010).

45. Imada, M., Fujimori, A. & Tokura, Y. Metal-insulator transitions. Rev. Mod. Phys. 70, 10391263 (1998).

46. Dagotto, E. Complexity in strongly correlated electronic systems. Science 309, 257262 (2005).

47. Filion, L. et al. Efcient method for predicting crystal structures at nite temperature: variable box shape simulations. Phys. Rev. Lett. 103, 188302 (2009).

48. de Graaf, J., Filion, L., Marechal, M., van Roij, R. & Dijkstra, M. Crystal-structure prediction via the Floppy-Box Monte Carlo algorithm: method and application to hard (non) convex particles. J. Chem. Phys. 137, 214101 (2012).

Acknowledgements

We thank F. Romano for helpful discussions. This study was partly supported by Grants-in-Aid for Scientic Research (S) (Grant No. 21224011) and Specially Promoted Research (Grant No. 25000002) from the Japan Society for the Promotion of Science (JSPS).

Author contributions

J.R. and H.T. proposed and supervised the study, K.A. performed the numerical simulations. All authors discussed the results and contributed to the writing of the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Akahane, K. et al. A possible four-phase coexistence in a single-component system. Nat. Commun. 7:12599 doi: 10.1038/ncomms12599 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

8 NATURE COMMUNICATIONS | 7:12599 | DOI: 10.1038/ncomms12599 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications