An Equation of State for Metals at High Temperature and Pressure in Compressed and Expanded Volume Regions

Abstract

A simple equation of state model for metals at high temperature and pressure is described. The model consists of zero-temperature isotherm, thermal ionic components, and thermal electronic components, and is applicable in compressed as well as expanded volume regions. The three components of the model, together with appropriate correction terms, are described in detail using Cu as a prototype example. Shock wave Hugoniot, critical point parameters, liquid–vapor phase diagram, isobaric expansion, etc., are evaluated and compared with experimental data for Cu. The semianalytical model is expected to be useful to prepare extended tables for use in hydrodynamics calculations in high-energy-density physics.

1. Introduction

Equation of State (EOS) of materials is an inevitable ingredient in several fields of solid state science like geophysics, hydrodynamic applications for the analysis of inertial confinement fusion systems, stellar structures, nuclear weapons, etc. Other applications include fast reactor accident analysis and study of weapon effects in various media. Euler equations of hydrodynamics, which expresses conservation laws of mass, momentum, and energy, are routinely used to describe the dynamical behavior of materials [1]. However, these equations describe the space-time evolution of four thermodynamic variables—viz., mass density (or specific volume), material velocity, specific internal energy, and pressure. The system of these equations is then closed with the addition of equation of state (EOS), which provides pressure when specific internal energy and density are given. The Mie–Grüneisen EOS [1] with an empirical specification for Grüneisen parameter is the most commonly used EOS of this type. There is also Tillotson’s EOS [2] which has a larger range of validity. However, a more complete EOS is specified by providing pressure and specific internal energy as functions of density and temperature. This temperature corresponds to thermodynamic equilibrium in the material, and can be eliminated from the expression for pressure in favor of specific internal energy to obtain the above mentioned relation between pressure, specific internal, energy, and density. An EOS of similar class, which treats pressure as independent variable, was proposed by Rice and Walsh [3] to model water. Here, specific volume is expressed in terms of enthalpy using the enthalpy–parameter which depends on pressure. This class of EOS, generally called enthalpy-based EOS, has been developed to model shock compression of porous materials [4], including explicit accounting of electronic effects [5].

In the following sections, we describe a semi-analytical EOS model of the general type, where volume and temperature are independent variables. The different components of the EOS model, including correction terms, are discussed in detail. Experimental data for Cu on Shock-Hugoniot, critical point parameters, liquid–vapor phase diagram, and isobaric expansion are used to test the model. Good agreement obtained shows that the model can be employed to prepare extended EOS tables for use in hydrodynamics calculations.

2. Three Component EOS Model

In general, EOS models consist of three components [6], which describe (i) the zero-temperature (or cold) isotherm, (ii) thermal ionic effects, and (iii) thermal electronic effects. Pressure and specific internal energy are then expressed as functions of volume $\left(V\right)$ and temperature $\left(T\right)$:

$$\begin{array}{c}\hfill P(V,T)=P_c\left(V\right)+P_{ti}(V,T)+P_{te}(V,T),\end{array}$$

(1)

$$\begin{array}{c}\hfill E(V,T)=E_c\left(V\right)+E_{ti}(V,T)+E_{te}(V,T),\end{array}$$

(2)

where the terms denote, respectively, the three components mentioned above. The subscripts c, $ti$, and $te$ stand for the terms ‘cold’, ‘thermal-ion’, and ‘thermal-electron’, respectively. There are interaction effects between ionic and electronic motion, however, these effects contribute only a few percent to pressure and energy and so may be neglected in high-pressure physics applications [6]. Models for the different terms in Equations (1) and (2) of varying degrees of sophistication are currently available in the literature [7].

First, principle methods using Density-functional theory (DFT) for $E_c$ and quasi-harmonic approximation (QHA), based on density-functional perturbation theory, for lattice vibration contributions to $E_{ti}$ are now quite common [8]. Such methods have proven to be extremely useful in thermodynamic studies of compounds of interest in Earth sciences [9]. Further, electronic density of states determined from DFT calculations yield accurate estimates of $E_{te}$ for lower ranges of temperature [10]. Availability of efficient computer implementation [11] of QHA, which uses DFT-generated data on $E_c$ and volume-dependent vibration frequencies, make QHA the method of choice for detailed studies of thermodynamic properties of materials in the solid phase. However, for developing global EOS models, which deal with very high temperatures (∼keV) and pressures (∼ tens of megabar), it is necessary to take only the relevant information from DFT computations and supplement it with other models [10] to incorporate effects of melting, extreme pressure and thermal ionization, expanded volume states, etc. The general approach is to use empirical fits to the cold isotherm [12], Debye–Grüneisen model for lattice thermal motion [13] and Thomas–Fermi model for thermal electron excitation [14]. Such extended models are essential for some of the hydrodynamic applications mentioned in the beginning.

The global EOS model we describe below, which is applicable even at very high temperatures and pressures in the compressed as well as expanded volume states, indeed uses different parameters obtained from DFT analyses, particularly when accurate experimental data on these are unavailable.

3. Zero-Temperature Isotherm

The zero-temperature isotherm is a manifestation of the Fermi-pressure developed in degenerate electron systems, and is a quantum effect just like zero-point vibration energy. This contributes significantly to the total pressure in compressed solids, and becomes the dominant contribution at extreme compression. A variety of approximate expressions to describe it quantitatively are available in the literature [15]. Computations using DFT, mentioned above, are now routinely used to determine energy versus specific volume (or volume per atom) tables, and thereafter the zero-temperature pressure–volume relation. Results of such analyses are then used in semi-empirical expressions. We propose to use a four-parameter model, developed by Li et al. [16], which is expressed as:

$$\begin{array}{ccc}\hfill E_{Li}\left(V\right)& =& -E_{coh}(1+a+\delta a^3)e^{-a},\eta =(\frac{9B_0{V}_{c0}}{E_{coh}}{)}^{1/2},x=(\frac{V}{V_{c0}}{)}^{1/3},\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill P_{Li}\left(V\right)& =& 3B_0\frac{(1-x)}{x^2}(1-3\delta a+\delta a^2)e^{-a},a=\eta (x-1),\delta =\frac{B_0^{\prime }-1}{2\eta }-\frac{1}{3}.\hfill \end{array}$$

(4)

The four parameters in the model are the specific volume $V_{c0}$, the bulk modulus $B_0$, its pressure derivative $B_0^{\prime }$, and the cohesive energy $E_{coh}$ at zero temperature. These parameters occur in terms of dimensionless quantities $\eta$ and $\delta$, and a is related to the dimensionless length variable x. Furthermore, if energy $E_{Li}$ is scaled with $E_{coh}$ and pressure $P_{Li}$ with $B_0$, then, these expressions are totally dimensionless—however, defined in terms of two parameters $\eta$ and $\delta$. The specific volume $V_{c0}$ is slightly lower than the volume at ambient conditions ($T_0=300$ K and $P_0=1$ bar). In our approach, we adjust the value of $V_{c0}$ such that the zero-temperature pressure together with thermal pressure of ions and electrons is just one bar at 300 K. The four-parameter model is a refinement over Rose equation [17] and Vinet equation [18], and is found to provide quite accurate descriptions of the zero-temperature energy and pressure over compressed volume up to ∼$V_{c0}/2$, which corresponds to about 100–150 GPa pressure, for about forty metals [16]. It also provides accurate representation for energy and pressure in the expanded volume up to about ∼$2V_{c0}$. However, the formulation is inadequate in the region of extreme compression, as is evident from Equation (3), which shows that $E_{Li}$ saturates as $V\to 0$. Theoretically, the zero-temperature energy and pressure should approach those of electron gas.

To rectify this problem, we use a procedure [12] to smoothly go over from the four-parameter model to the the quantum statistical model (QSM) [19], which is known to provide accurate descriptions of pressure and energy of electrons above few hundred GPa pressure. The QSM accounts for exchange and correlation effects in addition to corrections for electron density gradients [20]. Electron pressure in a compressed atom within the QSM model is expressed as:

$$\begin{array}{ccc}\hfill P_{QSM}\left(V\right)& =& \frac{e^2}{5}\frac{a_0}{V^{5/3}}{\left(3{\pi }^2\right)}^{2/3}{Z}_n^{5/3}exp[-\alpha -\beta ],\hfill \\ \hfill \alpha & =& 0.3225R_w{Z}_n^{[0.495-0.039{log}_{10}{Z}_n]},\hfill \\ \hfill \beta & =& \frac{5}{3}{R}_w^2[0.068+0.078{log}_{10}{Z}_n-0.086{\left({log}_{10}{Z}_n\right)}^2].\hfill \end{array}$$

(5)

Here, e is electron charge, $Z_n$ is atomic number, $a_0$ is the Bohr radius, and $R_w$ is the Wigner-Seitz cell radius in units of $a_0$. Specific internal energy $E_{QSM}$ is obtained from pressure by integrating the thermodynamic relation $P=-dE/dV$ from a suitable initial volume, say $V_{c0}$.

Now, choose a value of V, say $V_m$, such that the four-parameter model $E_{Li}\left(V\right)$ is accurate for $V\ge V_m$. That is, we assume that the zero-temperature isotherm $E_{cold}\left(V\right)=E_{Li}\left(V\right)$ and $P_{cold}\left(V\right)=P_{Li}\left(V\right)$ for $V\ge V_m$. Then, for lower values of V, these are defined as

$$\begin{array}{ccc}\hfill E_{cold}\left(V\right)& =& (E_{QSM}\left(V\right)-E_{QSM}\left(V_m\right))B_{int}\left(V\right)+E_{Li}\left(V_m\right),V\le V_m,\hfill \\ \hfill P_{cold}\left(V\right)& =& P_{QSM}\left(V\right)B_{int}\left(V\right)+(E_{QSM}\left(V\right)-E_{QSM}\left(V_m\right))B_{int}^{\prime }\left(V\right),V\le V_m,\hfill \\ \hfill B_{int}\left(V\right)& =& (1+b_{1}V+b_2{V}^{4/3}+b_3{V}^{5/3}).\hfill \end{array}$$

(6)

Here, $B_{int}\left(V\right)$ is a suitable interpolating function. Note that, by definition, $E_{cold}\left(V\right)$ is continuous at $V_m$. Now, the parameters $b_k$$(k=1,3)$ in $B_{int}\left(V\right)$ are chosen such that $P_{cold}\left(V\right)$ and its first two derivatives are also continuous at $V_m$ [12]. This procedure gives a smooth transition from the four-parameter model to the QSM. Plots of energy versus V for Cu using the two models are shown in Figure 1A with the choice $V_m=V_{c0}/1.4=0.07998$ cm${}^3$/g.

