
ABSTRACT

Title of Dissertation: PHASE TRANSITIONS AFFECTED BY
MOLECULAR INTERCONVERSION

Thomas Joseph Longo
Doctor of Philosophy, 2023

Dissertation Directed by: Professor Mikhail A. Anisimov
Institute for Physical Science and Technology &
Department of Chemical and Biomolecular Engineering

Typically, pure substances may be found with only one gaseous or liquid state, while their

solid state may exist in various polymorphic states. The existence of two distinct liquid forms in

a single component substance is more unusual since liquids lack the long-range order common

to crystals. Yet, the existence of multiple amorphous states in a single component substance, a

phenomenon known as “liquid polyamorphism,” has been observed or predicted in a wide variety

of substances. In contrast to standard phase transitions, it has been suggested that polyamor-

phic liquid-liquid transitions are caused by the interconversion of molecular or supramolecular

states. To investigate this phenomenon, a nonequilibrium thermodynamic model was developed

to quantitatively describe the interplay between the dynamics of molecular interconversion and

fluid-phase separation. The theory has been compared to a variety of interconverting systems,


and has demonstrated a quantitative agreement with the results of Monte Carlo and Molecular

Dynamics simulations.

In this thesis, it is shown that there are two major effects of molecular interconversion on

the thermodynamics and the kinetics of fluid-phase separation: if the system evolves to an equi-

librium state, then the growth of one of the alternative phases may result in the destruction of

phase coexistence - a phenomenon referred to as “phase amplification.” It is demonstrated that

depending on the experimental or simulation conditions, either phase separation or phase am-

plification would be observed. Previous studies of polyamorphic substances report conflicting

observations of phase formation, which may be explained by the possibility of phase amplifica-

tion occurring. Alternatively, if the system evolves to a nonequilibrium steady state, the phase

domain growth could be restricted at a mesoscopic length scale. This phenomenon (referred to as

“microphase separation”) is one of the simplest examples of steady-state dissipative structures,

and may be applicable to active matter systems, hydrodynamic instabilities, and bifurcations in

chemical reactions, in which the nonequilibrium conditions could be imposed by an external flux

of matter or energy.


PHASE TRANSITIONS AFFECTED BY
MOLECULAR INTERCONVERSION

by

Thomas J. Longo

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2023

Advisory Committee:
Professor Mikhail Anisimov, Chair
Professor Christopher Jarzynski
Professor Jeffery Klauda
Professor Pratyush Tiwary
Professor John Weeks


© Copyright by
Thomas J. Longo

2023


Preface

This thesis consists of results from my original research performed in collaboration with

my co-authors. The contents of this thesis have been adopted with permission from our publica-

tions, of which I am either the first or second author. The investigation of the fluid-fluid phase

transition in high-pressure hydrogen (Section 2.2), high-density sulfur (Section 2.3), and fluids

exhibiting water-like anomalies (Section 2.4) was conducted in Refs. [1–3], respectively. The

phenomenological theory that describes the dynamics of phase transitions affected by molec-

ular interconversion (Chapter 3) was developed in Refs. [4, 5]. This theory was verified via

Monte Carlo (MC) and molecular dynamics (MD) simulations of three physically distinct micro-

scopic models of mixtures with species interconversion (Chapter 4), which were investigated in

Refs. [6–8].

Throughout my Ph.D. project, I developed the theory, participated in the computer simula-

tions of the microscopic models, performed the data analysis, made the figures (except Figs. 4.10-

4.13, which were created by Betül Uralcan), and was a major contributor to the text of the pub-

lished manuscripts. I would like to ascribe appropriate credit to my research advisor, who con-

ceived the project and supervised the research, and to our collaborators, whose contribution was

crucial for the success of the project.

The two-state thermodynamic approach (Section 2.1) was originally developed by my advi-

sor, Mikhail Anisimov, and his co-workers [9]. The equation-of-state of high-pressure hydrogen

(Section 2.2) was developed in collaboration with Nathaniel Fried [1]. The maximum-valence

ii


model (Section 2.3), the HL model (Section 4.1), and the HCS model (Section 4.3) were de-

veloped and simulated by Sergey Buldyrev and Nikolay Shumovskyi [2, 5, 6, 8]. The blinking-

checkers lattice model (Section 2.4) was originally developed by Frédéric Caupin and Mikhail

Anisimov [10]; I described the interfacial properties of this model [3]. My major contribution to

the research project was the generalization of the Cahn-Hilliard theory of spinodal decomposition

and the development of the generalized Cook-Binder theory of the time-evolution of the structure

factor in fluids with interconversion of species (Chapter 3) [4,5,8]. Numerical calculations of the

time evolution of the order parameter and structure factor (Section 3.3) were performed in col-

laboration with Salim Asadov, Nikolay Shumovskyi, and Sergey Buldyrev [5]. The simulations

of the CM model (Section 4.2) were performed by Betül Uralcan, Frank H. Stillinger, and Pablo

Debenedetti [7, 8].

Thomas J. Longo

April 15, 2023

iii


Acknowledgments

I would like to express my sincere gratitude to my research advisor, Professor Mikhail

Anisimov. His constant tutelage, endless support (often working, literally, every day with me),

and insightful suggestions not only benefited the success of this project, but also enabled me to

grow as a researcher.

I am grateful for the collaborators that I had the privilege to work closely with, namely:

Salim Asadov (University of Maryland, College Park and Azerbaijan National Academy of Sci-

ences), Sergey Buldyrev (Boston University and Yeshiva University), Frédéric Caupin (Université

Claude Bernard Lyon 1 and Institut Universitaire de France), Pablo Debenedetti (Princeton Uni-

versity), Nathaniel Fried (University of Maryland, College Park), Nikolay Shumovskyi (Boston

University), and Betül Uralcan (Princeton University and Bogazici University). This endeavor

would not be possible without their contributions.

I am also thankful to the Chemical Physics Program at the Institute for Physical Science

and Technology as well as the Department of Chemical and Biomolecular Engineering for the

continuous encouragements and financial support during my time at the University of Maryland,

College Park.

This dissertation is a product of the research collaboration between the University of Mary-

land (the Lead institution), Princeton University, Boston University, and Arizona State University,

supported by the National Science Foundation (NSF). I acknowledge the financial support of the

iv


NSF (award No. CHE-1856479) as well as the partial support of the Chateaubriand Fellowship

of the Office for Science & Technology of the Embassy of France in the United States.

I would be remiss if I did not mention my family, in particular, my parents, my siblings, and

my fiancé. Thank you all for your endless support and encouragements throughout my studies.

v


Table of Contents

Preface ii

Acknowledgements iv

Table of Contents vi

List of Tables viii

List of Figures ix

List of Abbreviations xxiv

Chapter 1: Introduction 1

Chapter 2: Two-State Thermodynamics of Fluid Polyamorphism 8
2.1 Interconversion of Two States . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Fluid-Fluid Phase Transition in High-Pressure Hydrogen . . . . . . . . . . . . . 13

2.2.1 Phase Behavior of Hydrogen . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.2 Phase Separation Coupled with Dimerization . . . . . . . . . . . . . . . 15

2.3 Liquid-Liquid Phase Transition in High-Density Sulfur . . . . . . . . . . . . . . 23
2.3.1 Maximum-Valence Model . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.2 Phase Separation Coupled with Polymerization . . . . . . . . . . . . . . 28

2.4 Liquid-Liquid Phase Transition in Fluids Exhibiting Water-Like Anomalies . . . 32
2.4.1 Blinking-Checkers Lattice Model . . . . . . . . . . . . . . . . . . . . . 33
2.4.2 Virtual Critical Line in Interconverting Mixtures . . . . . . . . . . . . . . 37
2.4.3 Anomalies of Interfacial Properties . . . . . . . . . . . . . . . . . . . . . 42

2.5 Conclusion of Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

Chapter 3: Phase Formation Affected by Species Interconversion 63
3.1 Generalizing Cahn-Hilliard Theory . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.1.1 Conserved and Nonconserved Order-Parameter Dynamics . . . . . . . . 64
3.1.2 Phase Amplification vs. Microphase Separation . . . . . . . . . . . . . . 73

3.2 Temporal Evolution of the Structure Factor . . . . . . . . . . . . . . . . . . . . . 78
3.2.1 Generalized Cook-Binder Theory . . . . . . . . . . . . . . . . . . . . . 78
3.2.2 Characteristic Length Scales . . . . . . . . . . . . . . . . . . . . . . . . 81

3.3 Temporal Evolution of the Order Parameter . . . . . . . . . . . . . . . . . . . . 85
3.4 Effects of Heat and Volume Change of Interconversion . . . . . . . . . . . . . . 89

vi


3.5 Effects of Critical Order-Parameter Fluctuations . . . . . . . . . . . . . . . . . . 93
3.6 Conclusion of Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Chapter 4: Application of Theory to Microscopic Models 99
4.1 Hybrid Ising / Binary Lattice (HL Model) . . . . . . . . . . . . . . . . . . . . . 99

4.1.1 HL Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.1.2 Phase Amplification in the HL Model . . . . . . . . . . . . . . . . . . . 103
4.1.3 Nonequilibrium Spatially-Modulated Stripes . . . . . . . . . . . . . . . . 115

4.2 Chiral-Mixture of Interconverting Enantiomers (CM Model) . . . . . . . . . . . 122
4.2.1 CM Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.2.2 Phase Amplification in the CM Model . . . . . . . . . . . . . . . . . . . 127
4.2.3 Formation of Dissipative Structures . . . . . . . . . . . . . . . . . . . . 130
4.2.4 Characteristic Time Scales in the CM Model . . . . . . . . . . . . . . . . 137
4.2.5 Finite-Size Restrictions in the CM Model . . . . . . . . . . . . . . . . . 139

4.3 Coarse-Grained Hard-Core-Shoulder Model (HCS Model) . . . . . . . . . . . . 144
4.3.1 HCS Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . 144
4.3.2 Nonequilibrium Bicontinuous Microemulsions . . . . . . . . . . . . . . 147

4.4 Sources of Microphase Separation . . . . . . . . . . . . . . . . . . . . . . . . . 150
4.4.1 Comparison of the Microscopic Models . . . . . . . . . . . . . . . . . . 150
4.4.2 Onset and Termination of Microphase Separation . . . . . . . . . . . . . 152

4.5 Conclusion of Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

Chapter 5: Overall Conclusion 158

Appendix A: Comparison of Exact Solution with Phenomenological Ansatzes 160

Bibliography 165

vii


List of Tables

2.1 The suggested locations of the FFCP and the SFF-TP in hydrogen. . . . . . . . . 15
2.2 Liquid-vapor critical points of interconverting systems, referred to as “actual”

