Generalization of the Unified Analytic Melt-Shear Model to Multi-Phase Materials: Molybdenum as an Example

Abstract

The unified analytic melt-shear model that we introduced a decade ago is generalized to multi-phase materials. A new scheme for calculating the values of the model parameters for both the cold $(T=0)$ shear modulus $\left(G\right)$ and the melting temperature at all densities $\left(\rho \right)$ is developed. The generalized melt-shear model is applied to molybdenum, a multi-phase material with a body-centered cubic (bcc) structure at low $\rho$ which loses its dynamical stability with increasing pressure (P) and is therefore replaced by another (dynamically stable) solid structure at high $\rho$. One of the candidates for the high-$\rho$ structure of Mo is face-centered cubic (fcc). The model is compared to (i) our ab initio results on the cold shear modulus of both bcc-Mo and fcc-Mo as a function of $\rho$, and (ii) the available theoretical results on the melting of bcc-Mo and our own quantum molecular dynamics (QMD) simulations of one melting point of fcc-Mo. Our generalized model of $G(\rho ,T)$ is used to calculate the shear modulus of bcc-Mo along its principal Hugoniot. It predicts that G of bcc-Mo increases with P up to ∼240 GPa and then decreases at higher P. This behavior is intrinsic to bcc-Mo and does not require the introduction of another solid phase such as Phase II suggested by Errandonea et al. Generalized melt-shear models for Ta and W also predict an increase in G followed by a decrease along the principal Hugoniot, hence this behavior may be typical for transition metals with ambient bcc structure that dynamically destabilize at high P. Thus, we concur with the conclusion reached in several recent papers (Nguyen et al., Zhang et al., Wang et al.) that no solid-solid phase transition can be definitively inferred on the basis of sound velocity data from shock experiments on Mo. Finally, our QMD simulations support the validity of the phase diagram of Mo suggested by Zeng et al.

1. Introduction: Unified Analytic Melt-Shear Model

A reliable model of the adiabatic (isentropic) shear modulus $\left(G\right)$ of a polycrystalline solid at all temperatures from zero to close to $T_m,$ the melting temperature, and up to pressures $\left(P\right)$ of order 100 GPa is needed for many applications, including the modeling of plastic deformation at extremes of pressure and temperature, numerical calculations of elastic and shock wave propagation, and even calculations of the oscillations of low-mass astrophysical objects. In 2003–2004 we developed a unified analytic model for the melting temperature as a function of density, $T_m\left(\rho \right)$, and the density-temperature dependence of the shear modulus, $G(\rho ,T)$, at all densities and $0\le T\le T_m$ [1,2]. The model for $T_m\left(\rho \right)$ is based in part on our analytic model for the Grüneisen parameter,

$$\gamma \left(\rho \right)=\frac{1}{2}+\frac{{\gamma }_1}{{\rho }^{1/3}}+\frac{{\gamma }_2}{{\rho }^q},$$

(1)

where the parameters ${\gamma }_1$, ${\gamma }_2$, and q can be determined from ambient data by solving a system of three algebraic equations [1,2]. Equation (1) in conjunction with the Lindemann melting criterion [3] yields

$$T_m\left(\rho \right)=T_m\left({\rho }_{ref}\right){\left(\frac{\rho }{{\rho }_{ref}}\right)}^{1/3}exp\left\{6{\gamma }_1\left(\frac{1}{{\rho }_{ref}^{1/3}}-\frac{1}{{\rho }^{1/3}}\right)+\frac{2{\gamma }_2}{q}\left(\frac{1}{{\rho }_{ref}^q}-\frac{1}{{\rho }^q}\right)\right\},$$

(2)

where ${\rho }_{ref}$ is a conveniently chosen reference density. The melting temperature and shear modulus along the solidus approximately satisfy the relation [1,4,5]

$$\frac{G(\rho ,T_m\left(\rho \right))}{\rho T_m\left(\rho \right)}=\frac{G({\rho }_{ref},T_m\left({\rho }_{ref}\right))}{{\rho }_{ref}{T}_m\left({\rho }_{ref}\right)}.$$

(3)

Preston and Wallace [6] modeled $G(\rho ,T)$ as a linear function of the reduced temperature $T/T_m$ with the correct value $G(\rho ,0)$ at $T=0:$

$$G(\rho ,T)=G(\rho ,0)\left(1-\beta \frac{T}{T_m\left(\rho \right)}\right).$$

(4)

The thermoelastic softening parameter $\beta$ may be density dependent. Combining Equations (2)–(4) results in

$$G(\rho ,0)=G({\rho }_{ref},0){\left(\frac{\rho }{{\rho }_{ref}}\right)}^{4/3}exp\left\{6{\gamma }_1\left(\frac{1}{{\rho }_{ref}^{1/3}}-\frac{1}{{\rho }^{1/3}}\right)+\frac{2{\gamma }_2}{q}\left(\frac{1}{{\rho }_{ref}^q}-\frac{1}{{\rho }^q}\right)\right\}.$$

(5)

Equations (2), (4), and (5) constitute our unified model of $T_m\left(\rho \right)$ and $G(\rho ,T)$ for single-phase metals.

For polymorphic (multi-phase) metals, the ${\gamma }_1$-${\gamma }_2$-q values for the ambient solid structure may differ significantly from those for the high-P structure(s), and therefore the model may fail in its description of those high-P structure(s). Here we generalize our unified melt-shear model to polymorphic materials. Let $G_i(\rho ,0)$ be the cold shear modulus of phase i, and similarly let $T_{mi}\left(\rho \right)$ denote the melting curve of phase i. In our multi-phase model both $G(\rho ,0)$ and $T_m\left(\rho \right)$ are accurate envelopes of the $G_i(\rho ,0)$ and $T_{mi}\left(\rho \right)$. We shall explain how such envelopes are constructed, and demonstrate the fidelity of the generalized model using molybdenum as an example.

2. Molybdenum: An Example of a Multi-Phase Material

