The re-entrant transition from the molecular to atomic phases of dense fluids: The case of hydrogen

A simple phenomenological thermodynamic model is developed to describe the chemical bonding and unbonding in homonuclear diatomic systems. This model describes the entire phase diagram of dimer-forming systems and shows a transition from monomers to dimers, with monomers favored at both very low and very high pressures, as well as at high temperatures. In the context of hydrogen, the former region corresponds to hydrogen present in most interstellar gas clouds, while the latter is associated with the long sought-after fluid metallic phase. The model predicts a molecular to atomic fluid transition in dense deuterium, which is in agreement with recently reported experimental measurements.


I. INTRODUCTION
The study of diatomic molecules has a long history, dating back to the foundation of modern chemistry.Molecules comprising only two atoms of the same kind (homonuclear) are in some sense the first step from physics into chemistry, having covalent bonds between the atoms and one more constraint than monoatomic species (the bond length).Interest in them dates back to Avogadro, who first proposed their existence as a means to reconcile measurements of the volumes of hydrogen and oxygen needed to form water with atomic theory. 1 The central question addressed here is "under what conditions do molecules break into their constituent atoms?"In general, this is a competition between the entropy of the atomic state and the energy of the molecular state, with the density also playing a role at high pressures.
6][37] In recent years, there has been a shift in the way we regard the formation of high pressure hydrogen, with density functional theory suggesting that the favored mechanism for metallization comes via molecular metals and an open crystal structure rather than full ionization, as traditionally expected.This significantly reduces the temperatures needed to dissociate and hence metallize fluid hydrogen.The most widely used tool for the theoretical investigation of high pressure behavior is first-principles calculations so that the vast majority of studies to date have employed density functional theory (DFT) or quantum Monte Carlo methods.Such calculations indicate that the The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp9][40] A detailed study along the same lines by Geng et al. 41 reinforced the idea of a critical point but also highlighted the strong dependencies of the predictions on the exchange-correlation functional used, 42 finite-size effects, as well as the miscibility of the atomic and molecular fluids.These studies base the claim of an atomic-molecular phase boundary on numerical indications of a discontinuity in density of pressurization.Contrariwise, a recent study using machine-learned potentials and thermodynamics 43 claimed that the transition is continuous, with no critical point.
Fried and co-workers develop a thermodynamic model to describe the transition. 44This approach uses a convex free energy function for non-ideal mixing to obtain a phase separation in a way similar to the Bragg-Williams model 45 in alloys.They incorporate a six-parameter reaction energy expanded to second order in pressure and temperature and a pressure-and temperature-dependent mixing enthalpy.The parameters are fitted to the experimental metallization curve and to DFT calculations for hydrogen, obtaining a model that exhibits a fluid-fluid critical point associated with the atomic to molecular transition.
In this article, we take an alternate, more phenomenological approach to give insights into the high pressure behavior of systems that form diatomic molecules, using a simple microscopic picture of the underlying microscopic processes to develop a thermodynamic model to describe the chemical bonding and repulsive interactions between atoms in dimer-forming species.Our aim is to map out the possible conditions for the atomic-molecular transition, building a general model that can be parameterized for particular materials.Such a model will enable us to not only study the dense phases, but also link continuously to extreme low density, such as the interstellar medium, 46 which would be completely intractable to DFT.
The key physical mechanisms in these systems that we want to qualitatively capture are the tendency of pairs of atoms to reversibly form (and break) chemical bonds, the repulsive interactions between atoms that tend to keep them from overlapping, and the difference of these repulsive interactions for atoms that are bonded or not bonded.
In the absence of non-bonded interactions between atoms, increasing temperature will tend to decrease the number of diatomic molecules, with the system favoring the increase in entropy due to the increase in mobility of unbonded atoms over bonded atoms over the lower energy of the bonded molecules.Applying pressure to the system will tend to promote bond formation, as well as increase the system density, as the balance shifts toward the bond energy rather than entropy.
Adding repulsive interactions between atoms will tend to enhance the effect of pressure, although the results will be qualitatively similar to that in the absence of interactions.However, if the effective volume occupied by bonded atoms is larger than that occupied by non-bonded atoms (i.e., the repulsion between bonded atoms is greater than that between non-bonded atoms), then increasing pressure will tend to favor the breaking of bonds.
Using different "sizes" for bonded and non-bonded atoms can lead to a qualitative change in the effect of pressure on diatomic molecule formation.The frequency at which bonds are broken or formed with changing pressure can be adjusted by making the repulsive interactions non-additive.The addition of positive nonadditivity, where the repulsion between dislike species is greater than that between like species, leads to a tendency for the two species to demix.Interestingly, non-additive hard mixtures have been shown [47][48][49] to exhibit a fluid-fluid phase separation.The larger the non-additivity, the greater the propensity for the two species to demix.Non-additive repulsive interactions can also arise from differences in the electrostatic screening behavior of various species in mixtures, 50,51 which is thought to lead to phase separation in mixtures of hydrogen and helium.
In this work, we develop a simple, microscopically motivated thermodynamic model to describe the high-pressure behavior of dimer-forming species, using a mean-field approach to account for the non-additive repulsion interactions between atoms.This approach leads to a theory that resembles the two-state model approach. 52,53The interplay between these elements provides a theory that captures the main features of the behavior experimentally observed through indirect measurements and predicted by density functional calculations for hydrogen at high pressures.We note in passing that the incorporation of chemical reactions with free energy models is not new and, in fact, is fundamental in the thermodynamic analysis of chemical reaction equilibria.This is typically done using an assumption of ideality (e.g., ideal gas or ideal solution) in thermodynamics; however, there is much of work where chemical reactions have been analyzed in combination with non-ideal free energy models. 54These are important for reactions under extreme conditions, such as at high temperatures and/or pressures, and processes where the non-ideality of the system plays a significant role, such as reactive distillation. 55,56This approach has also been used to construct models to account for systems where the molecules can strongly associate with one another, such as through hydrogen bonds, and include statistical association fluid theory (commonly referred to as SAFT) [57][58][59][60][61] and cubic plus association models. 62e want the theory to capture the very low pressure gas region, which is primarily diatomic at low temperatures and becomes increasingly monoatomic with increasing temperature.More importantly, however, we want to describe the relatively sharp transition of the system from consisting primarily of diatomic molecules to primarily unbonded atoms at high pressures and relatively moderate temperatures.
In addition, this transition has an associated critical point, where below the critical temperature, the transition is first order and discontinuous, and above the critical temperature, it is gradual.We explore the general topologies of phase diagrams the model produces and highlight key qualitative and quantitative changes arising from changes in a small number of parameters that have a major impact on the general aspect of a phase diagram.We show that by making recourse to existing well-measured data for H 2 at ambient pressure, we are able to qualitatively reproduce the highpressure molecular-to-atomic transition in hydrogen including its critical point.By a modest amount of fine-tuning, our model also becomes able to quantitatively reproduce the most recent experimental measurements of this behavior for deuterium in the region of ∼1000 K by Knudson et al. using the Z-machine 63 and Celliers et al. at the National Ignition Facility. 64Due to its simplicity, our model can easily be transferred to other chemical species by fitting to their known parameters, as well as extended to quantitatively The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcpdetermine parameters of interest, such as density jumps, the bulk modulus, and the isothermal heat capacity.In future, more sophisticated approximations can be used to better relate the parameters of the thermodynamic model to an underlying molecular interaction model.The remainder of this article is structured as follows.We first present the basic mean field model in Sec.II.This model is then formulated in terms of the Gibbs free energy in Sec.III.In Sec.IV, the qualitative behavior of the model is analyzed and compared to experimental systems.Finally, the main findings of this work are summarized in Sec.V.