critical points, for the seven systems considered in this section (with εAA = 1.6,
εBB = 2.0, ê = 3, and ŝ = 4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.3 Liquid-liquid critical points (LLCP) for the three systems exhibiting liquid polyamor-
phism and a LLCP (with εAA = 1.6, εBB = 2.0, ê = 3, and ŝ = 4). . . . . . . . . 40

2.4 Asymptotic amplitudes of the liquid-liquid interfacial tension and liquid-liquid
correlation length of concentration fluctuations for the three systems exhibiting
liquid polyamorphism and a liquid-liquid critical point. The asymptotic mean-
field behavior is illustrated in Fig. 2.21. . . . . . . . . . . . . . . . . . . . . . . 58

2.5 The surface tension σ, normalized interfacial thickness ζ̂ = ζ/`, and normalized
shift δ̂ = δ/` for the three coexisting phases (liquid A, liquid B, and vapor) for
the system with εBA = 1.00 at the triple point temperature (T = 0.6843). . . . . . 58

3.1 Limiting cases of the interplay between diffusion, natural interconversion, and
forced interconversion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

3.2 Conditions for phase amplification and microphase separation as illustrated in
Fig. 3.1. The left column corresponds to the solid lines and the right column
corresponds to the dashed lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.1 Fit parameters of the theoretical prediction in Fig. 4.8a (top, pr = 1/32) and
Fig. 4.8b (bottom, pr = 1/256) of the inverse steady-state domain size, q−. . . . . 118

4.2 Parameters for the CM model in real and reduced units . . . . . . . . . . . . . . 124
4.3 Fit parameters used for the theoretical predictions in Fig. 4.14a of the inverse

steady-state domain size, q−. The mobility is given by the Stokes-Einstein rela-
tion with amplitude b0 and characteristic temperature T0 = 1.2. Two additional
parameters, the amplitude, A, and the characteristic size of the inhomogeneities
(q0 ∼ 1/R0) are included. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

4.4 Fit parameters used for the theoretical prediction illustrated in Fig. 4.25a for the
HL model, with pr = 1/128. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

A.1 The degree of asymmetry for each profile of the system with εBA = 1.04 at T =
0.6943, determined from the symmetric ansatz (sym), the Fisher-Wortis ansatz
(FW), or the exact solution as calculated from Eq. (A.3). . . . . . . . . . . . . . 163

viii


List of Figures

1.1 Hypothesized global pressure-temperature phase diagram for water. The black
curves indicate the solid-solid and solid-liquid-gas phase transitions experimen-
tally reported for water [11], while the blue curve indicates the hypothesized
liquid-liquid phase transition. The dotted black curve represents the line of spon-
taneous ice nucleation [12]. The red circle is the liquid-vapor critical point (LVCP)
and the red star is the hypothesized liquid-liquid critical point (LLCP) around
220 K and 60 MPa as suggested in Ref. [13]. . . . . . . . . . . . . . . . . . . . . 2

2.1 A symmetric binary mixture is quenched along the critical molecular fraction
(where the order parameter is ϕ = 0) from a point, T̂1 in the one phase region
to a point, T̂2 below the critical temperature of demixing in the unstable region.
The solid curve indicates the liquid-liquid coexistence (binodal), while the dashed
curve indicates the limit of stability (spinodal). . . . . . . . . . . . . . . . . . . . 9

2.2 The global pressure-temperature phase diagram for hydrogen. (a) The full range
from low to extreme pressures in logarithmic scale. The crosses indicate the
experimental data for the solid-liquid melting transition presented in Diatschenko
et al. [14] (blue), Datchi et al. [15] (cyan), Gregoryanz et al. [16] (pink), and Zha
et al. [17] (purple). The solid black curves at low pressure (P ≤ 0.1 GPa) are
the liquid-gas-solid phase transitions [18], while the solid black curve at high
pressure (P > 0.1 GPa) is the Kechin equation [19] as reported in ref. [17].
The black dashed curve is the predicted continuation of the melting line based
on experimental and computational evidence [17, 20, 21], while the dotted lines
represent the highly-debated prediction [22–28] of the domain of solid metallic
hydrogen [22,23,28–33]. The red line is the first-order fluid-fluid phase transition
adopted in this section. (b) The phase diagram of hydrogen at extreme conditions,
in the area of the box in (a). The open circles are experimental data presented in
Zaghoo et al. [34–36] (dark brown), McWilliams et al. [37] (light brown), and
Ohta et al. [38] (orange). Simulation results [21, 39–46] are spread within the
grey area and shown in detail in Fig. 2.3. The fluid-fluid phase transition (solid
red) and Widom line (dotted red) are represented by Eq. (2.18). The red star is
the location of the fluid-fluid critical point (FFCP) as adopted in this section. . . . 14

ix


2.3 Unifying the different simulation results with experimental data of hydrogen by
the generalized law of corresponding states. (a) Experimental and simulation
results for the fluid-fluid phase transition (FFPT). (b) Unified representation of
the FFPT by reducing pressure, P̂ = P/Pc, temperature, T̂ = T/Tc, and the
critical value of the entropy, Ŝ = S/Sc. In (a) and (b), the open circles are
the experimental data of Zaghoo et al. [34–36] (dark brown), McWilliams et al.
[37] (light brown), and Ohta et al. [38] (orange). The computational results are
indicated by the triangles: blue tints correspond to the Density Functional Theory
(DFT) simulations of Bonev et al. [47] (dark blue), Morales et al. [48] (blue),
Hinz et al. [44] (sky blue), and Karasiev et al. [45] (light blue). Meanwhile, green
tints correspond to the Quantum Monte Carlo (QMC) simulations of Morales et
al. [48] (dark sea green), Lorenzen et al. [39] (green), Perlioni et al. [41] (dark
green), Mozzola et al. [42] (lime green), and Tirelli et al. [46] (yellow green).
The colored stars correspond to the reported (or adopted in this section) critical
points for each data set. The solid black curve is the solid-fluid phase transition
line as discussed in Fig. 2.2, and the red solid line is the FFPT predicted in this
section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.4 The components of the Gibbs energy (per atom) for hydrogen in the vicinity of
the fluid-fluid phase transition. (a) The Gibbs energy of reaction, GBA, as given
by Eq. (2.9). The isotherms are T = 0.5Tc (orange), T = 0.75Tc (blue), T = Tc

(green), T = 1.25Tc (red), and T = 1.5Tc (purple). (b) The Gibbs energy of
mixing, Gmix, as given by Eq. (2.15). Gmix is shown as a function of the fraction
of hydrogen atoms in the diatomic state, x, for isotherms T = 0.5Tc (blue),
T = 0.75Tc (green), and T = Tc (red) at P = Pc. The solid curve corresponds to
the fluid-fluid coexistence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.5 Equilibrium fraction of hydrogen atoms in the diatomic state, xe. (a) Equilibrium
fraction-pressure diagram for T = 0.5Tc (orange), T = 0.75Tc (blue), T = Tc

(green), and T = 1.25Tc (red). (b) Equilibrium fraction-temperature diagram for
P = 0.75Pc (blue), P = Pc (green), P = 1.25Pc (red), P = 1.5Pc (purple). The
solid and dashed black curves are, respectively, the fluid-fluid coexistence and the
limit of thermodynamic stability (spinodal). . . . . . . . . . . . . . . . . . . . . 20

2.6 The pressure-density phase diagram of hydrogen based on the equation of state
developed in this section. The open circles correspond to the QMC simulations
of Tirelli et al. [46]. Isotherms are T̂ = 0.67 (orange), T̂ = 0.73 (red), T̂ = 0.8
(brown), T̂ = 0.87 (purple), T̂ = 0.93 (green), and T̂ = 1.0 (blue). The fluid-
fluid coexistence is shown by the solid black curve. The red star is the FFCP
adopted in this section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.7 A schematic representation of the phase diagram of sulfur. The solid black curves
represent the solid-liquid-gas phase transition in sulfur [49–51], while the blue
curve represents the experimentally verified liquid-liquid phase transition [51].
The dashed black line is the predicted continuation of the melting line based on
experimental evidence [49]. The black circles are the solid-liquid-liquid triple
point (TP-SLL) and the solid-liquid-gas triple point (TP-SLG), while the red cir-
cle is the liquid-gas critical point (LGCP) and the red star is the liquid-liquid
critical point (LLCP). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

x


2.8 Reactions and interactions in the maximum-valence model. (a) The three types
of covalent bond-forming reversible chemical reactions that may occur in the
system. If two atoms without bonds (S0) collide with each other, they may form
a bond and become S1 atoms. If a S0 and S1 atom collide, they may form a bond
and become S1 and S2 atoms, respectively. If two S1 atoms collide with each
other, they form an additional bond and become S2 atoms. (b-d) The three major
interactions between atoms, in which each atom is composed of a core and shell,
both with a radius σ and mass m. U(r) is the pair potential energy and r is the
distance from the center of an atom. (b) The cores of each atom interact with an
attractive square well of depth ε = 1 and width w = 0.4. (c) The shells may react
to form covalent bonds that consist of a narrow well with depth εb = 1 and width
wb = 0.02. (d) Phase segregation is coupled to polymerization via the additional
attractive interactions between atoms in state S2, described by a square well of
depth ε22 = 0.5 and width w22 = 0.3. . . . . . . . . . . . . . . . . . . . . . . . 27

2.9 Phase diagrams for the maximum-valence model (with ε22 = 0.5 and εb = 1.0)
obtained in an NVT ensemble after t = 106 time units. (a) The isotherms in
the P -ρ plane are T = 0.96 − 1.20 (red-purple) in steps ∆T = 0.02. (b) The
liquid-gas and liquid-liquid critical isochores in the P -T plane are ρLG

c = 0.35
and ρLL

c = 0.81 as indicated by the lower and upper dashed lines, respectively.
In both figures, the liquid-gas and liquid-liquid coexistence curves are calculated
via the Maxwell construction and indicated by the solid curves. The liquid-gas
(T LG

c = 1.023, P LG
c = 0.0922, ρLG

c = 0.35) and liquid-liquid (T LL
c = 1.187,

P LL
c = 2.28, ρLL

c = 0.81) critical points are indicated by the red open circles,
while the triple point (P TP = 0.0738, T TP = 0.995) is indicated by the black
open circles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.10 (a) T -ρ phase diagram for the maximum-valence model (with ε22 = 0.5 and
εb = 1.0) obtained in an NVT ensemble after t = 106 time units. The temperature
dependence of the fraction of atoms with two bonds, x2, (b) and the average
chain-length, 〈n〉, (c) in two coexisting liquid phases. The simulation data in (b)
is fit to a second order polynomial. . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.11 (a) The density correlation function g(r) and (b) the structure factor S(q) across
the liquid-liquid transition at T = 1.00 for densities of ρ = 0.65 (blue), ρ = 0.70
(orange), ρ = 0.75 (green), ρ = 0.80 (red), ρ = 0.85 (purple), and ρ = 0.90
(black). In (a), the sharp peak, around r = 1 (in units of σ), corresponds to
the length of the covalent bond, which increases upon increasing density. Simul-
taneously, in (b), the maximum of the structure factor (the first peak) shifts to
larger wavenumbers upon increasing density, while the second peak acquires a
characteristic bump, similar to what was recently observed in sulfur [51]. The di-
vergence of the structure factor at q = 0 indicates the divergence of the isothermal
compressibility in the vicinity of the LLCP. The insets (dashed boxes) highlight
the behavior of the maximum of the correlation function and second peak of the
structure factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

xi


2.12 Density-temperature and pressure-temperature phase diagrams for the blinking-
checkers model with εAA = 1.6, εBB = 2.0, εAB = 1.08, ê = 3, and ŝ = 4.
The liquid-vapor coexistence (blue curves) terminates at the liquid-vapor critical
point, LVCP. The liquid-liquid coexistence (red curves) terminates at the liquid-
liquid critical point, LLCP. The limits of thermodynamic stability (spinodals) are
given by the dashed curves. In (a,b), the dotted line corresponds to the condition,
x = 1/2, which qualitatively separates the regions enriched either by species A
(at low temperatures and low densities) or species B (at high temperatures and
high densities). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.13 Liquid-vapor (a,b) and liquid-liquid (c,d) coexistence curves for the seven sys-
tems with εAA = 1.6, εBB = 2.0, e = 3, s = 4, and with various values of εBA:
εBA = 1.00 (blue), εBA = 1.04 (orange), εBA = 1.08 (green), εBA = 1.12 (red),
εBA = 1.16 (purple), εBA = 1.20 (pink), and εBA = 1.24 (gray). The critical
points, indicated by the stars in (a,b) and open circles in (c,d), are the unique
liquid-vapor (LVCP) or liquid-liquid (LLCP) critical points in the interconvert-
ing system, referred to as “actual” critical points. As also indicated in Fig. 2.12,
species A is enriched in the low-density, low-temperature region, while species
B is enriched in the high-density, high-temperature region. For the system with
εBA = 1.00, the dashed blue curves in (c,d) indicate the liquid-vapor coexistence,
while in (a-d), the dotted blue line indicates the discontinuity at the triple point. . 38

2.14 Illustration of the thermodynamic path selected by the interconversion reaction
for ê = 3 and ŝ = 4 at constant volume, V , and at constant number of occupied
lattice sites,N1+N2, represented through (a) the liquid branch of the liquid-vapor
temperature-concentration coexistence (see Fig. 2.13b), and (b) the activity, a =

1/[1+e−GBA/T̂ ], for the systems with different interaction parameters: εBA = 1.00
(blue), εBA = 1.04 (orange), εBA = 1.08 (green), εBA = 1.12 (red), εBA = 1.16
(purple), εBA = 1.20 (pink), and εBA = 1.24 (gray). For each system, the liquid-
vapor critical line (LVCL) is shown by the dashed curves, while the collapsed
coexistence, in (a), is illustrated by the black curve. The insets show the LV
critical points for each scenario. Note that for the system with εBA = 1.00, the
thermodynamic path crosses through the triple point, indicated by the dotted line
in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.15 Liquid-vapor (a,b) and liquid-liquid (c,d) diameters of the density, Eq. (2.60), and
the concentration, Eq. (2.61), as a function of the distance to the virtual LVCL (or
LLCL) along the thermodynamic path selected by interconversion for systems
with εAA = 1.6, εBB = 2.0, ê = 3, ŝ = 4, and with various values of εBA:
εBA = 1.00 (blue), εBA = 1.04 (orange), εBA = 1.08 (green), εBA = 1.12 (red),
εBA = 1.16 (purple), εBA = 1.20 (pink), and εBA = 1.24 (gray). In (a-d), the
dotted blue lines indicate the discontinuity for the system with εBA = 1.00 at the
triple point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

xii


2.16 The liquid-vapor interfacial tension as a function of temperature (a), and also pre-
sented in reduced units (b), for the system with εBA = 1.00 (blue), εBA = 1.04
(orange), εBA = 1.08 (green), εBA = 1.12 (red), εBA = 1.16 (purple), εBA = 1.20
(pink), and εBA = 1.24 (gray). In (b), the critical temperature is given by the “vir-
tual” critical point (for the non-reactive binary mixture) for each concentration
along the thermodynamic path selected by the interconversion reaction. The blue
arrows indicate the direction of warming. (c) The reduced interfacial thickness,
ζ̂ = ζ/`, and (d) the reduced relative distance between the concentration and
density profiles, δ̂ = δ/`. In (a-d), the dotted lines indicate the discontinuity of
the interfacial properties for the system with εBA = 1.00 at the triple point, shown
by the vertical bars in (c,d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

2.17 (a) The behavior of the liquid-vapor interfacial tension follows the power law, σ =
σ0|∆T̂ |3/2, where the amplitude was found to be σ0 = 0.71, asymptotically close
to the actual liquid-vapor critical temperature (see Table 2.2). (b) The behavior
of the reduced liquid-vapor interfacial thickness, ζ̂ = ζ/`, follows the power law,
ζ̂ = ζ̂0|∆T̂ |−0.38, where the amplitude was found to be ζ̂0 = 1.50 asymptotically
close to the actual critical point. . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

2.18 Asymptotic behavior of the liquid vapor coexistence for systems with εAA = 1.6,
εBB = 2.0, ê = 3, ŝ = 4, and with various values of εBA: εBA = 1.00 (blue), εBA =
1.04 (orange), εBA = 1.08 (green), εBA = 1.12 (red), εBA = 1.16 (purple), εBA =
1.20 (pink), and εBA = 1.24 (gray). (a) The temperature-density LV coexistence
follows the meanfield power law, ∆T̂ ∼ |∆ρ̂|2, where ∆T̂act = 1 − T/T act

c

and T act
c is the actual critical temperature selected by the interconverting path.

Likewise, ∆ρ̂act = 1 − ρ/ρact
c , where ρact

c is the actual critical density. (b) The
temperature-average concentration LV coexistence follows the meanfield power
law, ∆T̂ ∼ |∆x̄act|2, where ∆x̄act = 1 − x̄/xact

c , in which x̄ = (xL + xV)/2
and xact

c is the actual critical concentration. (c,d) Illustrate, as an example, the
asymptotic behavior of the system with εBA = 1.08, in which (d) shows the
asymptotic behavior of each branch of the concentration coexistence. . . . . . . . 53

2.19 Comparison between the liquid-liquid (dashed curves) and liquid-vapor (solid
curves) interfacial tensions as a function of temperature for the system with εBA =
1.00 (blue), εBA = 1.04 (orange), εBA = 1.08 (green), εBA = 1.12 (red). The
dotted blue line indicates the discontinuity in the liquid-vapor interfacial tension. 56

2.20 Liquid-liquid interfacial properties of the systems exhibiting liquid polyamor-
phism with εAA = 1.6, εBB = 2.0, ê = 3, ŝ = 4, and with various values of
εBA: εBA = 1.00 (blue), εBA = 1.04 (orange), εBA = 1.08 (green), εBA = 1.12
(red). (a) the reduced thickness, ζ̂ = ζ/` of the liquid-liquid interface, and (b)
the reduced shift, δ̂ = δ/`, between the density and concentration liquid-liquid
profiles. In (a,b) the thickness and shift reach a finite value (marked with a blue
circle) at the triple point temperature. . . . . . . . . . . . . . . . . . . . . . . . . 56

xiii


2.21 (a) The behavior of the liquid-liquid interfacial tension follows the meanfield
power law, σ = σ0|∆T̂ |3/2, (dashed lines) asymptotically close to the actual
liquid-liquid critical temperature (see Table 2.3). (b) The behavior of the reduced
liquid-liquid interfacial thickness, ζ̂ = ζ/`, follows the meanfield power law, ζ̂ =
ζ̂0|∆T̂ |−1/2, (dashed lines). In (a,b) the systems exhibiting liquid polyamorphism
and a liquid-liquid critical point are shown: εBA = 1.04 (orange), εBA = 1.08
(green), εBA = 1.12 (red), and the amplitudes of the asymptotic meanfield power
laws are provided in Table 2.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

2.22 Normalized density and concentration liquid-vapor profiles as a function of the
coordinate perpendicular to the planar interface, ẑ = z/`, given by Eqs. (2.54)
and (2.55) for the system with εAA = 1.6, εBB = 2.0, εBA = 1.08, ê = 3, and
ŝ = 4 at the two temperatures (a,b) that correspond to the two extrema of the
liquid-vapor interfacial tension (shown in Fig. 2.16). Normalized (c) density and
(d) concentration profiles for three-phase coexistence at the triple point, TTP =
0.6843, for the system with εAA = 1.6, εBB = 2.0, and εBA = 1.00. . . . . . . . . 59

2.23 Interfacial profiles of species B, ρB = ρx, in the blinking-checkers model demon-
strate surface enrichment near the TP temperature, TTP = 0.68429. (a) Surface
enrichment of species B for the system with εBA = 1.00. The colored curves indi-
cate temperatures from T = 0.68989 to T = TTP in steps of ∆T̂ = −0.0008 in or-
der of purple to red. The black curves are T = 0.68389 (dashed) and T = 0.68309
(solid). (b) Surface enrichment of species B for the system with εBA = 1.04. The
curves are T = 0.6882 to T = 0.6826 in steps of ∆T = −0.0008 (blue to pink).
In (a,b), the black arrows indicate the direction of decreasing temperature. Note
that while the transition of a surface enriched profile (T > TTP) to a smooth pro-
file (T < TTP) is discontinuous in the system with εBA = 1.00, it is continuous in
the system with εBA = 1.04. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

3.1 The characteristic growth rate, Eq. (3.14), affected by the competition between
diffusion, natural interconversion, and forced interconversion at ∆T̂ = −0.5.
Complete phase separation (as predicted by Cahn-Hilliard’s theory for L = 0
and K = 0) is illustrated by the red curve for M = 100. Phase amplification is
illustrated by the purple curves for restricted (M = 100 - solid) and unrestricted
(M = 10 - dashed) cases, in which L = 10 and K = 2. Microphase separation
for M = 100, L = 1, and K = 2 is illustrated by the solid green curve. When
the growth rate is always negative, as illustrated by the green dashed curve (for
M = 1, L = 1, and K = 2), there is no phase domain growth corresponding
to a homogeneous steady state. The green circles indicate the three characteristic
wavenumbers of the amplification factor: the maximum, qm, the lower cut-off,
q−, and the upper cut-off, q+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

xiv


3.2 Steady-state phase domain morphology for different magnitudes of forced inter-
conversion (after ∼ 105 time steps) numerically computed from the time evolu-
tion of the order parameter, Eq. (3.7), with M = 1, L = 1/127, ∆T̂ = −0.1,
` = 64, σi = 0.1, and η = 10−5, as discussed in Section 3.3. Morphologies
are shown for the middle slice of the three-dimensional system at (a) K = 0, (b)
K = 5× 10−4, (c) K = 15× 10−4, and (d) K = 25× 10−4. The red regions
correspond to where the value of the normalized order parameter is ϕ/ϕmax = 1,
the purple regions correspond to where the value of the normalized order param-
eter is ϕ/ϕmax = −1, and the other colors depict the interface between these two
regions. The image in (a) depicts a metastable structure toward phase amplifica-
tion [6], while the images in (b-d) are modulated steady-state structures with a
characterize size, 1/q−. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

3.3 a) The amplification factor, ω(q), given by Eq. (3.15) with κ = 1, ∆T̂ = −0.1,
M = 1, L = 1/127, and K = 1.3× 10−3. b) The time evolution of the structure
factor, given by Eq. (3.23) in the presence of natural and forced interconversion.
The black dotted line depicts the evolution of the maximum of the structure factor.
Due to the external source of forced interconversion, the maximum of the struc-
ture factor is interrupted at the wavenumber q−, while for complete phase separa-
tion and phase amplification, the maximum of the structure factor will evolve to
q = 0 for an infinite-sized system. . . . . . . . . . . . . . . . . . . . . . . . . . 80

