Evaluation of the grand-canonical partition function using expanded Wang-Landau simulations. V. Impact of an electric field on the thermodynamic properties and ideality contours of water

The properties of water subjected to an electric field are of key importance for many applications, ranging from nanofluidics,^1–4 biochemistry,^5,6 to chemical processes.^7,8 This is, for instance, the case in living organisms, where electric fields control the orientation of the water molecules and thus the proton transfer direction driving ATP synthesis.⁶ Similarly, in atmospheric sciences, ions are well-known to have a direct impact on the nucleation process by, depending upon the conditions, promoting or preventing the transition from vapor to liquid through the nucleation of droplets.^9–11 Alternatively, an electric field can also be used as a great tool to control the behavior, trigger phase transitions, and tune the properties of liquids and materials.^12–15 

However, the response of a polar liquid subjected to very strong fields still remains poorly understood. A recent work shows that the critical temperature of water subjected to very large electric fields is lower than that of the bulk,^16,17 while others indicate either the opposite behavior^18–21 or, in the case of nanoconfined water, a complex dependence of the critical temperature on the extent of the nanoconfinement and the field strength.²² There are also few studies available on the impact of the electric field on the Gibbs free energy and entropy of water. Furthermore, in the recent years, the elucidation of the thermodynamics of supercritical water has drawn considerable interest in geochemistry²³ and as a green solvent for the catalytic conversion of biomass into fuel,²⁴ and it is unknown how strong fields change the thermodynamics and structure of water in the supercritical region of the phase diagram. The goal of this work is twofold: (i) we extend the recently developed Expanded Wang-Landau (EWL) simulation method^25–28 to determine the thermodynamic properties of systems subjected to an external field and (ii) we apply the resulting method to analyze the impact of the electric field on the properties of water over a wide range of conditions, i.e., at the vapor-liquid phase boundary, for compressed liquids as well as under supercritical conditions. In particular, in the supercritical region of the phase diagram, we focus on analyzing the effect of the field on the ideality contours,^29–35 known as the Zeno line and the curve of ideal enthalpy, to establish a correspondence between the results obtained for different fields. These ideality contours have recently emerged as a new way to bridge the gap in our understanding of supercritical fluids^36–39 and have paved the way for the development of new similarity laws³⁰ and maps of the supercritical region of the phase diagram.^40,41

This paper is organized as follows. In Sec. II, we discuss how we extend the EWL method to a polar fluid subjected to an electric field. We then detail the simulation model as well as the technical details before presenting the results obtained in this work. We start by discussing how the electric field impacts the grand-canonical partition function for water, and from there, we examine the effect of the field on the thermodynamics of phase transition, on the location of the critical point and on the response of compressed liquids of water and supercritical water to the field. We particularly focus on the interplay between the changes in density and structure and the thermodynamic properties of water under an electric field. We finally discuss the impact of the field on the ideality contours of water before drawing the main conclusions of this work in Sec. IV.

The goal of the recently developed Expanded Wang-Landau (EWL) simulation method consists in determining a high-accuracy estimate for the grand-canonical partition function, which, in turn, yields all thermodynamic properties of the system through the statistical mechanics formalism.⁴² The first papers of the series have discussed developments of the method as applied to single-component systems^25,26 and mixtures²⁷ modeled with classical and quantum (tight-binding)²⁸ force fields. Here we extend this approach to the case of systems subjected to an external field and consider the case of water subjected to an electric field. The grand-canonical partition function⁴² is given by

where $β = 1 / k B T$⁠, N the number of molecules, $μ$ the chemical potential, and Q(N, V, T) is the canonical partition function written as

where $q t r a n s = [ 2 π M k B T h 2 ] 3 / 2 V = V Λ 3$ is the translational partition function for an individual water molecule of mass M, q_int is the intramolecular partition function, and $𝚪$ denotes a specific configuration of the system. The equation for q_int depends on the assumptions underlying the force field used for water as discussed in Sec. II B.

In Eq. (2), $U ( 𝚪 )$ denotes the potential energy of the system for a given configuration. U is calculated as the sum of the interaction energy between water molecules and of the energy resulting from the interaction with the electric field. For a system of N molecules, the energy due to the field is obtained through

where m_i is the dipole moment of water molecule i and E is the applied electric field.