As an application of the zero-temperature energy $E_{cold}\left(x\right)$, we use the lattice inversion method [21] and obtain an effective inter-particle potential between Cu atoms in the solid. We may imagine that the lattice is formed by assembling shells successively around a central atom. Then, the zero-temperature energy $E_{cold}\left(x\right)$ per atom, where $x={(V/V_{c0})}^{1/3}$ is the scaled nearest-neighbor distance, can be readily written as a lattice sum involving the inter-particle potential $U\left(x\right)$. The inversion method provides a similar formula [22] for $U\left(x\right)$ in terms of $E_{cold}\left(x\right)$. For an FCC lattice, the direct formula for $E_{cold}\left(x\right)$ and the inverse formula for potential $U\left(x\right)$ are given by

$$\begin{array}{c}\hfill E_{cold}\left(x\right)=\frac{1}{2}\sum _{n=1}^{N_s}{r}_{n}U\left(b_{n}x\right),\\ \hfill U\left(x\right)=2\sum _{n=1}^{N_s}{w}_n{E}_{cold}\left(b_{n}x\right).\end{array}$$

(7)

Here, n is the shell index, $r_n$ is the number of atoms in the shell, $b_n$ is the normalized radius of the shell, and $w_n$ is the weight factor for the shell. This is a truncated formula, and the total number $\left(N_s\right)$ of shells considered should be sufficiently large for convergence. The total potential (curve-1) and its repulsive (curve-2) and attractive (curve-3) components, as per the Weeks–Chandler–Andersen prescription [23], are shown in Figure 1B. For the sake of completeness, we have listed the constants $r_n,b_n$ and $w_n$$(1\le n\le N_s)$ for FCC lattice in

Table 1. Constants of inversion formula for FCC lattice.

n	${\mathit{r}}_{\mathit{n}}$	${\mathit{b}}_{\mathit{n}}$1	${\mathit{w}}_{\mathit{n}}$	n	${\mathit{r}}_{\mathit{n}}$	${\mathit{b}}_{\mathit{n}}$1	${\mathit{w}}_{\mathit{n}}$	n	${\mathit{r}}_{\mathit{n}}$	${\mathit{b}}_{\mathit{n}}$1	${\mathit{w}}_{\mathit{n}}$
1	12	$\sqrt{1}$	1/12	6	8	$\sqrt{6}$	1/9	11	24	$\sqrt{11}$	−1/6
2	6	$\sqrt{2}$	−1/24	7	48	$\sqrt{7}$	−1/3	12	24	$\sqrt{12}$	7/12
3	24	$\sqrt{3}$	−1/6	8	6	$\sqrt{8}$	1/32	13	72	$\sqrt{13}$	−1/2
4	12	$\sqrt{4}$	−1/16	9	36	$\sqrt{9}$	1/12	14	0	$\sqrt{14}$	−1/3
5	24	$\sqrt{5}$	−1/6	10	24	$\sqrt{10}$	0	15	48	$\sqrt{15}$	1/3

¹ For FCC lattice $b_n=\sqrt{n}$. These are different for other lattices [21].

. The factors $1/2$ and 2 in Equation (7) arise because $U\left(x\right)$ is the energy for two atoms; while $E_{cold}\left(x\right)$ is the energy per atom.

4. Ionic Thermal Component

The treatment of the ionic contribution to the EOS is best described by considering the compressed region and the expanded volume regions separately. In the former, we start with the low-temperature solid phase and go over to the melted fluid phase on increasing temperature. However, in the latter region, the material is always in the fluid phase, which also encompasses what is called the warm dense matter region. We consider different models for describing the EOS in these regions.

4.1. Compressed Region—Johnson’s Model

The model for ionic thermal energy should describe the low-temperature properties of solids, fluid phase for temperatures above melting, and ideal gas behavior at higher temperatures. Thus, for a specified volume, the constant-volume molar specific heat of ions—denoted by $C_{Vi}$—must vary from the low temperature $T^3$ law to $3R$ above Debye’s temperature and finally to $3R/2$. This feature is essential since shock compression of materials produce high temperatures and the solid melts and becomes a fluid after shock traversal. The mean field model used by Wang does not possess this crucial feature [24]. The parameters in the model used by Kormer et al. [25] need adjustments for every material.

So, we propose to use Johnson’s ionic model [26] in the compressed region. In addition to the general constraints on $C_{Vi}$ mentioned above, the model adds an extra contribution (3RT/Tm), linearly varying with T in the interval $T_m$ to $1.2T_m$ to account for the heat of fusion. This corresponds to an increment of $0.6R$ in entropy, which has been determined from studies of several materials. Furthermore, the typical decrease of $C_{Vi}$ from its value $3R$ at $T_m$ to $9R/4$ at $5T_m$—and thereafter a linear variation in $ln\left(T\right)$ to the ideal gas value $3R/2$—are built in to the model. Specific internal energy and pressure within the model (version-I) are thus given by [26]:

$$\begin{array}{ccc}{E}_{ti}(V,T)\hfill & =& E_D+N{k}_{B}T(E_0+{ϵ}_{\psi }).\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill P_{ti}(V,T)& =& \frac{{\mathsf{\Gamma }}_i}{V}{E}_D+\frac{1}{V}(2{\mathsf{\Gamma }}_i-2/3)N{k}_{B}T(E_0+{ϵ}_{\psi }),\hfill \\ \hfill E_D(V,T)& =& N{k}_{B}T[\frac{9}{8}\frac{{\theta }_D}{T}+9(\frac{T}{{\theta }_D}{)}^3\underset{0}{\overset{{\theta }_D/T}{{\displaystyle \int }}}\frac{z^3}{e^z-1}dz].\hfill \end{array}$$

(9)