3.4 Temporal evolution of the structure factor: a) for a system undergoing diffusion
dynamics (M = 1) toward an equilibrium state in the absence of natural intercon-
version (L = 0) and forced interconversion (K = 0); b) for a system undergoing
a hybrid of diffusion (M = 1) and natural interconversion (L = 0.01) dynam-
ics in the presence of forced interconversion (K = 1.5× 10−3) toward a steady
state. The structure factor, given by Eq. (3.23), exhibits a crossover from spinodal
decomposition to the nucleation regime. The dashed-black curves indicate the de-
velopment of the maximum of the structure factor. The characteristic crossover
time is defined in Eq. (3.24) and adopted as τ = 100. In (a) the evolution of the
maximum of the structure factor moves to q = 0 for infinite-size system and satu-
rates at Sm(q = 0, t→∞) = 1/(2ξ2) = 5 for ∆T̂ = −0.1. In contrast, in (b) the
evolution of the maximum is interrupted at a characteristic cut-off wavenumber
predicted by the characteristic phase domain growth rate, qs

m(t→∞) ∝ q−, and
it saturates at Sm(q−, t→∞) = 417. . . . . . . . . . . . . . . . . . . . . . . . . 82

3.5 The time evolution of the wavenumber corresponding to the maximum of the
structure factor, given by Eq. (3.23), during the crossover from the early stage
of spinodal decomposition, q ∝ t1/4 (green - dashed), to the nucleation regime,
q ∝ t1/3 (orange - dashed) for a system undergoing diffusion dynamics in the
absence of natural interconversion (L = 0) and forced interconversion (K = 0)
under conditions: M = 1, ∆T̂ = −0.1, τ = 100. . . . . . . . . . . . . . . . . . 83

xv


3.6 The effect of increasing interconversion force on the phase domain growth rate
for M = 100, L = 1, and ∆T̂ = −0.5. The red dashed line corresponds to the
inverse maximum size of the phase domain on the length scale of the simulation
box, q∗. When q− > q∗ microphase domains will form. Alternatively, when
q− < q∗, the size of the simulation box will cut-off the growing phase domains.
The conditions where ω(q) < 0 (dashed-dot portions of the curves) corresponds
to non-growing wavenumbers. As the rate of forced interconversion increases,
the growth rate is shifted down from the onset of phase separation where q− = q∗

(red, K = 1), to the microphase region (green, K = 3.75), to the termination
point of domain growth (blue, K = 6.5) where q− = qm = q+ = q∗∗, and to the
no growth regime for any wavenumber (orange, K = 9.25). . . . . . . . . . . . 84