Molybdenum is one of the most important technological materials. Although its phase diagram is not known in detail, there are several features of the phase diagram that have been firmly established. First, melting on the shock Hugoniot of Mo occurs at ∼400 GPa [7,8,9], and presumably at $T\sim$ 10,000 K [7]. Second, the stability of body-centered cubic (bcc) Mo has been confirmed experimentally to 560 GPa under room-T isothermal compression [10], and to 1000 GPa under ramp-wave compression [11]. There is presently no experimental evidence for a solid-solid (s-s) transition in Mo. However, there is compelling theoretical evidence for a s-s phase transition at high pressures [12,13,14,15,16]. Our calculations of the phonon spectrum of bcc-Mo presented below (see Figure 1) clearly show its dynamic instability at 1000 GPa and low T, and previous simulations indicated that bcc-Mo dynamically destabilizes at ∼700 GPa [17,18]. Hence, a s-s phase transformation to another solid phase that is dynamically stable is expected to occur at ∼700 GPa at low T. Two candidates for the high-P structure of Mo are face-centered cubic (fcc) [17,19,20] and double hexagonal close-packed (dhcp) [12,18,21]. Simulations show that the bcc-fcc (bcc-dhcp) phase boundary has positive slope at $P\ge 700$ GPa [12,17]. We emphasize that this theoretical phase boundary is consistent with the aforementioned experimental results for room-T isothermal compression to 560 GPa, and with ramp-wave compression to 1000 GPa. The ramp-wave trajectory (RWT) in the P-T plane has a positive slope and lies at sufficiently high temperatures that the bcc-fcc (bcc-dhcp) transition boundary lies below the RWT; e.g., $T\sim 5000$ K on the RWT at 700 GPa [11], and at 800 GPa the bcc-fcc transition is predicted to occur at ∼5000 K [17] or ∼4500 K [12], while $T\sim 6500$ K on the RWT [11]. Therefore at $P\le 700$ GPa the crystal structure is bcc from 0 to $T_m$, and is also bcc at high T to at least 1000 GPa.

For our present purposes of demonstrating the generalization of the unified melt-shear model, it is not essential to definitively identify the high-P structure (fcc or dhcp). Numerous examples show that the elastic properties of different hexagonal arrangements of the same material (AB for hcp, ABC for fcc, ABAC for dhcp, etc., where A, B and C denote the layer stacking in the unit cell) are very close to each other; e.g., the results of [22] reveal that the shear moduli of both hcp-Be and fcc-Be are virtually identical, and so are the corresponding melting curves [23]. Hence, we expect both the shear moduli and the melting curves of fcc-Mo and dhcp-Mo to be close to each other. Although our QMD simulations strongly suggest that the high-density phase of Mo is dhcp, simply because fcc elastic constants are generally easier to calculate than the dhcp ones, for the purpose of comparing our ab initio results to the new generalized melt-shear model here we assume that the high-P structure of Mo is fcc.

An analysis of longitudinal sound speed data along the principal Hugoniot of Mo by Errandonea et al. [24] shows a decrease in the shear modulus with increasing pressure. They concluded that such behavior is indicative of a transformation of bcc-Mo into another solid phase. However, the conclusion of earlier theoretical work by Cazorla et al. [25] and Zeng et al. [17] is that there is no s-s phase transformation on the Mo Hugoniot. More recent calculations of $C_L$ along the Hugoniot in bcc-Mo by Lukinov et al. [16] again suggests that a s-s transition occurs at ∼200 GPa. However, we have carefully analyzed the calculations of Lukinov et al., and we find that their conclusion is incorrect. We maintain that a s-s transition along the Mo Hugoniot is ruled out by both experiment and theory. A detailed discussion of the experimental and theoretical/computational investigations of a possible s-s transition on the Mo Hugoniot is presented in Appendix A.

In the following we construct a thermoelasticity model for bcc-Mo that predicts a decrease in G with P on the principal Hugoniot at pressures above ∼240 GPa, hence the decrease in G with P is not evidence for another solid phase of Mo. Thus our model supports the conclusion of References [8,9,26] that no s-s phase transition should be inferred from the sound velocity data on Mo.

3. Generalization of the Unified Melt-Shear Model to Multi-Phase Materials

As mentioned above, the unified melt-shear model [1,2] is based on the simplifying approximation that both $G(\rho ,0)$ and $T_m\left(\rho \right)$ are related to a common Grüneisen parameter. In the general case, this approximation may fail. Indeed, since $G(\rho ,0)\sim {\Theta }_D^2(-3,\rho ){\rho }^{1/3}$ [27,28],

$$\frac{dlnG(\rho ,0)}{dln\rho }\approx 2\gamma (-3,\rho )+\frac{1}{3},\gamma (-3,\rho )\equiv \frac{dln{\Theta }_D(-3,\rho )}{dln\rho }.$$

(6)

Similarly, since in the Lindemann formulation $T_m\left(\rho \right)\sim {\Theta }_D^2(-2,\rho ){\rho }^{-2/3}$ [27,28],

$$\frac{dln{T}_m\left(\rho \right)}{dln\rho }\approx 2\gamma (-2,\rho )-\frac{2}{3},\gamma (-2,\rho )\equiv \frac{dln{\Theta }_D(-2,\rho )}{dln\rho }.$$

(7)

In the above equations, ${\Theta }_D(n,\rho )\sim \hslash {\langle {\omega }_i^n\rangle }^{1/n}/k_B,n\ge -3$ are the Debye characteristic temperatures [28]. As noted in [27], Wallace argued that ${\Theta }_D(0,\rho )$ should be used in the Lindemann melting formula, Equation (7), instead of ${\Theta }_D(-2,\rho ).$ The two gammas appearing in Equations (6) and (7) must asymptote to 1/2 as $\rho \to \infty$ [1,2,28]. For single-phase materials they are approximately equal at all densities and are accurately modeled as $\gamma \left(\rho \right)=1/2+{\gamma }_1/{\rho }^{1/3}+{\gamma }_2/{\rho }^q$ [29]. However, this approximate equality of $\gamma (-3,\rho )$ and $\gamma (-2,\rho )$ most likely fails for multi-phase materials having solid phases that dynamically stabilize (or destabilize) at elevated $P$. For example, consider a high-P solid structure stable at low T but unstable at high $T,$ and assume that the s-s phase boundary has a positive slope. Then there is a range of P in which the cold shear modulus exists and is described by Equation (6), but $T_m$ does not exist for this high-P phase since the material melts from another stable solid phase. The same Grüneisen gamma cannot be used for both $G(\rho ,0)$ and $T_m\left(\rho \right)$. If the solid-solid boundary has negative slope then there is a range of P for which $T_m$ exists and is described by Equation (7) but the cold shear modulus does not exist, and again a single Grüneisen gamma cannot be used for both cold G and $T_m$.

In such cases we assume that both $G(\rho ,0)$ and $T_m\left(\rho \right)$ can be constructed as the envelopes of the $G_i(\rho ,0)$ and $T_{mi}\left(\rho \right)$ for a sequence of solid phases, but they are not necessarily coupled to each other in terms of a common set of parameters. We introduce the corresponding two Grüneisen gammas, ${\gamma }_G\left(\rho \right)$ and ${\gamma }_{T_m}\left(\rho \right),$ so that Equations (6) and (7) are generalized to