Here, $E_D$ is Debye’s specific internal energy, $T_M\left(V\right)$ is melting temperature, ${\theta }_D\left(V\right)$ is Debye’s temperature, $k_B$ is Boltzmann’s constant, N is number of atoms per gram, and $\psi =T/T_M$ is scaled temperature. Further, ${\mathsf{\Gamma }}_i$(V) is Grüneisen parameter for ions, to be defined below. The energy parameters $E_0$ and ${ϵ}_{\psi }$ are fitted functions of $\psi$ in order to account for the constraints on $C_{Vi}$ mentioned above, and are given by

$$\begin{array}{ccc}\hfill E_0& =& -\frac{3}{2}+\frac{3}{2}{a}_4({\psi }^{-3/2}-{\psi }^{-2}/2)+a_2\frac{(a_{3}y+{\psi }^{1-y})}{\left({\psi }^y{(a_3+{\psi }^{1-y})}^2\right)},1\le \psi <\infty ,\hfill \\ \hfill {ϵ}_{\psi }& =& \frac{3}{2}(\psi -1/\psi ),1\le \psi \le 1.2,\hfill \\ \hfill {ϵ}_{\psi }& =& 0.66/\psi ,1.2\le \psi <\infty .\hfill \end{array}$$

(10)

Thus, Johnson starts with Debye’s model in the region $T\le T_m$. The contribution to specific internal energy due to heat of fusion is given by ${ϵ}_{\psi }$; and $E_0$ describes the variation of $C_{Vi}$ after melting. Note that $E_0$ varies from zero at $T_m$ to $-3/2$ in the high temperature limit. An equally important feature is that the factor $(2{\mathsf{\Gamma }}_i-2/3)$ facilitates correct approach of the effective ${\mathsf{\Gamma }}_i$ to its ideal gas limit (2/3). The constants $a_2,a_3,a_4$, which take care of the constraints, are given by

$$\begin{array}{ccc}\hfill a_1& =& -5.7-\sum _i{n}_{i}ln[{(A_i/\overline{A})}^{3/2}/n_i],a_3=200,\hfill \\ \hfill y& =& (201{[1600a_1^2+2398(4a_1+5)]}^{1/2}-40(5-197a_1))/\left[3980(4a_1+5)\right],\hfill \\ \hfill a_2& =& \frac{3}{2}{(1+a_3)}^3/\left[a_3(1-y)(a_{3}y+2-y)\right],\hfill \\ \hfill a_4& =& -\frac{8}{5}(a_1+a_2/(1+a_3)),\hfill \end{array}$$

(11)

where, $n_i$, $A_i$, and $\overline{A}$ denote the number fraction, mass number, and its average, respectively. The summation in $a_1$, which contributes only for mixtures, accounts for ideal entropy mixing. Thus, the model is applicable to the case of compounds as well.

In a more elaborate method (version-II), region $1\le \psi <\infty$ is divided into three segments, viz., $1\le \psi \le 5$, $5\le \psi \le \zeta$, and $\zeta \le \psi <\infty$. In the first region, $C_{Vi}$ varies linearly in T; while in the the second, the variation is linear in $ln\left(T\right)$. The value of $\zeta$ is determined so that entropy approaches ideal gas limit. Heat of fusion is added, as in the first version. Finally, specific energy and pressure are expressed as

$$\begin{array}{ccc}\hfill E_{ti}(V,T)& =& E_D+N{k}_{B}T(E_{\psi }+{ϵ}_{\psi }),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill P_{ti}(V,T)& =& \frac{{\mathsf{\Gamma }}_i}{V}{E}_D+\frac{1}{V}(2{\mathsf{\Gamma }}_i-2/3)N{k}_{B}T(E_{\psi }+{ϵ}_{\psi }),\hfill \end{array}$$

(13)

where ${ϵ}_{\psi }$ is the same as that given in Equation (11), while $E_{\psi }$ is given by

$$\begin{array}{ccc}\hfill E_{\psi }& =& 3/16-3\psi /32-3/\left(32\psi \right),1\le \psi \le 5,\hfill \\ \hfill E_{\psi }& =& -(3/4)+bln(\psi /5)-b+5(b+9/20)/\psi ,5\le \psi \le \zeta ,\hfill \\ \hfill E_{\psi }& =& -(3/2)+5(b+9/20)/\psi -\zeta b/\psi ,\zeta \le \psi <\infty .\hfill \end{array}$$

(14)

The new constants b and $\zeta$ are given by

$$\begin{array}{ccc}\hfill b& =& 9/\left(32\right(a_1+3/4+27ln\left(5\right)/16\left)\right),\hfill \\ \hfill \zeta & =& 5exp[-3/(4b\left)\right].\hfill \end{array}$$

For illustration, we show in Figure 2 the variation of specific heat of Cu with temperature at normal volume $V_0$, using first version in graph-A and second version in graph-B. Both versions, thus, produce almost identical results. From the discussion summarized above, it is clear that Johnson’s model is to be used in the compressed volume region. So, we next discuss a suitable model in the expanded volume region.

4.2. Expanded Region—Modified Soft-Sphere