3.7 Time evolution of the structure factor computed from the Fast Fourier trans-
form (FFT) of Eq. (3.7) for M = 1, L = 1/127, ∆T̂ = −0.1, ` = 64,
σi = 0.1, η = 10−5 depicted at times t = 6× 103 (green), t = 1.2× 104 (blue),
t = 2.4× 104 (orange), t = 5× 104 (red), t = 1× 105 (pink), and t = 2× 105

(black). The open circles in (a-d) depict the computed structure factors for the
four selected magnitudes of forced interconversion averaged over N = 100 real-
izations with 95% confidence interval error bars, while the solid lines illustrate the
behavior of the structure factors assuming a Gaussian distribution. The wavenum-
ber is normalized by the size of the system, such that q = 1 corresponds to phase
domains with a characteristic size of half the simulation box, `/2. . . . . . . . . . 86

3.8 The temporal evolution of the symmetry of phase separation. a) The time evolu-
tion of the average order parameter, calculated by first averaging over all space
and second by averaging the absolute value over N = 100 realizations, for
M = 1, L = 1/127, ∆T = −0.1, σi = 0.1, η = 1.0× 10−5, and various mag-
nitudes of forced interconversion, K. b) The time evolution of the N -averaged
standard deviation of the averaged order parameter, calculated by first determin-
ing the standard deviation of the spatially averaged order parameter and second
by averaging over N = 100 realizations. This method of averaging highlights the
behavioral deviation from an equal concentration of species A and B, ϕ = 0. . . . 88

3.9 Three hypothesized binary mixture systems exhibiting interconversion of species
and liquid-liquid phase separation quenched from high temperature to low tem-
perature (without volume change). The black dashed curve corresponds to the
liquid-liquid phase coexistence in this system without interconversion and with
interaction energy, ε = 2. The open circle indicates the liquid-liquid critical point
(LLCP), while the crosses show the locations of T = TBA, the points correspond-
ing to 50:50 interconversion for different energy change of reaction. For a system
with TBA = 1.05Tc (green), no liquid-liquid phase transition will be observed
upon quenching. For a system with TBA = Tc (red), the quenching process passes
through the critical point. For a system with TBA = 0.95Tc (purple), there are two
equilibrium solutions for the fraction of interconversion, such that upon quench-
ing to the cross, phase amplification occurs with equal probability of forming an
A-rich or B-rich phase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

xvi


3.10 The phase diagram (suggested in ref. [13]) for supercooled water that exhibits
a liquid-liquid phase transition (global phase diagram illustrated in Fig. 1.1). A
hypothesized quenching process by compression for supercooled water is shown
from the one phase region at P1 = 20 MPa (orange) to the two phase region at
P2 = 120 MPa (green) along the Widom line (dashed black) which corresponds
to a line of constant fraction of interconversion, lnK = 0. Two additional isobars
are shown for reference at P = 40 MPa (blue) and P = 80 MPa (purple) along
with the liquid-liquid coexistence (black). Phase amplification would only be
possible in a system where the number of molecules changes to compensate the
volume change of the interconversion reaction. . . . . . . . . . . . . . . . . . . . 92

3.11 a) Characteristic phase domain growth rate in the vicinity of the critical point
(∆T̂ = −0.001) forM0 = 1, L = 0.002, andK = 2.25×10−5 calculated through
Eq. (3.13), red curve, with use of the diverging molecular mobility, Eq. (3.27),
and scaling inverse susceptibility in the first order epsilon expansion, χ̂−1

q=0 ∼
|∆T̂ |−γ with γ = 1 + ε/6 (ε = 4 − d). The meanfield approximation is shown
by the green curve, Eq. (3.13). b) The onset, red solid curves, Eq. (3.28), and
termination, red dashed curves, Eq. (3.29), of microphase separation affected by
critical fluctuations for M = 1, L = 0.01, ` = 100, ν = 1/2 + ε/12. The
meanfield approximation is shown by the green curves. . . . . . . . . . . . . . . 96

4.1 The spontaneous equilibrium order parameter (ϕ = ϕ0) in the lattice gas / lattice
binary mixture along the liquid-vapor phase coexistence (red domain). One of
the two alternative magnetizations (ϕ0 > 0 and ϕ0 < 0) in the Ising ferromagnet
in zero field are shown in the red domain with blue arrows. The solid curve is the
crossover from meanfield behavior (dashed) to the asymptotic scaling power law
ϕ ∝ ∆T β with β = 0.326 [52, 53], while the crosses correspond to simulation
results of the HL model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

4.2 Phase amplification in the HL model for systems with T = 4.0, ` = 100, and
interconversion probabilities pr: (a) 100%, (b) 10%, (c) 10−4%, and (d) 10−7%.
The inset in (d) shows 10−7% interconversion probability at a shorter time scale.
Each curve represents one of the 100 different realizations of the MC simulations,
while the solid zero line depicts the behavior of a system with 0% interconversion.
The simulation snapshots in (a,d) correspond to the end of the simulation time. . . 104

4.3 The evolution of the order parameter during phase amplification. (a) The RMS
of the distribution of the growth rates for different probabilities captured at the
same time, t = 300. The solid curve is the crossover between σ ∝ √pr and
σ ∝ pr, approximated as σ = a

√
pr(1 + bpr)/(1 +

√
pr). (b-d) The growth of the

order parameter for different (b) probabilities at ∆T̂ = −0.11 and ` = 100, (c)
system sizes at pr = 1.0 and ∆T̂ = −0.11, and (d) distances to the critical point
at pr = 1.0 and ` = 100; The inset shows the power law for the initial growth of
the reduced order parameter, ϕ/ϕ0 ∝ t3/4. . . . . . . . . . . . . . . . . . . . . . 106

xvii


4.4 Scaling properties of the growth of the reduced order parameter, ϕ̂ = 〈ϕ〉/ϕ0. (a)
The order parameter growth with time rescaled by system size. The size depen-
dence of the rescaling parameter, τ(`), is shown in the inset; in the log-log scale
with a slope of 1. The colors are the same as in Fig. 4.3c. (b) The order pa-
rameter growth with time rescaled by probability; the rescaling parameter τ0(pr),
inversely proportional to the probability, is shown in the inset. The colors are the
same as in Fig. 4.3b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

4.5 Topological characteristics of the time dependence of phase amplification. (a) For
spherical domains, the reduced deviation from the equilibrium order parameter,
∆ϕ̂ = 1 − ϕ/ϕ0 scales as t3/2. This is shown for temperatures, ∆T̂ , as: −0.11
(cyan), −0.025 (purple), −0.014 (orange), and −0.009 (brown). The inset shows
the effect for a cylindrical domain at ∆T̂ = −0.11. (b) A zero curvature Schwarz-
P interface is initially formed by simulating a system with Kawasaki dynamics
(pr = 0) for a long time. At t − t0 = 0, the system obtains Glauber dynamics
(pr = 1) and the collapse of one of the phases is shown. Amplification transitions
from random-walk behavior,

√
t, at short times (see inset) to exponential behavior

before saturation shown by the straight line. . . . . . . . . . . . . . . . . . . . . 110
4.6 The time evolution of the concentration for the HL model with pure Ising dy-

namics (pr = 1), given by Eq. (4.5) with the time-dependent susceptibility in the
form of Eq. (4.4). For ∆T̂ = −0.32 (green), τ = 0.2, ϕ∞ = 0.91, a = 1.9 and
b = 0.3, and for ∆T̂ = −0.11 (blue) τ = 0.2, ϕ∞ = 0.73, a = 4.6, and b = 0.4.
The open circles are the computational data presented in Fig. 4.3d. . . . . . . . . 113

4.7 Effect of forced interconversion on domain size, R, normalized by the system
size, `, in the HL model. (a) The time dependence of the domain growth for en-
ergy source E = 5 and interconversion probability pr = 1/128 at T/Tc = 0.24
(green), T/Tc = 0.27 (blue), and T/Tc = 0.40 (red), where Tc = 4.511 [54]. The
horizontal dashed line indicates the size of the simulation box, R = `. (b) Tem-
perature dependence of the steady-state domain size for E = 5 and pr = 1/128.
The vertical dashed line indicates the onset temperature, T ∗/Tc. (c) Depen-
dence of the steady-state domain size on the energy of forced interconversion
for pr = 1/256 and T/Tc = 0.31. The vertical dashed line denotes the onset
source energy, E∗. In (a-c), the system is simulated on a 3-dimensional lattice of
size ` = 100. The open circles are the results of MC simulations, the images are
snapshots of the system for selected conditions, and the curves are the theoretical
predictions. In (a-c), black denotes up-spins and white denotes down-spins. . . . 116

4.8 Domain size as a function of temperature in the HL model. (a) pr = 1/32 and (b)
pr = 1/256, for a system of size ` = 100 and energies: E = 1 (orange), E = 2
(green), E = 3 (red), E = 4 (purple), E = 5 (brown), E = 6 (pink), E = 7
(gray), E = 8 (yellow), E = 9 (cyan), and E = 10 (dark blue). The solid curves
are the theoretical predictions of q−, the inverse steady-state domain size, where
the fit parameters are provided in Table 4.1. . . . . . . . . . . . . . . . . . . . . 118

xviii


4.9 a) Steady-state structure factors computed for the HL (open circles) and the pre-
diction given by Eq. (4.7) (solid lines) for selected external energy sources (E)
at ∆T̂ = −0.4, M = 1, L = 1/127, ` = 100, and averaged over N = 60 re-
alizations with 95% confidence interval error bars. The insets show steady-state
(t ∼ 3× 105) domain morphologies observed in the HL model for the selected
energies. b) The dependence of forced interconversion on the wavenumber cor-
responding to the maximum of the structure factor, qs

m, in the steady-state limit.
The open circles are numerical computations of structure factors determined from
FFTs of the time evolution of the order parameter, given by Eq. (3.7), in the
steady-state limit (t ∼ 105) for M = 1, L = 1/127, σi = 0.1, and η = 10−5,
averaged over N = 100 realizations. The triangles correspond to the predic-
tions of K determined from fits of Eq. (4.7) to the structure factor for the HL
model, like those presented in (a). The curves illustrate the theoretical predic-
tion qm(t→∞) ∝ q−, given by the full expression for q−, found from evaluating
ω(q, 0) = 0 using Eq. (3.15). The colors correspond to temperatures: ∆T̂ = −0.1
(blue), ∆T̂ = −0.2 (green), ∆T̂ = −0.3 (red), and ∆T̂ = −0.4 (purple). The
inset shows the relationship between K and E. . . . . . . . . . . . . . . . . . . . 121

4.10 Molecular representation and geometrical features of tetramer molecules. L-
enantiomers (green), D-enantiomers (blue) and achiral transition states (cis- or
trans-, red) are shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

4.11 Time dependence of the instantaneous dihedral angle of a typical tetramer in a
racemic mixture at P = 0.1 for the conservative formulation of the CM model, at
several values of T and kd. a) T = 0.6, b) T = 1.7, and c) T = 2.3, with kd = 5
(green), kd = 9.86 (orange) and kd = 19.86 (purple). d) Behavior of the dihedral
angle at a very low value of the dihedral constant, kd = 0.001. . . . . . . . . . . 128

4.12 The phase diagram showing phase amplification in the CM model with conser-
vative intermolecular forces, heterochiral bias parameter λ = 0.5, and rigidity
spring constant kd = 0.001. The circles on the solid curve are the computational
data for the critical temperature of equilibrium phase separation and the curve
is the fit of Eq. (4.12). The snapshots depict the equilibrium states for the pres-
sures P = 0.1, P = 1, P = 5, and P = 10 below the critical temperature and
at P = 1 above the critical temperature. The triangles show the prediction of
the critical temperature from the extrapolation of the CM model with dissipative
intermolecular forces to the limit kd →∞. . . . . . . . . . . . . . . . . . . . . . 129

4.13 The change of compositional heterogeneity with chiral interconversion kinetics
at T = 1.7 and P = 0.1 in the dissipative-force formulation. a) Steady-state
snapshots of chiral liquid systems at various dihedral force constants (kd). b) The
steady-state domain size as a function of interconversion rate, 1/τINC. The solid
line is the approximation given by 1/τINC = a1/R

2
∞ + a2/R

4
∞, which follows

from Eq. (4.15), where a1 = 4.6 × 10−3 and a2 = 3.8 × 10−4. The inset shows
the linear correlation betweenR∞ and kd. The colored points highlight the results
corresponding to the three dihedral force constants for which the domain growth
is shown in Fig. 4.15. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