II. FORMULATION IN TERMS OF THE HELMHOLTZ FREE ENERGY
We consider a system of atoms that can exist either as unbonded monomers (labeled as a) or bonded together in pairs as dimeric molecules (labeled as m).This system of atoms can be considered to be a binary mixture of monomers and dimers.A fraction x of the atoms are bonded to another atom (i.e., within a molecule), while the remaining fraction 1 − x of these atoms are unbonded (i.e., remain in monomeric form).Consequently, for a system of N atoms, the number of molecules is Nm/2 = Nx/2, and the number of monomers is Na = N(1 − x).It follows that N = Na + Nm is independent of x.
The Helmholtz free energy atom F of the system is a function of the absolute temperature T, the volume per atom v, and the amount of each species.It can be divided into an ideal gas contribution F ig and a residual contribution F res as follows: where T is the absolute temperature.The ideal gas contribution is given by the following standard expression: where p 0 = 1 bar is the reference pressure, β = (kBT) −1 , kB is the Boltzmann constant, μ ○,ig a (T, p 0 ) is the formation free energy of the unbonded atomic species, μ ○,ig m (T, p 0 ) is the formation free energy of the diatomic molecular species, and ρ = N/V is the total number density of atoms (either bonded or unbonded) in the system.The terms μ ○,ig a (T, p 0 ) and μ ○,ig m (T, p 0 ) account for the contribution of the internal degrees of freedom (e.g., rotation, vibration, or electronic energy states) of the atoms and molecules, respectively, to the free energy.
The residual free energy [Eq.(5)] accounts for interactions between atoms that are not directly related to chemical bonding.To obtain an approximate expression for this term, we propose a simplified picture for the system of atoms, which consists of a mixture of positively charged nuclei and negatively charged electrons.Taking a coarse-grained perspective, this can be viewed more simplistically as a collection of particles, composed of a nucleus and partially bound electron(s), that have an effective charge.The electrons can be considered to be either "bound" to a specific nucleus (making its effective charge significantly lower than that of the bare nucleus), participating in chemically bonding two nuclei together (forming a diatomic molecule) or as "free" particles that screen the electrostatic potential in the system.The "free" electrons are treated as a neutralizing continuum background through which the particles interact.
There are two types of particles in our model: atoms that are monomeric, which are denoted by the subscript a, and atoms that are part of a diatomic molecule, which are denoted by the subscript m.Note that each of these types of particles can, in principle, have a distinct effective charge.The interaction energy u αα ′ (r) between these effective particles can be taken to be of the Yukawa form as follows: where r is the distance between particles, e 0 is the fundamental unit of charge, ε 0 is the permittivity of free space, κ −1 is a screening length, Z is a nominal effective charge of an atom, and Mαα takes into account any differences in the effective charge due to the state of the atom (e.g., unbound or chemically bonded).
The parameters Maa and Mmm characterize the relative magnitudes of the effective charges of the atoms.The quantity Maa characterizes the size of unbonded atoms.In the present work, this is taken as a reference and Maa = 1, so an unbonded atom has an effective charge of Z.The effective charge of an atom in a diatomic molecule is Z √ Mmm.This may be expected to differ from that of a monomeric atom, as the electrons that are shared between nuclei in a diatomic molecule may screen their charge in a different manner.
The quantity Mam characterizes the interaction between a nonbonded atom and an atom within a molecule.Typically, one would expect that this would be equal to the geometric mean of Maa and Mmm.However, there may be deviations from an expected behavior, as the electrons within molecules may be correlated in a different manner than those in monomeric atoms.Consequently, we write where δ is a parameter that characterizes the deviation from the expected mixing rule. 51For δ < 0, two similar atoms will be repelled from each other more than they will be repelled by a different type of atom.This will tend to cause unbonded atoms to mix with molecules.For δ > 0, two similar atoms will repel each other more that they will repel a different atom.This will tend to cause different types of atoms to de-mix (i.e., unbonded atoms will tend to separate from molecules).The screening length κ −1 provides a length scale for the theory, and the combination ε = Z 2 e 2 0 κ 4πε 0 provides an energy scale.Given values for these parameters, the nominal effective charge Z can be determined as Z = √ 4πε 0 εκ −1 /e 0 .The interaction potential is more conveniently written as The Journal of Chemical Physics