The expanded volume region of the material, which covers the warm dense region and liquid–vapor transition, is important in several hydrodynamic applications. For instance, the material undergoing expansion induced via rarefaction wave is in this region. Similarly, highly porous materials (e.g., copper with porosity more than fifty percent) reach this state after shock-compression. It is necessary to have a separate model for this region as the physics here is mainly determined with excluded volume effect and a weak van der Waals type attractive interaction [27]. Young developed a soft-sphere model for liquid metals [28] using Monte Carlo simulations data for the thermal properties of particles interacting via soft-sphere potential, $ϵ{(\sigma /r)}^n$. Here, $ϵ$, $\sigma$, and n define the parameters of the inverse power law potential. Together with the van der Waals attractive interaction, five parameters in the model were fitted to liquid–vapor co-existence data. This procedure is an extension of van der Waals theory of fluids, which treats attractive part as a zero temperature component of the EOS. We developed a modified soft-sphere model [29], wherein the use of simulation data was retained, however a generalized Lennard–Jones model was used as the attractive component. Instead of numerical fitting, all the parameters of the modified model were determined in terms of $V_{c0}$, $E_{coh}$, $B_0$, and $B_0^{\prime }$ at normal conditions.

To discuss the formulation briefly, we start with a generalized Lennard–Jones pair potential between the atoms:

$$U_{LJ}\left(R\right)={ϵ}_{rep}(\frac{\sigma }{R}{)}^n-{ϵ}_{att}(\frac{\sigma }{R}{)}^m.$$

(15)

Here, respectively, ${ϵ}_{rep}$ and ${ϵ}_{att}$ denote repulsive and attractive energy parameters, and $\sigma$ is a spatial scale factor. This particular form is chosen because, as already mentioned above, extensive simulation data is available [27] for the repulsive component in Equation (15). Further, a proper choice of the exponent m will lead to an attractive part similar to that in van der Waals theory. The zero-temperature ionic energy per gram, $E_{LJ}=N{\mathcal{E}}_{LJ}$, where ${\mathcal{E}}_{LJ}$ is the energy per atom, can be expressed as

$$E_{LJ}\left(V\right)=N\frac{{ϵ}_{rep}}{2}{C}_n(\frac{V_{c0}}{V}{)}^{n/3}-N\frac{{ϵ}_{att}}{2}{C}_m(\frac{V_{c0}}{V}{)}^{m/3},$$

(16)

where $V_{c0}=N{\sigma }^3/\sqrt{2}$ is the specific volume at zero-temperature, and $C_n$ and $C_m$ are lattice sums [30]. Now, on imposing the conditions that the zero-temperature pressure $P_{LJ}=-d{E}_{LJ}/dV$ vanishes at $V_{c0}$ and the corresponding energy $E_{LJ}\left(V_{c0}\right)=E_{coh}$, we readily find that

$${ϵ}_{rep}=\frac{E_{coh}}{N}\frac{2}{C_n}\frac{m}{n-m},{ϵ}_{att}=\frac{E_{coh}}{N}\frac{2}{C_m}\frac{n}{n-m}.$$

(17)

The expressions for $E_{LJ}$ and $P_{LJ}$ can be re-written as

$$\begin{array}{ccc}\hfill E_{LJ}\left(V\right)& =& \frac{E_{coh}}{n-m}\left[m\right(\frac{V_{c0}}{V}{)}^{n/3}-n\left(\frac{V_{c0}}{V}{)}^{m/3}\right],\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill P_{LJ}\left(V\right)& =& \frac{E_{coh}}{n-m}\frac{nm}{3V_{c0}}\left[\right(\frac{V_{c0}}{V}{)}^{1+n/3}-\left(\frac{V_{c0}}{V}{)}^{1+m/3}\right].\hfill \end{array}$$

(19)

The exponents n and m are yet to be determined, however, bulk modulus $B_0$ and its pressure derivative $B_0^{\prime }$ can be computed from Equation (19) as

$$B_0=E_{coh}\frac{nm}{9V_{c0}},B_0^{\prime }=2+\frac{n}{3}+\frac{m}{3}.$$

(20)

Even though these relations imply that n and m are the roots of the quadratic equation $x^2-3(B_0^{\prime }-2)x+9B_0{V}_{c0}/E_{coh}=0$, experimental parameters generally lead to complex roots [31]. Therefore, in such situations, we can use only three parameters. Jiuxun has shown [32] that the spinodal condition $\left|B\right|\sim \sqrt{(P-P^{*})}$, where $P^{*}$ is the pressure corresponding to $B\left(P^{*}\right)=0$, follows if we use the relation $m=(n-3)/2$. Then, the first relation in Equation (20) shows that n is determined from the quadratic equation $n^2-3n-18B_0{V}_{c0}/E_{coh}=0$. The positive root gives $n=8.803$ and $m=2.901$ for Cu. The accuracy of the generalized Lennard–Jones model can be assessed by comparing it with the four-parameter model [16] discussed earlier. In Figure 3A, we compare the repulsive component of the interparticle potential for Cu. The energy parameter ${ϵ}_{rep}$ and the exponent n in the repulsive component are to be used in the model for ionic thermal energy. Figure 3B shows the generalized Lennard–Jones energy and pressure, which compare quite well.

As explained in the beginning of this section, the ionic thermal energy of the modified soft-sphere model is precisely contributed by the repulsive part of the inter-particle potential. Monte Carlo simulation data of thermal pressure of soft particles, interacting via the potential ${ϵ}_{rep}{(\sigma /r)}^n$, are accurately expressed in the parametric form [27]:

$$P=\frac{N{k}_{B}T}{V}[1+\frac{n}{3}\frac{C_n}{2}{\left(\frac{{\rho }_{sc}}{\sqrt{2}}\right)}^{n/3}+\frac{1}{18}n(n+4){\left(\frac{{\rho }_{sc}}{\sqrt{2}}\right)}^{n/9}],{\rho }_{sc}=\frac{N}{V}{\sigma }^3[\frac{{ϵ}_{rep}}{k_{B}T}{]}^{3/n}.$$