xix


4.14 Steady-state domain size, R, normalized by the size of the system, `, in the CM
model. (a) The time evolution of the domain size for different interconversion
rates, tuned by the rigidity parameter, kd, as kd = 3 (purple), kd = 5 (green),
and kd = 9.86 (red) at the reduced pressure P = 0.1 and T/Tc = 0.35, where
Tc(P = 0.1) = 2.32, as indicated on Fig. 4.12. (b) The normalized steady-state
domain size as a function of temperature at P = 0.1 and kd = 5. The verti-
cal dashed line indicates the onset temperature, T ∗/Tc. In (a) and (b), the open
circles correspond to simulation results, the curves correspond to the theoretical
predictions (see Table 4.3), and the images show snapshots of the system simu-
lated at the indicated conditions. In (a-b), dark/clear spheres correspond to the
L-/D-configuration of a chiral tetramer (spheres are located a tetramer’s center of
mass). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

4.15 Phase domain growth for the dissipative-force formulation of the CM model. The
open circles represent computational data [7], and the curves illustrates predic-
tions of the time evolution of the domain size, R(t) = 1/qs

m, from Eq. (3.23) for
dihedral force constants: kd = 5 (green), kd = 9.86 (red), kd = 19.86 (blue), and
kd → ∞ (black) for T = 1.8, Tc = 2.3, τ = 2, and M = 0.0022. The domain
size is normalized by the size of the computational box, `. The dashed curves
represent the predictions of the domain growth if it is not restricted by the finite
size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

4.16 The mutual diffusion coefficient affected by interconversion, given by Eq. (4.17),
as a function of the normalized steady-state domain size R∞, Eq. (4.18) in the
CM model. The open circles are computational data for the dissipative-force
formulation of the CM model at three different dihedral force constants: kd = 5
(green), kd = 9.86 (red), and kd = 19.86 (blue) [7]. The molecular mobility in
the limit kd → ∞ is approximated as M = b0T/η, where b0 = 0.94 and η is the
viscosity approximated by the Arrhenius equation η ∼ eT0/T . The characteristic
temperature T0 is 3.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136

4.17 Temperature dependence of the characteristic time scales in the dissipative-force
formulation of the CM model for dihedral force constants: kd = 5 (green), kd =
9.86 (orange), and kd = 19.86 (purple) at P = 0.1. a) Characteristic time for
liquid-liquid phase separation at the length scale of the simulation box, q = q∗.
The curves are τLLPS ∝ 1/ω̃(q∗) where ω̃(q) is given by Eq. (4.16) where it was
found that q∗ = 0.15 and Tc = 2.35. The condition τLLPS → ∞ corresponds to
T = T ∗. b) Characteristic times for chiral interconversion. τINC ∝ 1/L, and the
curves are given by Eq. (4.13). c) Characteristic self-diffusion times, where the
curves are given by Eq. (4.19). . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

xx


4.18 Dihedral force constant dependence of the onset temperature of phase separation
on the length scale of the simulation box in the dissipative-force formulation CM
model for P = 1.0 (red circles), P = 0.5 (blue circles), and P = 0.1 (black
circles). The curves are numerically calculated from the first solution of ω̃ = 0,
given by Eq. (4.16), when q− = q∗ = 0.11 ≈ 1/Rmax ∼ 1/`, T = T ∗, and
Tc(P = 1.0) = 3.45, Tc(P = 0.5) = 2.91, and Tc(P = 0.1) = 2.3 for different
pressures (a) and in rescaled coordinates (b). The triangles, shown in (b), are
obtained from the asymptotic limits of the time of liquid-liquid phase separation,
τLLPS →∞, as shown in Fig. 4.17a for q∗ = 0.15 and Tc = 2.35. . . . . . . . . . 141

4.19 The conditions for microphase separation at a length scale, R∞ = 1/q−. The
solid curve depicts the onset temperature of microphase separation, T = T ∗,
at the scale of the computational box, Eq. (4.20), for R∞ = 1/q∗ = 7.1 and
Tc = 2.4. The two dashed curves depict the lines of steady-state domain size
that are smaller than the size of the computational box, R∞ = 3.6 (lower) and
R∞ = 2.0 (upper), while the dotted curve depicts the growth-termination tem-
perature, T = T ∗∗, Eq. (4.21). The symbols are computational data (from the
observed size of the phase domain) in the dissipative-force formulation of the
CM model for pressure P = 0.1. The symbols indicate the onset of phase sep-
aration (open circles), the microphase region (R∞ = 1/q∗, closed circles), and
the two phase region (R∞ > 1/q∗, triangle). The two simulation snapshots illus-
trate the behavior of the system composed of A-rich (green), B-rich (blue), and
intermediate (red) states of the CM model. [7] . . . . . . . . . . . . . . . . . . . 143

4.20 Self-diffusion (molecular mobility) of the HCS model. Simulations were con-
duced for a system with b = 0 (no interconversion, either natural or forced).
The open circles indicate simulation results, while the curve is the fit of M =
0.056

√
T + 0.015T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

4.21 Critical temperature of the HCS model. Internal energy, ∆U , of the HCS model
as a function of temperature (blue) for the case b = 0 (no interconversion, either
natural or forced). The blue circles indicate the results of simulations, while the
red circles correspond to the calculated isochoric heat capacity, CV = ∂U/∂T |V .
The solid red curve is an empirical fit ofCV , where the temperature corresponding
to the maximum of CV is the critical temperature for the HCS model, Tc = 3.6±
0.05, in units of ε0/kB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

4.22 Steady-state domain size, R, normalized by the size of the system, `, in the HCS
model. (a) The temperature-dependence of the normalized steady-state domain
size for b = 0.005 (blue), b = 0.050 (black), and b = 0.075 (green) at ε =
10. (b) The normalized steady-state domain size as a function of interconversion
probability, b, for the energy source ε = 12 and T/Tc = 0.22. The vertical dashed
line indicates the onset interconversion probability, b∗. In (a) and (b), the open
circles correspond to simulation of 64, 000 particles, the curves correspond to the
theoretical predictions, and the images show snapshots of the system simulated
at the indicated conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

xxi


4.23 Structure factor, S(q, t), in the HCS model. The formation of a steady-state bi-
continuous microemulsions observed in MD simulations of the HCS model for
64,000 particles at interconversion energy source ε = 12, T/Tc = 0.1, and inter-
conversion probability b = 0.05 (a) and b = 0.15 (b). The images show snapshots
of the steady-state system simulated at the specified conditions. . . . . . . . . . . 149

4.24 Domain size as a function of temperature for the HCS model. (a) ε = 6 and
(b) ε = 8, for system size ` = 40σ, averaged over N = 10 realizations, for
interconversion probabilities, b = 0.005 (blue), b = 0.05 (black), and b = 0.075
(green). It should be noted that for low temperatures, T/Tc ≤ 0.8, the simulated
systems freeze, such that the characteristic domain size is always less than the
predicted size (as indicated by the triangles). . . . . . . . . . . . . . . . . . . . . 150

4.25 The onset and termination of microphase separation. (a) The steady-state domain
size in the HL model for pr = 1/128 and different forced interconversion ener-
gies from E = 1 (orange) to E = 10 (dark blue) in steps of ∆E = 1, just as
in Fig. 4.8 (Fit parameters provided in Table 4.4). For E > E∗∗ = 7 (the ter-
mination energy), the data collapse into a single line (black), indicating that the
characteristic steady-state domain size remains on the order of the microscopic
length scale R0(T ), which corresponds to homogeneous steady-state systems for
all temperatures. For E ≤ 7, the onset of microphase separation is observed at
T = T ∗(E∗), where E∗ is the onset energy. For T < T ∗, the steady-state domain
size is equal to the size of the system, R = `. (Caption continues on the next page.)154