$$\frac{dlnG(\rho ,0)}{dln\rho }\equiv 2{\gamma }_G\left(\rho \right)+\frac{1}{3},$$

(8)

$$\frac{dln{T}_m\left(\rho \right)}{dln\rho }\equiv 2{\gamma }_{T_m}\left(\rho \right)-\frac{2}{3}.$$

(9)

Using the approach developed in [2], it is easy to show that both gammas must be analytic functions of ${\rho }^{-1/3},$ essentially the interatomic distance. Since both gammas asymptote to 1/2 as $\rho \to \infty ,$

$${\gamma }_G\left(\rho \right)=\frac{1}{2}+\sum _{i=1}^{\infty }{g}_{Gi}{\left({\rho }^{-1/3}\right)}^i\equiv \frac{1}{2}+\frac{{\gamma }_{1,G}}{{\rho }^{1/3}}+\frac{{\gamma }_2}{{\rho }^{q_2}},$$

(10)

$${\gamma }_{T_m}\left(\rho \right)=\frac{1}{2}+\sum _{i=1}^{\infty }{g}_{T_{m}i}{\left({\rho }^{-1/3}\right)}^i\equiv \frac{1}{2}+\frac{{\gamma }_{1,T_m}}{{\rho }^{1/3}}+\frac{{\gamma }_3}{{\rho }^{q_3}}.$$

(11)

Below we show that ${\gamma }_1$ is common to both G and $T_m,$ and that

$${\gamma }_{1,G}\approx {\gamma }_{1,T_m}\equiv {\gamma }_1=\frac{7}{40}{Z}^{2/3},$$

(12)

where Z is the atomic number of the material. In Equations (10) and (11) the terms ${\gamma }_2/{\rho }^{q_2}$ and ${\gamma }_3/{\rho }^{q_3}$ represent the contributions of the quadratic and higher-order terms. Note that since these equations are meant to describe the envelopes of the $G_i(\rho ,0)$ and $T_{mi}\left(\rho \right)$ for each of the stable solid phases of the material, the signs in front of the terms $g_{Gi}/{\rho }^{i/3}$ and $g_{T_{m}i}/{\rho }^{i/3}$ can in principle be arbitrary, and so may be the signs of both ${\gamma }_2$ and ${\gamma }_3.$ Based on general grounds, however, we expect that both $q_2$ and $q_3$ are positive.

With Equations (10)–(12) it follows from Equations (8) and (9) that

$$G(\rho ,0)=G({\rho }_0,0){\left(\frac{\rho }{{\rho }_0}\right)}^{4/3}exp\left\{6{\gamma }_1\left(\frac{1}{{\rho }_0^{1/3}}-\frac{1}{{\rho }^{1/3}}\right)+\frac{2{\gamma }_2}{q_2}\left(\frac{1}{{\rho }_0^{q_2}}-\frac{1}{{\rho }^{q_2}}\right)\right\},$$

(13)

$$T_m\left(\rho \right)=T_m\left({\rho }_m\right){\left(\frac{\rho }{{\rho }_m}\right)}^{1/3}exp\left\{6{\gamma }_1\left(\frac{1}{{\rho }_m^{1/3}}-\frac{1}{{\rho }^{1/3}}\right)+\frac{2{\gamma }_3}{q_3}\left(\frac{1}{{\rho }_m^{q_3}}-\frac{1}{{\rho }^{q_3}}\right)\right\},$$

(14)

where ${\rho }_0$ and ${\rho }_m$ are the corresponding initial densities, usually taken as those at ambient P.

The Value of ${\gamma }_1$: Ultrahigh Density Limit

The value of ${\gamma }_1$ is found in the limit of ultrahigh densities. In this limit, it follows from Thomas-Fermi theory that $P/Z^{10/3},$ $B/Z^{10/3},$ and $G/Z^{10/3}$ are functions of $ZV,$ where Z is the atomic number and V is the molar volume [30]. The theoretical limiting forms of $P$, B and G as $V\to 0$ imply that $P/Z^{10/3},B/Z^{10/3}\sim {\left(ZV\right)}^{-5/3}exp\{-a{\left(ZV\right)}^{1/3}\}$ (P and B have similar $ZV$ dependences since $B\equiv -VdP/dV$), but $G/Z^{10/3}\sim {\left(ZV\right)}^{-4/3}exp\{-a{\left(ZV\right)}^{1/3}\}.$ The limiting form of P includes an exponential screening correction to the EOS of a free electron gas [2], and this correction also enters the limiting form of G in view of $G\sim (B-2/3tP),$ where $t\to 5/2$ as $V\to 0$ [2]. The exact expressions for P and G in the $V\to 0$ limit are [31,32]

$$P\left(V\right)=1003.6Z^{5/3}{V}^{-5/3}exp\{-a{\left(ZV\right)}^{1/3}\}\mathrm{GPa}·{\mathrm{cm}}^5,$$

(15)

$$G\left(V\right)=190Z^2{V}^{-4/3}exp\{-a{\left(ZV\right)}^{1/3}\}\mathrm{GPa}·{\mathrm{cm}}^4.$$

(16)

The corresponding limiting form of $T_m$ is $({\mathsf{\Gamma }}_m=175)$ [2]

$$T_m\left(V\right)={\left(\frac{4\pi }{3}\right)}^{1/3}\frac{Z^2{e}^2}{k_B{\mathsf{\Gamma }}_m}{V}^{-1/3}exp\{-b{\left(ZV\right)}^{1/3}\}=1300Z^2{V}^{-1/3}exp\{-b{\left(ZV\right)}^{1/3}\}\mathrm{K}·\mathrm{cm}.$$

(17)

Below we show that $b\approx a.$

We now convert the above functions of V into functions of $\rho$ to directly compare with Equations (13) and (14) and extract the value of ${\gamma }_1.$ This is done by using $V=M/\rho ,$ where M is the molar mass, so that Equations (16) and (17) reduce to

$$G\left(\rho \right)=190\frac{Z^2}{M^{4/3}}{\rho }^{4/3}exp\left\{-a\frac{Z^{1/3}{M}^{1/3}}{{\rho }^{1/3}}\right\}\mathrm{GPa}·{\mathrm{cm}}^4$$

(18)

and

$$T_m\left(\rho \right)=1300\frac{Z^2}{M^{1/3}}{\rho }^{1/3}exp\left\{-b\frac{Z^{1/3}{M}^{1/3}}{{\rho }^{1/3}}\right\}\mathrm{K}·\mathrm{cm}.$$

(19)