EWL simulations take advantage of an efficient scheme for the insertion/deletion of molecules, known as the expanded ensemble approach,^43–53 which splits the insertion and deletion of a full molecule into M stages. The implementation of efficient schemes for the insertions/deletion steps, e.g., in expanded ensemble-transition matrix Monte Carlo methods⁵⁴ or with the continuous fraction component methods,^55,56 has been shown to yield very accurate results. Here, the combination of a Wang-Landau sampling with the expanded ensemble approach ensures an efficient sampling of all possible N values and results in highly accurate predictions for the thermodynamic properties in the low temperature high density regime.^25,57 In the EWL method, the simulated system is composed of N full molecules and of a fractional molecule at stage l (with $0 ≤ l ≤ M − 1$⁠). The insertion/deletions steps are handled through changes in the value of l. If, during the simulation, l is increased beyond M, the fractional molecule becomes a full molecule and a new fractional molecule at stage l $−$ M is created. This results in the insertion of a new full molecule as the system now contains N + 1 full molecule and a new fractional atom at stage l $−$ M. The deletion of a molecule is similarly achieved though a decrease of the value of l for the fractional molecule. Finally, when l = 0, the fractional molecule is void and the system contains N full molecules. The resulting simplified expanded grand-canonical ensemble (SEGC) partition function²⁵ for this system is

in which Q(N, V, T, l) is the canonical partition function for a system of N full atoms and a fractional atom at stage $l > 0$⁠, given by

Here, the mass, as well as the moments of inertia, for the fractional molecule is chosen to be the same as that of a full molecule, leading to identical translational and rotational partition functions, q_l,trans = q_trans and q_l,rot = q_rot.

The Q(N, V, T, l) functions are then determined numerically during the EWL simulations, through the iterative evaluation of the biased distribution $p b i a s ( 𝚪 , N , l )$⁠.^25–28 For Wang-Landau simulations,^57–69 the Metropolis criterion for a move from an old state (⁠$𝚪 𝐨 , N o , l o$⁠) to a new state (⁠$𝚪 𝐧 , N n , l n$⁠) is given by

Following the derivation for p_bias for the EWL simulations,^25–28 we obtain for the case of water the following equation for p_bias:

This leads to the numerical determination of p_bias and thus of Q(N, V, T, l). A specific advantage of carrying out a Wang-Landau sampling in the grand-canonical ensemble is that the variable sampled, i.e., N the number of molecules, is a discrete function, which circumvents the discretization of the energy range that would be required in other ensembles,⁷⁰ and allows the simulations to obtain accurate free energy measurements. The grand-canonical partition function $Θ ( μ , V , T )$ can then be calculated for any value of $μ$ from Eq. (1) using the numerical values of Q(N, V, T, l = 0) (we note this function as Q(N, V, T) from now on and drop the l = 0 specification). The conditions for vapor-liquid coexistence are determined from the number distribution p(N) as follows. We start from p(N) defined as

We then solve numerically the equation below to determine the chemical potential at coexistence,

where N_b is the point at which the function p(N) reaches its minimum, and the left hand side and the right hand side of the equation correspond to the probability of the vapor and of the liquid phase, respectively. The other thermodynamic properties can then be determined through the usual statistical mechanics relations.^25,42

We use the SPC/E force field⁷¹ to model H₂O. Each molecule is described as a distribution of three Lennard-Jones (LJ) sites and three point charges (one on each atom) with the interaction between atoms (i, j) as

The parameters for the SPC/E model are taken from the work of Berendsen et al.⁷¹ To model the interaction between an atom of the fractional molecule with an atom of a full molecule, we scale the interaction parameters $ϵ i j , σ i j$ and the product q_iq_j by (l/M)^1/3, (l/M)^1/4 and (l/M)^1/3, respectively. We scale the size of the fractional molecule (i.e., the bond length between H and O) by (l/M)^1/4.

The SPC/E model⁷¹ treats water as a rigid molecule. This implies that the effect of the intramolecular vibrations on the thermodynamics of water is not taken into account, and that, in this case, q_int is equal to the rotational partition function of the water molecule. This gives

where $σ$ is the symmetry number (equal to 2 in the case of water), and I_A, I_B, and I_C are the 3 principal moments of inertia.

The electric field E is applied along the x-axis. This means that the interaction between the dipole moment m_i of water molecule i can be calculated as follows:

where E is the norm of the field E, $q α$ is the point charge associated with atom $α$ of molecule i, and $x i , α$ is the coordinates of atom $α$ of molecule i along the x-axis.

The use of a classical force field, like the SPC/E potential, to model water under an electric field^18,19,72 implies that the dissociation of water molecules is not taken into account in this work. Recent ab initio molecular dynamics (MD) simulations⁷³ have identified a dissociation threshold for water molecules for an electric field of about 0.35 V/Å. This is in agreement with the findings from previous experimental work^74,75 that reported the onset of water dissociation under external fields of about 0.32 V/Å–0.44 V/Å, and from previous simulation results,⁷⁶ that showed that external fields of 0.3 V/Å–0.6 V/Å enhance water ionization. In this work, we consider fields E of 0.05 V/Å, 0.1 V/Å, and 0.2 V/Å, i.e., well below the dissociation threshold and for which the dissociation of water molecules is extremely rare and short-lived.⁷³ We also include, for comparison, results for a larger field of 0.5 V/Å, a field for which ab initio MD indicates that 8% of water molecules are ionized.⁷³ Very interestingly, these results show that such a strong field yields to a qualitatively different behavior of water, even if water dissociation is not taken into account. Finally, we add that the use of a rigid model amounts to neglecting the elongation of the intramolecular O–H bonds as a result of the electric field. However, according to the ab initio MD results,⁷³ the range of electric fields studied here have a very moderate impact on the O–H bond length. Saitta et al.’s work reveal a non-monotonic dependence of the bond length on the electric field, with an O–H bond length varying between 1 Å and 1.01 Å.

We use 3 different types of Monte Carlo (MC) moves during the EWL simulations of water under an electric field. The first type of MC move (37.5% of the total number of moves) corresponds to the translation of a single water molecule (randomly chosen among the N + 1 molecules, i.e., the N full water molecules plus the fractional molecule). The second type of MC move (37.5% of the total number of moves) involves the rotation of a single molecule (randomly selected as one of the N full molecules or the fractional molecule). The 25% remaining moves are changes in (N, l) values for the system, resulting in the sampling of the range of N, and hence densities, studied in this work. Here we carry out simulations from N = 0 to up to N = 300 molecules in cubic cells with box lengths of L = 20 Å for all systems, with the usual periodic boundary conditions.⁷⁷ For the LJ part of the potential, we use a cutoff distance set to half the box length and apply tail corrections beyond that cutoff distance.⁷⁷ The long-range electrostatic interactions are handled using Ewald sums, with the screening parameter for the charge Gaussian distribution set to 5.6/L and the reciprocal cutoff vector set to $k m a x = 6 ( 2 π ) / L$⁠. The parameters for the EWL simulations are the same as in prior work.²⁷ The number of stages M is set to 100, the starting value for the convergence factor f in the iterative Wang-Landau scheme to e, and its final value to 10⁻⁸, with each (N, l) value being visited at least 1000 times for a given value of f.

We start by analyzing the results obtained for the grand-canonical partition function for water subjected to an electric field. Figure 1 shows that the behavior for this function at $T = 575 K$ for increasing values of the field. We also include in Fig. 1 the results obtained in the absence of field for comparison purposes. Figure 1 shows that increasing the field strength has two main effects on $ln Θ ( μ , V , T )$⁠. First, it leads to a shift in the chemical potential for which $ln Θ$ exhibits a steep increase. Since this sharp increase corresponds to the vapor $→$ liquid transition,⁷⁸ this means that $μ c o e x$⁠, the chemical potential at coexistence, becomes lower and lower with the field. At $T = 575 K$⁠, and taking as reference $μ c o e x$ in the absence of field, applying a field of $0.05 V / Å$ results in a decrease of 0.7% in $μ c o e x$⁠, while a field of $0.1 V / Å$ yields a decrease in $μ c o e x$ of 1.6%. This decrease becomes even more significant as the field is further increased with a value of $μ c o e x$ lower by 6.3% for a field of $0.2 V / Å$ and by 25% for a field of $0.5 V / Å$⁠. This result can be directly connected to the increase in the $ln Q ( N , V , T )$ function, shown in the inset of Fig. 1, since the slope of $ln Q ( N , V , T )$ is proportional to the chemical potential. Second, the steep increase in $ln Θ ( μ , V , T )$ becomes sharper and sharper with the field strength. This implies that the difference between the vapor and liquid phases becomes more and more important as the field increases. A corollary to this result is the fact that for a given temperature, applying an electric field will result in making the two phases more and more different or, in other words, postponing (in terms of temperature) the onset of criticality, for which the difference between the properties of the two phases vanishes. This suggests that the critical point increases as a result of the increase in field strength. This point will be further studied inSec. III B.