(21)

We show the accuracy of this fit in Figure 4, where the scaled pressure $Z=PV/\left(N{k}_{B}T\right)$ is plotted versus the scaled density, denoted as ${\rho }_{sc}$, for values of exponent $n=4,5,6$. Similar results for the exponent $n=7,8,9,12$ are shown in Figure 5. Impressive agreement between the data and the parametric fit is evident in these figures. On adding this contribution to $E_{LJ}$ and $P_{LJ}$, the ionic energy and pressure within the modified soft-sphere model are given by

$$\begin{array}{ccc}\hfill E_{soft}(V,T)& =& N{k}_{B}T[\frac{3}{2}+\frac{1}{6}(n+4)(\frac{V_{c0}}{V}{)}^{n/9}\left(\frac{{ϵ}_{rep}}{N{k}_{B}T}{)}^{1/3}\right]+E_{LJ}\left(V\right),\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill P_{soft}(V,T)& =& \frac{N{k}_{B}T}{V}[1+\frac{1}{18}n(n+4)(\frac{V_{c0}}{V}{)}^{n/9}\left(\frac{{ϵ}_{rep}}{N{k}_{B}T}{)}^{1/3}\right]+P_{LJ}\left(V\right),\hfill \end{array}$$

(23)

where we have used the expression for zero-temperature reference volume, $V_{c0}=N{\sigma }^3/\sqrt{2}$, and the subscript $soft$ to denote the term ‘soft-sphere’.

As mentioned earlier, the zero-temperature isotherms of the solid and the correction to ideal gas thermal energy (second term in brackets), determined from Monte Carlo simulation data, define the modified soft-sphere model. All the parameters in this model are determineda priori, so there is no need to fit its parameters.

4.3. Grüneisen Parameter

Ionic Grüneisen parameter ${\mathsf{\Gamma }}_i\left(V\right)$, Debye temperature ${\theta }_D\left(V\right)$, and melting temperature $T_M\left(V\right)$ are needed in modeling ionic energy and pressure in the compressed volume region. While there are several empirical fits [13] for ${\mathsf{\Gamma }}_i\left(V\right)$, the one due to Preston et al. [33], which has the correct asymptotic behavior for strong compression, is expressed as

$${\mathsf{\Gamma }}_i\left(V\right)=\frac{1}{2}+c_1{V}^{1/3}+c_2{V}^q.$$

(24)

The parameters $c_1,c_2$, and q are determined using experimental data on ${\mathsf{\Gamma }}_i$ at $T=300$ K, at zero-pressure melting point and asymptotic ($V\to 0$) approximation for free electron states. These constraints define three nonlinear equations to determine the parameters. Expressions for ${\theta }_D\left(V\right)$ and $T_M\left(V\right)$, which follow from Debye–Grüneisen law and Lindemann’s law, respectively, are given by

$$\begin{array}{ccc}\hfill {\theta }_D\left(V\right)& =& {\theta }_0(\frac{V_0}{V}{)}^{1/2}exp[3c_1(V_0^{1/3}-V^{1/3})+\frac{c_2}{q}(V_0^q-V^q)].\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill T_M\left(V\right)& =& T_{M0}(\frac{V_r}{V}{)}^{1/3}exp[6c_1(V_r^{1/3}-V^{1/3})+\frac{2c_2}{q}(V_r^q-V^q)].\hfill \end{array}$$

(26)

Values of ${\theta }_0$, $T_{M0}$, and $V_r$, which reference Debye temperature, melting temperature, and melting volume at $P=0$, respectively, are tabulated for a variety of materials [33]. Further, the expression for $T_M\left(V\right)$ is found to compare well with experimental data, as shown below.

As an application of the models discussed in this section, we show the Grüneisen parameter for Cu in the compressed ($V\le V_0$) as well as the expanded ($V\ge V_0$) volume regions for 300 K (curve-1) and 5000 K (curve-2) in Figure 6A. The region to the left of the (vertical) dashed line in the figure is the compressed phase while the expanded region is on its right side. The thermodynamic definition for ionic Grüneisen parameter, ${\mathsf{\Gamma }}_i=(1/V){(\partial P/\partial E)}_V=(1/V){(\partial P/\partial T)}_V(1/C_{Vi})$, can be readily used with analytical expressions for P and E, which are described in both volume regions (Johnson’s model for $V\le V_0$; and the modified soft-sphere model for $V\ge V_0$). Experimental value 2.19 at 300 K is also shown (filled circle) in the figure. For a specified temperature, say 5000 K, the material goes over to the melted fluid region as volume is increased, and hence ${\mathsf{\Gamma }}_i\left(V_0\right)$ is reduced to 1.96 from its value at 300 K. For the expanded volume region, the modified soft-sphere model shows explicit temperature dependent of ${\mathsf{\Gamma }}_i\left(V\right)$.

In Figure 6B, we have also sketched ${\theta }_D\left(V\right)$ and $T_M\left(V\right)$, given in Equations (25) and (26). Experimental data on $T_M$ taken from Preston et al. [33], and the value of ${\theta }_D\left(V_0\right)$ are also shown. Good comparison of $T_M$ data, which are not used in fixing the parameters in ${\mathsf{\Gamma }}_i$ ($c_1,c_2$, and q in Equation (24)), shows the internal consistency of the procedure. In addition, we have shown ${\theta }_D$ (dashed line) derived from Slater’s formula, ${\mathsf{\Gamma }}_i=(1/6)-(1/2)(d{B}_S/dP)$ [11], where $B_S$ is the zero-temperature bulk modulus obtained from Equation (4). Again, the agreement between ${\theta }_D$ obtained from two totally independent sources (Preston’s and Slater’s ${\mathsf{\Gamma }}_i$) demonstrates its accuracy.