Comparison with Equations (13) and (14) shows that the following two conditions must be fulfilled: ${\gamma }_{1,G}=a{Z}^{1/3}{M}^{1/3}/6$ and ${\gamma }_{1,T_m}=b{Z}^{1/3}{M}^{1/3}/6.$ Fitting the functional form ${\left(ZV\right)}^{-5/3}exp\{-a{\left(ZV\right)}^{1/3}\}$ to the three sets of data on $P/Z^{10/3},$ $B/Z^{10/3}$ and $G/Z^{10/3}$ [30] leads to $a\approx 0.8$ cm${}^{-1}$.

Now consider the ratio $A/Z$ as a function of Z; here $A=M$ g${}^{-1}$ is the atomic mass. Its values for the first six elements (H through C) are 1.0, 2.0, 2.3, 2.25, 2.16, and 2.0. Starting with N, this ratio is a monotonically increasing function of Z with a nearly constant slope, such that $A/Z\approx 2$ for N and O and reaches ∼2.6 by the end of the table [33]; hence $A/Z=2.3\pm 0.3.$ We expect $A/Z=2.3$ to be a reliable approximation for the vast majority of elements, specifically for those in the middle of the table (Mo with Z = 42 is one of them).

With $a=0.8$ cm${}^{-1}$ and $A=2.3Z,$ we get

$${\gamma }_{1,G}=\frac{0.8\times {2.3}^{1/3}}{6}{Z}^{2/3}\approx 0.176Z^{2/3}$$

(20)

in units of (g/cc)${}^{1/3}$. Our independent analysis of the numerical sets of parameters for 32 melting curves (of 28 elemental solids and 4 alloys) [29] revealed that fitting $Z^{2/3}$ to the set of 32 values of the corresponding ${\gamma }_{1,T_m}$ leads to

$${\gamma }_{1,T_m}=\frac{4}{23}{Z}^{2/3}\approx 0.174Z^{2/3}$$

(21)

in units of (g/cc)${}^{1/3}$. Equations (20) and (21) show that $b\approx a$ in Equations (16) and (17) and that

$${\gamma }_{1,G}\approx {\gamma }_{1,T_m}\equiv {\gamma }_1=0.175Z^{2/3}=\frac{7}{40}{Z}^{2/3}.$$

(22)

Relation (22) is expected to be valid for all the elements except perhaps (i) both low-Z and high-Z elements; (ii) the alkali and alkaline-earth metals which are not described by the above equation of state with an exponential screening correction in the Thomas-Fermi regime [2]; and (iii) rare-gas elements for which the corresponding values of ${\gamma }_{1,T_m}$s are different from those predicted by $\frac{4}{23}{Z}^{2/3}$ [29].

By equating Equations (13) and (14) to (18) and (19), respectively, and taking the $P\to \infty$ limit we obtain the following two systems of equations for the two sets of parameters:

$$\frac{6{\gamma }_1}{{\rho }_0^{1/3}}+\frac{2{\gamma }_2}{q_2{\rho }_0^{q_2}}=ln\frac{190Z^2{\rho }_0^{4/3}}{A^{4/3}G\left({\rho }_0\right)},$$

(23)

$$\frac{1}{2}+\frac{{\gamma }_1}{{\rho }_0^{1/3}}+\frac{{\gamma }_2}{{\rho }_0^{q_2}}={\gamma }_G\left({\rho }_0\right),$$

(24)

and

$$\frac{6{\gamma }_1}{{\rho }_m^{1/3}}+\frac{2{\gamma }_3}{q_3{\rho }_m^{q_3}}=ln\frac{1300Z^2{\rho }_m^{1/3}}{A^{1/3}{T}_m\left({\rho }_m\right)},$$

(25)

$$\frac{1}{2}+\frac{{\gamma }_1}{{\rho }_m^{1/3}}+\frac{{\gamma }_3}{{\rho }_m^{q_3}}={\gamma }_{T_m}\left({\rho }_m\right),$$

(26)

where G is in GPa, $T_m$ in degrees K, and both ${\rho }_0$ and ${\rho }_m$ are in g/cc. With ${\gamma }_1=\frac{7}{40}{Z}^{2/3},$ the solution to this system of equations is

$$q_2=2\frac{{\gamma }_G\left({\rho }_0\right)-\frac{1}{2}-\frac{7}{40}\frac{Z^{2/3}}{{\rho }_0^{1/3}}}{ln\frac{190Z^2{\rho }_0^{4/3}}{A^{4/3}G\left({\rho }_0\right)}-\frac{21}{20}\frac{Z^{2/3}}{{\rho }_0^{1/3}}},$$

(27)

$$q_3=2\frac{{\gamma }_{T_m}\left({\rho }_m\right)-\frac{1}{2}-\frac{7}{40}\frac{Z^{2/3}}{{\rho }_m^{1/3}}}{ln\frac{1300Z^2{\rho }_m^{1/3}}{A^{1/3}{T}_m\left({\rho }_m\right)}-\frac{21}{20}\frac{Z^{2/3}}{{\rho }_m^{1/3}}};$$

(28)

the corresponding values of ${\gamma }_2$ and ${\gamma }_3$ are then found from (24) and (26). The values of ${\gamma }_G\left({\rho }_0\right)$ and ${\gamma }_{T_m}\left({\rho }_m\right)$ are determined by Equations (8) and (9) with $\rho ={\rho }_0$ and $\rho ={\rho }_m,$ respectively:

$${\gamma }_G\left({\rho }_0\right)=\frac{1}{2}\left(\frac{dlnG\left(\rho \right)}{dln\rho }{|}_{\rho ={\rho }_0}-\frac{1}{3}\right)=\frac{1}{2}\frac{{\rho }_0}{G\left({\rho }_0\right)}\frac{dG\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_0}-\frac{1}{6},$$

(29)

$${\gamma }_{T_m}\left({\rho }_m\right)=\frac{1}{2}\left(\frac{dln{T}_m\left(\rho \right)}{dln\rho }{|}_{\rho ={\rho }_m}+\frac{2}{3}\right)=\frac{1}{2}\frac{{\rho }_m}{T_m\left({\rho }_m\right)}\frac{d{T}_m\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_m}+\frac{1}{3};$$

(30)

both $\frac{dG\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_0}$ and $\frac{d{T}_m\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_m}$ come from experimental and/or theoretical data.

4. Generalized Melt-Shear Model for Molybdenum

In the remainder of this paper we apply the generalization of the unified melt-shear model developed in the previous section to molybdenum, and compare it to the available experimental and theoretical results.

We first determine the values of the two sets of parameters for G and $T_m.$

Model Parameters for Mo