We now focus on the effect of temperature for a fixed field strength. Figure 2 shows the temperature dependence of $ln Θ ( μ , V , T )$ for a field of $0.2 V / Å$⁠. The plots exhibit two qualitatively different behaviors. A sharp increase, corresponding to the fluid undergoing the vapor $→$ liquid transition, is observed for the two lower temperatures. On the other hand, a much smoother increase in $ln Θ ( μ , V , T )$⁠, typical of a supercritical fluid for which there is no longer a transition, is found for the two higher temperatures. Increasing the temperature has a direct effect on $μ c o e x$ as shown by the shift observed between the results for $T = 300 K$ and for $T = 550 K$⁠. At $T = 550 K$⁠, $μ c o e x$ is 38% lower than $μ c o e x$ at $T = 300 K$⁠. As in Fig. 1, this shift in $μ c o e x$ can be correlated with the variations for the slope of $ln Q ( N , V T )$⁠, shown in the inset of Fig. 2, which is found to decrease as the temperature increases from $300 K$ to $550 K$⁠. We also find a markedly different behavior in $ln Q ( N , V , T )$ between the results for subcritical and supercritical fluids. While $ln Q ( N , V , T )$ is a monotonic function over the whole range of densities for subcritical fluids, we observe that, for a supercritical fluid, $ln Q ( N , V , T )$ is non-monotonic anymore and exhibits a maximum. Furthermore, $ln Q ( N , V , T )$ becomes less and less dependent on the temperature in the supercritical regime, which can be attributed to the lesser impact of the intermolecular interactions at very high temperatures.

Using the results obtained for the partition function, we calculate the number probability distribution p(N) to determine the densities at coexistence for the vapor-liquid transition²⁵ using Eqs. (8) and (9). We show in Fig. 3 the densities at coexistence obtained in the absence of field and for field strengths ranging from $0.05 V / Å$ to $0.5 V / Å$⁠. We also plot in Fig. 3 the critical point for each field value. The critical parameters are obtained as follows. We use a scaling law for the critical temperature,

where B is a fitting parameter, $β$ is the 3D Ising critical exponent adjusted for real substances (⁠$β = 0.326$⁠), and $ρ l$ and $ρ v$ are the densities for the liquid and vapor phases at coexistence given by the EWL simulations. The critical density is calculated from the law of rectilinear diameters

where $ρ c$ and A are two fitting parameters and T_c is the estimate for the critical temperature obtained from Eq. (13). The critical points so obtained are plotted in Fig. 3 and their numerical values are given in Table I.

The electric field has two main effects on the phase diagram of water. First, the whole phase diagram is shifted towards the higher temperatures. This can best be seen through the increase in T_c which starts from $641 K$ in the absence of field, and increases by 0.3% for $E = 0.05 V / Å$ and by 3 % for $E = 0.1 V / Å$⁠. This increase becomes even more pronounced as the field further increases with a T_c for $E = 0.2 V / Å$ greater by 6.6 % than in the absence of field, and a T_c for $E = 0.5 V / Å$ that is 10.1% above the value of T_c in the absence of field. Second, applying an electric field results in tilting the phase diagram, with a decrease in the critical density for very strong fields. More specifically, taking as the reference the critical density in the absence of field, we find that $ρ c$ decreases by 1.9% for $E = 0.05 V / Å$⁠, by 2.3% for $E = 0.1 V / Å$⁠, by 3.9% for $E = 0.2 V / Å$⁠, and by 6.1% for $E = 0.5 V / Å$⁠. This combined increase in T_c and decrease in $ρ c$ is consistent with the findings from Gibbs Ensemble Monte Carlo simulations on polar compounds.¹⁸ Comparing the results obtained for $T = 575 K$⁠, we find that the density of the liquid at coexistence increases with the field, whereas the density of the vapor at coexistence decreases with the field. Overall, the results confirm that fields of up to $0.05 V$Å have a limited impact on the densities at coexistence¹⁹ and on the critical parameters. This is consistent with the moderate changes in the partition functions for such fields, as shown in Fig. 1 and discussed in Sec. III A. We therefore focus in the rest of the paper on fields of the order of $0.1 V / Å$ and above.