ARTICLE pubs.aip.org/aip/jcp
Within the mean field approximation, the average interaction energy between a particle of type α and particles of type α ′ surrounding it is where v is the specific volume.
Summing over the combinations of α and α ′ gives the expression for the extensive residual Helmholtz free energy, where the parameter a is defined as Note that a is an intensive variable that only depends on the fraction x of atoms that are within molecules.It has no explicit dependence on the temperature or volume of the system.The final expression for the Helmholtz free energy is This is the free energy of an ideal gas of atoms, an ideal gas of molecules plus the rotational/vibrational energy of the molecules [contained within μ ○,ig m (T, p 0 )], an entropy of mixing, and a van der Waals-like interaction term.The molecules and atoms are implicitly treated as distinguishable.Later, we will include chemical bonding.

A. Thermodynamic properties
All thermodynamic properties of the system can be determined from the free energy.The pressure of the system can be determined from the equilibrium free energy as The first term is simply the ideal gas pressure for Na + Nm objects.
The second term has the same form as the interactions in the van der Waals equation.Thus, with fixed x, the model reduces to the van der Waals equation with no excluded volume correction.In the Yukawa picture, the parameter a describes the magnitude of the repulsive interaction.
The chemical potential of monomeric atoms in the system is given by The chemical potential of atoms within a molecule is given by

B. Other thermodynamic parameters in relation to the Helmholtz free energy
The entropy of the system is where S ○,ig a (T, p 0 ) and S ○,ig m (T, p 0 ) are the entropies of isolated individual atoms and molecules, respectively, due to their internal degrees of freedom (e.g., those associated with different electronic states, rotations, and vibrations) The isochoric heat capacity is The bulk modulus is

C. Chemical bonding
Previously, we developed an approximate free energy for a mixture of atomic and molecular particles, based on a simple mean field interaction model.A key parameter in the theory is the free energy of bond formation ΔG bond (T), which is given by the difference in the formation Gibbs free energies of a dimer/molecular species and two individual atoms, The Journal of Chemical Physics where the subscripts "m" and "a" refer to the molecular and atomic species, respectively.This contains the free energy of the covalent bonding, oscillation and rotations of the molecule, and also the difference in ideal gas free energy from halving the number of particles.(Following convention, ΔG is an intensive quantity, defined per molecule not per atom.) In terms of the bonding free energy, the Helmholtz free energy can be simplified as To get a feel for the order of magnitudes to expect for the bond formation term, this quantity is shown in Fig. 1 for hydrogen, deuterium, and oxygen.The slope of the line is the bond entropy ΔS bond (T), while the curvature of the line is related to the bond heat capacity ΔC p,bond (T).The bond enthalpies for hydrogen, deuterium, and oxygen are shown in the inset of Fig. 1.They slightly decrease as the temperature increases.Note that the slopes of these curves are the "bond heat capacities," which include a contribution from reduced translational kinetic energy of one molecule vs two atoms.For hydrogen at 25 ○ C, the bond enthalpy is ∼436 kJ mol −1 (4.52 eV), for deuterium, it is 443 kJ mol −1 (4.60 eV), while for oxygen, it is about 498 kJ mol −1 (5.16 eV).