The values of ${\gamma }_1,$ ${\gamma }_2,q_2\mathrm{and}{\gamma }_3,q_3$ are calculated using the following input from the literature: $Z=42,$ $A=95.95,$ ${\rho }_0=10.25{\mathrm{g}/\mathrm{cm}}^3,$ ${\rho }_m=9.56{\mathrm{g}/\mathrm{cm}}^3,$ $G\left({\rho }_0\right)=128.2\mathrm{GPa},$ and $T_m\left({\rho }_m\right)=2896\mathrm{K}.$

Our ab initio study of the thermoelasticity of Mo described below resulted in a $T=0$ shear modulus for bcc-Mo that is described by the following equation (G in GPa, $\rho$ in g/cm${}^3$):

$$G(\rho ,0)=113.0+33.5(\rho -10.25)+4.65{(\rho -10.25)}^2-0.68{(\rho -10.25)}^3+0.018{(\rho -10.25)}^4.$$

(31)

Hence $\frac{dG\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_0}=33.5,$ and Equation (24) gives

$${\gamma }_G\left({\rho }_0\right)=1.173.$$

(32)

To determine the corresponding value of ${\gamma }_{T_m}\left({\rho }_m\right)$ we use the theoretical bcc-Mo melting curve of Alfe et al. [34] based on ab initio quantum molecular dynamics (QMD) simulations. The experimental melting curve of Errandonea et al. [35,36] measured using laser-heated diamond anvil cells (DAC) is highly controversial: it is “flat” and lies considerably below the shock-melting datapoint at ∼(400 GPa, 10,000 K), now firmly established, as discussed above (Figure 2 and Figure 3). The very recent experimental study by Hrubiak et al. [37] demonstrates that, on increasing $T,$ compressed Mo undergoes a transformation that results in a texture (microstructure) change: large Mo grains become unstable at high T due to high atom mobility and reorganize into smaller crystalline grains. This transformation occurs below melting, and the pressure dependence of the transformation temperature is consistent with the previous DAC melting curve by Errandonea et al. [35,36]. On the other hand, the experimental melting curve of Mo determined in [37] is consistent with the theoretical melting curve of Alfe et al. [34]. Hence, most likely, Errandonea’s curve is an intermediate DAC transition boundary, while Alfe’s curve is the actual melting curve of Mo. Therefore, it is indeed Alfe’s theoretical melting curve that should be used for the determination of the generalized melt-shear model parameters for Mo.

The theoretical data of Alfe et al. [34] are consistent with $T_m=2896$ K at ${\rho }_m=9.933$ g/cm${}^3$. The density shift of 3.9%, compared to the experimental value of ${\rho }_m=9.56$ g/cm${}^3,$ comes from the specific implementation of density functional theory (DFT) in their calculations. For the same reason, a density shift of similar magnitude is present in our calculations as well: as follows from (31), the experimental $G=128.2$ GPa corresponds to ${\rho }_0=10.67$ g/cm${}^3$ which is 4.1% higher than the experimental ${\rho }_0=10.25$g/cm${}^3.$ The same density shift is responsible for a small offset of the theoretical curve relative to the experimental points of Reference [38], as seen in Figure 4.

The melting data of Alfe et al. [34] are described by the following equation:

$$T_m\left(\rho \right)=2896+695(\rho -9.933)-5.5{(\rho -9.933)}^2+21{(\rho -9.933)}^3-1.4{(\rho -9.933)}^4.$$

(33)

with the density shift taken into account, this curve should model the actual melting curve of Mo with $\frac{d{T}_m\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_m}=695;$ hence, via (30),

$${\gamma }_{T_m}\left({\rho }_m\right)=1.48.$$

(34)

we note that $\frac{d{T}_m\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_m}\sim 700$ also follows from an alternative approach: $\frac{d{T}_m\left(\rho \right)}{d\rho }=\frac{B_m}{{\rho }_m}\frac{d{T}_m}{dP}{|}_{P=0},$ where $B_m\equiv B({\rho }_m,T_m\left({\rho }_m\right))$ is the bulk modulus at the ambient melting point. There are no data on $B_m$ for solid Mo, but there are data on $B_m$ for liquid Mo, and we assume that the $B_m$s are close to each other so that $B_m$ for liquid Mo can be used in the above formula. Its value can be found from the data on the density dependence of the sound velocity in liquid Mo: $C=C\left(\rho \right)=0.8573\rho -3.3241$ [39]. Using $B_m={\rho }_m{C}^2$ with ${\rho }_m=9.33$ g/cm${}^3$ for liquid Mo [39], we obtain $B_m=201$ GPa, and then, with $d{T}_m/dP{|}_{P=0}=34\pm 6$ K/GPa [40], we get $\frac{d{T}_m\left(\rho \right)}{d\rho }{|}_{\rho ={\rho }_m}=715\pm 26$ K/(g cm${}^{-3}$).

With ${\gamma }_1=\frac{7}{40}{42}^{2/3}=2.1145,$ the solution to (27) and (28) is

$$q_2=0.63,{\gamma }_2=-1.31,\mathrm{and}{q}_3=0.45,{\gamma }_3=-0.045.$$

(35)

5. Comparison to Data

In this section we compare our generalized melt-shear model for Mo to (i) our ab initio calculations of the $T=0$ shear moduli of both bcc-Mo and fcc-Mo and (ii) our ab initio melting datapoint on fcc-Mo, as well as the theoretical melting data on bcc-Mo from the previous section. We also incorporate our generalized analytic model for $T_m\left(\rho \right)$ in a model for $G\left(P\right)$ along the principal Hugoniot and find excellent agreement with data.

5.1. Cold Shear Modulus of bcc-Mo and fcc-Mo from VASP

In the following we present the results of our calculations of the $T=0$ shear moduli of both bcc-Mo and fcc-Mo, as well as the finite-T shear modulus of bcc-Mo using the DFT code VASP (Vienna Ab initio Simulation Package). We use the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional. We model Mo using two core-valence representations, namely [${}^{30}$Zn$]4p^{6}4d^{5}5s^1$ and [${}^{18}$Ar$3d^{10}]4s^{2}4p^{6}4d^{5}5s^1$, i.e., we assign, respectively, the 12 or 14 outermost electrons of Mo to the valence. The valence electrons are represented with a plane-wave basis set with a cutoff energy of 400 eV for the 12-electron representation, and 550 eV for the 14-electron one, while the core electrons are represented by projector augmented-wave (PAW) pseudopotentials. The core radii (the largest values of RCUTs among those for each of the quantum orbitals) of the two pseudopotentials are, respectively, 2.5 a.u., or 1.323 Å and 2.3 a.u., or 1.217 Å. Numerical errors in the calculations using VASP will remain almost negligible until the nearest neighbor distance reaches 2 × RCUT$/(1.25\pm 0.05)$ [41]. Hence, with these pseudopotentials, one can study systems with densities up to ∼35 g/cm${}^3$ and ∼45 g/cm${}^3.$ At higher densities one should expect the overlap of interatomic spheres starting to develop. But since this overlap is generally allowed and handled well by VASP, in contrast to some other electronic structure codes that give huge errors when the atomic spheres overlap considerably [41], numerical errors in the VASP calculations should remain small at densities much higher than the above two values, specifically at $T=0$ when atoms remain in the equilibrium positions. Hence, we expect our calculations of the cold elastic constants to be fairly accurate to the highest density considered.