We now turn to the results obtained for the thermodynamic properties at coexistence. Figure 4 shows the vapor pressure obtained from the grand-canonical partition function through

As for the densities at coexistence, the electric field essentially shifts the curve for the vapor pressure. For instance, in the absence of field, the vapor pressure reaches $50 b a r s$ around $550 K$⁠. For $E = 0.1 V / Å$⁠, the vapor pressure reaches this value at a temperature of about $565 K$⁠, i.e., 3% above. Similarly, this value for the vapor pressure is achieved at temperatures larger by 7% and by 12%, for fields of $0.2 V / Å$ and $0.5 V / Å$⁠, respectively, than in the absence of field. If we now look at the effect of the field at fixed temperature, we see that the vapor pressure is decreased by 23% if a field $E = 0.1 V / Å$ is applied, by 47% for $E = 0.2 V / Å$ and by 65% for $E = 0.5 V / Å$⁠. This decrease in vapor pressure with the field intensity can be directly related to the decrease in the density of the vapor phase at coexistence observed for strong fields.

Figure 4 also shows the variation of the chemical potential at coexistence with temperature. This plot highlights the decrease in $μ c o e x$ with the applied field. For instance, at $T = 575 K$ and using $μ c o e x$ in the absence of the field as the reference, we see that $μ c o e x$ is decreased by 3% for $E = 0.1 V / Å$⁠, by 7% for $E = 0.2 V / Å$⁠, and by 21% for $E = 0.5 V / Å$⁠. This shift towards the lower values of $μ$ directly stems from the results obtained for the partition function (see Fig. 1), which revealed that $ln Θ$ was increasingly shifted towards the lower end of the $μ$ range as the applied field increased. The bottom panel in Fig. 4 shows the dependence of the entropy of vaporization $Δ S$ on the field. We find that the curve for $Δ S$ against T is shifted towards the top of the temperature range for strong fields. Using the results obtained in the absence of field at $T = 575 K$ as the reference, we find an increase in $Δ S$ of 22% for $E = 0.1 V / Å$⁠, of 35% for $E = 0.2 V / Å$⁠, and of 49% for $E = 0.5 V / Å$⁠. This can be attributed to the greater difference (in terms of density) between the two phases at coexistence for strong fields, as evidenced by the larger value for $( ρ l − ρ v )$ shown in the phase diagrams (see Fig. 3). However, another possible cause for this large $Δ S$ may also be the structural changes induced by strong electric fields, a point that we aim to elucidate in the next paragraph.

We focus here our analysis on the structure of the more organized phase, the saturated liquid, since the structure of vapor phases of water is known to be much less sensitive.¹⁸ We start by analyzing the radial distribution functions for the saturated liquid at $T = 575 K$⁠. We show in Fig. 5 a comparison between the radial distribution functions g_OO(r) obtained for different fields. The applied field only has a mild effect on this distribution function, which is consistent with prior results on low temperature liquids of water, which only revealed notable changes in radial distribution functions^79–82 for fields in excess of $0.5 V / Å$⁠. We are, however, able to notice a slight decrease in the height of the first peak in g_OO(r) for high field strengths, which is reminiscent of the general effect of other external fields (e.g., shear) on these functions.^83,84 We show in Fig. 5 the results for the g_OH(r) distribution functions, which confirm the moderate effect of the field on the radial distribution functions^79–82 for fields up to $E = 0.5 V / Å$⁠. This set of results would a priori indicate that the structure of the fluid has not been dramatically changed by the applied field. To ascertain this point, we refine our analysis by computing the average values for different order parameters for the saturated liquid.

We calculate the three following order parameters: q_t, q_l, and N_H. q_t is tetrahedral order parameter which characterizes the arrangement of the four nearest neighbors of each water molecule.^85,86 q_t is calculated by averaging over all water molecules i the local quantity defined as follows:

In this equation, j and k are two water molecules chosen among the four nearest neighbors of molecule i, and $ϕ j i k$ is the angle between the vectors connecting i and its two neighbors j and k. q_t is equal to 1 if the environment around i is perfectly tetrahedral.

q_l measures the alignment of water molecules in the direction of the electric field. It is calculated by averaging over all water molecules i the following quantity:

where $θ ( i )$ is the angle between the electric field E and m_i, the dipole moment of molecule i. q_l reaches a value of 1 when all dipoles are perfectly aligned with the field, a value of 0 when there is no preferred orientation of the dipole, and a value of −0.5 when the dipoles are perpendicular to the field.

N_H is the number of hydrogen bonds within the saturated liquid. N_H is obtained by applying the following 3 criteria:^87–90 (i) the OO distance between two water molecules must be less than 3.5 Å, (ii) the $O ⋯ H$ distance between the O of the first molecule and the H of the second molecule involved in the hydrogen bond is less than 2.5 Å, and (iii) the $O ⋯ H – O$ angle along the hydrogen bond is less than $30 °$⁠.

We plot in Fig. 6 the variations along the coexistence line of the 3 order parameters for the saturated liquid. Unlike the radial distribution functions that very weakly depended on the field, the order parameters all reveal that the amount of structural order in the liquid increases with the field. Figure 6 shows that an electric field of $0.1 V / Å$ increases q_t by 3% and N_H by 4%. Applying a field of $0.2 V / Å$ leads to an increase in q_t of 7% and of N_H by 8%. Similarly, a field of $0.5 V / Å$ results in an increase by 15% for both q_t and N_H. This increase in both q_t and N_H values concomitantly occurs with an increase in the alignment order parameter q_l which goes from 0 (no preferred orientation) in the absence of field to up to more than 0.6 when the applied field is of $0.5 V / Å$⁠. This shows that the dramatic increase in the entropy of vaporization for strong fields (up to 49% for $E = 0.5 V / Å$⁠) therefore results from the cumulative effects of the greater density difference between the two coexisting phases and of the greater structural organization in the saturated liquid.

The results for the partition functions can also be used to shed light on the thermodynamics of compressed liquids. For this purpose, we vary the chemical potential $μ$ and calculate the resulting number distribution p(N) through Eq. (8) and the thermodynamic properties from p(N).^25,27 We show in Fig. 7 the results for the density $ρ$⁠, the Gibbs free energy G, and the entropy S of compressed liquids at $T = 575 K$⁠. In this plot, the variations of $ρ$⁠, G, and S with pressure are reported in the absence of field and for fields ranging from $0.1 V / Å$ to $0.5 V / Å$⁠. At fixed pressure, the dependence of the liquid density on the field is as follows: $ρ$ (⁠$0 V / Å$⁠) $< ρ ( 0.1 V / Å ) < ρ ( 0.2 V / Å ) < ρ ( 0.5 V / Å$⁠). For instance at $P = 500 b a r s$⁠, the density for $E = 0.1 V / Å$ is 5% greater than in the absence of field. Similarly, for $E = 0.2 V / Å$⁠, the density is greater by 10% than in the absence of field, while for $E = 0.5 V / Å$⁠, the density is 15% greater than in the absence of field. The increased density as a result of the applied field has two main consequences on liquid properties as evidenced by the behavior observed for G and S of the liquid. For $P = 500 b a r s$ and using G in the absence of field as the reference, G is found to decrease by 2% for $E = 0.1 V / Å$⁠, 7% for $E = 0.2 V / Å$ and 21% for $E = 0.5 V / Å$⁠. This decrease in G can be attributed to the decrease in potential energy due to the larger number of attractive water-water interactions per unit volume and to the stronger interaction of water with the field, the latter becoming increasingly significant as the field gets stronger. We now focus on the variations of S with the field. At $P = 500 b a r s$⁠, we find that S decreases by 4% when a field of $0.1 V / Å$ is applied, by 6% when a field of $0.2 V / Å$⁠, and by 10% when a field of $0.5 V / Å$⁠. This decrease in S can be attributed to the increase in the density of the fluid, but also to the greater organization within the compressed liquid arising from the applied field, as discussed below.

We plot in Fig. 8 the results obtained for the three order parameters q_t, q_l, and N_H at $P = 500 b a r s$⁠. All 3 order parameters point to a greater organization within the liquid as the field gets stronger. q_t is found to increase by 2% for $E = 0.1 V / Å$⁠, by 5% for $E = 0.2 V / Å$⁠, and by 9% for $E = 0.5 V / Å$⁠. The degree of alignment also steadily increases with the field and shows that q_l goes from 0 in the absence of field to up to 0.66 for $E = 0.5 V / Å$⁠. The number of hydrogen bonds within the compressed liquid is also found to increase with the field. With respect to the liquid at $P = 500 b a r s$ in the absence of field, we find that N_H is greater by 3% for $E = 0.1 V / Å$⁠, by 7% for $E = 0.2 V / Å$⁠, and by 11% for $E = 0.5 V / c$⁠. These results, together with the increase in density with the field, account for the observed decrease in entropy for strong fields.