III. FORMULATION IN TERMS OF THE GIBBS FREE ENERGY
Here, we relate the Helmholtz free energy model developed in Sec.II to a two-state model. 52In this approach, the problem formulated within the Gibbs ensemble, where the independent variables of the system are temperature, pressure, and the total number of atoms in the system.The corresponding free energy, in this case, is the Gibbs free energy, which can be determined from the usual Legendre transformation of the Helmholtz free energy, Equation ( 8) is a quadratic in 1/v, which can be solved to give an expression for the volume v in terms of pressure, Note that this reduces to the ideal gas law in the limit a → 0. The specific Gibbs free energy is given by This is the key result in this paper, which we will later use to analyze the behavior of dimerizing systems.The right side contains five contributions to the free energy.The first term is the free energy of an ideal gas of monomeric atoms at fixed pressure p 0 and absolute temperature T. The second term is the free energy of forming a diatomic molecule, which incorporates covalent bond formation, molecular vibrations and rotations, and the change in free energy from converting x/2 atoms to molecules.The third term accounts for the ideal entropy of mixing.The fourth term represents the confinement entropy (arising from volume/pressure changes) of the ideal gas.The final term arises from the non-bonding intermolecular interactions [see Eq. ( 5)].

A. Thermodynamic properties
Knowledge of the Gibbs free energy allows for the calculation of all the thermodynamic properties of the system.For example, the volume v per atom in the system is given by which leads to the expression given in Eq. ( 17).
The specific enthalpy per atom in the system H can be determined from the Gibbs-Helmholtz relation, a (T, p 0 ))/∂β is the enthalpy of atom in the ideal gas state and ΔH bond (T) = ∂(βΔG bond (T))/∂β is the bond enthalpy.The first two terms represent the enthalpy associated with the internal degrees of freedom of the atoms and molecules in the system.The third term is the per-atom specific ideal gas enthalpy of (1 − x/2)N particles from N atoms.The final term is the contribution of non-bonded interactions between atoms [Eq.( 5)].
The chemical potential of unbonded atoms is and the chemical potential of molecules (which is twice the chemical potential of an atom in a molecule) is

B. Reaction equilibrium
In Sec.III A, we inherently assumed that the monomer and dimer forms of the atoms did not interconvert; however, we can easily relax this assumption.At equilibrium, the Gibbs free energy of a closed system held at constant temperature and pressure is minimized with respect to any internal degrees of freedom, which in this case is the fraction x of atoms that are within molecules.In other words, the equilibrium Gibbs free energy is where x * is the value of x that minimizes the Gibbs free energy at fixed T, p, and N. Note that it is the equilibrium fraction of atoms that are within molecules.
If the Gibbs free energy is convex in the variable x, then there is a unique minimum in the NPT ensemble corresponding to This equation can be used to solve for x * .Note that for our approximate Gibbs free energy, given in Eq. ( 18), we find where If we note that the number of monomeric atoms is Na = (1 − x)N and the number of molecules is Nm = Nx/2, then So, we see that the condition for equilibrium is given by which is the standard condition for chemical reaction equilibrium.This is equivalent to the condition that the chemical potential of unbonded atoms is the same as that for atoms within molecules.

C. Equilibrium properties
The equilibrium Gibbs free energy is simply equal to the Gibbs free energy evaluated at x * .As for systems with "fixed" composition, properties that are related to the first derivatives of the Gibbs free energy, such as the volume, enthalpy, or species chemical potentials, are simply these properties evaluated at x * , v eq (T, p) = v(T, p, x * ), H eq (T, p) = H(T, p, x * ), μ eq m (T, p) = 2μ eq a (T, p) = μm(T, p, x * ).