We used a $5\times 5\times 5$ (250-atom) supercell for bcc-Mo and a $4\times 4\times 4$ (256-atom) supercell for fcc-Mo with a $5\times 5\times 5$ k-point mesh in both cases. Full energy convergence (to ≲0.5 meV/atom) was checked in each simulation. We followed the scheme for the calculation of the shear modulus suggested by Hill, namely, averaging the Voigt $\left(V\right)$ and Reuss $\left(R\right)$ upper and lower bounds on the isotropic shear modulus, which for a cubic system are

$$G_V=\frac{2C^{\prime }+3C_{44}}{5},G_R=5{\left(\frac{2}{C^{\prime }}+\frac{3}{C_{44}}\right)}^{-1}.$$

(36)

The Voigt-Reuss-Hill average shear modulus obtained in this way is

$$G_{VRH}=\frac{G_V+G_R}{2}=\frac{3C^{\prime }{}^2+3C_{44}^2+19C^{\prime }{C}_{44}}{15C^{\prime }+10C_{44}}.$$

(37)

Here $C^{\prime }\equiv (C_{11}-C_{12})/2$ and $C_{44}$ are the single-crystal elastic constants. A methodology for calculating $C^{\prime }$ and $C_{44}$ was developed by Söderlind et al. [42]. To evaluate $C^{\prime },$ the cubic lattice is deformed by the (volume conserving) transformation

$$\left(\begin{array}{ccc}1+\delta & 0& 0\\ 0& 1+\delta & 0\\ 0& 0& \frac{1}{{(1+\delta )}^2}\end{array}\right).$$

The resulting energy change is

$$\frac{\Delta E}{V_{\mathrm{eq}}}=6C^{\prime }{\delta }^2+O\left({\delta }^3\right),$$

(38)

where $V_{\mathrm{eq}}$ is the equilibrium volume of the system. Similarly, $C_{44}$ is obtained by applying the (volume conserving) transformation

$$\left(\begin{array}{ccc}1& \delta & 0\\ \delta & 1& 0\\ 0& 0& \frac{1}{1-{\delta }^2}\end{array}\right),$$

resulting in an energy change

$$\frac{\Delta E}{V_{\mathrm{eq}}}=2C_{44}{\delta }^2+O\left({\delta }^4\right).$$

(39)

In our calculations, we have evaluated the energy as a function of $\delta$ at intervals of 0.01 up to $\delta =0.05$ for the $C_{44}$ case. For the $C^{\prime }$ case, the energy is not an even function of $\delta ,$ and so negative values of $\delta$ were used as well $(-0.03\le \delta \le 0.03).$ The resulting $E\left(\delta \right)$ were fit to third- and fourth-degree polynomials for $C^{\prime }$ and $C_{44}$, respectively, and the quadratic coefficient was used to evaluate the elastic constant from Equation (38) or Equation (39).

In Figure 5 our multi-phase analytic model for $G(\rho ,0)$ is compared with our ab initio results for bcc-Mo and fcc-Mo.

5.2. bcc-fcc Phase Transition in Mo

Figure 4 and Figure 5 show that with increasing density the cold shear modulus of bcc-Mo first increases but then decreases and eventually turns negative at $\rho \gtrsim 30$ g/cm${}^3.$ This behavior is due to a dynamic instablity which sets in at ∼20 g/cm${}^3$ ($V\sim 8$ Å${}^3$/atom, or $P\sim 700$ GPa), and which is confirmed by the DFT-based phonon spectra shown in Figure 1. At the same time, at $T=0$ fcc-Mo dynamically stabilizes at ∼14.5 g/cm${}^3$ ($P\sim 200$ GPa). The $T=0$ free energy calculations demonstrate that fcc-Mo becomes thermodynamically more stable than bcc-Mo at $\rho \sim 20$ g/cm${}^3$ ($P\sim 700$ GPa); see Figure 6. The accurate calculation of the two enthalpies gives the following bcc-fcc transition parameters:

$$P_{\mathrm{bcc}-\mathrm{fcc}}=680.0\mathrm{GPa},\rho {\left(\mathrm{bcc}\right)}_{\mathrm{bcc}-\mathrm{fcc}}=20.033{\mathrm{g}/\mathrm{cm}}^3,\rho {\left(\mathrm{fcc}\right)}_{\mathrm{bcc}-\mathrm{fcc}}=20.404{\mathrm{g}/\mathrm{cm}}^3,$$

(40)

that is, a volume change of ∼2%. At lower P fcc-Mo dynamically stabilizes at nonzero T (the P-T stability regions of fcc-Mo are discussed in [17]) but, as we noted above, full free energy calculations rule out the possibility of fcc-Mo being thermodynamically more stable than bcc-Mo below 700 GPa.

5.3. Melting Curve of Mo

Figure 7 compares our generalized melt-shear model for Mo to the available theoretical data on the melting of Mo [34,43,44], including our own ab initio QMD result for the melting temperature of fcc-Mo at a density of 23.6 g/cm${}^3$. Our simulations were carried out using the Z method implemented with VASP which is described in detail in Reference [45]. We used a 500-atom ($5\times 5\times 5)$ fcc-Mo supercell. The simulations were carried out with a single $\mathsf{\Gamma }$-point (full energy convergence, to ≲1 meV/atom, was achieved with such a large supercell), having 14 electrons of Mo in the valence. We chose to generate a melting datapoint for fcc-Mo to check that the unified model describes this phase as well as it describes bcc-Mo on which the model is based. Since bcc-Mo is dynamically unstable at a density of 23.6 g/cm${}^3$ (a pressure of ∼1300 GPa), it is reasonable to assume that Mo melts from fcc at this density.

We carried out seventeen QMD runs on fcc-Mo at $\rho =23.6$ g/cm${}^3$ (a lattice constant of 3.0 Å) with initial temperatures from 42,500 K to 47,500 K separated by an increment of 312.5 K. Figure 8 and Figure 9 illustrate the time evolution of T and P in three of these seventeen QMD runs, with initial Ts of 44,687.5, 45,000, and 45,312.5 K. The 45,000 K run is the melting run [45] for which $(P_m,T_m)\sim (1300\mathrm{GPa},18,500\mathrm{K}).$