We now examine the impact of the field on the thermodynamic properties of supercritical water. Figure 9 shows plots of $ρ$⁠, G, and S against pressure. The impact of the field on the density of supercritical fluids is more limited than for liquids. For instance, for $P = 1000 b a r s$⁠, the increase in density only starts to become noticeable for $E = 0.2 V / Å$⁠, which results in a density increase of 3% with respect to the density in the absence of field. Increasing further the field to $0.5 V / Å$ yields a density 10% greater than in the absence of field. The impact on G and S follows a similar trend. For G, increasing the field from $0 V / Å$ to $0.2 V / Å$ decreases the G by about 1.5% and an increase in the field from $0 V / Å$ to $0.5 V / Å$ decreases G by 6%. As for the liquid, we attribute this decrease in G to the greater number of attractive water-water interactions and to the greater interaction energy with the field. For S, decreases in entropy with respect to the value in the absence of field are of 1% and 4% for fields of $0.2 V / Å$ and $0.5 V / Å$⁠, respectively. This stems from the increase in density with the field and from the increase in structural order within the fluid (see Fig. 10). Figure 10 shows that $E = 0.5 V / Å$ for q_t is, e.g., greater by 10% than in the absence of field, that q_l reaches up to 0.3, and that N_H increases by 24%.

We continue our characterization of the thermodynamics of supercritical water with the determination of the ideality contours for water. These contours^29,30 have recently emerged as a means to bridge the gap in our understanding of the supercritical region of the phase diagram,^36–39 since they provide a way to map these regions and to establish a correspondence between the supercritical states of different fluids.⁴¹ Here we focus on two of these contours. We start by determining the locus for the Zeno line, defined as the contour along which supercritical water behaves as an ideal gas from the standpoint of the ideal gas law. This locus is obtained by varying numerically $μ$ such that the following condition is satisfied:

where $V ¯$ is the reciprocal density and $N ¯ = ∑ N p ( N )$ is the average number of particles in the system. Previous work in the field has focused on establishing the shape of this contour for the van der Waals equation,²⁹ model systems,⁹¹ Argon,³⁰ metals,^33,35,40,92 and a few molecular fluids including water.^31,41 Remarkably, it has been shown that the Zeno line is a straight line that extends over several hundred degrees. However, the effect of an electric field on the shape of the Zeno line has not been studied so far, and it remains to be seen how the field impacts the shape of this contour.

The second contour we study in this work is the H line, defined as the curve of ideal enthalpy.^29,30 We obtain the locus for this contour by varying numerically $μ$ such that the condition written below is obeyed,

Similar to the Zeno line, the locus for this contour has been shown to correspond to a straight line for the van der Waals equation,²⁹ Argon,³⁰ and several molecular fluids.⁴¹ It has, however, been studied much less extensively than the Zeno line and, to our knowledge, the impact of an electric field on the H line has yet to be investigated.