IV. RESULTS AND DISCUSSIONS
Here, we examine the types of qualitative thermodynamic behavior that can be exhibited by the model described above.The key physical parameters of the model relate to the repulsive interactions between the atoms in the system and the strength of the chemical bond.The parameter ε sets the unit of energy, and κ −1 sets the unit of length.As mentioned previously, we set Maa = 1.
To simplify calculations, we assume that the free energy of bonding ΔG bond (T) is linear in temperature, where the parameter A is related to the bond entropy (i.e., ΔS bond /kB = A and the parameter B is related to the bond enthalpy.The results of fitting Eq. ( 22) to the data in Fig. 1 for different temperature ranges are given in Table I.
In the analysis that follows, we assume that A and B are known parameters (fixed from experiment) and treat the remaining parameters of the model as fitting parameters with which we can explore possible phase diagrams.Those free parameters are Mmm the relative repulsion of an atom in a molecule in comparison with an unbonded atom, δ the non-additive mixing parameter between monomers and dimers, A, and B.  The first step in evaluating the performance of our model is to investigate the default behavior, in the absence of tight fits to thermodynamic data available for real chemical systems, inclusion of known physical changes at the atomic level when molecule formation occurs or potential immiscibility.In practical terms, keeping our focus on hydrogen for which we have empirical data on the free energy of bond formation (βΔG bond ), this means our choices for model parameters are: Maa = 1.0 = Mmm, δ = 0, ε = 3.627 eV, κ −1 = 2 Å, A = 14.1, and B = −4.65 eV (for βΔG bond ).Note that for these parameter values, Z = 0.710, and the effective charge of both unbonded and bonded atoms is 0.710e 0 .
The phase diagram produced by this choice of parameters is depicted in Fig. 2 for temperatures ranging from 0 to 30 000 K and pressures ranging from ambient to 1.7 TPa (17 Mbar).The model does not display any distinct pressure-induced molecular-to-atomic transition.The main feature is the entropy-driven breaking of molecule at increased temperature.There is also a weak preference for molecules at pressure.
In the following, we examine the influence of each of the parameters in turn.First, we explore the influence of the temperature dependence of the bond free energy.Then, we examine the influence Mmm, which characterizes the relative repulsion between two unbonded atoms and two bonded atoms in different molecules.Next, we study the influence of the parameter δ.Finally, we relate the model to hydrogen.
A. Limits: Non-interacting atoms, high pressure, and zero temperature Before engaging in any exploration of our model, it is worth checking the behavior of certain limiting scenarios.The most relevant ones are the non-interacting limit, which essentially should recover an ideal gas; the high pressure limit, to ensure a physically reasonable T(P) curve without divergences; and the zero temperature limit.
Let us now examine the first of these cases.In the limit that a → 0 [Eq.( 6)], the system reduces to an ideal gas with N(1 − x/2) particles.In this situation, the condition for equilibrium given by Eq. ( 19) simplifies, and we recover the standard law of mass action for ideal gases, where (x/2)/(1 − x/2) is the fraction of dimers and ) is the fraction of dissociated atoms.This can be rearranged, which yields the following expression for the temperature: We now turn our attention to the limit of extremely high pressures.In this situation, where p → ∞, from Eq. ( 17), we find that the specific volume varies as The specific volume is dominated by the interatomic repulsion in this regime.This relation allows us to simplify the derivative of the Gibbs free energy [see Eq. ( 20)] to In addition, if we take again βΔG bond (T) = A + Bβ as the functional from for Eq. ( 14) and let x = 1/2, we obtain According to our definitions, we have The In the final equation for βε obtained above, let us denote the square bracket by χ and the right-hand side of the equation by Q.The equation can then be rewritten as where W k is the Lambert W function.For χ < 0, the solution is given by W −1 , while for χ > 0, the solution is W 0 .A graphical depiction of these behaviors (ideal gas and high pressure limits) and the effect of various approximations on them is presented in Fig. 3. Finally, to ensure the physicality of our model, we consider the zero temperature limit, where the Gibbs free energy becomes the enthalpy, The first term is independent of x, so we have a competition between molecular bonding, interparticle interactions, and pv.Typically, ΔG bond (0) will favor molecules.The interparticle term becomes more important with pressure and will favor whichever species is favored by a (i.e., the smaller of Maa vs Mmm).For δ = 0, Maa = Mmm, the parameter a becomes independent of x, and this term has no effect.Finally, the explicit pv term favors whichever is smaller out of two atoms or one molecule.Typically, we can only expect a zero-T transformation if the atomic form is denser than the molecular one, which here means Maa < Mmm."default" case with the exception of Mmm, which will now be greater than Maa.Mma is linked to Mmm via δ, so also increases.
The fraction of atoms in molecular form in the system at different temperatures and pressures is presented in Fig. 4 for Mmm = 1.2 (where the effective charge of a bonded atom is 0.778e 0 ) in the left panel and for Mmm = 1.4 (where the effective charge of a bonded atom is 0.840e 0 ) in the right panel.When Mmm > Maa, there is a gradual transition from molecules back to atoms at high enough pressures, and the pressure at which the system transforms from consisting mainly of dimers to mainly of monomers decreases with increasing Mmm.In addition, there is a large region where monoatomic and molecular species are mixed.As temperature increases, the fraction of atoms in atomic form also increases, similar to the "default" case.
In Fig. 5, we show the variation of the fraction atoms x in molecules with pressure at a temperature of 30 000 K. Increasing Mmm from 1 to 2 (where the effective charge of bonded atoms is 1.00e 0 ) causes the molecules to turn back into atoms at higher pressures.The exact atomic phase fraction of the system at any higher pressure is proportional to the value of Mmm.
In Fig. 6, we show the variation of the fraction of atoms in molecules with temperature at a constant pressure 5 GPa for different values of Mmm.A progressive lowering of the required temperatures for atomization is readily visible for smaller values of Mmm, which switches to an increase in the temperature needed for molecular formation when the atomic phase is stable at low temperatures, for higher values of Mmm.
The influence of Mmm on the Gibbs free energy of the system is shown in Fig. 7.As Mmm increases, the interaction between atoms The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcp that are with molecules become more repulsive, and, as a consequence, we see that the Gibbs free energy of the system increases when there are more molecules in the system.In all cases, the Gibbs free energy is convex, and the filled circles denote the minima of free energy.We see that the value of x at the minimum decreases as Mmm increases; physically, this means that there are fewer molecules in the system as the atoms within molecules become more repulsive.The transition from molecules to monomers with Mmm becomes sharper as the temperature decreases.In Fig. 8, we examine the influence of Mmm on the variation of x and density with pressure.The left plots are for a temperature of 6000 K, and the right plots are for a temperature of 1000 K.At the higher temperature, the fraction of atoms within molecules increases monotonically with pressure for values of Mmm near 1.For larger values of Mmm, there is a pressure at which x reaches a maximum, and at higher pressures, x continuously decreases.As expected, increasing the repulsion between molecules decreases the density of the system.
At the lower temperature of 1000 K, we see a similar behavior of the variation of x with pressure as at 6000 K, although in this case, the decrease in the fraction of molecules with increasing pressure is much sharper.In addition, the pressure at which this transition occurs decreases as Mmm increases.