As Figure 7 clearly demonstrates, the analytic melting curve of our unified model is in excellent agreement with all the theoretical data on bcc-Mo, including our QMD datapoint on fcc-Mo.

5.4. Thermoelastic Softening Parameter

In the Preston-Wallace (PW) thermoelasticity model [6] the thermoelastic softening parameter $\beta$ can in principle be a function of density, so the most general analytic form of the PW model is

$$G(\rho ,T)=G(\rho ,0)\left(1-\beta \frac{T}{T_m\left(\rho \right)}\right),$$

(41)

instead of Equation (4). It was determined in [6] that for Mo at $P=0,$ $\beta =0.23.$ In order to ascertain if there is any density dependence of $\beta ,$ we carried out QMD simulations of elastic constants at a set of four nonzero Ts, specifically, ∼$T_m/4,$ $T_m/2,$ $3T_m/4$ and $T_m,$ at three different densities. The slope of a straight line fitted to the four nonzero T values combined with the $T=0$ value calculated earlier gives the value of $\beta$ at the corresponding density. We followed the same strategy as for the $T=0$ simulations described above. We employed a Nosé-Hoover thermostat and used the Nosé algorithm with controlled temperature oscillations around the initial T so that their frequency is the same as the characteristic phonon frequency of the atomic system at the P-T conditions of the QMD run. The finite-T simulations typically require between 5000 and 10,000 time steps of 1 fs to achieve full energy convergence and to produce sufficiently long output for the extraction of reliable averages for the values of energy that are used in Equations (38) and (39) to calculate the corresponding values of the elastic constants.

Our QMD simulations of G as a function of T at three different densities, shown in Figure 10, demonstrate that, at least for bcc-Mo, $\beta$ does depend on density. The best fit to the three datapoints subject to the constraint that $\beta (\rho =10.25)=0.23$ ($\rho$ in g/cm${}^3$) is

$$\beta =\beta \left(\rho \right)=0.23+0.0212{(\rho -10.25)}^{1.11}.$$

(42)

6. Thermoelasticity Model of Molybdenum

The unified analytic melt-shear model for a multi-phase material developed in this paper gives the cold shear modulus and melting temperature at all densities across all solid phases. For example, in Figure 5 the model curve covers bcc-Mo at densities up to ∼20 g/cc and smoothly transitions to fcc-Mo at higher $\rho$. In order to have a model of the thermoelasticity of any given solid phase of a multi-phase material, the corresponding unified model must be supplemented with the information on (i) the density interval in which the solid phase is dynamically stable; and (ii) the density dependence of the thermoelastic softening parameter over this density interval. In the following, we consider bcc-Mo as an example of such a complete thermoelasticity model.

Equation (41) with $G(\rho ,0)$ and $T_m\left(\rho \right)$ given by Equations (13) and (14), respectively, and $\beta \left(\rho \right)$ as in Equation (42) define the complete thermoelasticity model of bcc-Mo over the range of densities in which it is dynamically stable. To illustrate its fidelity, we now consider shear modulus along the principal Hugoniot of Mo. This is one of the applications of our model that can be directly compared to experiment.

Shear Modulus and Sound Velocities Along the Principal Hugoniot

Using the cold G from Equation (31) in the PW elasticity model, along with (i) $T_m\left(\rho \right)$ from the new analytic model; (ii) $\beta =\beta \left(\rho \right)$ from Equation (42); and (iii) the temperature along the Hugoniot gives the shear modulus along the Hugoniot.

We found two different sources of data on $T\left(\rho \right)$ along the Hugoniot: [46,47]. The best fits to the two sets of data are, respectively,

$$T\left(\rho \right)=293+58.98{(\rho -10.206)}^{2.78586},$$

(43)

and

$$T\left(\rho \right)=293+52.2891{(\rho -10.215)}^{2.87282}.$$

(44)

The intersection of the melting curve (14) with the parameter values from (35) with either of the above $T\left(\rho \right)$ along the Hugoniot predicts the corresponding Hugoniot melting point, and it turns out to be $(\rho ,T)=(16.19,8900)$ and $(16.10,8800),$ respectively, in good agreement with the theoretical estimate $(16.5,9300)$ by Hixson et al. [7] based on their experimental results. Here we take Equation (43) as that for $T\left(\rho \right)$ along the principal Hugoniot of Mo.

We use the experimental $P\left(\rho \right)$ (those of [46,47] are virtually identical) to calculate G as a function of pressure along the bcc-Mo Hugoniot; our results are shown in Figure 2 and are in good agreement with experiment. The experimental data are consistent with bcc-Mo along the entire Hugoniot; there is no evidence for any other crystal structure such as Phase II [24]. Our $G=G\left(P\right)$ model curve can be approximated by

$$G\left(P\right)=122.0+0.94P-3.0\times {10}^{-3}{P}^2+1.2\times {10}^{-5}{P}^3-2.8\times {10}^{-8}{P}^4-2.4\times {10}^{-11}{P}^5+8.3\times {10}^{-14}{P}^6.$$

(45)

We also calculated sound velocities $C_L$ and $C_B$ along the Hugoniot; they are shown in Figure 3. We used $B=\rho dP/d\rho$ calculated from the experimental $P\left(\rho \right).$

Our results on the shear modulus and sound velocities in bcc-Mo along its principal Hugoniot are compared to experimental data in Figure 2 and Figure 3. The shear modulus of bcc-Mo increases with P up to ∼240 GPa but then decreases. This P dependence is not due to a high-P phase such as Phase II proposed by Errandonea et al. [24] for its explanation but is rather a property of bcc-Mo.

In Appendix B we consider two additional examples of transition metals with ambient bcc structure that dynamically destabilizes with increasing pressure, namely, tantalum and tungsten. We find the same behavior of G along their Hugoniots as that of Mo. These two additional examples strongly suggest that this behavior may be typical for transition metals with ambient bcc structure that becomes dynamically unstable at higher P. The parameters for the thermoelasticity models of Mo, Ta, and W are given in

Table 1. Numerical values of the parameters for the thermoelasticity model, Equation (41) (in the corresponding units). In Equation (41) $G(\rho ,0)$ and $T_m\left(\rho \right)$ are given by Equations (13) and (14), respectively. The last three columns of the table apply to the bcc phases of the three transition metals for which $\beta =\beta \left(\rho \right)=\beta \left({\rho }_0\right)+a{(\rho -{\rho }_0)}^b.$

	Z	${\mathit{\gamma }}_1$	${\mathit{\gamma }}_2$	${\mathit{q}}_2$	${\mathit{\gamma }}_3$	${\mathit{q}}_3$	${\mathit{\rho }}_0$	${\mathit{\rho }}_{\mathit{m}}$	$\mathit{G}({\mathit{\rho }}_0,0)$	${\mathit{T}}_{\mathit{m}}\left({\mathit{\rho }}_{\mathit{m}}\right)$	$\mathit{\beta }\left({\mathit{\rho }}_0\right)$	a	b
