Scaling laws are often used to determine the critical temperature T c of a system [1][2][3][4][5]. For instance, T c can be obtained by fitting an Ising scaling law for an order parameter defined as the difference between the liquid and vapor densities at coexistence (ρ l -ρ v ). This approach has been used extensively and has been shown to hold for a wide range of atomic and molecular systems. This has proved to be extremely useful in the case of metals since critical temperatures are generally extremely high [6], require advanced experimental techniques [7] and often exhibit large uncertainties with for instance, for Al, estimates ranging from around T c = 5500 K to T c = 9600 K [8]. On the other hand, the determination of other critical properties, their determination often relies on the use of empirical laws such as, e.g., the law of rectilinear diameter [9,10] for the critical density. This law, which assumes a linear behavior for diameter ρ m of the liquid-vapor coexistence curve as a function of temperature, with ρ m = 0.5(ρ l + ρ v ). However, several systems [11][12][13] exhibit strong departures from this law. For instance, for alkali metals like Cs and Rb, experiments showed that the two branches of the coexistence curves were strongly asymmetric and the law of rectilinear diameter was found to break down over a large temperature range [14]. As a result, this law breaks down for metals [6,14,15] as a result of many-body effects [16]. Alternative methods have been developed in recent years to find ρ c via extrapolation, through power series law for the diameter [6] or novel symmetrized equations for the vapor-liquid coexistence curve [17][18][19].
In recent years, entropy has become increasingly key to our understanding of phase transitions, self-assembly and of the pathways underlying such processes. This has sparked a number of studies aiming at the determination of entropy and at a deeper understanding of its behavior along phase transition curves. Scaling theory [20] provides an equation for the variation of the entropy density S = s/V along the vapor-liquid coexistence curve as the system approaches criticality. Scaling theory predicts that the entropy density of vaporization ∆S, as well as the density change upon the phase transition ∆ρ = (ρ l -ρ v ), vary with the temperature T as (T c -T ) β , in which β = 0.326 is the corresponding 3D-Ising critical exponent. Here we focus on entropy, rather than on entropy density, and identify how entropy varies along the coexistence curve, how data for the entropies of the two coexisting phases can be used to estimate the critical entropy and on whether entropy exhibits the same behavior for ideal systems and metallic systems. To address these questions, we carry out flat histogram simulations to determine the thermodynamic locus for coexistence, as well as the entropy for the two coexisting phases. We apply this approach to model systems, using an effective potential that encompasses both ideal, Lennard-Jones-like, behavior and metallic behavior, and to real systems, modeled with a many-body potential. More specifically, for the effective potential, the onset of the first Friedel oscillation [21][22][23] is controlled through a switching parameter. This allows us to assess systematically the impact of the increase in metallic character in the system on entropy and its behavior along the vapor-liquid transition. For real systems, we use the quantum corrected Sutton Chen Embedded Atoms Model (qSC-EAM) and analyze the impact of the many-body interactions on entropy.
The paper is organized as follows. We first present the models and simulation methods used in this work. In particular, we discuss how our simulation approaches allow for the determination of the conditions for coexistence and of the thermodynamic properties along the vapor-liquid equilibria, including the entropy of the two coexisting phases. We then analyze the symmetry of the coexistence curve in the T -S and P -S planes and characterize the impact of an increasing metallic character, focusing first on model systems and then on the realistic models for metals. We finally draw the main conclusion of this work in the last section.
To understand the impact of the metallic nature of the system, we use two different types of force fields. First, we model the interactions between metal atoms via an effective pair potential that account for the onset of the first Friedel oscillation that arises in metals.
We achieve this through the combination of a Lennard-Jones (LJ) functional form and of a Dzugutov (DZ) pair potential [22,24], that has been shown to favor the formation of quasicrystalline [25][26][27][28][29] and crystalline σ-phases [30]. While previous work has focused on the determination of phase boundaries for the liquid-solid [22,23] and vapor-liquid equilibria [31], there has not been, to our knowledge, any study on the behavior of entropy along the phase envelope. The resulting effective potential is given by
with
and φ DZ (r) = φ 1 (r) + φ 2 (r)
in which X is a switching parameter (0 ≤ X ≤ 1), m = 16, A = 5.82, C = 1.1, a = 1.87, B = 1.28, d = 0.27 and b = 1.94. The switching parameter controls the transition from an ideal-like system (X = 1) to a metallic system (X < 1), and provide insight into the transition from an ideal, corresponding-states, behavior to a metallic behavior.
Second, to take into account many-body effects, we use a many-body force field, known as the quantum-corrected Sutton-Chen embedded atoms model (qSC-EAM), to model the interactions between metal atoms. In this case, the potential energy for a system of N metal atoms is equal to
in which r ij is the distance between two atoms i and j and the density term ρ i is given by
We give in Table I the qSC-EAM potential parameters for the metals considered in this work. As in previous work [32][33][34][35], the cutoff distance to evaluate the interactions is set to twice the parameter a. The qSC-EAM force field has been shown to perform well for a broad range of thermodynamic [32,36] and transport properties [37,38] of metals. For the solid phase, the qSC-EAM potential has been shown to provide results in very good agreement with the experimental data for elastic constants, cohesive energy and surface energy [39], as well as for the melting points of pure metals [32] and alloys [40]. In the liquid phase, this model performs very well for the density and viscosity of liquid metals [37,38,41]. Finally, the qSC-EAM force field has also been shown to predict accurately the vapor pressure and critical properties [42] and the boiling points for transition metals [43,44]. To identify the conditions for coexistence and the properties of the coexisting phases, we perform flat-histogram sampling simulation approach known as Expanded Wang-Landau (EWL) simulations [45][46][47]. We briefly describe the simulation method here (more details may be found in refs. [45][46][47]). EWL simulations rely on the Wang-Landau (WL) scheme to sample extensively all possible configurations of the system [48][49][50][51][52][53][54]. The EWL approach is Monte Carlo (MC) method implemented within the grand-canonical (µ, V, T ) ensemble, meaning that, in this case, the WL sampling is applied to sample evenly all possible number of atoms and determine Q(N, V, T ) for all N values and the grand-canonical partition function of the system as
As shown in previous work [45], an accurate determination of the partition function hinges on high acceptance rates of the MC steps leading to changes in the number of atoms. Since, for instance, the random insertion of additional metal atoms in dense liquid phases is often associated with a low acceptance rate, we combine the WL sampling with an expanded ensemble approach [55][56][57][58][59][60][61][62]. Thus, we split the insertion (or deletion) of entire atoms into M stages and gradually grow (shrink) the extra atom to be inserted (deleted). In other words, at all times, the system contains N atoms and a fractional particle at stage l (with 0 ≤ l < M ). Throughout the EWL simulations, histograms for Q(N, V, T, l) are continuously updated every time a configuration with a given set of (N, l) is visited. When the simulation has converged, we gather the canonical partition functions obtained for all systems with a void fractional particle (l = 0) and calculate the grand-canonical partition function through Eq. 6. Once the partition functions have been determined, all thermodynamic properties, including the entropy, can be determined for a wide range of conditions.
We perform EWL simulations for systems modeled with the effective pair potential (LJ+DZ) and for Ag, Pd, Ni and Ir using the qSC-EAM many-body force field. More specifically, we implement the EWL approach within a MC framework, with 75% of attempted moves allocated to random translations of a randomly chosen particle (either an atom or a fractional particle) and the remaining 25% of attempted moves allocated to changes in (N, l)
values. The number of stages is set to M = 100, the initial value for the convergence factor f is set to e and the minimum number of visits for each (N, l) set is fixed to 1000. Once all possible (N, l) sets have been visited at least 1000 times, the convergence factor is set to f → √ f and we initialize the histogram for the number of visits to 0. This is repeated until f < 10 -8 . Throughout the simulations, the histograms for the canonical partition functions are collected allowing for the determination of the thermodynamic properties for the system.
Throughout the manuscript, the results are given in units reduced with respect to the LJ parameters, i.e. using σ as the unit of length, as the unit of energy and by setting the De Broglie wavelength Λ = 1 for the property calculations, for the (LJ+DZ) potential and in real units for Ag, Pd, Ni and Ir. In addition, results obtained for the entropy at coexistence for Cu using the grand-isobaric adiabatic ensemble [63] are also presented. The interactions are calculated using a spherical cutoff (r c = 3σ) for the (LJ+DZ) potential, with the usual tail corrections applied beyond the cutoff distance [64], and with a spherical cutoff of 2a for the qSC-EAM potential.
We first examine the behavior of entropy along the vapor-liquid coexistence curve in an ideal system using, as a reference, the available experimental data on Argon [65]. Starting from the results of scaling theory [20], the density of the two coexisting phases ρ ± is predicted to behave as
in which ρ c is the critical density, t = 1 -T /T c , β is a critical exponent equal to 0.326 and α, θ, θ 5 , A 2β , A 1-α , b θ and b 2β are additional critical exponents and parameters (see Kim et al. for more details [20]). Turning to the difference between the two densities at coexistence ∆ρ = ρ + -ρ -, most terms cancel out and the following result is obtained
The top panel in Fig. 1 confirms that the experimental data on Argon can be fitted accurately by the functional form ∆ρ = 2ρ c B|t| β with β = 0.326, showing that the series of terms between brackets in Eq. 8 can be estimated to be equal to 1 and leading to a value of 2ρ c B = 1800.55 kg/m 3 .
Similarly, for the entropy density S = s/V of the two coexisting phases, scaling theory gives
This leads to a difference between the two entropy densities at coexistence ∆S = S + -S -
The middle panel in Fig. 1 indicates that ∆S = 2ρ c k B B S |t| β models accurately the behavior of ∆S. The fit to the experimental data also shows that the series of terms between brackets is close to 1 and provides a value for the parameter 2ρ c k B B S = 3227.05 kJ/m 3 /K.
We now focus on the behavior of the entropy of vaporization ∆(S = s/N ) = ∆(S/ρ), rather than of the entropy density ∆S. Scaling theory provides a path towards its determination, provided that a few approximations are made. Keeping only the (critical) constant term and the leading order term in β in the equations for both the density and the entropy density (Eqs. 7 and 9), one can obtain
and, hence, to the leading order, to 1/ρ ± = (1/ρ c ){1 ∓ B|t| β }. In turn, this gives the following equation for the entropy of vaporization
The critical entropy density is given by S c = ρ c k B k 0 , which leads to
Using the parameters B and B s obtained from the two previous fits, as well as the experimental data for Argon ρ c = 536 kg/m 3 and S c = 2.261 kJ/kg/K, Eq. 13 provides a relation to the first order in β for the dependence of ∆S. As shown in the bottom panel of Fig. 1, the resulting equation does not perform nearly as well as the previous two equations for the density and entropy density difference over the same temperature interval. While power laws with an exponent of 0.326 model well both the density difference and entropy density difference under these conditions, it turns out that the experimental data on Argon for ∆S is more accurately modeled with a power law of exponent of 0.47.
We now turn to the systems simulated in this work and discuss on the example of the X = 0.9 (LJ+DZ) system how we determine the entropy along coexistence from the simulation results. The plots are given in Figs. 2 and3. The canonical partition functions collected over the EWL simulations are shown for T = 1 (inset of Fig. 2(a)) and exhibit the expected increase with N for subcritical conditions [45]. Applying Eq. 6 leads to the calculation of the grand-canonical partition function across the range of chemical potentials. Most notably, we see in Fig. 2(a) a sharp increase in Θ(µ, V, T ) for chemical potentials greater than µ = -3.2, corresponding to the transition from the vapor (low µ) to the liquid phase (high µ). To determine accurately the conditions of coexistence, we evaluate the number distribution p(N ) as
and obtain numerically the chemical potential at coexistence, where the two phases are equally probable, by solving
in which N b corresponds to the minimum for p(N ), and the two sides of the equation are the total probability for the vapor (left) and for the liquid phase (right). The plot for p(N ) at coexistence is shown in Fig. 2(b), and exhibits two peaks of equal area corresponding to two different average values for < N >, thus establishing that we have correctly identified the conditions for coexistence. This means that we now have access to µ coex at coexistence and to the average number of particles in the two coexisting phases, < N vap > and < N liq >, or equivalently the number densities for the two coexisting phases < ρ liq > and < ρ vap >. Next, to determine the energy of the two phases, we collect, during the EWL simulations, an histogram for the potential energy E(N ) as a function of the number of particles N (shown in the bottom panel of Fig. 2(b) for X = 0.9 and T = 1). This, in turn, provides access to the average potential energy for the two coexisting phases, noted as < E liq > and < E vap >, which will be key for the determination of the entropies of coexistence as discussed below.
We now turn to the analysis of the fluid properties as the conditions approach criticality. To compare the results for the coexisting entropies along the phase boundary, we need to be able to obtain scaled entropy-temperature and entropy-pressure plots for all systems, with a scaling performed with respect to the critical parameters for each system. For this purpose, we first focus on determining the critical properties (T c , P c , S c ) Results are first shown for the difference between the densities of the two coexisting phases as a function of temperature in Fig. 3(a) for X = 0.9. A power-law with an Ising-type exponent of 0.326 fits very well the simulation results, providing a reliable estimate for the critical temperature T c (critical properties are listed for all systems studied in this work in Table II). Moving on to the dependence of the pressure at coexistence upon temperature, we calculate P coex from the usual statistical mechanical relation We report in Fig. 3(b) the results obtained for P coex when X = 0.9. We fit the simulation results with the following exponential function P = A exp B T +C (A, B and C are fitting parameters). Fig. 3(b) shows that the simulation results for the pressure at coexistence are very well accounted for by this functional form. This fit, also known as Antoine's law, allows to determine the dependence of P coex upon temperature and to obtain an estimate for the critical pressure P c by plugging in the value previously obtained for T c (P c estimates are gathered in Table II for all systems).
The coexisting entropies are calculated from the EWL simulation results as
FIG. 4: Entropy for the two phases at coexistence for X ranging from 0.7 to 1 (LJ system), X = 0.9. In (a), the entropy scaled with respect to the critical entropy S c is plotted against the scaled pressure (with respect to P c , while (b) shows a scaled entropy-temperature plot.
in which the average energy for each phase is calculated from E(N ), shown in Fig. 2(b), weighted by the number distribution p(N ) over each of the two peaks. This gives, e.g., for the vapor phase, < E vap >= N b N =0 p(N )E(N ). We finally estimate the critical entropy S c by averaging the entropies of the two coexisting phases at the highest temperature available from the simulations (see Table II for numerical values). This is motivated by the observation that the average entropy exhibits little dependence upon temperature close to the critical point. For instance, when X = 0.9, it varies by less than 5% for T ≥ 1. We show in Fig. 5 the entropies at coexistence along the phase coexistence for X ranging from 0.7 to 1 (LJ system), using as a scaling factor the value of S c for each system. The results are plotted in the scaled entropy-pressure plot in Fig. 5(a) and in the scaled entropytemperature plot in Fig. 5(b). Despite strong changes in the effective pair potential and the gradual onset of the first Friedel oscillaton as X decreases, we find that the behavior of the scaled entropy for the two coexisting phases remains remarkably the same for all systems.
Indeed, the simulation results almost fall onto the same plot for all systems over the range of temperatures studied here. Our results also show that the same conclusion applies to both plots, i.e. when the coexisting entropies are plotted either against pressure, or against temperature.
The next step consists in applying the same analysis to systems modeled with a manybody force field. We start by determining the critical properties from the EWL simulation results and give in Table II the numerical values we use to scale the temperature, pressure and entropy. We show in Fig. 5(a) the resulting scaled entropy-pressure plot and in Fig. 5(b) the scaled entropy-temperature plot for a series of transition metals (Ag, Ir, Ni and Pd).
The plots show that the four metals exhibit very similar behaviors, with the entropies for the two coexisting phases aligning well for the various systems. This conclusion holds over the range covered by the different potential parameters, but also for the different values taken by the exponents used in the repulsive two-body part (see Eq. 4) and in the attractive many-body part (see Eq. 5). This means that we will be able to analyze and characterize as a group, i.e. the many-body group, the variations of entropy for the two coexisting phases along the phase boundary.
Turning to the entropy difference between the two phases, we compare in Fig. 6 the results obtained for the effective pair potential to those found for the many-body potential for all four metals. Fig. 6 shows the emergence of two distinct groups, which establishes the qualitatively different behavior exhibited by both classes of models. Perhaps surprisingly, the onset of the first Friedel oscillation does not significantly impact how the entropy difference varies as a function of the scaled pressure. Indeed, the behavior found for all X values studied here remains remarkably similar to that found for the ideal (LJ) system, simulated for X = 1. On the other hand, the qSC-EAM many-body potential yields results that are markedly different. Another interesting finding from this plot is the similarity found between entropy differences plot across the range of metals considered here. This points towards the predominant role played by many-body effects in such systems.
The dependence of ∆S and ∆S is shown on an example for the two classes of systems studied here, i.e. X = 1 and Cu for the qSC-EAM potential (See Fig. 7). The results for the entropy density difference at coexistence show that, for both classes of systems, a power law with the expected exponent of 0.326 models well the results along the coexistence curve.
However, when turning to the entropy difference at coexistence, two very different behaviors are observed, with the many-body potential for Cu leading to a power law fit with a higher exponent (0.67) than for X = 1 (0.47).s To characterize this further, we examine the dependence of the difference between the scaled entropies of the two coexisting phases on the scaled temperature. Fig. 8 shows the onset of two different groups, the effective pair potential group, which conforms to the ideal behavior as shown with the LJ system (X = 1), and the many-body potential group. We then set on to extract from the simulation data a scaling law for the entropy difference. For this purpose, we fit the two groups of results with the following power-law
in which A and s are two fitting parameters. The resulting fits are also plotted in Fig. 7 and shows that this functional form performs well for both groups. The value found for the exponent s takes distinct values for the two groups. Specifically, s is found to be equal to be close to 0.5 for the ideal behavior, with a value of 0.495 obtained for the effective s pair potential group. On the other hand, the exponent of 0.67, identified above for for Cu, holds over the entire group modeled with the qSC-EAM many-body-potential. This result departs notably from the ideal critical exponent of 0.5 found for the ideal group. As shown in Table I, the parameters for the qSC-EAM model undergo significant changes from one metal to another, most notably for c that varies between ∼ 84 -225 and n that varies between 10 -13 in the cases studied here. However, we observe a remarkable agreement between the exponents obtained for the metals considered in this work and attribute this finding to the general functional form of the qSC-EAM force field. Future work will focus on extending this conclusion to a broader range of metals and models.
In this work, we characterize the behavior of entropy along the vapor-liquid phase coexistence and identify a scaling equation for the difference in entropy between the two coexisting phases. For this purpose, we carry out series of flat-histogram simulations, based on the Expanded Wang-Landau method, to determine the partition function, coexistence loci, critical parameters and entropies along the binodal. From the simulation results, we are able to obtain two seldom explored projections of the phase boundary in the S -P and S -T planes and to determine estimates for the critical entropy. Then, we characterize the impact of the onset of the first Friedel oscillation in a series of effective pair potentials and of many-body effects in a series of transition metals, modeled with an embedded-atoms force field. We find two classes of behaviors. First, the effective pair potential yields results that, when shown on scaled entropy-pressure plots (using the critical entropy S c and the critical pressure P c as scaling factors) and on scaled entropy-temperature plots (using the critical temperature T c as a scaling factor in this case), exhibit an ideal behavior similar to the Lennard-Jones system. The results also show that the onset of the first Friedel oscillation in the effective pair potential does not result in any qualitative change, leading to a critical exponent of 0.5 for the difference in entropy between the two phases at coexistence as a function of temperature. Second, the many-body potential results in an entropy behavior that is qualitatively the same for all transition metals. However, in this case, the presence of a many-body term in the interactions leads to markedly different coexistence curve in the S -P and S -T planes, when compared to ideal systems, as shown by the critical exponent of about 0.83 for the difference in coexisting entropies as a function of temperature. Future work will focus on elucidating the relation between the nature of the interatomic and intermolecular interactions and the scaling of entropy as conditions approach criticality.