C. Influence of non-additivity/partial miscibility
We now introduce the final element of realistic physics into our model, turning on the non-additivity of the interactions (i.e., δ ≠ 0).This parameter influences the degree of mixing between the monomer atomic species and dimer molecular species.For a positive non-additivity, the individual species prefer to interact with themselves, rather than with each other.This will produce a tendency for unlike species to demix; when this effect becomes sufficiently large, it can lead to phase separation between the atoms and molecules.Below, we examine the influence of δ > 0 on.For the following analysis, we fix A = 14.1 and B = −4.65 eV.
In Fig. 9(a), we show the variation of the Gibbs free energy with the fraction x of bonded atoms for increasing values of the non-additivity δ, from δ = 0 (purple) to δ = 0.1 (red), for a system with Mmm = 1.2 at T = 3000 K and P = 100 GPa.The filled circles locate the minimum of the Gibbs free energy for each curve.At δ = 0, the Gibbs free energy is convex for all values of x.In fact, for δ ≤ 0, the Gibbs free energy will remain convex, and there will be no demixing.As the value of δ increases, the curvature of the Gibbs free energy gradually decreases until eventually there are regions where the free energy is concave.In this situation, there can be multiple local minimum in the free energy, which implies that the system can possess multiple equilibrium fractions of molecules, which may possibly be metastable.This leads to the potential for the demixing of the monomers and the molecules into two different phases.
In the bottom of Fig. 9, we plot the variation of the Gibbs free energy with x along different isobars for systems with δ = 0.05.At low pressures, the Gibbs free energy is everywhere convex, and there is only a single minimum, which is located near x ≈ 1. Increasing the pressure, the Gibbs free energy develops a concave region for intermediate values of x, and there are two minima: one is stable, near x ≈ 1, and the other is metastable, near x ≈ 0, and a local minimum.At higher pressures, the Gibbs free energy of the minimum near x ≈ 0 decreases until it becomes equal to that near x ≈ 1.At this point, the system splits into two phases, with one rich in monomers and the other rich in dimers.Increasing the pressure further, the Gibbs free energy of the minimum near x ≈ 0 becomes lower than near x ≈ 1, and, so, the system transitions to a single phase predominantly composed of monomers.
In Fig. 11, we examine the influence of the parameter δ on the pressure dependence of the fraction of atoms within molecules and of the density.In these calculations, the temperature of the system is 6000 K, Mmm = 1.12 (where the effective charge of bonded atoms is 0.751e 0 ), and the parameter δ is varied from 0.0 to 0.05.At extremely low pressures, the atoms prefer to be in the form of monomers; as the pressure increases, the atoms increasingly form dimeric molecules.The density continuously increases with increasing pressure, suggesting that no phase transition occurs.In the low pressure region, the non-additivity parameter has little effect on the properties of the system, as the non-bonded atoms do not interact significantly with each other.At higher pressures, the fraction of atoms within molecules increases with pressure, eventually reaching a maximum, and then decreasing afterward.Non-additivity has a relatively significant effect at these higher pressures; the larger the value of δ, the greater the tendency for the atoms to form diatomic molecules.For small values of δ, the decrease in x with increasing pressure is gradual and continuous; however, for larger values, there is a discontinuous drop in x with pressure, which becomes wider as δ increases.This discontinuity reflects the formation of two coexisting phases: one rich in monomers and the other rich in dimers.
In Fig. 10(a), we plot the fraction of atoms within a molecule as a function of pressure, along several isotherms.At high temperatures, these curves are continuous, with the atoms existing mainly as molecules at low pressures and gradually transitioning to predominantly monomers at high pressures, due mainly to the atoms in the molecules being more repulsive (i.e., Mmm > Maa).At low temperatures, we see the same general trend with pressure, but now there is a discontinuous drop of in x.At this discontinuity, there is a coexistence between two phases: one rich in monomers and the other rich in diatomic molecules.The black dashed line denotes the boundaries of these discontinuities with pressure, and the gray shaded region represents the state points in which the system will separate into these two coexisting phases.The variation magnitude of the composition difference between the phases with temperature and pressure is shown in Figs.10(b) and 10(c), respectively.
The density of the system increases monotonically with pressure, as shown in the bottom plot of Fig. 11.The non-additivity parameter also has a greater influence on the density of the system at high pressures.Interestingly, at lower pressures, increasing the value of δ will decrease the density of the system, while at higher pressures, increasing δ will increase the density (see the inset in the bottom of Fig. 11).The discontinuity in x with pressure at sufficiently large values of δ is associated with a discontinuity in the system density.
In the case of a discontinuous molecular-to-atomic phase transition, from 100% atoms bound in molecules to 100% free atoms, the , where f = 2πε/κ 3 and p is the pressure of the transition at the given temperature.We can see from here that the density jump arising from a discontinuous molecular to atomic transition is proportional to the value of Mmm.
It is readily visible that non-additivity introduced a critical point in the system, between a discontinuous molecular-atomic transition and a continuous one.The former is evidenced by an associated discontinuity in the fraction of dimers, is also mirrored in the density of the system either side of the transition.
The critical point, which is characterized by the critical temperature Tc, pressure p c , and composition xc, is given by the solution of the following equation: With the chosen model parameters Mmm = 1.12 and δ = 0.01, we find that the critical point is located at T ≈ 2655 K, P ≈ 235 GPa, and xc ≈ 0.3944.The Gibbs free energy of the system at a temperature below (T = 1000 K) and above (T = 4000 K) the critical temperature for a selection of pressures is presented in Fig. 12.For T < Tc, the Gibbs free energy has a convex region, and each curve has two local minima: one near x = 0 and the other near x − 1, while a single minimum that moves with increasing pressure is visible for T > Tc.
In Fig. 13, we examine the influence of the non-additivity parameter δ on the discontinuous transition between the molecular and monomer phases.The solid lines denote the coexistence curve for FIG.12. Gibbs free energy of a system with Maa = 1.0,Mmm = 1.12, and δ = 0.01 at temperatures below (T = 1000 K, upper panel) and above (T = 4000 K, lower panel) the critical point.The filled circles denote the minimum of the Gibbs free energy.
The Journal of Chemical Physics ARTICLE pubs.aip.org/aip/jcphydrogen, while the filled circles locate the critical point of the transition.At lower pressures, the system consists mainly of diatomic molecules, while at higher pressures, it consists mainly of unbonded, monomeric atoms.Upon crossing the coexistence curve, there is a discontinuous transition between these two phases.As the value of δ increases, the critical point moves to higher temperatures and lower pressures.
To examine the isotope effect, we also plot the coexistence curve and critical point of deuterium as the dashed lines and crosses, respectively, in Fig. 13.The ε and κ parameters for deuterium are the same as for hydrogen, but the bond free energy differs (A = 14.3 and B = −4.71eV).
The critical temperatures of both species are similar.At lower temperatures, the monomer to molecular transition for deuterium is shifted to higher pressures, as compared to hydrogen; however, the coexistence curves approach each other at higher temperatures.Interestingly, for higher values of the non-additivity parameter δ, the coexistence curves cross, and the critical pressure of deuterium becomes lower than that of hydrogen.
In Fig. 14, we show the fraction of atoms bound in diatomic molecules for systems with different values of δ.At low temperatures and low pressures, the system tends to prefer molecules.As temperature increases, the system favors monomers, which have a higher entropy.At high pressures, monomers are again favored, as they have less repulsive interactions, due to the fact that we have Mmm > Maa.For the case δ = 0, the transition from mainly monomers to mainly molecules is continuous.