Mo	42	2.11450	0.63	−1.31	0.45	−0.045	10.25	9.56	128.2	2896	0.23	0.0212	1.11
Ta	73	3.05670	1.0	−6.6	3.0	−3000	16.74	15.30	72.2	3290	0.23	0.0210	1.09
W	74	3.08456	0.84	−5.8	1.5	−11.1	19.31	17.96	163.4	3695	0.23	0.0211	1.07

.

7. The Phase Diagram of Molybdenum

The results obtained in this work suggest that the phase diagram of Mo cannot be extended beyond Figure 1 of Reference [11]. Indeed, if there are other solid phase fields on the phase diagram of Mo, in addition to bcc which is well established experimentally, and perhaps some other solid structure (fcc, hcp, dhcp, etc.) which is so far only predicted theoretically, neither Hugoniot nor ramp-wave compression curves cross into them. Hence, further refinement of the phase diagram of Mo is not yet feasible. However, we have confirmed that the topology of the phase diagram of Mo is essentially the same as in Figure 1 of Reference [11] and is represented in Figure 11. To this end, we carried out two sets of independent inverse Z runs [48] to solidify liquid Mo and to check whether there is any solid-solid phase boundary so that liquid Mo solidifies into bcc on one side of this boundary and into another solid structure on the other side. We used a computational cell of 512 atoms prepared by melting a $8\times 8\times 8$ solid simple cubic (sc) supercell which would eliminate any bias towards solidification into bcc or any other solid structure (fcc, hcp, dhcp, etc.). We used sc unit cells of 2.225 and 1.985 Å; the dimensions of bcc unit cells having the same volume as the sc ones are 2.80 and 2.50 Å, respectively, which corresponds to ∼200 and 700 GPa.

We carried out $NVT$ simulations using the Nosé-Hoover thermostat with a timestep of 1 fs. Complete solidification typically required from 15 to 25 ps, or 15,000–25,000 timesteps. The inverse Z runs indicate that liquid Mo only solidifies into bcc at ∼200 GPa in the whole temperature range from 0 to essentially the corresponding $T_m$. However, at ∼700 GPa it solidifies into bcc above the transition boundary in Figure 11, while below this boundary it solidifies into another solid structure. The radial distribution functions (RDFs) of the final solid states are noisy; upon fast quenching of the three structures to low $T,$ where RDFs are more discriminating, and by comparing them to the RDFs of fcc, hcp and dhcp, we conclude that dhcp is the closest strucure to those that liquid Mo solidies into below the transition boundary.

The RDFs of the solidified states at ∼700 GPa above the transition boundary are shown in Figure 12, and of those solidified below the transition boundary in Figure 13. The 2500 K state virtually lies on the boundary. We tentatively assign it to bcc, because it definitely has features of bcc (RDF peaks at $R\sim 55$, 65 and 90, small peak at $R\sim 80$, etc.). At the same time it certainly has some features that are both uncharacteristic of bcc (e.g., the disappearance of the bcc split-peak at $R\sim 40$) and characteristic of dhcp (a small peak at $R\sim 50$, etc.). The RDF of this 2500 K state is shown in Figure 13 also, for comparison to other dhcp states, including the real dhcp-Mo (at 2500 K). Most likely, this 2500 K state is a mixture of both bcc and dhcp.

A few more comments are in order. The 8750 K state at ∼700 GPa did not solidify, most likely for the reason of not being supercooled enough to intiate the solidification process [48]. Indeed, 8750 K constitutes ∼0.8 of the corresponding $T_m$ (20% of supercooling) while for the other set of points at ∼200 GPa, the highest solidification T of 5000 K constitutes ∼0.75 of the corresponding $T_m$ (25% of supercooling) which apparently allows for the solidification process to go through in this case. The 625 K state at ∼700 GPa is amorphous, with a characteristic RDF shoulder at $R\sim$ 60–70, for the reason of the corresponding solidification T being too low, most likely below the corresponding T of vitrification.

8. Concluding Remarks

We have generalized the unified analytic melt-shear model that we introduced a decade ago to multi-phase materials. We have described a new scheme for calculating the values of the model parameters for both G and $T_m$ at all densities. The new model was applied to molybdenum, a multi-phase material having body-centered cubic (bcc) structure at low $\rho$ and face-centered cubic (fcc) structure at high $\rho .$ The new model was compared to (i) our ab initio results on the $T=0$ shear modulus of both bcc-Mo and fcc-Mo as a function of $\rho$; and (ii) the available theoretical results on the melting of bcc-Mo and our own quantum molecular dynamics simulations of one high-P melting point of fcc-Mo.

We have constructed a model for $G(\rho ,T)$ of bcc-Mo at all $\rho$ for which bcc-Mo remains dynamically stable, and all T from 0 to $T_m.$ The model is described by Equation (41) in which $G(\rho ,0)$ and $\beta \left(\rho \right)$ are given by Equations (31) and (42), respectively, and $T_m\left(\rho \right)$ is given by Equation (14) with ${\gamma }_1=2.1145$ and ${\gamma }_3$ and $q_3$ from Equation (35). We then applied this thermoelasticity model to the calculation of the shear modulus of bcc-Mo along its principal Hugoniot. Our results were compared to experimental data and showed that G of bcc-Mo increases with P up to ∼240 GPa but then decreases. We emphasize that this P dependence is not due to a high-P phase such as Phase II proposed by Errandonea et al. [24] for its explanation. Similar behavior of G along the Hugoniot is also found in the transition metals Ta and W, thereby suggesting that it may be typical for transition metals with ambient bcc structure that become dynamically unstable at higher $P.$ Thus, we agree with the common conclusion of a number of recent papers [8,9,26] that no s-s phase transition can be definitively inferred from the sound velocity data on Mo. If there are other solid phase fields on the phase diagram of Mo, in addition to bcc, neither Hugoniot nor ramp-wave compression curves cross into them. In the absence of additional experimental and/or theoretical information on other solid phases, our current knowledge of the phase diagram of Mo cannot be extended beyond Figure 1 of Reference [11].

Finally, our thermoelasticity model for bcc-Mo predicts a decrease in G with P on the principal Hugoniot at pressures above ∼204 GPa, hence the decrease in G with P is not evidence for a s-s transition on the Hugoniot.