5. Electronic Thermal Component

Electronic thermal component of energy and pressure is significant at temperatures reached in shock compression of porous materials. The EOS codes mentioned earlier uses the well-known Thomas–Fermi (TF) model [14] to describe electronic properties. However, the pressure and energy resulting from the approximations involved in this model need to be corrected in the low-temperature range. For example, the low-temperature specific heat constant predicted by the model differs from experimental values. Further, it is necessary to employ pressure and energy tables since in-line solutions of the TF equation are time consuming. Therefore, analytical fits [34] to results of Latter’s calculations are sometimes employed [25] in high-pressure studies. However, it is important to note that this fit is valid in the compressed volume region, even though TF model as such may be applied even to an isolated atom.

So, we propose a somewhat different approach. First of all, following Atzeni et al [35], we assume that the Fermi gas model can be used to compute the thermal energy and pressure of electrons with a suitable average ionization degree $Z^{*}$ of atoms, which depends on density and temperature. We use an excellent analytical fit for $Z^{*}$ given by More [36] using results of Thomas–Fermi model:

$$\begin{array}{ccc}\hfill R_5& =& \rho (g/c{m}^3)/\left(Z_{n}A\right),tc=T\left(ev\right)/A^{4/3},t=tc/(1+tc),\hfill \\ \hfill a_5& =& 3.323\times {10}^{-3}\times t{c}^{0.971832}+9.26148\times {10}^{-5}\times t{c}^{3.10165},\hfill \\ \hfill b_5& =& -exp[-1.763+1.43175\times t+0.315463\times t^2],\hfill \\ \hfill c_5& =& -0.36667\times t+0.98333,q_6=a_5\times R_5^{b_5},\hfill \\ \hfill q_5& =& {(\mathbf{factor}\times R_5^{c_5}+q_6^{c_5})}^{1/c_5},x_5=14.3139\times q_5^{0.6624},\hfill \\ \hfill Z^{*}& =& Z_n\times x_5/(1+x_5+{[1+2x_5]}^{1/2}),\hfill \end{array}$$

(27)

where $Z_n$ and A are atomic number and mass number, respectively. We have introduced a multiplicative correction term, termed ‘factor’ in the equation above, which needs to be adjusted so that $Z^{*}$ agrees with the experimental value of average ionization degree at ${\rho }_0$ and $T_0$. It takes value 0.079 for Cu, and the corrected and uncorrected variation of $Z^{*}$ with temperature is shown in Figure 7A. Note that $Z^{*}$ is corrected only in the lower ranges of temperature, and it correctly saturates to $Z_n$ at high temperature. The insert in this figure shows specific heat variation for three densities: (1) ${\rho }_0/10$, (2) ${\rho }_0$, and (3) $10{\rho }_0$. Similar correction factors for some other metals are given in

Table 2. $Z_0^{*}$ using More’s formula [36].

Element	${\mathit{\rho }}_0(\mathit{g}/\mathit{c}\mathit{c})$	${\mathit{Z}}_{0\mathit{U}\mathit{C}}^{*}$1	${\mathit{Z}}_{0\mathit{E}}^{*}$1	factor	${\mathit{Z}}_{0\mathit{C}}^{*}$1
Cu	8.93	4.38	1.01	0.079	1.02
Al	2.74	2.46	1.0	0.210	1.09
Fe	7.89	5.85	2.01	0.270	2.06
W	19.41	4.15	1.34	0.095	1.36

¹$Z_{0UC}^{*}$ (uncorrected), $Z_{0E}^{*}$ (desired), $Z_{0C}^{*}$ (corrected).

.

To test this approach, we have plotted in Figure 7B the zero-temperature electron pressure $P=C_p{(N{Z}^{*}/V)}^{5/3})$, where $C_p={(3/\pi )}^{2/3}(h^2/20m_e)$ versus density $\left(\rho \right)$ using the uncorrected (curve-1) and corrected (curve-2) formula for $Z^{*}$. Numerical results [37] employing TF theory (filled circles) are also shown. It is interesting to note the TF results are accurately reproduced with the uncorrected $Z^{*}$, while the correction gives lower pressures. Finally, the thermal component of electron energy and pressure, within the free electron gas model, are given by

$$\begin{array}{ccc}\hfill \frac{N}{V}{Z}^{*}& =& C_0{\left(k_{B}T\right)}^{3/2}{I}_{1/2}[E_F/k_{B}T)].\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill E_e[V,T]& =& C_0{\left(k_{B}T\right)}^{5/2}V{I}_{3/2}[E_F/k_{B}T)]-\frac{3}{5}\left(N{Z}^{*}\right)E_F(V,0).\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill P_e[V,T]& =& \frac{2}{3}\frac{1}{V}{E}_e[V,T],C_0=\frac{4\pi }{h^3}{\left(2m_e\right)}^{3/2}.\hfill \end{array}$$

(30)

Here, $E_F(V,T)$ is the Fermi energy, which is implicitly defined via Equation (28); and $I_n$ are Fermi–Dirac integrals. Very accurate rational approximations for these integrals are now available [38]. The next level of improvement to the electron EOS is to add Coulomb interaction, exchange, and correlation energies, thereby accounting for all the terms in the uniform electron gas model [39].