We report in Fig. 11 the sets of $( T , ρ )$ satisfying the conditions given in Eqs. (18) and (19) for water in the absence of field and for fields ranging from $0.1 V / Å$ to $0.5 V / Å$⁠. We then fit the EWL results to linear laws (also shown in Fig. 11). Both the Zeno and H line remain remarkably straight lines regardless of the field strength. To establish further this point, we evaluate the departure of the Zeno lines from straight lines by calculating the Average Absolute Relative Deviation (AARD) error (⁠$A A R D ( % ) = 1 N ∑ i | T F i t − T s i m T s i m |$⁠). For all fields, the AARD remains small (between 1% and 3%) and of the same order as in the absence of field (3.7%). This means that the shape of the Zeno line is not altered by the electric field and remains the same as in the absence of the field, thus opening the door for the application of the similarity laws based on the Zeno line⁹² to systems under an electric field. As shown in Fig. 11, only high temperature results are available for the H line, since the low temperature/high density domain for the H line lie within the domain of stability of the solid. These linear fits allow us to determine the Boyle and H parameters, which are key input parameters in the similarity laws of Apfelbaum and Vorob’ev.⁹³ The Boyle parameters are obtained as the intercept of the Zeno line with the temperature axis (this gives the Boyle temperature T_B) and with the density axis (this provides the Boyle density $ρ B$⁠). The results for the Boyle parameters are summarized in Table I. While the Boyle density remains essentially constant throughout the range of fields considered in this work, the Boyle temperature steadily increases with the field strength. T_B is greater by 6% for $E = 0.1 V / Å$ than in the absence of field, by 7.5% for $E = 0.2 V / Å$⁠, and by 18.5% for $E = 0.5 V / Å$⁠. Similarly, the intercepts of the H line with the temperature and density axes provide the two parameters T_H and $ρ H$ (given in Table I). In line with the results obtained for T_c and T_B, T_H is shown to increase with the field, with a 7% increase for $E = 0.1 V / Å$⁠, a 13% increase for $E = 0.2 V / Å$⁠, and a 61% increase for $E = 0.5 V / Å$⁠. As observed for $ρ c$⁠, $ρ H$ is found to decrease with the field with a 4% decrease for $E = 0.1 V / Å$ and a 13% increase for $E = 0.5 V / Å$ (both with respect to the value for T_H in the absence of field).

To provide a comparison between the binodal curve and ideality contours obtained with and without the field, we rescale the phase diagram of water by the critical temperature and critical density found for each value of the field. Figure 12 shows that the binodals obtained for all fields can all be superimposed on top of one another. This shows that the behavior of subcritical water remains qualitatively the same for very strong fields (up to $0.5 V / Å$⁠). Similarly, the loci obtained for the Zeno line in the scaled temperature-density plane are in good agreement for all fields, with maximum deviations of less than 7% for the scaled T_B and $ρ B$ for all fields. Similar conclusions only apply for the H line for fields of up to $0.2 V / Å$⁠. This can best be seen on the scaled T_H for $E = 0.5 V / Å$⁠. For this field, the scaled T_H is 46% greater than in the absence of field and exhibits a markedly larger value than for all other conditions. This shows that while a correspondence between the results obtained for different fields can be made on the basis of the scaled phase diagrams for fields of up to $0.2 V / Å$⁠, the behavior of water subjected to fields of $0.5 V / Å$ and beyond is qualitatively different. This discrepancy observed in the H line for the strongest field is likely due to the predominant contribution of the potential energy due to the interaction of water molecules with the applied field, and the strong alignment of water molecules with the field, along the H line in these low density-supercritical states.

In this work, we extend the expanded Wang-Landau simulation method to determine the impact of an electric field on the phase diagram of water, on the thermodynamic properties of the vapor-liquid transition, of compressed liquids, and of supercritical phases of water, and on the loci for the ideality contours of water. Through staged insertions/deletions of water molecules, the EWL method allows us to calculate the grand-canonical partition function of water under an electric field and to determine its properties using the statistical mechanics formalism. Our results show that the impact of the electric field on the partition function becomes significant for fields greater than 0.05 V/Å and that it steadily increases with the field. This, in turn, has a number of important consequences on the phase diagram with a shift of the binodal towards the lower densities/higher temperatures as the field gets stronger. This also leads to a decrease in the chemical potential at coexistence (by up to 21%) and to a significant increase in the entropy of vaporization (by up to 49%). This result is attributed to the greater difference between the two phases at coexistence, both in terms of density and structural order, as shown by the analysis of the dependence on the field of the tetrahedral order parameter, of the extent of the alignment of water molecules with the field and of the number of hydrogen bonds. These conclusions extend to the thermodynamics of compressed liquids and of supercritical phases of water. Finally, a correspondence between the results obtained for different field strength is carried out through the analysis of the ideality contours. The results show that the Zeno and H lines remain straight over the entire range of electric fields studied in this work. This correspondence, however, starts to break down for a field of 0.5 V/Å, which lead to markedly different H parameters as a result of the predominance of the field-water interaction in these low density-supercritical states. Further work using more sophisticated models for water, including, e.g., polarizable models, quantum effects,⁹⁴ and allowing for water dissociation,^73,76,95,96 will allow for a refinement of these findings.

Partial funding for this research was provided by NSF through CAREER Award No. DMR-1052808.