4.25 (b) The onset energy E∗ (black circles and curve) for the HL model as a function
of temperature for pr = 1/128. Colored open circles and dashed curves corre-
spond to steady-state domain sizes: R = 0.143 (blue), R = 0.095 (green), and
R = 0.074 (red). (c) The inverse onset rigidity parameter 1/k∗d (black circles
and curve) for the CM model at P = 0.1. Colored circles and dashed curves
correspond to steady-state domain sizes: R = 0.32 (blue), R = 0.22 (green), and
R = 0.18 (red). In (b) and (c), the blue area corresponds to phase separation on
the scale of the simulation box, the white area corresponds to microphase separa-
tion, and the yellow area corresponds to homogeneous steady states. The images
in (b) and (c) correspond to the different final states of the systems below E∗(T ∗)
and 1/k∗d(T

∗) where the size of the phase domain is on the scale of the simulation
box (q∗ ∼ 1/`). In (b), phase amplification is observed since for pr = 1/128 nat-
ural interconversion is relatively fast, while in (c), where natural interconversion
is relatively slow for the simulated range of kd, complete two-phase separation,
on the scale of the simulation box, is found. . . . . . . . . . . . . . . . . . . . . 155

A.1 Comparison between the Fisher-Wortis (FW) and an alternative symmetric ansatz
(sym). obtained for the liquid-vapor coexistence, for two systems with εBA =
1.04 (a) and εBA = 1.12 (b) with εAA = 1.6, εBB = 2.0, ê = 3, and ŝ = 4. (c) The
relative deviation between the symmetric and FW ansatzes. . . . . . . . . . . . . 161

A.2 Numerical calculations of the liquid-vapor surface tension for three temperatures
in the system with εBA = 1.04, εAA = 1.60, εBB = 2.00, ê = 3, and ŝ = 4. The
FW ansatz is given at iteration 0. The temperatures were chosen based on surface
tensions with similar values, as predicted by the FW ansatz. . . . . . . . . . . . . 162

xxii


A.3 Comparison of the liquid-vapor density (a) and concentration (b) interfacial pro-
files, as obtained from numerical calculations after eight iterations, ρV4 and xV4,
and for the FW ansatz, ρV0 and xV0, for the system with εBA = 1.04, εAA = 1.60,
εBB = 2.00, ê = 3, and ŝ = 4 at temperature, T = 0.6943 (red curve in Fig. A.2).
(c,d) The difference between the exact solution for the density and concentration
profiles and the FW and symmetric ansatzes. In (a-d), the different profiles were
aligned along their Gibbs dividing surface, such that the excess density is zero for
all profiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

xxiii


List of Abbreviations and Key Notations

Acronyms
FFCP Fluid-fluid critical point
FFPT Fluid-fluid phase transition
FFT Fast Fourier transform
FW Fisher-Wortis
HDL High-density liquid
INC Characteristic interconversion
LDL Low-density liquid
LLCP Liquid-liquid critical point
LLPS Liquid-liquid phase separation
LLPT Liquid-liquid phase transition
LVCL Liquid-vapor critical line
LVCP Liquid-vapor critical point
MC Monte Carlo
MD Molecular dynamics
DMD Discrete molecular dynamics
QMC Quantum Monte Carlo
RG Renormalization Group
RMS Root mean square
SFF-TP Solid-fluid-fluid triple point
SLL-TP Solid-liquid-liquid triple point
SLG-TP Solid-liquid-gas triple point
TP Triple point

Simulated Models
HL Hybrid Ising / binary lattice
CM Chiral mixture of interconverting enantiomers
HCS Coarse-grained hard-core-shoulder

Key Notation
b Interconversion probability in the HCS model
ε Non-ideality interaction parameter
` Size of a system
Ft Total free energy (Helmholtz energy)
f = Ft/N Free energy per molecule (or per lattice site)
Gt Total Gibbs energy

xxiv


G = Gt/N Gibbs energy per molecule (or per lattice site)
Gmix Gibbs energy of mixing
GBA Gibbs energy of reaction
κ The square of the range of molecular interactions
kB Boltzmann’s constant
kd Force constant of dihedral motion in the CM model
K Kinetic coefficient of forced interconversion
L Kinetic coefficient of natural interconversion
M Kinetic coefficient of diffusion (Mobility)
µ Chemical potential
µth Thermodynamic “bulk” part of the chemical potential
N Number of molecules (or lattice sites)
P Pressure
pr Probability (rate) of interconversion in the HL model
ϕ Order parameter
R Phase domain size
Rid Ideal gas constant
ρ Density
S Structure Factor
T Temperature
q Wavenumber
qm Wavenumber corresponding to the maximum of the amplification factor
qs

m Wavenumber corresponding to the maximum of the structure factor
q− Wavenumber corresponding to the first root of the amplification factor
q∗ Wavenumber corresponding to the onset of microphase separation
q∗∗ Wavenumber corresponding to the termination of microphase separation
ω Amplification factor (phase domain growth rate)
Ωt Total grand thermodynamic potential
Ω = Ωt/N Grand thermodynamic potential per lattice site
x Molecular fraction
ξ Correlation length of order-parameter fluctuations

xxv


Chapter 1: Introduction

Typically, pure substances may be found with only one gaseous or liquid state, while their

solid state may exist in various polymorphic crystalline states. The existence of multiple dis-

tinct liquid forms in a single component substance is more unusual since liquids lack the long-

range order common to crystals. Yet, the existence of multiple amorphous liquid states in a

single component substance, a phenomenon known as “liquid polyamorphism,” [9, 55, 56] has

been observed or predicted in a wide variety of substances, such as superfluid helium [57, 58],

high-pressure hydrogen [1, 35, 37, 38, 59], high-density sulfur [2, 51], phosphorous [60, 61], car-

bon [62], silicon [63–66], silica [67,68], selenium and tellurium [69,70], and cerium [71]. Liquid

polyamorphism is also highly plausible in deeply supercooled liquid water [9, 13, 55, 56, 72–80].

A substance may be found to be polyamorphic by experimentally or computationally de-

tecting a liquid-liquid phase transition (LLPT), which terminates at a liquid-liquid critical point

(LLCP). For example, liquid polyamorphism, via the existence of two alternative supramolec-

ular structures, has been hypothesized to explain the remarkable anomalies in the thermody-

namic properties of supercooled water, namely a maximum in the temperature dependence of

its density and its isothermal compressibility [81, 82], a maximum in the isobaric heat capac-

ity [83], and an inflection point in its surface tension [84–87]. Simulations of water-like mod-

els [74–77, 79, 88–91], apparently supported by experiment [92, 93], have demonstrated the hy-

1


Figure 1.1: Hypothesized global pressure-temperature phase diagram for water. The black curves
indicate the solid-solid and solid-liquid-gas phase transitions experimentally reported for wa-
ter [11], while the blue curve indicates the hypothesized liquid-liquid phase transition. The dotted
black curve represents the line of spontaneous ice nucleation [12]. The red circle is the liquid-
vapor critical point (LVCP) and the red star is the hypothesized liquid-liquid critical point (LLCP)
around 220 K and 60 MPa as suggested in Ref. [13].

pothesized LLPT in supercooled water, see Fig. 1.1.

Fluid polyamorphism ultimately originates from the complex interactions between molecules

or supramolecular structures. However, it may be phenomenologically modeled through the re-

versible interconversion of two alternative molecular or supramolecular states [4, 9, 10]. The

application of this “two-state” approach to the variety of polyamorphic substances could be just

as a useful phenomenology or it may reflect the true microscopic origin of fluid polyamorphism.

Indeed, there are a few substances, such as hydrogen, sulfur, phosphorous, and carbon, where

the existence of alternative liquid or dense-fluid states is explicitly induced by a reversible chem-

ical reaction: dimerization in hydrogen [1] or polymerization in sulfur, phosphorus, and car-

bon [2]. While equilibrium phase transitions in simple fluids have been well-studied, the descrip-

2


tion of phase transformations in the presence of interconversion between alternative molecular or

supramolecular states, and in fluids far from equilibrium, is much less investigated [94–96].

The interconversion of two-states significantly effects both the thermodynamics and dy-

namics of fluid mixtures. In particular, molecular interconversion imposes an additional ther-

modynamic constraint, the chemical-reaction equilibrium condition, which reduces the number

of thermodynamic degrees of freedom. In this case, the concentration of a binary mixture is no

longer an independent variable, being a function of temperature and pressure. Consequently, a bi-

nary mixture with interconversion of species follows the Gibbs phase rule for a single-component

substance. In the absence of interconversion, the binary mixture exhibits a critical point for each

concentration. Collectively, these critical points make up a critical locus. The interconversion

of species selects a single path through the planes of phase coexistence at fixed concentration,

crossing the critical locus of the binary mixture at a unique point. Moreover, this path could cross

the critical line, the line of triple points, or any other unique line on the phase diagram more

than once. The resulting phase diagram can be viewed as the phase diagram of a polyamorphic

single-component fluid with multiple fluid-fluid critical points.

The thermodynamic properties, being state functions, under interconversion equilibrium

conditions are the same as those in the non-reacting binary mixture at the corresponding equi-

librium temperature, pressure, and overall concentration. Therefore, at any point along the ther-

modynamic path selected by interconversion, the properties are affected by the proximity to the

critical line of the non-reacting binary mixture. It should be noted that the evolution of the ther-

modynamic properties along the path (where the equilibrium concentration changes) is generally

very different from the evolution of the same property at constant composition, which would

typically be studied in non-reacting binary mixtures. This line in the interconverting mixture is

3


referred to as the “virtual critical line,” which may be used to describe the anomalous behavior

of the thermodynamic properties along this path.

The interconversion of species also effects the dynamics of phase formation. Generally, a

mixture of two liquids will phase separate if their interactions are not favorable for mutual mix-

ing. The most recognizable example is almost complete separation of water and oil. Another

example is a possible demixing of structural isomers, such as enantiomers with opposite molec-

ular chirality [7, 97–99]. In contrast, solid ferromagnetic (such as iron) or ferroelectric (such as

barium titanate) materials, in the absence of a magnetic or electric field, do not establish equilib-

rium coexistence between phases with alternative magnetizations or polarizations [94, 100, 101].

However, if liquids exhibit interconversion of species, similar to the flipping of magnetic spins or

electric dipoles, the thermodynamics and dynamics of phase separation will dramatically change.

In this thesis, it is shown that molecular interconversion may destroy or restrict liquid-liquid

phase separation. After a binary mixture, initially containing equal amounts of the alternative

molecules, is quenched from the one-phase homogeneous region at a high temperature into the

unstable region below the critical temperature of demixing at constant pressure, the species will

separate through a process known as spinodal decomposition [102]. However, if the two species

with the same density may rapidly interconvert (like in a mixture of enantiomers [7]), then to

avoid the formation of an energetically unfavorable interface, the phases will compete with each

other until one of them is eliminated [6]. This is the phenomenon of phase amplification, the

result of the competition between the diffusive dynamics of phase separation and the “flipping”

dynamics of interconversion. In magnetic and ferroelectric materials, away from the Curie tem-

perature and in the absence of a magnetic or electric field, the phenomenon of phase amplification

occurs naturally. In such materials, there is no restriction on the direction of magnetization or

4


polarization, meaning that there is no conservation of the number of magnetic spins or electric

dipoles with a particular orientation. Similarly, in fluids with fast molecular interconversion,

the conservation of the number of alternative molecules is broken and phase amplification could

occur.

In interconverting fluids, the growth of one phase at the expense of another stable phase

occurs to avoid the formation of an energetically unfavorable interface, and consequently, the

coexistence of the alternative phases is destroyed [6]. Since the interface is crucial to phase am-

plification, in macroscopic systems where the interfacial energy is much smaller than the bulk

energy, the system may enter a metastable state in which the energetic benefit of phase amplifi-

cation is negligible [4]. Thus, the size of the system plays a crucial role in phase amplification.

As discussed in Chapter 3, it is also shown that, in addition to smaller system sizes, a faster rate

of the interconversion reaction and closer proximity to the LLCP of demixing are necessary for

one to observe phase amplification [4, 6].

In a nonequilibrium system, if an external force causes the alternative molecules to stay

in equal numbers, the striking phenomenon of steady-state, restricted (“microphase”) separation

into mesoscale domains may be observed. In equilibrium, examples of mesoscale structures are

present in bicontinuous or spatially modulated microemulsions [103, 104] and microphase sepa-

ration of diblock or polyelectrolyte copolymers [105,106], where these mesoscale patterns are the

result of the minimization of the equilibrium free energy [9,107]. As discussed in Chapter 4, it is

found that the structure of the phase domains formed from nonequilibrium microphase separation

resembles modulated or bicontinuous microemulsion structures [9,103–107]. However, contrary

to the patterns formed in equilibrium or in metastable (“frozen”) conditions [108,109] (like those

commonly observed in glasses [110,111]), these nonequilibrium structures persist in steady-state

5


due to the continuous energy supply. Thus, steady-state microphase separation is one of the sim-

plest examples of dissipative structures in condensed matter. The characteristic length scale of

this dissipative structure emerges as a result of the competition between forced interconversion

and phase growth. If the source of forced interconversion is not sufficiently strong to overcome

the natural interconversion of alternative species, then the phenomenon of phase amplification,

the growth of one stable phase at the expense of another phase, is observed [4, 6, 112]. An ex-

ternal racemizing energy source can be achieved in physical systems through the interactions of

energy-carrying particles, such as photons, that may break intramolecular bonds [111]. It can

also be seen in biological cells through a flux of energy produced by ATP [113, 114] or it could

be achieved chemically through an external flux of matter or heat [115–117].

The size of the system, rate of forced interconversion, and proximity to the LLCP are also

important to observe microphase separation. In this thesis, it is found that there are two key

conditions for microphase separation to be observed. First, if the characteristic size of the meso-

scopic steady-state microphase domains is comparable to half the size of the system, then the

system will produce the same two alternative phases that would be observed without intercon-

version [4, 5, 8]. As a result, the size of the system may “cut off” the system’s ability to phase

separate into microdomains. Second, if the rate of forced interconversion is much faster than

the natural interconversion or diffusion rate, then the external force dominates the systems’ ki-

netics and no phase formation is possible. Consequently, it is shown that the size of the system

and the rate of the forced interconversion reaction are crucial components to observe microphase

separation [4, 5, 8].

This thesis is organized as follows: In Chapter 2, the two-state thermodynamic approach

to model fluid polyamorphism is discussed. The application of this approach is considered for

6


fluid-fluid transitions in hydrogen, sulfur, and substances exhibiting water-like anomalies. In

Chapter 3, the Cahn-Hilliard model of spinodal decomposition (phase separation process that

occurs after a binary mixture is quenched into the unstable region), is generalized to include

equilibrium (natural) and nonequilibrium (forced) interconversion of species. In Chapter 4, the

generalized Cahn-Hilliard theory is applied to describe the behavior of microscopic models of

mixtures, simulated through molecular dynamics (MD) or Monte Carlo (MC) methods, which

exhibit both natural and forced molecular interconversion. General concluding remarks and sug-

gestions for future research are presented in Chapter 5.

7


Chapter 2: Two-State Thermodynamics of Fluid Polyamorphism

In this chapter, the two-state thermodynamic approach for fluid polyamorphism is intro-

duced (Sec. 2.1). The approach is applied to describe the fluid-fluid transition in high-pressure

hydrogen (Sec. 2.2), high-density sulfur (Sec. 2.3), and fluids exhibiting water-like anomalies

(Sec. 2.4). In addition to the phase behavior of these substances, the approach is also used to

describe the anomalies in the interfacial properties of fluids exhibiting interconversion between

species.

2.1 Interconversion of Two States

Consider a symmetric binary mixture of two species A and B with molecular fractions,

xB = x and xA = 1 − x. Initially, it may be assumed that both species have the same densities

(ρ = 1), viscosities, and molecular weights. This system can be described by a Landau-Ginzburg

free-energy functional with a single order parameter uniquely linked to the fraction of species B

as, ϕ = x/xc − 1, where xc is the critical fraction. This functional reads as

F [{ϕ}] =
1

ρ

∫
V

(
Ĝ(ϕ, T, P ) +

1

2
κ|∇ϕ|2

)
dV (2.1)

8


Figure 2.1: A symmetric binary mixture is quenched along the critical molecular fraction (where
the order parameter is ϕ = 0) from a point, T̂1 in the one phase region to a point, T̂2 below the
critical temperature of demixing in the unstable region. The solid curve indicates the liquid-liquid
coexistence (binodal), while the dashed curve indicates the limit of stability (spinodal).

where the first term represents the thermodynamic “bulk” free energy and the second term is

included to describe the contribution to the free energy due to inhomogeneities within the system.

For an isotropic system, the coefficient κ is the square of the range of intermolecular interactions,

on the order of the square of the molecular size.

The bulk free energy density for the system, Ĝ, is the reduced Gibbs energy, Ĝ = G/kBTc,

where kB is Boltzmann’s constant and the critical temperature for liquid-liquid phase separation

is Tc. The bulk free energy density may be expressed in the general form,

Ĝ = ĜA + xĜBA + Ĝmix (2.2)

where ĜBA = µ̂0
B− µ̂0

A is the reduced difference between the Gibbs energies (chemical potentials)

of pure species A and B, referred to as the Gibbs energy change of reaction. For the symmetric

9


binary-lattice (“regular solution”) model, Ĝmix is formulated through the order parameter, ϕ, and

in the meanfield approximation, it reads as

Ĝmix = T̂

[(
1 + ϕ

2

)
ln

(
1 + ϕ

2

)
+

(
1− ϕ

2

)
ln

(
1− ϕ

2

)]
+
ε

4
(1− ϕ2) (2.3)

where T̂ = T/Tc is the reduced temperature and ε is a non-ideality interaction parameter, which

generally depends on temperature and pressure. The conditions for liquid-liquid phase equilib-

rium is

∂Ĝmix

∂ϕ

∣∣∣∣
T,P

=
T̂

2
ln

(
1 + ϕ

1− ϕ

)
− ε

2
ϕ = 0 (2.4)

The critical molecular fraction, xc, and critical temperature, Tc, are obtained from the thermody-

namic stability conditions

∂2Ĝ

∂ϕ2

∣∣∣∣
T,P

= 0 and
∂3Ĝ

∂ϕ3

∣∣∣∣
T,P

= 0 (2.5)

as

xc = 0.5 and Tc = ε/(2kB) (2.6)

Consider the interconversion between molecular states A and B through a reversible chem-

ical reaction of the form

A
k1−−⇀↽−−
k2

B (2.7)

where k1 and k2 are the forward and reverse reaction rates, respectively. Since, for this system, the

chemical reaction coordinate is the order parameter, the chemical reaction equilibrium condition

10


reads

∂Ĝ

∂ϕ

∣∣∣∣
T,P

=
∂Ĝmix

∂ϕ

∣∣∣∣
T,P

− ĜBA = 0 (2.8)

The reaction-equilibrium condition constrains the number of thermodynamic degrees of freedom

for the system. Consequently, the fraction of interconversion, given through the order parameter

ϕ, is no longer an independent thermodynamic variable, but instead, becomes a function of tem-

perature and pressure. The reduced Gibbs energy change of the reaction, ĜBA, can be expressed

through the reaction equilibrium constant, K(T̂ , P̂ ) = k1/k2, and is given, up to second-order in

temperature and pressure, as

ĜBA = −T̂ lnK(T̂ , P̂ ) = ê− ŝT̂ + υ̂P̂ + γT̂ P̂ +
1

2
δT̂ 2 − 1

2
κP̂ 2 (2.9)

where ê = e/(kBTc), ŝ = s/kB, and υ̂ = υ/(kBTc) are the reduced energy, entropy, and volume

changes of the reaction, while γ, δ, and κ are proportional to the volumetric expansivity, isobaric

heat capacity, and isothermal compressibility changes of the reaction, respectively [9, 10].

The reduced chemical potential for such a system undergoing spinodal decomposition to-

wards both chemical-reaction and phase equilibrium, µ̂ = µ/kBTc, is the reduced deviation of the

chemical-potential difference in solution (µ = µA − µB) from its equilibrium value, µ = 0. The

reduced time-dependent chemical potential is found from the functional derivative of Eq. (2.1) as

µ̂ = µ̂th − κ∇2ϕ (2.10)

where the order parameter depends on space and time, ϕ = ϕ(r, t). In this form, the chemical po-

tential is comprised of a thermodynamic potential as µ̂th = ∂Ĝ/∂ϕ|T,P , and a spatial-dependent

11


part, κ∇2ϕ. The thermodynamic component [118] is

µ̂th =
T̂

2
ln

(
1 + ϕ

1− ϕ

)
− ε

2
ϕ− ĜBA (2.11)

Expanding the logarithmic term to first order in ϕ via a Taylor series around ϕ = 0 (the value

of the order parameter at the initial time t = 0), gives ∂Ĝ/∂ϕ ≈ χ̂−1
q=0ϕ, where the inverse ther-

modynamic susceptibility in the limit of zero wavenumber is χ̂−1
q=0 = ∂2Ĝ/∂ϕ2. In the meanfield

regular-solution model, the inverse susceptibility scales as χ̂−1
q=0 ∼ ∆T̂ , where ∆T̂ = T/Tc − 1

is the reduced distance to the critical temperature. Therefore, the chemical potential defined in

Eq. (2.10) in the first-order approximation becomes

µ̂ ≈ ∆T̂ϕ− ĜBA − κ∇2ϕ (2.12)

It should be emphasized that the reduced chemical potential, given by Eq. (2.12), when the term

related to interconversion is absent, is the same as obtained in the classical Cahn-Hilliard the-

ory [102].

In addition to the chemical potential, the density and entropy are obtained from derivatives

of the Gibbs energy, Eq. (2.2), as

ρ̂(P̂ , T̂ ) =

(
∂Ĝ

∂P̂

)−1

ϕ,T̂

=

[
1

ρ̂A
+ xe

∂ĜBA

∂P̂

∣∣∣∣
ϕ,T̂

+
∂Ĝmix

∂P̂

∣∣∣∣
ϕ,T̂

]−1

(2.13)

S(P̂ , T̂ ) = −∂Ĝ
∂T̂

∣∣∣∣
ϕ,P̂

= SA − xe
∂ĜBA

∂T̂

∣∣∣∣
ϕ,P̂

− ∂Ĝmix

∂T̂

∣∣∣∣
ϕ,P̂

(2.14)

where ρ̂A = ρ̂A(P, T ) and SA = SA(P, T ) are the density and entropy (per molecule) of state A

12


and xe = x(P̂ , T̂ ) is the equilibrium fraction of species. Thus, from the appropriate specification

of state A and the equilibrium reaction constant, K, one may obtain all other thermodynamic

properties, such as the isothermal compressibility and heat capacity, as well as the global phase

diagram that includes both the vapor-liquid and the liquid-liquid transitions. In the next section,

this is demonstrated by applying the two-state approach to the fluid-fluid phase transition in high-

pressure hydrogen.

2.2 Fluid-Fluid Phase Transition in High-Pressure Hydrogen

Experiments and simulations have discovered that at extremely high pressures, highly-

dense fluid (dimeric) hydrogen dissociates into atomistic fluid hydrogen [21, 34–45, 48, 59, 119–

123]. In this section1, the two-state thermodynamic approach (Sec. 2.1) and the generalized

law of corresponding states (obtained from reducing the temperature, pressure, and entropy by

their critical values), in combination with the available experimental data and with the results of

computations [21, 39–46, 48, 124], are utilized to predict the equation of state of hydrogen near

the fluid-fluid phase transition (FFPT).

2.2.1 Phase Behavior of Hydrogen

The suggested global phase diagram of hydrogen, based only on the available experimental

evidence for the fluid-fluid phase transition [34–38], the solid-liquid melting transition [14–17,

125–127], and the location of solid-metallic hydrogen [22, 29–33], which is supported by the

most recent computational studies [21, 39–45, 47, 48, 124, 128–130], is shown in Fig. 2.2a. It

1This section was reproduced from Nathaniel R. Fried, Thomas J. Longo, and Mikhail A. Anisimov, J. Chem.
Phys., 157, 101101 (2022); https://doi.org/10.1063/5.0107043, with the permission of AIP Publishing.

13


Figure 2.2: The global pressure-temperature phase diagram for hydrogen. (a) The full range
from low to extreme pressures in logarithmic scale. The crosses indicate the experimental data
for the solid-liquid melting transition presented in Diatschenko et al. [14] (blue), Datchi et al. [15]
(cyan), Gregoryanz et al. [16] (pink), and Zha et al. [17] (purple). The solid black curves at low
pressure (P ≤ 0.1 GPa) are the liquid-gas-solid phase transitions [18], while the solid black
curve at high pressure (P > 0.1 GPa) is the Kechin equation [19] as reported in ref. [17]. The
black dashed curve is the predicted continuation of the melting line based on experimental and
computational evidence [17,20,21], while the dotted lines represent the highly-debated prediction
[22–28] of the domain of solid metallic hydrogen [22, 23, 28–33]. The red line is the first-order
fluid-fluid phase transition adopted in this section. (b) The phase diagram of hydrogen at extreme
conditions, in the area of the box in (a). The open circles are experimental data presented in
Zaghoo et al. [34–36] (dark brown), McWilliams et al. [37] (light brown), and Ohta et al. [38]
(orange). Simulation results [21, 39–46] are spread within the grey area and shown in detail in
Fig. 2.3. The fluid-fluid phase transition (solid red) and Widom line (dotted red) are represented
by Eq. (2.18). The red star is the location of the fluid-fluid critical point (FFCP) as adopted in
this section.

illustrates the fact that a huge pressure gap separates the liquid-gas [18] and fluid-fluid phase

transitions in hydrogen.

The adopted locations of the FFCP and the solid-fluid-fluid triple point (SFF-TP) in this

section are based on the available experimental data [34–38] and on discussions present in the

literature [42–44] (Table 2.1). It should be noted that the exact location of the FFCP is uncertain,

as the interpretation of both the Zaghoo et al. [35] and Ohta et al. [38] experimental data have

been highly debated [43, 131–133]. Most authors suggest that the experimental data of Ohta et

al. [38], on the anomalies of the heating efficiency, are obtained in the supercritical region [43].

14


Table 2.1: The suggested locations of the FFCP and the SFF-TP in hydrogen.

P [GPa] T [K] ρ [g/cm3]

FFCP 105 1900 0.8
SFF-TP 250 600 -

The results observed by Ohta et al. [38] have been interpreted [1] as being the anomalies of

the heating efficiency along the “Widom line”, the line corresponding to the maximum of the

fluctuations of the order parameter, which emanates from the critical point [9, 78, 134].

The significant discrepancy between the results of different computational models makes

it impossible to utilize these results for a single equation of state. However, presenting the same

results in reduced variables, as suggested by the law of corresponding states, allows the com-

putational results to be used along with the experimental data for thermodynamic modeling. In

Fig. 2.3, all of the available computational and experimental data on the fluid-fluid phase transi-

tion are presented in real units of pressure and temperature (Fig. 2.3a) and in reduced variables

(Fig. 2.3b), P̂ = P/Pc and T̂ = T/Tc, where Pc and Tc are the critical pressures and temper-

atures obtained (or adopted) from different works. It has been found that the simulation data

based on the Quantum Monte Carlo (QMC) approach could also be collapsed into the universal

phase diagram by reducing the entropy by its critical value, Ŝ = S/Sc [1]. In classical thermo-

dynamics, the reference value for the entropy is arbitrary. Commonly, the value, Sc, is adopted

as Sc = dP̂ /dT̂ |T=Tc [135–138], which was found to be Sc = 0.8 for all QMC simulations.

2.2.2 Phase Separation Coupled with Dimerization

Thermodynamically, the fluid-fluid transition in hydrogen can be modeled through the two-

state approach (Sec. 2.1), in which state A represents the free atoms of hydrogen and state B

15


represents dimerized hydrogen atoms. The total Gibbs energy per hydrogen atom (reduced by

RidTc, where Rid is the ideal-gas constant) is given by Eq. (2.2), where µ̂◦A and µ̂◦B are the Gibbs

energies of hydrogen in the monatomic or diatomic states, respectively, and x is the fraction of

hydrogen atoms in the diatomic state. Unlike the symmetric system considered in Sec. 2.1, the

mixing of free hydrogen atoms and diatomic hydrogen is not symmetric, and thus, Ĝmix may be

modeled as

Ĝmix(x) = T̂
[(x

2

)
ln
(x

2

)
+ (1− x) ln (1− x)

]
+ ε(T̂ , P̂ ) x(1− x) (2.15)

The dimensionless non-ideality parameter, ε = ε(T̂ , P̂ ), may be approximated up to first order

in ∆T̂ and ∆P̂ , as

ε(T, P ) = ε0 − εT∆T̂ + εP∆P̂ (2.16)

where ∆T̂ = T̂ − 1 and ∆P̂ = P̂ − 1.

The FFCP parameters are determined from the thermodynamic stability criteria (Eq. 2.5),

such that the critical fraction and critical temperature of hydrogen atoms is given by

xc =
√

2− 1 and Tc = 2(2−
√

2)2ε0 (2.17)

It should be noted that the first study to apply the two-state thermodynamic approach to high-

pressure hydrogen was presented by Cheng et al. [21]. While the predictions of Cheng et al.

for the FFPT are not in agreement with the results of all other simulations and experimental

studies [45], their study provides a reasonable idea for how the non-ideality parameter, ε, might

depend on pressure and temperature. Based on the suggested trend, εT and εP are optimized to

16


Figure 2.3: Unifying the different simulation results with experimental data of hydrogen by the
generalized law of corresponding states. (a) Experimental and simulation results for the fluid-
fluid phase transition (FFPT). (b) Unified representation of the FFPT by reducing pressure, P̂ =
P/Pc, temperature, T̂ = T/Tc, and the critical value of the entropy, Ŝ = S/Sc. In (a) and (b),
the open circles are the experimental data of Zaghoo et al. [34–36] (dark brown), McWilliams
et al. [37] (light brown), and Ohta et al. [38] (orange). The computational results are indicated
by the triangles: blue tints correspond to the Density Functional Theory (DFT) simulations of
Bonev et al. [47] (dark blue), Morales et al. [48] (blue), Hinz et al. [44] (sky blue), and Karasiev
et al. [45] (light blue). Meanwhile, green tints correspond to the Quantum Monte Carlo (QMC)
simulations of Morales et al. [48] (dark sea green), Lorenzen et al. [39] (green), Perlioni et al. [41]
(dark green), Mozzola et al. [42] (lime green), and Tirelli et al. [46] (yellow green). The colored
stars correspond to the reported (or adopted in this section) critical points for each data set. The
solid black curve is the solid-fluid phase transition line as discussed in Fig. 2.2, and the red solid
line is the FFPT predicted in this section.

agree with the behavior of hydrogen from the available computational data [41–44, 46, 48], and

consequently, these parameters were adopted as εT = 2.062 and εP = −0.175. The asymmetric

Gibbs energy of mixing is illustrated in Fig. 2.4a along with the fluid-fluid coexistence, calculated

via the common tangent method, and the limit of absolute stability (spinodal), calculated via the

thermodynamic stability conditions.

The condition for chemical-reaction equilibrium is given by ∂Ĝ/∂x|T,P = 0, resulting in

the balance of the Gibbs energy of reaction, ĜBA, and the exchange chemical potential of mixing,

17


µ̂mix = ∂Ĝmix/∂x|T̂ ,P̂ , such that thermodynamic equilibrium follows from

ĜBA = −µ̂mix (2.18)

The Gibbs energy of reaction, ĜBA = ĜBA(T̂ , P̂ ), may be approximated up to second order in T̂

and P̂ , as in Eq. (2.9). To balance the Gibbs energy of reaction, Eq. (2.9), with the derivative of

the Gibbs energy of mixing, ĜBA is expressed as an expansion in ∆T̂ and ∆P̂ as

ĜBA = u− a∆T̂ + b∆P̂ + g∆T̂∆P̂ +
d

2
(∆T̂ )2 − k

2
(∆P̂ )2 (2.19)

where the modified coefficients of the thermodynamic balance, Eq. (2.19), are related to the

coefficients of reaction, Eq. (2.9), as:

ê = u+ a− b+ g +
d

2
− k

2

ŝ = a+ g + d

υ̂ = b− g + k

(2.20)

along with γ = g, δ = d, and κ = k.

If the Gibbs energy of mixing, Ĝmix, would be symmetric with respect to x, then ĜBA =

−µmix = 0, could describe the conditions for both reaction equilibrium and fluid-fluid phase

equilibrium [9]. However, since the monatomic and diatomic mixing is asymmetric, the condition

for the balance of phase and reaction equilibrium, Eq. (2.18), is given through

µ̂mix

T̂
= a2

(
ε(T, P )

T̂
− ε0

)2

+ a1

(
ε(T, P )

T̂
− ε0

)
+ a0 (2.21)

18


Figure 2.4: The components of the Gibbs energy (per atom) for hydrogen in the vicinity of the
fluid-fluid phase transition. (a) The Gibbs energy of reaction, GBA, as given by Eq. (2.9). The
isotherms are T = 0.5Tc (orange), T = 0.75Tc (blue), T = Tc (green), T = 1.25Tc (red), and
T = 1.5Tc (purple). (b) The Gibbs energy of mixing, Gmix, as given by Eq. (2.15). Gmix is shown
as a function of the fraction of hydrogen atoms in the diatomic state, x, for isotherms T = 0.5Tc

(blue), T = 0.75Tc (green), and T = Tc (red) at P = Pc. The solid curve corresponds to the
fluid-fluid coexistence.

where the coefficients a0 = −0.502, a1 = 0.166, and a2 = −0.071.

The developed equation of state is formulated through the Gibbs energy for the system

as a function of temperature and pressure. Due to the interconverting nature, the two-states of

hydrogen are thermodynamically equivalent to a single component system. Consequently, this

produces an equation of state in terms of the equilibrium fraction of dimerized atoms, xe =

xe(T, P ), and the density of the system, ρ = ρ(T, P ). The equation of state presented in this

section contains seven adjustable parameters: five from the Gibbs energy of reaction, GBA, (u,

a, b, g, and k), Eq. (2.9) and two from the non-ideality parameter in the Gibbs energy of mixing

(εT and εP ), Eq. (2.16). The number of adjustable parameters are reduced from the following

analysis of the available computational data on hydrogen in the vicinity of the fluid-fluid critical

point.

From the computational heat capacity data presented by Karasiev et al. [45], the heat-

19


Figure 2.5: Equilibrium fraction of hydrogen atoms in the diatomic state, xe. (a) Equilibrium
fraction-pressure diagram for T = 0.5Tc (orange), T = 0.75Tc (blue), T = Tc (green), and
T = 1.25Tc (red). (b) Equilibrium fraction-temperature diagram for P = 0.75Pc (blue), P =
Pc (green), P = 1.25Pc (red), P = 1.5Pc (purple). The solid and dashed black curves are,
respectively, the fluid-fluid coexistence and the limit of thermodynamic stability (spinodal).

capacity change of reaction is approximated to be δ ≈ 0, and from the computational isothermal-

compressibility data presented in the supplemental material of Pierleoni et al. [41], it is approx-

imated that κ ≈ 0.625 [mm3/GPa ·mol]. Additionally, e = −108 [kJ/mol] is adopted based on

the known value of the bond dissociation energy of H2 [139]. As discussed above, εT and εP are

assigned as εT = 2.062 and εP = −0.175.

From these findings, the number of free parameters have been reduced to three: a, b, and

g. The values of the remaining free parameters are determined as a = −4.95, b = 0.044,

and g = 0.0124 from the computational and experimental data utilizing the generalized law of

corresponding states (Fig. 2.3). Using the relations between these parameters and the physical

parameters in Eq. (2.9), the following parameters are estimated: the entropy change of the reac-

tion as s = −34.0 [J/K ·mol], the volume change of the reaction as υ = 393 [mm3/mol], and

the volume-expansivity change of the reaction as γ = 0.0677 [mm3/K ·mol]. The Gibbs energy

change of reaction is shown in Fig. 2.4b. It demonstrates that the pressure is the major factor in

20


the behavior of GBA.

Using the Gibbs energy of mixing, Eq. (2.15), the Gibbs energy of reaction, Eq (2.9), and

the variables determined from the universal phase diagram, the equilibrium fraction of hydrogen

atoms in the dimerized state, xe, is determined from Eq. (2.18). The corresponding equilibrium-

fraction phase diagrams are presented in Fig. 2.5(a,b). At higher temperatures and lower pres-

sures, the equilibrium composition changes from the dimeric state xe = 1 to the monomeric state

xe = 0.

The density of species is expressed through the equilibrium fraction via Eq. (2.13), and

with use of Eq. (2.15) for Ĝmix is given in the form [9]

ρ̂(P̂ , T̂ ) =

(
∂Ĝ

∂P̂

)−1

T̂

=

[
1

ρ̂A
+ xe

∂ĜBA

∂P̂
+
∂ε

∂P̂
xe (1− xe)

]−1

(2.22)

where ρ̂A = ρ̂A(P̂ , T̂ ) is the volume of the monatomic hydrogen state, and may be expressed up

to second-order in ∆T̂ and ∆P̂ as

ρ̂A = ρ̂c − ρ̂0∆T̂ + ρ̂1∆P̂ + ρ̂2∆T̂∆P̂ − ρ̂3

(
∆P̂
)2

+ ρ̂4

(
∆T̂
)2

(2.23)

Using the most recent QMC simulations presented in Tirelli et al. [46], ρ̂A is estimated by

Eq. (2.23) with coefficients: ρ̂c = 1.01, ρ̂0 = 0.25, ρ̂1 = 0.56, ρ̂2 = 0.56, ρ̂3 = 0.21, and

ρ̂4 = 0.12. The corresponding pressure-density phase diagram is presented in Fig. 2.6, and

demonstrates a good agreement with the computational data in the vicinity of the FFCP.

It is noted that the properties observed in experimental studies (e.g. conductivity, reflec-

tivity, thermal efficiency, etc.) could be indirectly related to the order parameter for the FFPT in

21


Figure 2.6: The pressure-density phase diagram of hydrogen based on the equation of state de-
veloped in this section. The open circles correspond to the QMC simulations of Tirelli et al. [46].
Isotherms are T̂ = 0.67 (orange), T̂ = 0.73 (red), T̂ = 0.8 (brown), T̂ = 0.87 (purple), T̂ = 0.93
(green), and T̂ = 1.0 (blue). The fluid-fluid coexistence is shown by the solid black curve. The
red star is the FFCP adopted in this section.

hydrogen, ϕ = x/xc−1. The measurable quantities (such as density or conductivity) are coupled

to the order parameter.

There is a remarkable analogy between the challenges in thermodynamic modeling of fluid

polyamorphism in hydrogen and that in supercooled water. In both cases, there is a reasonable

agreement among the scientific community on the shape and location of the first-order transition

line, while the position of the corresponding FFCP and LLCP are highly uncertain and subjects

of current debate in the literature [13, 34–36, 43, 78, 131–133]. This uncertainty, in both hydro-

gen and water, is due to the extreme conditions of the phenomena. In supercooled water, the

liquid-liquid transition is hidden below the temperature of spontaneous ice formation [75, 78],

22


while in hydrogen, the fluid-fluid transition occurs at immensely high pressures (millions of

atm) [43]. Consequently, it is not surprising that the available computational or experimental

data are scarce [40, 43, 44]. In this section, it has been shown that despite the uncertainty in de-

termining the location of the FFCP in hydrogen, thermodynamic modeling provides a principle

direction to predict the equation of state for the system. Remarkably, the law of corresponding

states can be utilized to reconcile the different computational models of hydrogen and experi-

ment [21, 34–45, 48, 124] into a unified equation of state. An additional parameter has been in-

cluded to generalize the law of corresponding states, the entropy at the critical point (Sc), which

provides the opportunity for further studies of hydrogen, both experimental and computational,

to be unified under the general approach presented in this section.

2.3 Liquid-Liquid Phase Transition in High-Density Sulfur

While the fluid-fluid phase transition in hydrogen was induced by the dimerization reac-

tion, a fluid-fluid phase transition may also be induced by more complex chemical reactions or

interactions, such as polymerization or gelation. An example of polymerization induced fluid

polyamorphism was recently discovered in high-density sulfur [51]. A polymerized state of sul-

fur exists above 433 K [140–144]). It has been experimentally demonstrated