6. Applications

In this section, we discuss three more applications of the EOS model. The first is its application to the calculation of the shock Hugoniot. Extensive data available in the shock wave database [40] is compared in Figure 8A with predictions of the model (curve-1). We find excellent agreement throughout the range of pressures obtained in the experiments. In the insert figure, we have shown the results without correcting the zero-temperature isotherm (curve-2), and without adding electronic terms (curve-3). Correction to the zero-temperature isotherm is found to be quite important as the four-parameter isotherm is valid up to 100–150 GPa. Addition of the electron component is found to improve the prediction even for pressures in the range of 200 GPa. Contrary to the common feeling, pressure on the Hugoniot is decreased when electron contribution is added. This is because the electron degrees of freedom reduces temperature, and consequently the pressure, for a specified volume on the Hugoniot. Temperature along the Hugoniot (curve-1) and the melting temperature (curve-2) displayed in Figure 8B show that melting occurs around 300 GPa pressure. Therefore, proper accounting of the melting transition is important even though the transition is not evident in the pressure–volume Hugoniot.

Next, we compare the liquid–vapor phase diagram of Cu employing the modified soft-sphere model, which was briefly considered in our earlier work [29]. It is well known that the attractive and thermal components of energy finely balance to produce the van der Waals loops in the isotherms in the vapor–liquid co-existence region. The critical point parameters we have obtained, via Maxwell’s construction, (${\rho }_c=2.246$ g/cm${}^3$, $T_c=8345$ K, and $P_c=0.8935$ GPa) are very well within the range quoted in the literature [41]. The phase diagram (curve-1) is shown in Figure 9A, and compared with simulation data (filled circles) [42]. These data were obtained via molecular dynamics simulations using an effective pair potential deduced from DFT calculations of energy–volume curve in the compressed and expanded volume regions. We find good comparison except in the liquid region of the phase diagram. The spinodal curve (curve-2), which is the locus of points where isothermal compressibility diverges; and the diameter (curve-3), which is the average of liquid and vapor phase densities on the phase diagram, are also shown in the figure.

In Figure 9B, we have compared results of the EOS model for the enthalpy (H) of Cu versus temperature, along the 1-bar isobar, with experimental data taken from Trainor et al. [43]. Starting from 50 K, the material expands from 0.1105 to 0.1268 (cm${}^3$/g), where it crosses the melting line around 1310 K. This results in a jump in enthalpy, and thereafter a smooth increase, as seen from the figure. The model for ${\mathsf{\Gamma }}_i$ employed in Johnson’s model is inappropriate in this volume region. Therefore, we have taken a constant value ${\mathsf{\Gamma }}_i=2.19$, which is the experimental value at normal conditions, and corresponding expressions for ${\theta }_D\left(V\right)$ and $T_M\left(V\right)$ in lieu of Equations (25) and (26). Furthermore, we find that the heat of fusion to be added at the melting point is about 225 kJ/kg, in good agreement with the experimental value of 205 kJ/kg.

Finally, we mention that the EOS model described above (except the correction term for the zero-temperature isotherm and the free electron gas model) has been applied very successfully for the development and application of enthalpy-based approach to describe shock wave propagation in porous materials. There, we started with the computation of enthalpy-parameter, defined as $\chi =P{(\partial V/\partial H)}_P$, and developed the method treating P as independent variable. First, we showed [5] (using analytical fits [34] to results of TF theory) the correct approach to incorporate electronic effects explicitly. Then, we developed the modified soft-sphere model to properly evaluate the Hugoniot of highly porous materials, as their final shocked states are in the fluid region [29]. In another development, we used the EOS to show that the enthalpy-based approach can be implemented [44] in hydrodynamic simulations to describe shock wave propagation in solid as well as porous materials.

Applications of the EOS model to other materials can be easily carried out using the experimental/theoretical values of the parameters [16,33] employed in the different components of the model. A list of these parameters for Cu, Al, Fe, and W is given in our earlier publication [44].

7. Summary

The main aim of this paper is to discuss the basic components of an EOS model for metals for high-pressure physics applications. Thus, we started with the division of the EOS in to three parts: Zero-temperature isotherm, thermal ionic component, and thermal electronic component. This division is convenient as there are theories of different levels of sophistication dealing with them. We mentioned that results of DFT-based electronic structure calculations can be suitably fitted into a functional form, as is done in the case of the four-parameter model. A method to correct this model—which is only valid in the lower pressure range, so as to approach the results of the quantum statistical model—was discussed. Similarly, in lieu of detailed lattice dynamic calculations for density of states of lattice vibration modes, we have used the simple Debye model with a suitable form for volume dependence of Debye temperature. However, we have stressed the need for incorporating the energy of fusion due to melting in the thermal ionic component. A model for this component, which incorporates continuous temperature dependence of constant volume-specific heat, was discussed. It is the valence electrons in the metal which contribute to the thermal electronic component, which we modeled using the results of the Fermi gas model. However, we have shown that the degree of ionization, and hence the free electron density, can be determined from the Thomas–Fermi model. A simple correction to get the experimental electron density at normal conditions was also discussed. Finally, for the purpose of demonstration, the model was applied to calculate the shock-Hugoniot, liquid-vapor phase diagram, and isobaric expansion of Cu. We hope that the model described here can be used to generate tables for hydrodynamic applications of impact experiments, shock wave studies, and above all, design and analysis of high-energy-density systems.