D. Application to hydrogen and deuterium
Since the stable phase is determined by minimizing the Gibbs free energy, all phase diagrams are based on the competition between internal energy U, entropy S, and density 1/v.Here, we consider the competition between molecules and atoms.At zero temperature, we need only consider the enthalpy.In the present case, the atomic form always has a higher entropy simply because of having more particles.Therefore, the atomic phase is always stable at high T.
An important mention is warranted at this point: the density jump in our model when crossing the discontinuous molecular-atomic transition at 1000 K is 4.92%, comparable to the one noted by Geng et al. 41 from a full quantum mechanical treatment of the system.
Hydrogen is also famous for displaying very strong isotopic effects, with the phase diagram of high pressure deuterium being noticeably skewed in comparison with H 2 .Keeping all the same interaction and additivity parameters for both species but allowing different empirical A and B parameters for ΔG bond , we obtain a curve for deuterium that mirrors that of hydrogen but has phase changes shifted to higher pressures at any given temperature.This situation reproduces the situation expected based on experimental investigations of solid phases, showing that close to ambient temperature the pressures necessary for the atomization of D 2 are higher than for H 2 .For example, at 1000 K, there is a 6 GPa difference in the pressure where the molecular-to-atomic transition occurs, at 290 GPa in H 2 and at 296 GPa in D 2 .This number is of course a likely underestimate of the pressure difference between H 2 and D 2 dissociation as the parameters ε, κ −1 , and most importantly Mmm should be different for the different isotopes.In Fig. 15, the quantitative agreement between our model (having the previously stated choice of parameters) and the most recent experimental data from shock compression of D 2 by Knudson et al. 63 and Celliers et al. 64 is readily visible.It is of course possible, by tuning the parameters detailed in the previous paragraph, to obtain quantitative agreement for H 2 as well, maintaining the same qualitative picture of the full phase diagram.However, given the main intention of the current work is to obtain a qualitatively accurate and order-of-magnitude only estimate of the phase behavior of diatomic systems, such fine-tuning was not attempted.Moreover, while experimental data for deuterium are very recent and shows good agreement between two distinct research groups and experimental machinery (Sandia Z-machine for Knudson and National Ignition Facility for Celliers, respectively), the reported data for H 2 are significantly more varied in the reported pressure and temperature conditions, and there is considerable discussion about the criteria employed for identifying the atomic phase, as discussed in detail in Norman and Saitov. 66

V. CONCLUSIONS
To summarize, we have built a thermodynamic model for a dimerizing fluid inspired by the behavior of hydrogen.The dimerization leads to dimers being favored at low temperatures thanks to their finite binding energy, while the atoms are preferred at high T due to their higher entropy.In the absence of intermolecular attraction, no condensed phase can form at zero pressure, so the atomic region extends to low temperature at extremely low densities.The primacy of the entropy can be easily understood since the entropy increases with ln V, while the dimer bonding is independent of V.
We then consider repulsive interactions between all objects, giving a density-dependent contribution to the free energy: this is based on a Yukawa model with scaling factors for atom-atom, atom-molecule, and molecule-molecule interactions.These scaling factors make the interaction both stronger and longer ranged.If the atom-atom repulsion is shorter ranged than the molecule-molecule, then the atomic phase is denser and favored at high pressures.This creates a situation where we have a dome of stability of molecules, stabilized by low energy from the bonds.The atomic phase is still stabilized at low densities by entropy and now also at high densities by the PV term.This creates a dome of stability for the molecular liquid.
In order to obtain a first-order atomic-molecular transition, it is necessary to have concavity in G as a function of x.Within our model, the non-additivity parameter δ > 0 can provide this, by making atom-molecular interactions relatively unfavorable.The first order transition ends at a critical point.δ appears in the interaction term a/v, so has most effect at high densities, and as we increase δ, we see that the critical point moves up the dome between molecular and atomic fluids.This produces a negative-sloped Clapeyron curve at high pressures.Surprisingly, as δ is further increased, the critical point moves to ever lower density and eventually crosses the top of the dome, beyond which the critical temperature drops.This leads to a positive slope in the Clapeyron curve.
With respect to hydrogen, the model implies that the high pressure metallic fluid and the low density interstellar medium are one and the same phase.One might expect a Mott transition when electrons become localized on the atoms; however, this is outside the scope of our model, which ignores the electrons-it also falls out with Kohn-Sham DFT, which always treats electrons as delocalized Bloch states.
The liquid-gas line within the molecular region is also missing in our model.It can easily be reproduced by using a van der Waals gas in place of the ideal gas.However, the van der Waals equation becomes incompressible at high densities, and if one fits the parameters to the hydrogen critical point, this unattainable density falls well below the metallization line.This transformation is far from the atomic-molecular dome, so our conclusions are unaffected by this failing.
Finally, there is a difference in ΔG bond between deuterium and hydrogen due to the zero point energy.This is easily incorporated into the model to show a pronounced isotope effect on the atomic-molecular transition, in accordance with experiment and path integral molecular dynamics calculations. 7,8THOR DECLARATIONS FIG. 1. Free energy of bond formation for hydrogen (black solid line), deuterium (black dashed line), and oxygen (red solid line).The intercept shows the strength of the chemical bond, and the slope is the entropy.The inset shows the enthalpy of bond formation ΔH bond .Data are from Ref. 65.

FIG. 5 .FIG. 4 .FIG. 6 .
FIG. 5. Fraction of atoms in molecules as a function of pressure along the isotherm T = 30 000 K for δ = 0.0 and various values of Mmm.

FIG. 7 .
FIG. 7. Gibbs free energy for a system with Maa = 1 and P = 100 GPa, considering different values of Mmm at temperature (top) T = 3000 K and (bottom) T = 4 K.The tendency for increasing values of Mmm to stabilize the atomic phase over the molecular one is readily visible under both conditions.

FIG. 8 .FIG. 9 .
FIG. 8. (left) Effect of Mmm on the formation of molecules and the system density with increasing pressure in the range 0-1 GPa at T = 6000 K.The changes at both low and high pressures are smooth and continuous and do not involve any discontinuous jumps in density.(right) Effect of Mmm on the breaking of molecules into atoms and the density with increasing pressure at T = 1000 K.

FIG. 10 .
FIG. 10.Behavior of a system with Mmm = 1.12 and δ = 0.01 around the critical point.(a) Variation of the fraction of atoms within molecules with pressure, (b) discontinuity in x with temperature along the coexistence line, and (c) discontinuity in x with pressure along the coexistence line.

FIG. 13 .
FIG. 13.Influence of the non-additivity parameter δ on the discontinuous transition between molecular-rich and monomer-rich phases.The solid lines denote the coexistence curve for hydrogen (A = 14.1 and B = −4.65 eV), and the filled circles mark the critical point.The dashed lines denote the coexistence curve for deuterium (A = 14.3 and B = −4.71eV), and the crosses mark the critical points.

TABLE I .
Parameters for the linear fit to βΔG bond [see Eq. (22)] for H 2 and D 2 over different temperature ranges.