*Article* **Symmetry-Adapted Finite Strain Landau Theory Applied to KMnF3**

#### **Andreas Tröster 1,\* , Wilfried Schranz <sup>1</sup> , Sohaib Ehsan 2, Kamal Belbase <sup>2</sup> and Peter Blaha <sup>2</sup>**


Received: 27 January 2020; Accepted: 11 February 2020; Published: 17 February 2020

**Abstract:** In recent years, finite strain Landau theory has been gradually developed as both a conceptual as well as a quantitative framework to study high pressure phase transitions of the group-subgroup type. In the current paper, we introduce a new version of this approach which is based on symmetry-adapted finite strains. This results in a substantial simplification of the original formulation. Moreover, it allows for replacing the clumsy use of truncated Taylor expansions by a convenient functional parametrization. Both the weaknesses of the traditional Landau approach based on infinitesimal strains as well as the major improvements made possible by our new parametrization are illustrated in great detail in an application to the ambient temperature high pressure transition of the perovskite KMnF3.

**Keywords:** high pressure; phase transitions; Landau theory; nonlinear elasticity theory; perovskites

#### **1. Introduction**

Through many decades, the Landau theory (LT) of phase transitions [1,2] (PTs) has proven to be one of the most valuable conceptual tools for understanding PTs of the group-subgroup type. In particular, the field of structural PTs abounds even with successful quantitative applications, and, for many classes of materials, complete collections of the corresponding coupling coefficients have been gathered in the literature (for ferroelectrics see, e.g., Appendix A of Ref. [3]). Effects of spontaneous strain that generally accompany temperature-driven structural PTs are sufficiently parameterized in terms of infinitesimal strain tensor components defined with respect to a baseline, which is obtained by extrapolating the generally small thermal expansion changes of the high symmetry reference phase. The corresponding Landau potential then involves only terms up to harmonic order, and any temperature dependence of the relevant parameters (high symmetry phase elastic constants and other coupling constants) can usually be completely neglected [4].

The situation changes drastically for high pressure phase transitions (HPPTs). Spontaneous strain components may still be numerically small, but they now must be defined with respect to a *P*-dependent base line. The total strain measured with respect to the ambient pressure reference state must then be calculated from a nonlinear superposition of finite background and spontaneous strain (see Equation (11) below). Furthermore, the Landau potential may be truncated beyond harmonic order only if calculated with respect to this *P*-dependent elastic background reference system. Therefore, neither the elastic constants nor the other elastic couplings can be assumed to be *P*-independent. In a high pressure context, clinging to the familiar infinitesimal strain Landau toolbox may result not only in a quantitatively, but also qualitatively completely erroneous description.

*Crystals* **2020**, *10*, 124; doi:10.3390/cryst10020124 www.mdpi.com/journal/crystals

It is not easy to construct a mathematically consistent and yet practically useful version of Landau theory taking into account the inherent nonlinearities and anharmonic effects that accompany HPPTs. In recent years, however, such a theory, for which we have coined the name finite strain Landau theory (FSLT), has been successfully developed. FSLT constitutes a careful extension of Landau theory beyond coupling to infinitesimal strain, fully taking into account the nonlinear elastic effects at finite strain. Its capabilities have been demonstrated in a number of applications to HPPTs [5–10]. However, as it stands, the numerical scheme underlying FSLT is still quite involved, and many practical workers in the field of HPPTs may be hesitant to go through the mathematical hardships it seems to pose. It is the purpose of this paper to show that FSLT is drastically simplified by switching from a formulation in terms of Cartesian Lagrangian strains to one in terms of symmetry-adapted finite strains. The enormous reduction of overall complexity of the approach as well as the vastly reduced numerical requirements of our new version of FSLT are demonstrated on the example of the HPPT in the perovskite KMnF3 (KMF).

#### **2. Review of Experimental Results on the Cubic-to-Tetragonal Transition of KMF**

In what follows, we focus on the antiferrodistortive high pressure phase transition of KMF from the cubic perovskite aristophase *Pm*3¯*m* to a tetragonal *I*4/*mcm* phase at room temperature which was experimentally investigated in Ref. [11] by X-ray diffraction up to 30 GPa. Since the ambient pressure Landau theory also provides a limiting reference frame for the description of the high pressure transition, we start our discussion with a detailed survey of the corresponding Landau theory.

According to Ref. [12], a similar transition observed at ambient pressure and temperature *Tc* = 186.5 K is weakly first order but close to critical, and the mechanism underlying these transitions is the same as in the well-studied *T* = 105 K cubic-to-tetragonal transition in strontium titanate, i.e., octahedral tilting with a critical wavevector at the *R*-point of the cubic Brillouin zone. Furthermore, in Ref. [12], it is argued that, even though a further structural transition to an orthorhombic phase related to further octahedral tilting at the *M*-point of the Brillouin zone at *TN* = 87 K is accompanied by antiferromagnetism, this coincidence between structural and magnetic transition temperatures may just be accidental, and actually there seems to be essentially no coupling between structural and magnetic order parameters (OPs) for these transitions. In passing we note that at *T* = 82 K there is a further transition to an orthorhombic canted ferromagnet [12].

As discussed in Ref. [12], the *Pm*3¯*<sup>m</sup>* <sup>→</sup> *<sup>I</sup>*4/*mcm* symmetry reduction corresponds to the isotropy subgroup of the three-dimensional irreducible representation *R*<sup>+</sup> <sup>4</sup> of *Pm*3¯*<sup>m</sup>* with respect to the order parameter direction (*Q*, 0, 0), the corresponding Landau expansion to sixth order being

$$\begin{aligned} F &= & \frac{A}{2} \left( Q\_1^2 + Q\_2^2 + Q\_3^2 \right) + \frac{B}{4} \left( Q\_1^4 + Q\_2^4 + Q\_3^4 \right) + \frac{B\_2}{4} \left( Q\_1^2 Q\_2^2 + Q\_2^2 Q\_3^2 + Q\_1^2 Q\_3^2 \right) \\ &+ \frac{C}{6} \left( Q\_1^6 + Q\_2^6 + Q\_3^6 \right) + \frac{C\_2}{6} \left( Q\_1^2 Q\_2^4 + Q\_2^2 Q\_3^4 + Q\_1^2 Q\_3^4 + Q\_2^2 Q\_1^4 + Q\_3^2 Q\_2^2 + Q\_3^2 Q\_1^4 \right) + \frac{C\_3}{6} Q\_1^2 Q\_2^2 Q\_3^2 \\ &+ \lambda\_4^{(0)} (Q\_1^2 + Q\_2^2 + Q\_3^2) \epsilon\_d + \lambda\_1^{(0)} (2 Q\_1^2 - Q\_2^2 - Q\_3^2) \epsilon\_l + \frac{K^{(0)}}{2} \epsilon\_a^2 + \frac{\mu^{(0)}}{2} \epsilon\_l^2 \end{aligned} \tag{1}$$

with volume and tetragonal symmetry-adapted strains *ea* = *e*<sup>11</sup> + *e*<sup>22</sup> + *e*33, *et* = <sup>2</sup>*e*33−*<sup>e</sup>* <sup>√</sup>11−*e*<sup>22</sup> <sup>3</sup> and bulk and longitudinal shear modulus related to the bare cubic elastic constants by *K*(0) = (*C*(0) <sup>11</sup> <sup>+</sup> *<sup>C</sup>*(0) <sup>12</sup> )/3, *<sup>μ</sup>*(0) <sup>=</sup> (*C*(0) <sup>11</sup> <sup>−</sup> *<sup>C</sup>*(0) <sup>12</sup> )/2. Targeting a transition into a single tetragonal domain where (*Q*1, *Q*2, *Q*3) ≡ (*Q*, 0, 0) and <sup>22</sup> = 33, we have

$$\begin{array}{rcl} F &=& \frac{A}{2}Q^2 + \frac{B}{4}Q^4 + \frac{C}{6}Q^6 + \lambda\_d^{(0)}Q^2\varepsilon\_d + 2\lambda\_l^{(0)}Q^2\varepsilon\_l\\ &+ \frac{K^{(0)}}{2}\varepsilon\_d^2 + \frac{\mu^{(0)}}{2}\varepsilon\_l^2 \end{array} \tag{2}$$

In a standard Landau approach, the coefficients *B*, *C*, *λ*(0) *<sup>a</sup>* , *<sup>λ</sup>*(0) *<sup>t</sup>* are assumed to be independent of temperature (and pressure), while for *A* the ansatz

$$A = A\_0(T - T\_0) \tag{3}$$

with *T*-independent coefficients *A*0, *T*<sup>0</sup> is made and quantum saturation has been neglected [13]. The elastic equilibrium conditions <sup>−</sup>*<sup>P</sup>* <sup>=</sup> <sup>1</sup> *V*0 *∂F ∂a* ¯*a*,¯*t* and 0 = <sup>1</sup> *V*0 *∂F ∂t* ¯*a*,¯*t* amount to

$$\mathfrak{E}\_d = \mathfrak{e}\_d(P) - \frac{\lambda\_d^{(0)}}{K^{(0)}} Q^2, \qquad \mathfrak{E}\_l = -\frac{2\lambda\_l^{(0)}}{\mu^{(0)}} Q^2 \tag{4}$$

with a background volume strain

$$\varepsilon\_4(P) = \mathfrak{Re}(P), \quad \varepsilon(P) = -\frac{P}{\mathfrak{K}^{(0)}} \tag{5}$$

Performing a Legendre transform yields the Gibbs potential

$$G\_{\phantom{\alpha}} = \frac{A\_R}{2}Q^2 + \frac{\check{B}}{4}Q^4 + \frac{\mathbb{C}}{6}Q^6 - \frac{P^2}{2K^{(0)}}\tag{6}$$

where

$$A\_R = \ \ A - 6\lambda\_a^{(0)}\varepsilon(P) = A - \frac{2\lambda\_a^{(0)}P}{K^{(0)}} \tag{7a}$$

$$\widetilde{B}^{(0)} = -B - \frac{2(\lambda\_a^{(0)})^2}{K^{(0)}} - \frac{8(\lambda\_t^{(0)})^2}{\mu^{(0)}} \tag{7b}$$

The values of coefficients *A*<sup>0</sup> = 63.118 kPa/K, *B* = −1.308 MPa and *C* = 13.032 MPa at *P* = 0, which imply a first order phase transition at *Tc* = 186.15 K, have been determined from caloric measurements by Salje and coworkers [13]. In principle, numerical values for the OP-strain coupling coefficients *λ*(0) *<sup>a</sup>* , *<sup>λ</sup>*(0) *<sup>t</sup>* may be extracted from experimental data on the temperature evolution of spontaneous strains ¯*a*, ¯*t*. Unfortunately, however, for KMF, this is easier said than done. At room temperature, the thermal expansion data *a*cubic(*T*) of Ref. [14] are observed to perfectly reproduce the value of the cubic lattice constant *a*cubic(*TR*) at ambient pressure as determined in Ref. [11]. However, the measurement data of thermal lattice parameter irregularities *a*(*T*), *c*(*T*) around *Tc* ≈ 186 K (Refs. [15–18]) available in the literature appear to be in rather poor mutual agreement. As Figure 1 illustrates, while a discontinuous behavior of the lattice parameters is clearly visible in all three data sets, the absolute values of the measured unit cell parameters differ considerably, yet none of the data sets seem to be compatible with extrapolating the thermal expansion data of Ref. [14].

Not unexpectedly, the agreement in relative splitting between the a-and c-axis in the tetragonal phase appears to be better, albeit far from perfect. Nevertheless, the cubic parts of the data of Refs. [15,18] exhibit slopes similar to the low temperature extrapolation of the thermal expansion data of Ref. [14]. In order to be able to collapse the data onto a common "master set", we thus shifted the data of Refs. [15,18] by constant absolute offsets to match the extrapolated baseline of Ref. [14] in an optimal way, treating the seemingly more precise measurements of Sakashita et al. (Refs. [16,17]) as an outlier. Figure 2 illustrates our results for a corresponding effort.

**Figure 1.** Compilation of experimental unit cell data from Refs. [15] (red) and [18] (green) (data range restricted to *T* > 155 K) [16,17] (blue) and [14] (gray). For comparison, the value (including error bars) of the room temperature (indicated by the vertical dashed line) lattice constant at ambient pressure as measured in Ref. [11] is illustrated by the gray horizontal area.

**Figure 2.** Collapse of data from Refs. [15] (red) and [18] (green), onto the low temperature expansion of the thermal expansion data of Ref. [14] (gray) by constant vertical shifts. The positive and negative branches of the data correspond to values *<sup>t</sup>* and *<sup>a</sup>* from the various references, respectively. As in Figure 1, the room temperature ambient pressure lattice constant of Ref. [11] is indicated by a horizontal gray bar for comparison. The shifted data from Ref. [16,17] (blue) clearly appear to be at odds with the other measurements.

With a meaningful baseline *a*cubic(*T*) for unit cell parameters *a*(*T*), *c*(*T*) in place, we calculate the spontaneous strain components <sup>1</sup> = *a*/*a*cubic − 1 and <sup>3</sup> = *c*/*a*cubic − 1 and thus the (infinitesimal) spontaneous volume and tetragonal strains

$$
\varepsilon\_d \quad = \ 2\varepsilon\_1 + \varepsilon\_3 = \frac{2a+c}{a\_{cubic}} - 3 \,, \tag{8}
$$

$$
\epsilon\_l \quad = \quad \frac{2(\epsilon\_3 - \epsilon\_1)}{\sqrt{3}} = \frac{2}{\sqrt{3}} \frac{c - a}{a\_{cubic}},
\tag{9}
$$

respectively. According to Equation (4), when plotted against *Q*2(*t*) at *P* = 0, *<sup>a</sup>* and *<sup>t</sup>* should resemble straight lines with slopes <sup>−</sup>*λ*(0) *<sup>a</sup>* /*K*(0) and <sup>−</sup>2*λ*(0) *<sup>t</sup>* /*μ*(0), respectively. Figure <sup>3</sup> illustrates corresponding fits with results

*λ*(0) *<sup>a</sup>* /*K*(0) <sup>=</sup> 0.002, *<sup>λ</sup>*(0) *<sup>t</sup>* /*μ*(0) <sup>=</sup> <sup>−</sup>0.005 (10)

**Figure 3.** Fits of *a*(*T*) and *t*(*T*) against *Q*2(*T*) according to Equations (4).

With the Landau theory of the temperature-driven transition at ambient pressure available, we are tempted to analyze the ambient temperature HPPT based on the same framework. In Ref. [11], the variation of the pseudo-/cubic lattice constants of KMF under pressure at room temperature was measured with X-ray scattering (Figure 4) and the cubic part of the data was fitted to a simple Murnaghan equation of state (EOS) with *K*<sup>0</sup> = 64 GPa and *V*<sup>0</sup> = 73.608 3. This provides a baseline to determine (Lagrangian) spontaneous strains *<sup>a</sup>* <sup>≈</sup> *a*, *<sup>t</sup>* <sup>≈</sup> *<sup>t</sup>* (Figure 5).

Comparing the thermal and pressure-induced spontaneous strains shown in Figures 2 and 5, respectively, we note that there is some spontaneous thermal volume strain *<sup>a</sup>* while practically *<sup>a</sup>* <sup>≈</sup> 0 for the pressure-induced case. Furthermore, even though it may be difficult to directly relate temperature and pressure scales, the pressure-induced tetragonal spontaneous strain *<sup>t</sup>* is observed to be roughly one order of magnitude larger than its thermal counterpart *t*. From the perspective of traditional Landau theory based on infinitesimal strain coupling, these findings are difficult to explain. Based on Equation (4), there are in principle two ways to alter the magnitude of spontaneous strains. One may either change

the value of the couplings *λ*(0) *<sup>a</sup>* /*K*(0) and *<sup>λ</sup>*(0) *<sup>t</sup>* /*μ*(0) or find a mechanism to increase the magnitude of *<sup>Q</sup>*¯ <sup>2</sup> which, of course, implicitly depends on all Landau coefficients. Since these are actually only known for *<sup>T</sup>* <sup>=</sup> 186.5 K, one could assume a thermal drift in *<sup>λ</sup>*(0) towards zero to be responsible for the vanishing of *<sup>a</sup>* at room temperature (more than counteracting against the thermal drift of *K*<sup>0</sup> which is generally expected to decrease with increasing *T*). For the tetragonal strain, a similar mechanism seems to be difficult to conceive, however. On the one hand, we would need to increase *λ*(0) *<sup>t</sup>* dramatically to explain the large values of *t*. On the other hand, Equation (7b) indicates that such an increase would send parameter *<sup>B</sup>* to negative values much larger than those found for the thermal transition, resulting in a pronounced first order character of the HTTP. This, however, is not observed. The only remedy therefore seems to find a way to increase *Q*¯ 2. Calculated from a standard 2–4 Landau potential, *Q*¯ <sup>2</sup> would be inversely proportional to 1/*B*. This may explain why advocates of an orthodox Landau description frequently resort to assuming HPPT's to be near a tricritical point, explicitly postulating some ad-hoc pressure dependence *B* = *B*(*P*) induced somehow by unspecified higher order coupling effects.

**Figure 4.** Pseudo-/cubic lattice constants of KMF under pressure at room temperature as measured in Ref. [11] with X-ray scattering. The transition pressure estimate *Ps* ≈ 3.4 GPa put forward in Ref. [11], which is indicated by the dashed vertical line, is clearly too low.

Further difficulties arise if we try to reconcile the observed value of *Pc* with the prediction of standard Landau theory. In Ref. [11], a brute force fit based on the above assumptions of a second order transition close to a tricitical point produced an estimate of *Pc* ≈ 3.4 GPa. In principle, for a second order or slightly first order phase transition, this pressure should coincide or be somewhat lower than the pressure value *P*0(*TR*) at which *AR*(*TR*, *P*0(*TR*)) vanishes at room temperature *TR*. Unfortunately, however, this is incompatible with extrapolating the Landau parametrization of Hayward et al. [19] to room temperature. In fact, inserting the parameter value (10) into Equation (7a) yields *P*0(*TR*) ≈ 1.7 GPa, which is completely at odds with *Pc* ≈ 3.4 GPa as reported in Ref. [11]. To reach this transition pressure at room temperature would require reducing our previous result *λ*(0) *<sup>a</sup>* /*K*(0) = 0.002 obtained at *T* = 186 K by a full factor of 2 (cf. Figure 6). However, even then, any unbiased reader should have a hard time believing that the bare data of Figure 5 should indicate a continuous transition at *Pc* = 3.4 GPa. In summary, we hope to have demonstrated that standard Landau theory is completely inadequate to describe the HPPT of KMF

unless one is willing to sacrifice any numerical meaning to Landau theory, leaving us with all coupling parameters as essentially unknown and with ad hoc pressure dependencies at room temperature.

**Figure 5.** Lagrangian strains *a*, *<sup>t</sup>* calculated from unit cell data and Murnaghan fit of cubic part according to Ref. [11].

**Figure 6.** *AR*(*<sup>T</sup>* <sup>=</sup> *TR*, *<sup>P</sup>*) as defined in Equation (10) for correct slope <sup>−</sup>2*λ*(0) *<sup>a</sup>* /*K*(0) <sup>=</sup> <sup>−</sup>0.004 in comparison to a slope of <sup>−</sup>2*λ*(0) *<sup>a</sup>* /*K*(0) ≈ −0.002 assumed in Ref. [11].

#### **3. A Quick Review of FSLT**

In a nutshell, in a generic high pressure experiment, a given crystal is observed to change under application of high hydrostatic pressure *P* from an ambient pressure "laboratory" state *X* to a deformed state to be denoted as *X* = *X*(*P*) with an accompanying total (Lagrangian) strain *ηij*. Frequently, a high pressure phase transition manifests itself in such an experiment through the observation of relatively small strain anomalies on top of a much larger "background strain" that in itself is unrelated to the actual transition. Recognizing that the concept of strain is always defined with respect to a chosen elastic

reference state, in Ref. [5], a corresponding background system *X* = *X*(*P*), defined as the (hypothetical) equilibrium state of the system with the primary OP clamped to remain zero, was introduced. Let *αik* and *eik* = <sup>1</sup> <sup>2</sup> (∑*<sup>n</sup> αniαnk* − *δik*) denote the deformation and Lagrangian strain tensors from *X* to *X*, respectively. The total strain *ηij* may then be disentangled as a nonlinear superposition

$$
\eta\_{ij} = \mathfrak{e}\_{ij} + \mathfrak{a}\_{kl}\widehat{\mathfrak{e}}\_{kl}\mathfrak{a}\_{lj} \tag{11}
$$

of the—generally large—Lagrangian background strain and a relatively small spontaneous strain *ij measured relative to the floating "background" reference state X*. Determining the proper background strain *eij* in the resulting reference scheme

$$\underbrace{\mathbf{x} \xleftarrow{\mathfrak{a}\mathfrak{z}} \hat{\mathbf{x}} \xleftarrow{\hat{\mathfrak{a}\mathfrak{z}}} \hat{\mathbf{x}}}\_{\mathfrak{a}\mathfrak{y}} \hat{\mathbf{x}} \tag{12}$$

is thus a crucial step in correctly identifying the actual spontaneous strain, which in turn is mandatory in a successful application of the concepts of Landau theory. Effectively "subtracting" the elastic baseline, the strategy of FSLT therefore consists of setting up Landau theory *within the background reference system X*. Based on the reasonable assumption that a harmonic expansion with pressure-dependent elastic constants *Cijkl*[*X*] suffices to capture the elastic energy originating exclusively from the relatively small spontaneous strain *kl*, one arrives at the Landau free energy

$$\frac{F(Q; \hat{\boldsymbol{\Sigma}}; \hat{\mathbf{X}})}{V[\hat{\mathbf{X}}]} = \Phi(Q; \hat{\mathbf{X}}) + \sum\_{\mu \ge 1} Q^{2n} d\_{ij}^{(2\mu)}[\hat{\mathbf{X}}] \hat{\boldsymbol{\varepsilon}}\_{ij} + \frac{F\_0(\hat{\boldsymbol{\Sigma}}; \hat{\mathbf{X}})}{V[\hat{\mathbf{X}}]} \tag{13}$$

where we have assumed a scalar OP *Q* for simplicity, and

$$\frac{F\_0(\hat{\mathbf{E}}; \hat{\mathbf{X}})}{V[\hat{\mathbf{X}}]} \approx \sum\_{ij} \tau\_{ij} [\hat{\mathbf{X}}] \hat{\varepsilon}\_{ij} + \frac{1}{2} \sum\_{ijkl} \mathbb{C}\_{ijkl} [\hat{\mathbf{X}}] \hat{\varepsilon}\_{ij} \hat{\varepsilon}\_{kl} \tag{14}$$

denotes the pure spontaneous strain-dependent elastic free energy at hydrostatic external stress *τij*[*X*] = −*δijP*. For the pure OP potential part, we assume the traditional Landau expansion

$$\Phi(Q; \hat{X}) = \frac{A[\hat{X}]}{2}Q^2 + \frac{B[\hat{X}]}{4}Q^4 + \frac{C[\hat{X}]}{6}Q^6 + \dots \tag{15}$$

with *P*-dependent coefficients yet to be determined. In Ref. [9], it was explicitly shown that the harmonic structure of *<sup>F</sup>*0(;*X*) *<sup>V</sup>*[*X*] yields the equilibrium spontaneous strain

$$\stackrel{\circ}{\hat{\mathfrak{c}}}\_{mn} = -\sum\_{\nu=1}^{\infty} \mathbb{Q}^{2\nu} \sum\_{ij} d\_{ij}^{(2\nu)} [\hat{\mathcal{X}}] \mathbb{S}\_{mnij} [\hat{\mathcal{X}}] \tag{16}$$

where the compliance tensor *Smnij*[*X*] is defined as the tensorial inverse of the so-called *Birch coefficients* [20,21]

$$B\_{ijkl}[\hat{\mathcal{K}}] = \mathbb{C}\_{ijkl}[\hat{\mathcal{K}}] + \frac{1}{2} \left( \tau\_{jk}[\hat{\mathcal{K}}]\delta\_{il} + \tau\_{ik}[\hat{\mathcal{K}}]\delta\_{jl} + \tau\_{jl}[\hat{\mathcal{K}}]\delta\_{lk} + \tau\_{il}[\hat{\mathcal{K}}]\delta\_{jk} - 2\tau\_{ij}[\hat{\mathcal{K}}]\delta\_{kl} \right) \tag{17}$$

of the background system *X* which effectively take over the role of the elastic constants at finite strain. Furthermore, if we eliminate ¯ *mn* from *<sup>F</sup>*(*Q*,;*X*) *<sup>V</sup>*[*X*] by this formula, we obtain the renormalized pure OP potential

$$\Phi\_R(Q; \hat{X}) \quad := \quad \Phi(Q; \hat{X}) - \sum\_{\mu, \nu = 1}^{\infty} \frac{2\mu}{2\mu + 2\nu} \left( \sum\_{ijkl} d\_{ij}^{(2\mu)}[\hat{X}] S\_{ijkl}[\hat{X}] d\_{kl}^{(2\nu)}[\hat{X}] \right) Q^{2(\mu + \nu)}$$

$$\equiv \quad \frac{A\_R[\hat{X}]}{2} Q^2 + \frac{B\_R[\hat{X}]}{4} Q^4 + \frac{C\_R[\hat{X}]}{6} Q^6 + \dots \tag{18}$$

from which the equilibrium OP *Q*¯ can be determined by minimization.

Information on the *P*-dependence of elastic constants *Cijkl*[*X*] is usually available from density functional theory (DFT) or may be extracted from experimental measurements. At this stage, it therefore remains to relate the potential coefficients *A*[*X*], *B*[*X*], ... and *d* (2*μ*) *ij* [*X*], which are still defined with respect to the reference system *X* i.e., *P*-dependent. In a generic application of FSLT, however, one starts from knowledge of an ambient pressure Landau potential, i.e., the lowest order coefficients of the free energy

$$\frac{F(Q, \eta; X)}{V(X)} = \Phi(Q; X) + \sum\_{\mu=1}^{\infty} Q^{2\mu} \left( \sum\_{ij} d\_{ij}^{(2\mu, 1)} \eta\_{ij} + \frac{1}{2!} \sum\_{ijkl} d\_{ijkl}^{(2\mu, 2)} \eta\_{ij} \eta\_{kl} + \frac{1}{3!} \sum\_{ijklmn} d\_{ijklmn}^{(2\mu, 3)} \eta\_{ij} \eta\_{kl} \eta\_{mn} + \dots \right)$$

$$+ \frac{1}{2!} \sum\_{ijkl} \Gamma\_{ijkl}^{(2)} \eta\_{ij} \eta\_{kl} + \frac{1}{3!} \sum\_{ijklmn} \Gamma\_{ijklmn}^{(3)} \eta\_{ij} \eta\_{kl} \eta\_{mn} + \dots \tag{19}$$

(defined at *τij*[*X*] ≡ 0) with *P*-*independent* coefficients are assumed to be known, which obviously places constraints on the possible *P*-dependence of the above coefficients *A*[*X*], *B*[*X*], ... and *d* (2*μ*) *ij* [*X*]. To explore the relations between the two set of coefficients defined in the *P*-dependent background system *X* and the laboratory system *X*, we insert the nonlinear superposition relation (11) into (19) and compare coefficients. Following Ref. [9], we content ourselves to just include OP-strain couplings of type *<sup>Q</sup>*2*ij* and obtain

$$A[\hat{\mathbf{X}}] \quad = \frac{1}{f(a)} \left[ A[\mathbf{X}] + 2 \left( \sum\_{ij} d\_{ij}^{(2,0)} \varepsilon\_{ij} + \frac{1}{2!} \sum\_{ijkl} d\_{ijkl}^{(2,1)} \varepsilon\_{ij} \varepsilon\_{kl} + \frac{1}{3!} \sum\_{ijklmn} d\_{ijklmn}^{(2,2)} \varepsilon\_{ij} \varepsilon\_{kl} \varepsilon\_{mn} + \dots \right) \right] \tag{20a}$$

$$d\_{st}^{(2)}[\hat{\mathbf{X}}] \quad = \frac{1}{J(\mathbf{a})} \left( \sum\_{ij} a\_{si} d\_{ij}^{(2,0)} a\_{tj} + \sum\_{ijkl} a\_{si} d\_{ijkl}^{(2,1)} a\_{tj} \varepsilon\_{kl} + \sum\_{ijklmn} a\_{si} d\_{ijklmn}^{(2,2)} a\_{tj} c\_{kl} \varepsilon\_{mn} + \dots \right) \tag{20b}$$

in addition to the trivial relations *B*[*X*] = *B*[*X*]/*J*(*α*), *C*[*X*] = *C*[*X*]/*J*(*α*), where *J*(*α*) = det(*αil*).

The above parametrization scheme, although mathematically correct, is certainly not very convenient for applications in which the background reference state *X* is of high symmetry. For the most important example of a cubic high-symmetry phase, the above equations simplify considerable, since the deformation tensor *<sup>α</sup>ij* <sup>≡</sup> *αδij* is diagonal, and so is the Lagrangian background strain *eij* <sup>≡</sup> *<sup>e</sup>δij* with *<sup>e</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> (*α*<sup>2</sup> − <sup>1</sup>), which yields

$$A[\hat{X}] \quad = \frac{1}{\mathfrak{a}^3} \left[ A[X] + 2 \left( e \sum\_{i} d\_{ii}^{(2,0)} + \frac{\mathfrak{c}^2}{2!} \sum\_{ik} d\_{iikk}^{(2,1)} + \frac{\mathfrak{c}^3}{3!} \sum\_{iku} d\_{iikmu}^{(2,2)} + \dots \right) \right] \tag{21a}$$

$$d\_{st}^{(2)}[\hat{X}] \quad = \ \frac{1}{a} \left( d\_{st}^{(2,0)} + e \sum\_{k} d\_{stk}^{(2,1)} + e^2 \sum\_{km} d\_{stkkmm}^{(2,2)} + \dots \right) \tag{21b}$$

These formulas are as far as we can get without committing to a specific set of irreducible representations that determine the symmetry-allowed couplings in the Landau potential. Since the background strains *<sup>e</sup>*(*P*) ∼ *<sup>P</sup>* + *<sup>O</sup>*(*P*2), they effectively represent a set of highly interrelated power series in powers of *P*. Unfortunately, this also comes with all the inherent drawbacks. On the one hand, going beyond the lowest order terms, which may be taken from an previously determined ambient pressure Landau potential, we are forced to truncate the above series at rather low order to limit the number of additional unknown parameters. Of course, such truncated series inevitable diverge for increasing values of *P*. In addition, if we want to employ the theory for the purpose of fitting experimental data, we would like to fix certain experimental observables, in particular the pressure *P*<sup>0</sup> at which *AR*[*X*] vanishes. Given the above parametrization, this is obviously difficult to do.

In Section 2, we have demonstrated the considerable structural simplification of a standard ambient pressure Landau approach upon replacing Cartesian strain tensors by symmetry-adapted strains. Remarkably, it turns out that, despite the additional complicated nonlinearities contained in formulas (21a) and (21b), similar manipulations may also be carried out in FSLT and are found to yield equally substantial structural simplifications. In the rest of the paper, the resulting scheme will be derived and illustrated by describing the HPPT of KMF.

#### **4. Symmetry-Adapted FSLT: The Cubic-to-Tetragonal HP Phase Transition of KMF**

For *ν* = 0, 1, 2 . . . , we set

$$
\lambda\_a^{(\nu)} \equiv \frac{1}{\mathfrak{Z}} \sum\_{k\_1 \ldots k\_{\nu}} \sum\_{i} d\_{iik\_1 k\_1 \ldots k\_{\nu} k\_{\nu}}^{(2,\nu)} \tag{22a}
$$

$$\lambda\_t^{(\nu)} \equiv \sum\_{k\_1 \dots k\_\nu} \frac{d\_{33k\_1k\_1\dots k\_\nu k\_\nu}^{(2,\nu+1)} - d\_{11k\_1k\_1\dots k\_\nu k\_\nu}^{(2,\nu)}}{2\sqrt{3}} \tag{22b}$$

in accordance with

$$
\lambda\_a[\hat{\mathbf{X}}] \quad \equiv \quad \lambda\_a(\mathbf{c}) \equiv \frac{1}{3} \sum\_{\mathbf{s}} d\_{\mathbf{s} \mathbf{s}}^{(2)}[\hat{\mathbf{X}}] \tag{23a}
$$

$$
\lambda\_l[\hat{\mathbf{X}}] \quad \equiv \quad \lambda\_l(\mathbf{e}) \equiv \frac{d\_{33}^{(2)}[\hat{\mathbf{X}}] - d\_{11}^{(2)}[\hat{\mathbf{X}}]}{2\sqrt{3}} \tag{23b}
$$

Equations (21a) and (21b) are then replaced by

$$\left[\boldsymbol{\mu}^{3}A[\hat{\mathbf{X}}]\right] = \left[A\left[\mathbf{X}\right] + 6\sum\_{\nu=0}^{\infty} \frac{\lambda\_{\boldsymbol{a}}^{\left(\nu\right)}}{\left(\nu+1\right)!} e^{\nu+1}\tag{24a}$$

$$\langle \boldsymbol{a}\lambda\_{\boldsymbol{a},t} \vert \hat{\boldsymbol{X}} \rangle \quad = \sum\_{\nu=0}^{\infty} \lambda\_{\boldsymbol{a},t}^{(\nu)} \mathbf{e}^{\nu} \tag{24b}$$

Furthermore, we recall from nonlinear elasticity theory [20] that for a cubic system the compliance tensor *Sijkl*[*X*] and the bulk modulus *K*[*X*]=(*B*1111[*X*] + 2*B*1122[*X*])/3 at finite pressure are related by ∑*<sup>i</sup> Siijj*[*X*] = *S*1111[*X*] + 2*S*1122[*X*] = 1/3*K*[*X*], while, for the longitudinal shear modulus, *μ*[*X*] = (*B*1111[*X*] − *B*1122[*X*])/2 the relation *S*1111[*X*] − *S*1122[*X*] = 1/2*μ*[*X*] holds. If we replace the Cartesian

equilibrium spontaneous strain components (16) by their symmetry-adapted volume and tetragonal counterparts using these relations, they acquire the representations

$$
\hat{\boldsymbol{\varepsilon}}\_{d} = -\frac{\lambda\_{d}[\check{\mathbf{X}}]}{\mathsf{K}[\hat{\mathbf{X}}]} \mathbb{Q}^{2}, \qquad \hat{\boldsymbol{\varepsilon}}\_{t} = -\frac{2\lambda\_{t}[\check{\mathbf{X}}]}{\mu[\hat{\mathbf{X}}]} \mathbb{Q}^{2} \tag{25}
$$

Furthermore, it is easy to check the identity

$$\sum\_{ijkl} d\_{ij}^{(2)}[\hat{\mathcal{K}}] S\_{ijkl}[\hat{\mathcal{K}}] d\_{kl}^{(2)}[\hat{\mathcal{K}}] = \frac{\lambda\_a^2[\hat{\mathcal{K}}]}{\mathcal{K}[\hat{\mathcal{K}}]} \tag{26}$$

which allows for similarly rewriting the renormalized potential Φ*R*(*Q*; *X*) as

$$\Phi\_{\mathcal{R}}(Q;\hat{\mathbf{X}}) = \Phi(Q;\hat{\mathbf{X}}) - \frac{Q^4}{4} \left( \frac{2\lambda\_a^2[\hat{\mathbf{X}}]}{\mathcal{K}[\hat{\mathbf{X}}]} + \frac{8\lambda\_t^2[\hat{\mathbf{X}}]}{\mu[\hat{\mathbf{X}}]} \right) \tag{27}$$

Combining these equations with Equations (24a) and (24b), we can summarize the symmetry-adapted parametrization of the coefficients of the renormalized Landau potential (18), whose minimization determines the equilibrium OP *Q*¯ with

$$A\_{\mathbb{R}}[\hat{X}] \quad = \quad \frac{A[X]}{a^3} + \frac{6\Delta\_{\mathbb{R}}[\hat{X}]}{a^3} \tag{28a}$$

$$B\_{\mathbb{R}}[\hat{X}] \quad = \begin{array}{c} \frac{B[\mathbf{X}]}{\mathfrak{a}^3} - \frac{2(\lambda\_a^{(2)}[\hat{X}])^2}{\mathcal{K}[\hat{X}]} - \frac{8\lambda\_t^2[\hat{X}]}{\mu[\hat{X}]} \\ \end{array} \tag{28b}$$

in addition to *CR*[*X*] = *C*[*X*], where we introduced the function

$$\Delta\_{\mathfrak{d}}[\hat{X}] \equiv \Delta\_{\mathfrak{d}}(\mathfrak{e}) \equiv \sum\_{\nu=0}^{\infty} \frac{\lambda\_{\mathfrak{a}}^{(\nu)}}{(\nu+1)!} \mathfrak{e}^{\nu+1} \tag{29}$$

Note the close formal similarity of Equations (28a), (28b) and (29) with their infinitesimal counterparts (7a) and (7b). In a generic application of this theory with cubic high symmetry, the lowest order Landau coefficients *A*[*X*], *B*[*X*], *C*[*X*] and the lowest order strain-OP coupling coefficients *λ*(0) *<sup>a</sup>*,*<sup>t</sup>* can be taken from an existing ambient pressure Landau theory. Furthermore, the (diagonal) background deformation components *α* = *α*(*P*) and the resulting Lagrangian strain *e* = *e*(*P*) may be determined by fitting a suitable EOS to the cubic unit cell volume data. Such a fit also immediately yields the pressure-dependent bulk modulus *K*[*X*] ≡ *K*(*P*). It is only for the pressure-dependence of the longitudinal shear modules *μ*[*X*] = *μ*(*P*) that we are truly forced to resort to additional input from DFT. To determine the EOS and the elastic constants of KMnF3, we have performed a series of fairly standard DFT calculations. We refer to Appendix A for further details of these simulations.

Since we do not need to maintain the highest possible precision in determining *μ*(*P*) at room temperature but can content ourselves with a reasonable approximation, we use the following heuristic strategy to promote the DFT result *μ*DFT(*P*) from *T* = 0 to ambient temperature. Figure 7 shows a comparison of the bulk moduli *K*DFT(*P*) calculated from DFT to the result for *K*(*P*) obtained from the Murnaghan fit of the data published in Ref. [11]. Numerically, one finds the fraction of these moduli stays

entirely within the narrow bounds 0.918 ≤ *K*(*P*)/*K*DFT(*P*) ≤ 0.92 over the whole interval 0 ≤ *P* ≤ 20 GPa. We therefore postulate a similar behavior for *μ*(*P*), setting

$$
\mu(P) \equiv \frac{\mathcal{K}(P)}{\mathcal{K}\_{\text{DFT}}(P)} \cdot \mu\_{\text{DFT}}(P) \tag{30}
$$

Figure 8 illustrates the resulting behavior of *μ*(*P*).

**Figure 7.** Comparison of pressure-dependent bulk modulus *K*(*P*) at room temperature *T* = *TR* extracted from the Murnaghan fit of the cubic part of the unit cell data of Ref. [11] to the *T* = 0 result obtained from DFT.

**Figure 8.** Pressure-dependent longitudinal shear modulus *μ*(*P*) at room temperature *T* = *TR* as determined by extrapolating the corresponding DFT result from *T* = 0 to *T* = *TR* based on Formula (30).

With all other ingredients in place, this leaves only the higher order coefficients *λ*(*ν*) *<sup>a</sup>*,*<sup>t</sup>* , *ν* > 0 as the remaining undetermined parameters of FSLT. It therefore seems that these parameters must be treated as unknown fit parameters in a practical application. In the next section, however, we will introduce a much more convenient and powerful approach.

#### **5. Efficient Parametrization**

Comparing the above symmetry-adapted system of equations to what we had before, a considerable structural simplification is obvious. However, the following drawbacks seem to persist:


In what follows, we propose a new scheme which finally allows for practically overcoming all of these problems in a single push. Let us start by taking a second look at Figure 6. Of course, the correct functional form of *AR*[*X*] ≡ *AR*(*TR*, *P*) that we are looking for should be that of a well-behaved function passing through zero around *Pc* ≈ 4 GPa. At *P* = 0, however, it should start out with roughly twice the initial slope of the purely linear yellow curve if we are to retain the ambient pressure Landau parameters. The pressure dependence of the correct function *AR*(*TR*, *P*) must therefore be far from linear. Experimentally, there is no indication of a re-entrant behavior below 25GPa, so we do not expect any second crossing point in this pressure range. Unfortunately, our numerical tests quickly revealed that, using a truncated version of Δ*a*(*e*(*P*)) as defined in (29), it seems virtually impossible to meet these requirements unless one is willing to go to prohibitively high truncation order, thus introducing a plethora of unknown fit parameters and the accompanying horrific numerical problems.

A way out is to propose a reasonable function Δ*a*(*e*(*P*)) in closed form that meets all of the above requirements while still containing some adjustable parameters that offer a certain amount of variational flexibility to allow improvement by least squares fitting. Of course, there are many ways to do this, and the choices are only limited by the reader's ingenuity. In fact, since the summand *ν* = 0 of Δ*a*(*x*) contributes the lowest order coefficient *xλ*(0) *<sup>a</sup>* , which is usually fixed from knowledge of an existing ambient pressure Landau theory, all candidate functions Δ*a*(*x*) that start out like

$$
\Delta\_a(\mathbf{x}) = \lambda\_a^{(0)}\mathbf{x} + O(\mathbf{x}^2) \tag{31}
$$

qualify as candidates for a meaningful function Δ*a*(*x*). In our current application to KMF, we parametrize

$$A\_R \equiv A\_R(T\_R, \varepsilon(P)) = \frac{A\_R(T\_R) + 6\Delta(\varepsilon(P); b, \mathfrak{c}, d)}{a^3(P)}\tag{32}$$

by introducing the function

$$
\Delta(e; b, c, d) := ce + 2bf\_d(e) \tag{33}
$$

with parameters *b*, *c*, *d* combining a linear part with slope *c* and a nonlinear contribution *fd*(*e*), which remains yet to be specified. Our general parametrization strategy is then as follows:


These steps eliminate parameters *b*, *c* in favor of the constants *α* and *e*0, leaving *d* as the single remaining free variational parameter. We still need to reconcile this parametrization of *AR* with that of the function *λa*[*X*] as defined in Equation (24b). We focus on the power series part

$$\Lambda(\mathbf{x}) \equiv \sum\_{\nu=0}^{\infty} \lambda\_a^{(\nu)} \mathbf{x}^{\nu} \tag{34}$$

which is based on the same set of coefficients *λ*(*ν*) *<sup>a</sup>* as Δ*a*(*x*) but lacks the accompanying factorials 1/*ν*!. These factorials can, however, easily be taken care of. Observe that

$$\int\_0^\infty dt \, e^{-t} \frac{(t\mathbf{x})^{\nu+1}}{(\nu+1)!} = \mathbf{x}^{\nu+1} \tag{35}$$

Therefore, using the *Borel transform*

$$(\mathcal{B}\Delta)(\mathbf{x}) \equiv \int\_0^\infty dt \, e^{-t\mathbf{x}} \Delta(t\mathbf{x}) \tag{36}$$

we may relate

$$
\Lambda(\mathbf{x}) = \frac{(\mathcal{B}\Lambda)(\mathbf{x})}{\mathbf{x}} \tag{37}
$$

It remains to specify a suitable function *fd*(*x*). Beyond producing a reasonable function *AR*(*TR*,*e*), the job profile for recruiting such a function includes at least two basic requirements:


For the present goal of understanding the HPPT in KMF, we came up with the choice

$$f\_d(\mathfrak{e}) := -1 + \sqrt{1 - d^2 \mathfrak{e}} \tag{38}$$

(note that *e*(*P*) < 0 for *P* > 0) which meets both of these minimal requirements, since, in this case, *<sup>b</sup>* = (*<sup>c</sup>* − *<sup>s</sup>*)/*d*<sup>2</sup> and

$$\varepsilon = \frac{s\left(\sqrt{1 - d^2 \varepsilon\_0} - 1\right) - ad^2 / 2}{d^2 \varepsilon\_0 / 2 + \sqrt{1 - d^2 \varepsilon\_0} - 1} \tag{39}$$

In this way, we have completely bypassed Taylor series expansions and their various accompanying drawbacks. Figure 9 illustrates the remaining variational freedom in our chosen parametrization.

**Figure 9.** Illustration of the remaining *d*-dependence of the parametrization of *AR*(*TR*,*e*(*P*)) by Equation (32) (upper panel) and the resulting function Λ(*e*(*P*)) as given by the Borel transform Equation (37) (lower panel). Parameters *<sup>b</sup>* and *<sup>c</sup>* were eliminated in favor of parameter *<sup>λ</sup>*(0) *<sup>a</sup>* = 0.128 GPa as prescribed from ambient pressure LT and a chosen pressure parameter *P*<sup>0</sup> = 6.0 GPa.

In summary, at this stage, the function Δ*a*[*X*] which governs the behavior of *AR*[*X*] and—through a Borel transform—also provides the coupling function *<sup>λ</sup>a*[*X*] between the spontaneous volume strain ¯ *<sup>a</sup>* and *Q*¯ <sup>2</sup> according to Equation (25) has been specified up to a single free parameter *d*. The remaining coupling function *λa*[*X*] ≡ *λa*(*P*) explicitly determines the proportionality between the spontaneous tetragonal strain ¯ *<sup>a</sup>* and *<sup>Q</sup>*¯ 2, However, both *<sup>λ</sup>a*[*X*] and *<sup>λ</sup>t*[*X*] also enter implicitly into the spontaneous strain via its implicit dependence on the quartic coefficient *BR*[*X*] of the renormalized Landau potential as given by Equation (28b), and, apart from the reasonable requirement lim*P*→<sup>0</sup> *<sup>λ</sup>t*[*X*] = *<sup>λ</sup>*(0) *<sup>t</sup>* , nothing is known in advance about its pressure dependence, such that introducing a truncated Taylor series and least squares

fitting seems unavoidable. However, we can actually do a lot better than this. Let us consider experimental high pressure spontaneous strain data in the form of *<sup>n</sup>* data points (*Pi*, *a*,exp(*Pi*), *a*,exp(*Pi*))*<sup>n</sup> <sup>i</sup>*=1. In fact, at a given prescribed pressure value *Pi* and with all other parameters in place, we can regard the room temperature values *<sup>a</sup>* <sup>=</sup> *a*(*λt*(*Pi*)) and *<sup>t</sup>* <sup>=</sup> *t*(*λt*(*Pi*)) as functions of the-unknown-function values *<sup>λ</sup>t*(*Pi*). The "best" value *<sup>λ</sup>t*(*Pi*) matching the data point (*Pi*, *a*,exp(*Pi*), *a*,exp(*Pi*)) may then be determined by numerically minimizing the function

$$\begin{split} s\_l(\lambda\_l(P\_l)) &:= \quad \quad w\_d \left[ \widehat{\mathfrak{e}}\_{d,\text{exp}}(P\_l) - \widehat{\mathfrak{e}}\_d(\lambda\_l(P\_l)) \right]^2 \\ &\quad \quad + w\_l \left[ \widehat{\mathfrak{e}}\_{l,\text{exp}}(P\_l) - \widehat{\mathfrak{e}}\_l(\lambda\_l(P\_l)) \right]^2 \end{split} \tag{40}$$

with weights *wa*, *wt* suitable adjusted to counterbalance size differences between *<sup>a</sup>* and *t*. Carried out for all *i* = 1, ... , *n*, this prescription results in a collection of *n* "optimal" values *λt*(*Pi*) from which one may hope to recover the full function *λt*(*P*) by interpolation. Figure 10 shows the result of our corresponding effort for KMF. Amazingly, we observe that all values *λt*(*Pi*) seem to accumulate on a straight line whose extrapolation *<sup>P</sup>* <sup>→</sup> 0 perfectly passes through the point (0, *<sup>λ</sup>*(0) *<sup>t</sup>* ), which is just the limiting value imposed by ambient pressure Landau theory. We believe that this behavior is not coincidental but a strong indication that the present parametrization is internally consistent and correct.

**Figure 10.** Results of minimizing the sum (40) using the relative weights (*wa*, *wt*)=(10, 1) and pressure parameters *<sup>P</sup>*<sup>0</sup> <sup>=</sup> 6GPa, *<sup>d</sup>* <sup>=</sup> 40. The dashed horizontal line indicates the limiting value *<sup>λ</sup>t*(*<sup>P</sup>* <sup>=</sup> <sup>0</sup>) <sup>≡</sup> *<sup>λ</sup>*(0) *<sup>t</sup>* .

A simple linear fit of *λt*(*P*) therefore completes our Landau parametrization. Our results for the pressure dependence of the couplings *AR*(*P*), *B*(*P*), *λa*,*t*(*P*) and the resulting pressure dependence of the equilibrium OP *Q*¯(*P*) are illustrated in Figure 11. Note that in the present description the transition appears to be of first rather than second order, with *P*<sup>0</sup> = 6 GPa yielding a transition pressure *Pc* ≈ 4.5 GPa. Finally, the resulting parametrization of the spontaneous strain is compared to experiment in Figure 12.

**Figure 11.** Upper left panel: pressure dependence of Landau parameters *AR*(*P*) as compared to the simple linear behavior of infinitesimal strain LT assumed in Equation (7a). Upper right panel: pressure dependence of *B*˜(*P*). The dashed horizontal line indicates the corresponding value from the parametrization of Hayward et al. [19] taken at *T* = 186.5 K, which is slightly displaced from our limiting value at *P* = 0 since we are taking into account thermal softening of *K*(*P*) and *μ*(*P*) for room temperature. Lower left panel: pressure dependence of coupling parameters *λa*(*P*)/*K*(*P*) and *λt*(*P*)/*μ*(*P*). Lower right panel: resulting behavior of the equilibrium order parameter *Q*¯(*P*) (right).

**Figure 12.** Parametrization of spontaneous strains *a*, *<sup>t</sup>* in comparison to experimental data from Ref. [11]. The thin vertical line indicates the pressure parameter *P*<sup>0</sup> = 6 GPa.

#### **6. Discussion**

For the high pressure community, our present paper contains some good news and some bad news—first, the bad news. We hope to have demonstrated convincingly that classical Landau theory is usually completely inadequate for "explaining" experimental findings in the field of high pressure phase transitions, and the conclusions drawn from it will often be misleading at best. Taking the example of KMF, both the first order character of the transition and the correct value of the transition pressure are completely obscured by sticking to Landau theory with infinitesimal strain coupling, even if one is willing to distort a pre-existing ambient pressure Landau parametrization beyond recognition. The good news, however, is that with the development of FSLT a mathematically consistent alternative incorporating nonlinear elasticity has recently become available. Up to now, however, FSLT has been rather complicated in structure, which quite likely scared off many potential users and thus did not lead to the widespread use that its developers were initially hoping for. The present paper, which exploits the enormous simplifications that arise by passing (i) from Cartesian to symmetry-adapted finite strains and (ii) by virtue of (i), from truncated Taylor expansions to a functional parametrization. These improvements should pave the way for routine use of our theory in successfully describing HPPTs. In particular, the present paper illustrates that, while the former version of FSLT involved delicate least-squares fitting procedures with a large number of unknown fit parameters and dealing with all the inadequacies implicit in the use of truncated Taylor expansions, in the present scheme, the unknown pressure dependencies can be systematically determined one-by-one in a step-wise, almost "deterministic" manner.

Admittedly, even the present parametrization of the HPPT in KMF is still less than perfect. For instance, the coupling function *λa*(*P*)/*K*(*P*) shown in the lower left panel of Figure 11 exhibits a steep initial decrease with increasing pressure. Since no spontaneous strain data are available for this pressure region, this does not change the physical values produced by the theory. However, it hints at a sub-optimal choice Equation (38) for the auxiliary function *fd*(*x*). The reader is invited to come up with an improved candidate function.

More importantly, in a full-blown application of FSLT, we should be able to predict e.g., the (*P*, *T*) cubic-to-tetragonal phase boundary of KMF and compute pressure—and temperature-dependent elastic constants. In principle, the ability to do so depends mainly on the successful construction of a pressure—and temperature-dependent baseline, i.e., the cubic EOS *V*cubic = *V*cubic(*P*, *T*). In a previous paper [10], this task has been successfully carried out for the perovskite PbTiO3 by combining zero temperature DFT calculations (see Appendix A for details) with the Debye approximation as implemented in the GIBBS2 package [22,23] to incorporate effects of thermal expansion (recently, we learned [24] that a similar approach also seems to work for MgSiO3). Unfortunately, our corresponding efforts to derive *V*cubic(*P*, *T*) for KMF along the same lines have failed so far, however. This failure manifests itself e.g., in the inability to simultaneously reproduce the thermal baseline at *P* = 0 measured experimentally and the ambient temperature EOS, even if we allowed for the introduction of a constant compensating background pressure which is frequently introduced in DFT calculations to compensate inadequacies of an employed exchange-correlation functional. Ref. [25] states that perovskites with octahedral tilting generally do not show an appreciable coupling between structural and magnetic order parameters. Nevertheless, this statement obviously does not exclude effects due to a coupling between magnetic degrees of freedom and the background volume strain. Since the Debye approximation is based exclusively on phonons, we may speculate that in KMF residual magnetic effects may be responsible for additional thermal energy consumption. An investigation of this problem is currently underway.

Finally, it may be argued that our current theory also does not "explain" the origins of the involved nonlinear *P*-dependencies on a fundamental level. In approaches based on infinitesimal strain couplings, similarly looking *P*-dependent couplings are also introduced, but in a more or less completely ad hoc

manner, blaming their existence on rather unspecific "higher order strain couplings". As a rule, such an approach results in a mathematically inconsistent theory. In contrast, the nonlinearities that arise in FSLT result for different but well-defined reasons. On the one hand, there are couplings between powers of the background strain *e*(*P*) and the Landau potential Φ*R*(*Q*; *X*) "floating" on this background strain, and there is a pressure-dependence of the elastic moduli *K*(*P*) and *μ*(*P*), which is in principle accessible e.g., to DFT calculations. This leaves the task of "explaining" the residual nonlinearities in the functions *λa*,*t*(*P*) describing the couplings between order parameter and spontaneous strains. In contrast to blaming their existence on the effects of some unspecified higher order couplings, our present theory provides a practical way to numerically determine these functions. In addition, we see no reason in principle as to why such *P*-dependent coupling constants could not eventually be extracted from DFT calculations along the general philosophy laid out in Refs. [26,27] and the subsequent follow-up literature.

**Author Contributions:** A.T.: concept, main theory, and writing; W.S.: additional theory, comparison to previous experimental data; S.E., K.B., and P.B.: DFT calculations. All authors have read and agreed to the published version of the manuscript.

**Funding:** A.T., S.E., and K.B. acknowledge support by the Austrian Science Fund (FWF) Project P27738-N28. W.S. acknowledges support by the Austrian Science Fund (FWF) Project P28672-N36. S.E. acknowledges support from H.E.C., Pakistan.

**Acknowledgments:** A.T. would like to thank Ronald Miletich-Pawliczek for useful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

AFN antiferromagnetic DFT density functional theory EOS equation of state FSLT finite strain Landau theory HPPT high pressure phase transition KMF *KMnF*<sup>3</sup> LDA local density approximation LT Landau theory OP order parameter PT phase transition NM nonmagnetic

#### **Appendix A. DFT Calculation Details**

The EOS and the elastic constants in the cubic phase of KMF were calculated using the WIEN2K DFT package, which is an all-electron code based on the (linearized) augmented plane-wave and local orbitals [(L)APW+lo] basis representation of the Kohn–Sham equations [28] of DFT. Here, we only content ourselves with a quick outline of the basic ideas and refer to Refs. [29,30] and the monograph Ref. [31] for more details.

In the (L)APW+lo method, the crystal's unit cell is partitioned into a set of atomic spheres surrounding the nuclei and a remaining interstitial region. Inside these atomic spheres, the wave functions are expanded into atomic-like basis functions, i.e., numerical radial functions times spherical harmonics while they are represented by plane waves throughout the interstitial region. These two regions are glued together by requiring continuity of the basis functions in value (and, depending on the flavor of the (L)APW+lo method, in radial slope) across the sphere boundaries. These LAPW calculations require making a couple

of choices regarding cell size, standard parameter values, etc. In detail, the calculations of the present work were done with *R*min MT *K*max = 9 and atomic sphere radii of 2.1, 2.0, and 1.6 bohr for K, Mn and F, respectively. The energy separation between core and valance states used was −6.0 Ry. Inside the sphere, the maximum angular momentum used in the spherical expansions was L*max* = 10, while the charge density in the interstitial was Fourier expanded up to a cutoff of *G*max = 14(*a*.*u*.)−1.

As to the use of exchange-correlation functionals, we have performed calculations using the standard local density approximation (LDA) [32] and three functionals of the generalized gradient approximation, namely, PBE from Perdew et al. [33,34], its solid-state optimized version PBEsol [35], and WC from Wu and Cohen [36]. For LDA and PBE, we observed the well-known tendency to underestimate and overestimate the lattice constants of solids, respectively, while PBEsol and WC produced more accurate results in between LDA and PBE [37,38].

To study the effect of magnetism on our results, we undertook calculations for a standard non-magnetic (NM) and ferromagnetic (FM) cubic perovskite unit cell with five atoms as well as an for a cubic supercell with 10 atoms in antiferromagnetic (AFM) structures of A-type, C-type and G-type [39]. After sufficient testing, we settled for a k-mesh sampling of 10 × 10 × 10 k-points for all types of structures. In FM and AFM structure, we observed a linearization error which is inherent to the basis functions inside the spheres. To overcome this problem, we use the second energy derivative of the radial part (HDLO) for d electrons for Mn atom (for more detail see Ref. [40]). Using this setup, we identified that G-type AFM structure to have lowest energy.

Cubic elastic elastic constants were calculated with the help of the WIEN2K add-on package by Charpin [41]. Results at *T* = 0 are compiled in Table A1. In particular, we conclude that FM and AFM structures give similar value of lattice parameters and bulk modulus, which suggests that the specific magnetic ordering is not overly important for these quantities. In passing, we note that non-magnetic KMF is found to be a metal (with all functionals), but all magnetic structures lead to insulators (with all functionals). Simulations with PBEsol+U with U = 4eV in the AFM G-type phase would even lead to slightly better agreement with experiment for lattice constants and bulk modulus, but overall the effect is small and we therefore used the PBEsol results in Figures 7 and 8.


**Table A1.** Lattice constant (Å) and bulk modulus (GPa) of cubic KMnF3 for different methods used in this work.

Estimating an additional softening due to finite temperature with various flavors of the Debye approximation which are implemented in the GIBBS2 software [22,23], the best agreement with the cubic part of the experimental data of Ref. [11] was reached for the combination of G-type antiferromagnetic structure and PBEsol functional.

#### **References**



© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Phase Transition and Dehydration in Epsomite under High Temperature and High Pressure**

**Linfei Yang 1,2, Lidong Dai 1,\* , Heping Li 1, Haiying Hu 1, Meiling Hong 1,2 and Xinyu Zhang 1,2**


Received: 15 January 2020; Accepted: 28 January 2020; Published: 30 January 2020

**Abstract:** The phase stability of epsomite under a high temperature and high pressure were explored through Raman spectroscopy and electrical conductivity measurements in a diamond anvil cell up to ~623 K and ~12.8 GPa. Our results verified that the epsomite underwent a pressure-induced phase transition at ~5.1 GPa and room temperature, which was well characterized by the change in the pressure dependence of Raman vibrational modes and electrical conductivity. The dehydration process of the epsomite under high pressure was monitored by the variation in the sulfate tetrahedra and hydroxyl modes. At a representative pressure point of ~1.3 GPa, it was found the epsomite (MgSO4·7H2O) started to dehydrate at ~343 K, by forming hexahydrite (MgSO4·6H2O), and then further transformed into magnesium sulfate trihydrate (MgSO4·3H2O) and anhydrous magnesium sulfate (MgSO4) at higher temperatures of 373 and 473 K, respectively. Furthermore, the established *P*-*T* phase diagram revealed a positive relationship between the dehydration temperature and the pressure for epsomite.

**Keywords:** epsomite; phase transition; dehydration reaction; Raman spectra; electrical conductivity; high pressure

#### **1. Introduction**

In recent decades, hydrated sulfates have attracted a large amount of interest due to their great importance in exploring the interior structure of icy satellites, such as Europa, Ganymede and Callisto. It was widely reported that hydrated sulfates might be dominant minerals in the interior of these icy satellites, which was proved by the infrared spectral result collected from the Galileo spacecraft [1,2]. In addition, the discovery that anhydrous sulfates occur in carbonaceous chondritic meteorites also provides evidence for the existence of hydrated sulfates in the icy mantle of satellites, and these hydrated sulfates can be formed during the accretion of icy satellites [3]. In consideration of the high-pressure and high-temperature environment in the interior of these icy satellites, it is possible that these hydrated sulfates undergo a series of pressure-induced phase transitions and dehydration reactions. More importantly, as illustrated by Fortes and Choukroun [4], these transformations would have great impacts on the internal structure and heat transport in the icy mantle. As a representative magnesium (Mg)-bearing hydrated sulfate for epsomite (MgSO4·7H2O), the investigation into its optical and electrical properties under a high pressure and temperature could help us to deeply understand the interior structure, composition and physical property the icy satellites.

Many previous works have been devoted to studying the phase stability of epsomite under a high pressure and room temperature, but their results showed great inconsistencies. Livshits et al. investigated the high-pressure behavior of epsomite and found four phase transitions, occurring at ~0.4,

~1.2, ~1.6 and ~2.5 GPa, respectively [5]. Gromnitskaya et al. employed ultrasonic and neutron scattering measurements on epsomite to reveal its elastic property under a high pressure of up to ~3.0 GPa, and three phase transitions were identified at pressures of ~1.4, ~1.6 and ~2.5 GPa, respectively [6]. Grasset et al. established the phase relation of MgSO4–H2O under a pressure of up to ~2.0 GPa in a diamond anvil cell (DAC), but they only determined one phase transition for epsomite at ~0.6 GPa [7]. The results from Nakamura et al. indicated that there were no any phase transitions in epsomite up to ~4.5 GPa [8]. Therefore, it is still unclear how epsomite behaves under a high pressure and room temperature, which requires more detailed experimental investigations. Raman spectroscopy and electrical conductivity measurements have been demonstrated to be efficient methods to characterize the pressure-induced structural variations in various hydrated sulfates. In view of this, systematic investigations on the optical and electrical properties of epsomite under a high pressure and room temperature will provide a good deal of insight into its high-pressure behaviors.

It is a common phenomenon for hydrated sulfates to undergo dehydration reactions to form lower hydrates under a high temperature and pressure. A great number of previous high-pressure studies were undertaken, mainly concerned with the dehydration process of gypsum, chalcanthite and blödite, and all of them were found to have a positive relationship between the dehydration temperature and pressure [9–13]. For example, Yang et al. found that the dehydration temperature of gypsum gradually increased with the rise in pressure, and the dehydration boundary between gypsum and bassanite was determined to be *P* (GPa) = –23.708 + 0.050*T* (K) [13]. As a similar hydrated sulfate, it is possible that the dehydration temperature of epsomite would also be significantly affected by pressure. However, to the best of our knowledge, there are no relative experimental reports on the issue of epsomite dehydration under a high pressure and temperature.

The present study clarified the phase transition and dehydration reaction for epsomite under a high pressure and temperature, using Raman scattering and electrical conductivity measurements in a diamond anvil cell (DAC). Furthermore, we determined a *P*-*T* phase diagram for epsomite, which is crucial to understanding the pressure effect on the dehydration reaction of hydrated sulfates, and also has great implications for modeling the internal structure of icy satellites.

#### **2. Materials and Methods**

Natural epsomite was used as the starting sample for all of the high-pressure experiments, which was obtained from a phosphorus-bearing rock series in Songlin town, Zunyi county, Guizhou province. The high pressure and room temperature experiments were implemented in a symmetric diamond anvil cell (DAC) with a 300 μm anvil culet. In the high-temperature and high-pressure experiments, we placed an external resistive heater around the diamond to achieve the high temperature conditions. A K-type thermocouple wire was directly pasted onto the diamond surface in order to measure the temperature in the pressure chamber. A classic ruby pressure calibration formula from Mao et al. was adopted to accurately calculate the pressure in the sample chamber, based on the shifts of R1 fluorescence lines [14]. Taking into account the inevitable influence of the high temperature on the pressure conditions in the sample chamber, the high temperature-corrected equation of pressure calibration was selected in the process of the high-temperature and high-pressure measurements [15]. The Raman spectra were collected by a micro-confocal Raman spectrometer combined with a 514.5 nm argon-ion laser as the excitation source. Helium was selected as the pressure medium to provide the hydrostatic environment in the sample chamber. The resolution of the Raman spectrometer and the repetition rate were 1 and 0.2 cm−1, respectively. Each Raman spectra was obtained in the wave number range of 900–1200 cm−<sup>1</sup> for the internal vibration modes of the sulfate and 3000–3800 cm−<sup>1</sup> for the water molecular stretching modes. The acquisition time for each Raman spectra was set to 120 s. As for electrical conductivity experiments, a Solartron-1260 impedance/gain phase analyzer was used in this study to acquire the AC impedance spectrum. No pressure medium was used in order to avoid the introduction of impurity substances. The laser drilling technique was used to drill a 200 μm hole in the pre-indented T-301 stainless steel gasket, with a thickness of 50 μm. A mixture of boron

nitride (BN) and epoxy was then compressed into the hole as the insulator, and another hole of 100 μm was drilled to provide an insulating sample chamber. More specific descriptions of the measurement procedure have been reported in previous works [16–18].

#### **3. Results and Discussion**

Under ambient conditions, three Raman peaks are observed to be located at positions of 983.3, 3322.6 and 3458.4 cm−1, respectively. The sharp peak at 983.3 cm−<sup>1</sup> is assigned to the ν<sup>1</sup> symmetric stretching mode of the sulfate tetrahedra, and the other two bands at 3322.6 and 3458.4 cm−<sup>1</sup> are ascribed to the ν<sup>1</sup> and ν<sup>3</sup> vibrational modes of water molecules, respectively. All of these acquired Raman modes for the epsomite are generally consistent with the previous data by Wang et al. [19].

Figure 1a displays the high-pressure Raman spectra of the epsomite in the SO4 vibrational mode range at pressures of 0–12.8 GPa and room temperature. The corresponding pressure shifts for the ν<sup>1</sup> (SO4) mode are plotted in Figure 1b, from which one obvious inflection point can be obtained, at a pressure of ~5.1 GPa. From ~0 to ~5.1 GPa, the vibrational mode of ν<sup>1</sup> (SO4) slightly shifts towards the higher wave numbers, with a slope of 0.16 cm−<sup>1</sup> GPa−1. However, upon further compression from ~5.1 to ~12.8 GPa, the ν<sup>1</sup> (SO4) mode exhibits a strong pressure dependence with a relatively high slope value of 4.46 cm−<sup>1</sup> GPa<sup>−</sup>1. In the meantime, it can be seen that the full width at half maximum (FWHM) of the ν<sup>1</sup> (SO4) peak also shows a discontinuous change at the same pressure point of ~5.1 GPa. As shown in Figure 1c, the FWHM of the ν<sup>1</sup> (SO4) band remains at an almost constant value of about 8.3 cm−<sup>1</sup> below ~5.1 GPa, but it gruadally increases with pressure above ~5.1 GPa. Figure 2 shows the Raman spectra of the epsomite in the OH-stretching mode range under a high pressure and room temperature. Obviously, both of the two OH-stretching peaks at 3322.6 and 3458.4 cm−<sup>1</sup> suddenly shift to lower frequencies above ~5.1 GPa. Upon decompression, the Raman spectra for the recovered epsomite exhibit a similar characterization to those of the starting sample.

**Figure 1.** (**a**) The Raman spectra of the epsomite under a high pressure, up to ~12.8 GPa, and at room temperature, collected in the frequency range of 900–1200 cm<sup>−</sup>1; (**b**) the pressure dependence of the

Raman shift for the ν<sup>1</sup> (SO4) mode; (**c**) the full width at half maximum (FWHM) of the ν<sup>1</sup> (SO4) peak as a function of the pressure.

**Figure 2.** The high-pressure Raman spectra of the epsomite at room temperature in the range of 3000–3800 cm<sup>−</sup>1.

In this study, we determined a pressure-induced structural phase transition in epsomite at ~5.1 GPa and at room temperature, based on the variation in Raman modes of SO4 tetrahedra and water molecules. The recovery of the Raman spectra after decompression reflects that this pressure-induced phase transformation is reversible. As for the nature of this structural phase transition, no evidence shows that it is associated with the dehydration process of epsomite, because the Raman spectra obtained in the whole pressure range are not in agreement with the dehydration product of the epsomite. We think that this is possibly related to the distortion of the sulfate tetrahedra and the strengthening of the hydrogen bonding. A recent work by Gonzalez–Platas et al. revealed that the chalcomenite underwent an isostructural phase transition at around ~4.0 GPa, which is very close to the transition point of epsomite that was revealed in this study (~5.1 GPa) [20]. The similarities between the transition pressures are closely associated with their crystal structure, because both of them have an orthorhombic structure with the same space group of *P*212121. In addition, other hydrated sulfates have also been revealed to undergo phase transitions under a high pressure and room temperature, such as gypsum, chalcanthite and mirabilite [12,21–23]. This implies that the pressure can play a very important role in tuning the structural and vibrational properties of most hydrated sulfates.

In order to further check the structural phase transition in epsomite, high-pressure electrical conductivity experiments were carried out in the pressure range of 0.5–12.4 GPa and under room temperature. Figure 3a presents the collected impedance spectra of epsomite in the frequency range from 10−<sup>1</sup> to 10<sup>7</sup> Hz and at a series of pressure points, and the horizontal and vertical axes represent the real and imaginary parts of the complex impedance, respectively. Each spectra datum is composed of two semicircular arcs, which stand for grain interior and boundary contributions, respectively. All of these displayed impedance spectra were fitted using the equivalent circuit method in the Z-View software to obtain the grain interior resistance of the sample. The equivalent circuit was composed of a resistor (R) and a constant phase element (CPE) in parallel. The electrical conductivity for the grain interior contribution of sample was then calculated by the formula as follows:

$$\sigma = L/SR$$

where σ stands for the electrical conductivity, *L* is the distance between two electrodes, *S* represents the electrode and the *R* denotes the fitting resistance. As shown in Figure 3b, according to the obvious change in the slope of electrical conductivity, this diagram can be divided into two distinct regions: (i) in the pressure range of 0.5–4.6 GPa, the electrical conductivity of the epsomite monotonously increases with pressure, and the corresponding pressure coefficient is determined to be 0.06 S cm−<sup>1</sup> GPa<sup>−</sup>1; (ii) from ~4.6 to ~12.4 GPa, its electrical conductivity still shows an increasing tendency but with a smaller slope value of ~0.04 S cm−<sup>1</sup> GPa<sup>−</sup>1. The abrupt change in resistivity indicates the occurrence of a phase transition [24]. This transition point at ~4.6 GPa obtained from electrical conductivity measurements is slightly lower than the data from our high-pressure Raman spectroscopic experiments, which is possibly caused by the different pressure environment in the sample chamber. In comparison with the hydrostatic environment for the Raman measurements, it is non-hydrostatic for electrical conductivity experiments. There exists relatively high deviatoric stress in the non-hydrostatic compression, which has been reported to significantly influence the high-pressure behavior of materials [25]. Our electrical conductivity results provide another critical clue to support the occurrence of a pressure-induced phase transition in epsomite.

**Figure 3.** (**a**) The typical impedance spectra for the epsomite, measured at 0.5–12.4 GPa and room temperature; (**b**) the calculated electrical conductivity of the epsomite with increasing pressure.

To reveal the effect of the high pressure on the dehydration process of epsomite, we in situ measured a series of Raman spectra of the epsomite under high temperatures of up to ~623 K and at constant pressure points of ~0.8, ~1.3, ~3.7 and ~6.4 GPa. Figure 4 presents the Raman spectra for the epsomite in the temperature range of 293–573 K and at a pressure of ~1.3 GPa. Characterized by the variation in the temperature dependence of the sulfate group and hydroxyl modes, three important temperature points were well-determined at ~343, ~373 and ~473 K, respectively. At ~343 K, it can be clearly observed from Figure 4a that the ν<sup>1</sup> (SO4) peak at 983.3 cm-1 abruptly shifted to a lower frequency of 982.4 cm<sup>−</sup>1. Moreover, we also detected some obvious variations in the two OH-stretching modes, which were changed into one mode at a position of ~3441.5 cm−<sup>1</sup> (Figure 4b). These observed new peaks at ~982.4 and ~3441.5 cm−<sup>1</sup> agreed well with the position of the characteristic Raman vibrational mode of hexahydrite (MgSO4·6H2O) and, hence, we think that the epsomite initially dehydrated to a new hydrous phase hexahydrite at ~343 K. At a higher temperature of ~373 K, the ν<sup>1</sup> (SO4) peaks suddenly broadened and moved to a higher wave number of 1023.2 cm−1, meanwhile, two new OH-stretching modes also started to appear at positions of 3153.03 and 3341.8 cm−1, respectively. The Raman spectra obtained at ~373 K exhibited similar peak characteristics to those of magnesium sulfate trihydrate (MgSO4·3H2O) [19]. Therefore, this provided robust evidence for the occurrence of another dehydration reaction, from hexahydrite to a new trihydrate phase. When the temperature was further enhanced, up to ~473 K, the complete dehydration of epsomite to form anhydrous magnesium sulfate (MgSO4) was finally observed by the appearance of a new peak at 1022.9 cm<sup>−</sup><sup>1</sup> and the vanishment of OH-stretching bands [19].

**Figure 4.** (**a**) The temperature-dependent Raman spectra of the epsomite at ~1.3 GPa for the sulfate internal vibration; (**b**) the high-temperature Raman spectra of the epsomite at ~1.3 GPa for the hydroxyl-stretching vibration.

Figure 5 shows the temperature-dependent Raman spectra of epsomite collected at conditions of 293–623 K and at a higher pressure of ~6.3 GPa. At temperatures below ~553 K, the mode of ν<sup>1</sup> (SO4) remained at an almost constant wave number, with an increasing temperature (Figure 5a). However, when the temperature was above ~553 K, this mode suddenly split into two new Raman peaks at 997.4 and 1020.0 cm<sup>−</sup>1. In addition, the Raman spectra for the OH-stretching modes were also greatly changed after ~553 K, as shown in Figure 5b. These discontinuous variations in the Raman modes of the sulfate tetrahedra and water molecules are associated with the dehydration reaction from the epsomite to magnesium sulfate trihydrate (MgSO4·3H2O). The high-temperature Raman spectra of the epsomite at another two pressure points of ~0.8 and ~3.7 GPa are displayed in Figure 6. At ~0.8 GPa, the epsomite was found to undergo thermal dehydration reactions at temperatures of ~313, ~353 and ~423 K, respectively. For the dehydration process at ~3.7 GPa, the epsomite initially dehydrated to hexahydrite at ~373 K, and then transformed to magnesium sulfate trihydrate at ~453 K.

**Figure 5.** (**a**) The Raman spectra of the epsomite at 293–623 K and ~6.4 GPa in the sulfate vibrational range; (**b**) the collected high-temperature Raman spectra of the epsomite at ~6.4 GPa in the hydroxyl-stretching vibrational range.

**Figure 6.** (**a**) The high-temperature Raman spectra of the epsomite at a pressure of ~0.8 GPa; (**b**) the Raman spectra of the epsomite with an increasing temperature at ~3.7 GPa.

On the basis of the data of dehydration temperatures obtained from the high-temperature and high-pressure Raman scattering experiments, a *P*-*T* phase diagram for epsomite was well established at a wide temperature range of 293–623 K and at pressures of 0.8–6.3 GPa. As shown in Figure 7, it is verified that the dehydration temperature of epsomite gradually increased with the rise in pressure. This positive relation between dehydration temperature and pressure was coincident with the result of other similar hydrated sulfates (gypsum, chalcanthite and blödite) [11–13]. In addition, we also found that the epsomite undergoes three-step dehydration reactions in the heating process: from epsomite (MgSO4·7H2O) to hexahydrite (MgSO4·6H2O) to magnesium sulfate trihydrate (MgSO4·3H2O) to anhydrous magnesium sulfate (MgSO4). The anhydrous MgSO4 has been studied and it is stable up to 17.5 GPa [26]. The dehydration sequence of the epsomite revealed in this study is different from the result reported by Brotton et al. [27]. They observed the direct formation of magnesium sulfate trihydrate without the occurrence of hexahydrite in the heated epsomite at room pressure. This discrepancy was possibly caused by some crucial factors, such as the heating rate, relative humidity and pressure condition, all of which have been demonstrated to significantly influence the dehydration process of hydrous minerals [9,28–30]. Furthermore, three dehydration boundaries were determined by the linear fitting of these data. The epsomite-hexahydrite dehydration boundary was determined to be *P* (GPa) = –14.645 + 0.048 *T* (K), and the dehydration boundary for the transformation from hexahydrite to magnesium sulfate trihydrate corresponded to *P* (GPa) = –9.172 + 0.028 *T* (K). As for the last dehydration reaction, from magnesium sulfate trihydrate to anhydrous magnesium sulfate, the corresponding boundary can be described by a linear equation: *P* (GPa) = –3.430 + 0.010 *T* (K).

**Figure 7.** The established *P*-*T* phase diagram for the epsomite in temperature ranges of 293-623 K and pressure ranges of 0.8-6.4 GPa.

#### **4. Conclusions**

In summary, we performed in situ Raman spectroscopic and electrical conductivity measurements on epsomite to investigate its phase transition and dehydration process under high pressure. The results from room-temperature and high-pressure experiments demonstrated a structural phase transition for the epsomite, occurring at ~5.1 GPa, and the hydrogen bonding in the structure of the epsomite is strengthened after this transformation. The high-temperature and high-pressure Raman measurements were conducted at four pressure points (~0.8, ~1.3, ~3.7 and ~6.4 GPa) for the purpose of revealing the influence of pressure on the dehydration process of the epsomite. In the heating process of the epsomite at ~1.3 GPa, we successively observed the appearance of hexahydrite, magnesium sulfate trihydrate and anhydrous magnesium sulfate at ~343, ~373 and ~473 K. At ~6.4 GPa, the epsomite dehydrated to magnesium sulfate trihydrate at ~553 K. Therefore, it can be concluded that the epsomite undergoes a three-step dehydration reaction under a high pressure, and its dehydration temperature gradually increases with the pressure.

**Author Contributions:** Conceptualization, L.D. and L.Y.; methodology, L.D., Y.H. and H.L.; software, L.D.; validation, L.Y., M.H. and X.Z.; formal analysis, L.Y. and L.D.; investigation, L.D.; data curation, L.Y.; writing—original draft preparation, L.Y. and L.D.; writing—review and editing, L.D.; visualization, L.Y.; supervision, L.D., Y.H. and H.L.; project administration, L.D.; funding acquisition, L.D. Methodology, supervison, H.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the strategic priority Research Program (B) of the Chinese Academy of Sciences (18010401), Key Research Program of Frontier Sciences of CAS (QYZDB-SSW-DQC009), Hundred Talents Program of CAS, NSF of China (41774099 and 41772042), Youth Innovation Promotion Association of CAS (2019390), Special Fund of the West Light Foundation of CAS, and Postdoctoral Science Foundation of China (2018M643532).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Topological Equivalence of the Phase Diagrams of Molybdenum and Tungsten**

**Samuel Baty 1,† , Leonid Burakovsky 2,\*,† and Dean Preston 3,†**


Received: 15 November 2019; Accepted: 26 December 2019; Published: 2 January 2020

**Abstract:** We demonstrate the topological equivalence of the phase diagrams of molybdenum (Mo) and tungsten (W), Group 6B partners in the periodic table. The phase digram of Mo to 800 GPa from our earlier work is now extended to 2000 GPa. The phase diagram of W to 2500 GPa is obtained using a comprehensive ab initio approach that includes (i) the calculation of the *T* = 0 free energies (enthalpies) of different solid structures, (ii) the quantum molecular dynamics simulation of the melting curves of different solid structures, (iii) the derivation of the analytic form for the solid–solid phase transition boundary, and (iv) the simulations of the solidification of liquid W into the final solid states on both sides of the solid–solid phase transition boundary in order to confirm the corresponding analytic form. For both Mo and W, there are two solid structures confirmed to be present on their phase diagrams, the ambient body-centered cubic (bcc) and the high-pressure double hexagonal close-packed (dhcp), such that at *T* = 0 the bcc–dhcp transition occurs at 660 GPa in Mo and 1060 GPa in W. In either case, the transition boundary has a positive slope *dT*/*dP*.

**Keywords:** phase diagram, quantum molecular dynamics, melting curve, Z methodology, multi-phase materials

#### **1. Introduction**

The high-pressure (HP) and high-temperature (HT) behavior of materials is of great importance for condensed matter physics and geophysics as well as for technological applications. In particular, the study of melting under extreme conditions is a subject of special interest. Among other topics, the HP melting of transition metals has been one of the main areas of research. Metals like tantalum (Ta), molybdenum (Mo), and nickel (Ni) have been studied both theoretically [1–10] and experimentally [11–25]. Research interest in HP-HT polymorphism in the body-centered cubic (bcc) transition metals has reemerged in connection with laser-heated diamond anvil cell (DAC) melting experiments in which melting curves with a small slope (*dT*/*dP*, where *T* stands for temperature and *P* stands for pressure) in the megabar pressure range have been determined [11–18]. These flat melting curves do not agree with the results of more recent experiments [19–25] and calculations [2,4,6,8–10]. We note that, in all the aforementioned experiments, *P* and *T* did not exceed 200 GPa and 5000 K, respectively. A pressure of ∼200 GPa was reached in the experimental study of Reference [21] and a temperature of ∼5000 K in those of both Reference [21] and Reference [24].

Several hypotheses have been proposed to explain apparent discrepancies between experiment and theory with regard to the flat melting curves. One of them suggests that the flat melting curves could in fact correspond to solid–solid (s–s) phase boundaries that occur under HP-HT below melting. In particular, HP-HT solid–solid transitions have been predicted to take place in Ta [3,5,9] and Mo [1]. In contrast to these metals, there are only a few studies on tungsten (W). The experimental melting curve of W is one of the flat ones [13]; it has the initial slope of ∼7 K/GPa and reaches ∼4000 K at ∼100 GPa, which contradicts both the initial slope of 44 K/GPa from isobaric expansion measurements [26] and the shock Hugoniot melting points of ∼12,000 K at ∼400 GPa [27]. The experimental melting curves of References [28–30] have even higher slopes of ∼90, 75, and 60 K/GPa, respectively. As noted in Reference [30], large experimental errors of Reference [31] lead to ∼50% error in the value of the corresponding slope, 70 ± 35 K/GPa. The other theoretical value of the slope is 28.7 K/GPa [32,33].

It is well known that some physical properties related to the phase diagrams of elements can be systematized across the corresponding groups or rows (or both) of the periodic table. Among those are the ambient melting point *Tm*(0); the initial slopes of melting curve *dTm*/*dP*(0); and the ambient volume change at melt, Δ*Vm*(0), and latent heat of melting, Δ*Hm*(0). For example, elements of group 7B, Fe, Ru, and Os, exhibit a very clear systematics with regard to *Tm*(0) and *dTm*/*dP*(0) : *Tm*(0) is, respectively, (in K) 1811, 2607, and ∼3400 [34], and *dTm*/*dP*(0) is (in K/GPa) 30.0 [21], 40.7 ± 2.2 [35], and 49.5 [36]. That is, *Tm*(0) is separated by ∼800 K and *dTm*/*dP*(0) is separated by ∼10 K/GPa. The analytic forms of the corresponding melting curves can also be systematized with regard to the numerical parameters involved: 1811 (<sup>1</sup> + *<sup>P</sup>*/30.2)0.50 [21], 2607 (<sup>1</sup> + *<sup>P</sup>*/(32.75 ± 1.75))0.51 [35], and 3370 (1 + *P*/36.1)0.53 [36]. The proper knowledge of such systematics and similarities in material properties allows one to make predictions and to validate new results. For instance, the melting curve systematics for the elements of the third row of the periodic table discussed in Reference [10] allowed to predict the yet unknown melting curve of Hf, which the preliminary theoretical results seem to confirm [10].

The phase diagrams themselves can be topologically similar (look-alike) or, in some cases, even topologically equivalent. The topological similarity of the phase diagrams of Ti, Zr, and Hf is well known (more detail will follow). The phase diagrams of Si and Ge are topologically equivalent at low *P* : both contain semiconducting diamond and metallic tetragonal *β*-Sn solid structures, and the slopes of the corresponding solid–solid and solid–liquid phase boundaries are almost identical [37]; in either case, *dTm*/*dP*(0) < 0. The proper knowledge of similarities of the phase diagrams can be useful in making predictions with regard to the phase diagram content and in offering suggestions as to what solid structures to look for in high-*PT* experiments.

The study of Frohberg [38] reveals that, for the elements of the three groups 4B, 5B, and 6B, namely, Ti, Zr, and Hf; V, Nb, and Ta; and Cr, Mo, and W, respectively, the corresponding *Tm*(0) and Δ*Hm*(0) obey very clear systematics. In addition, the phase diagrams of 4B's Ti, Zr, and Hf are topologically similar; their only difference is in the slopes of the solid–solid transition boundaries between the three solid structures found on each of the three phase diagrams: hexagonal close-packed (hcp), hexagonal omega (hex-omega), and bcc. Since the phase diagrams of 5B's Nb and Ta may also be topologically similar, as we discuss below, it is natural to assume the topological similarity of the phase diagrams of Mo and W since they both belong to 6B, the last of these three groups.

Here, we present an extensive theoretical study of the phase diagram of W to a pressure of 2500 GPa (25 Mbar). We demonstrate that the phase diagram of W is topologically equivalent to that of Mo; the latter was considered in our previous publication [39] and is now extended to 2000 GPa (20 Mbar). As it turns out, the critical features of the two phase diagrams that make them equivalent, that is, the location of the solid–solid–liquid triple point and the analytic form of the solid–solid phase boundary, are revealed at much higher pressure that goes far beyond the *P* and *T* ranges of ∼0–200 GPa and ∼0–5000 K, respectively, that all the experimental studies on transition metals mentioned above have covered. Specifically, the solid–solid phase boundary between the two solid phases that are identical in both cases lies in a *P* range of ∼650–1050 GPa for Mo and of ∼1050–1700 GPa for W. The corresponding *T* ranges are, respectively, ∼0–12,300 and ∼0–23,700 K. We therefore expect that, in the near future, our findings may not be verified in experiment; they will remain

theoretical predictions until pressures of ∼1000 GPa (10 Mbar) and temperatures of ∼10,000 K become experimentally attainable.

#### **2. The Phase Diagram of Mo**

The construction of the phase diagram of Mo using state-of-the-art theoretical techniques was discussed in detail in Reference [39], where we presented the phase diagram of Mo to 800 GPa. Here, we extend this phase diagram to 2000 GPa by calculating the melting curve of the double hexagonal close-packed (dhcp) structure of Mo, which is the stable solid phase of Mo at *P* ≥ 660 GPa [39], and by determining the solid–solid phase transition boundary between the body-centered cubic (bcc) structure, which is the stable solid phase of Mo at *P* ≤ 660 GPa, and the dhcp structure. Our results can be summarized as follows:

(i) the melting curve of bcc-Mo in the Simon–Glatzel form [40] (*Tm* in K, *P* in GPa):

$$T\_{\rm m}(P) = 2896 \left( 1 + \frac{P}{36.6} \right)^{0.43},\tag{1}$$

(ii) the melting curve of dhcp-Mo in the Simon–Glatzel form (*Tm* in K, *P* in GPa):

$$T\_m(P) = 1880 \left( 1 + \frac{P}{24.5} \right)^{0.50},\tag{2}$$

Our melting curve of bcc-Mo, Equation (1), is in excellent agreement with another theoretical calculation, *Tm*(*P*) = <sup>2894</sup>(<sup>1</sup> + *<sup>P</sup>*/37.2)0.433 [41]. For Equation (1), *dTm*/*dP*|*P*=<sup>0</sup> = 34.0 K/GPa, which is in excellent agreement with 32 K/GPa [26] or 34 ± 6 K/GPa [42], each from isobaric expansion measurements.

Figures 1 and 2 demonstrate the time evolution of *T* and *P*, respectively, in the *Z* method runs of the bcc-Mo melting point (*P* in GPa, *T* in K) (*P*, *T*) = (1185, 12,770), which is one of the melting points that we calculated in the course of the present study. Figures 3 and 4 demonstrate the same for the dhcp-Mo melting point (1186, 13,430). These two points are chosen as examples and are shown in Figure 5 as open blue and green circles, respectively.

**Figure 1.** Time evolution of temperature for body-centred cubic (bcc) Mo in three QMD runs with initial temperatures (T0s) separated by 625 K: The middle run is the melting run, during which *T* decreases from ∼17,000 K for the superheated state to ∼12,800 K for the liquid at the corresponding melting point.

The melting curves of bcc-Mo and dhcp-Mo cross each other at (*P* in GPa, *T* in K) (*P*, *T*) = (1045, 12,320) which is the bcc–dhcp-liquid triple point. The choice of *<sup>T</sup>*(*P*) = *<sup>a</sup>*(<sup>660</sup> − *<sup>P</sup>*)*b*, 0 < *<sup>b</sup>* < 1 as a functional form for the bcc–dhcp phase boundary leads to

(iii) the bcc–dhcp solid–solid phase transition boundary:

$$T(P) = 58.5 \left(P - 660\right)^{0.90} \text{\textdegree} \tag{3}$$

which crosses the triple point as well as the point (*P*, *T*)=(725, 2500) that comes from the solidification simulations using the inverse Z method; for more detail see Reference [39].

The phase diagram of Mo to 2000 GPa is shown in Figure 5. It includes the two melting curves, the solid–solid phase boundary, and the results of the solidification of liquid Mo into the final states of either solid bcc or solid dhcp using the inverse Z method. Figure 5 actually represents the extended version of Figure 11 of Reference [39] which covered a *P* range of 0–800 GPa.

**Figure 2.** The same as in Figure 1 for the time evolution of pressure (in kbar; 10 kbar = 1 GPa): During melting, *P* increases from ∼1180 GPa for the superheated state to ∼1185 GPa for the liquid state at the corresponding melting point.

**Figure 3.** The same as in Figure 1 for double hexagonal close-packed (dhcp) Mo.

**Figure 4.** The same as in Figure 2 for dhcp-Mo.

**Figure 5.** Ab initio phase diagram of molybdenum.

#### **3. The Phase Diagram of W**

Although the phase diagram of W has not been known in detail, there are several features of this phase diagram that have been firmly established. First, melting on the shock Hugoniot of W occurs at ∼400 GPa and presumably at *T* ∼12,000 K [27]; these values are extracted from the corresponding *P* = *P*(*Up*) and *T* = *T*(*Up*) dependences on the particle velocity *Up* as those at *Up* ∼ 2.5 km/s at which melting occurs on the shock Hugoniot of W [27]. Second, the stability of bcc-W has been confirmed experimentally under room-*T* isothermal compression to 420 GPa [43] and to 500 GPa [44]. There is presently no experimental evidence for an s–s transition in W. However, there is compelling theoretical evidence for an s–s phase transition at high *P*. Calculations of the phonon spectra of both the bcc and face-centered cubic (fcc) structures of W show that bcc-W becomes mechanically unstable at pressures above 12 Mbar while fcc-W being mechanically unstable at low *P* becomes fully stabilized above 4 Mbar [45]. A very recent theoretical study [46] demonstrates that fcc-W becomes mechanically stable at ∼450 GPa while bcc-W becomes mechanically unstable with increasing *P*; above 12 Mbar, bcc-W is unstable at low *T* but remains stable above 1000 K. Hence, an s–s phase transformation to another solid phase that is mechanically stable is expected to occur at *<sup>P</sup>* <sup>&</sup>lt; ∼ 12 Mbar at low *T* and the s–s transition boundary is expected to have positive slope (*dT*/*dP* > 0) since the bcc-W stability range widens with increasing *T*. Because fcc-W also becomes thermodynamically more favorable than bcc-W [46], fcc-W is one of the candidates for the high-*P* solid structure of W. In fact, the bcc-fcc s–s phase transition was predictied to occur at ∼12 Mbar and calculations show that the bcc-fcc phase boundary does have positive slope at *P* ≥ 12 Mbar [46]. However, in Reference [47], a transformation to a different solid structure, namely, dhcp, which is thermodynamically more stable than fcc, is predicted to occur, albeit at much lower pressure of 650 GPa. If a bcc–dhcp s–s phase transition does occur in W, just like in Mo, the two phase diagrams may look similar. As a matter of fact, as we demonstrate in what follows, the two phase diagrams are topologically equivalent.

We begin our theoretical study of the phase diagram of W with the calculation of the cold (*T* = 0) free energies (i.e., enthalpies) of a number of different solid structures of W as a function of *P*.

#### *3.1. Cold Enthalpies of Different Solid Structures of W*

The calculations were based on density-functional theory (DFT) with the projector-augmented wave (PAW) [48] implementation and the generalized gradient approximation (GGA) for exchange-correlation energy in the form known as Perdew–Burke–Ernzerhof (PBE) [49]. All the calculations were done using VASP (Vienna Ab initio Simulation Package). Since the simulations were performed at high-*PT* conditions, we used accurate pseudopotentials where the semi-core 5s and 5p states were treated as valence states. Specifically, W was modeled with 14 valence electrons per atom (5s, 5p, 5d, and 6s orbitals). The valence electrons were represented with a plane-wave basis set with a cutoff energy of 500 eV, while the core electrons were represented by projector-augmented wave (PAW) pseudopotentials. The core radius (the largest value of RCUTs among those for each of the quantum orbitals) of this W pseudopotential is 2.3 a.u. or 1.217 Å. Since numerical errors in the calculations using VASP will remain almost negligible until the nearest neighbor distance reaches 2 × RCUT/(1.25 ± 0.05) [39], with this pseudopotential, one can study systems with densities up to ∼80 g/cm3 (a pressure of ∼8000 GPa).

Cold enthalpies were calcualted using unit cells with very dense *k*-point meshes (e.g., 50 × 50 × 50 for bcc-W) for high accuracy. In all the non-cubic cases, we first relaxed the structure to determine its unit cell parameters at each volume. Tight convergence criteria for the total energy (10−<sup>5</sup> meV/atom) and structural relaxation (residual forces < 0.1 meV/Å and stresses ≤ 0.1 kbar) were employed.

Our calculated cold enthalpies of five different solid structures of W are shown in Figure 6. It is seen that, with inceasing *P*, a number of other solid structures become thermodynamically more favorable than bcc, but it is dhcp that does it first, at a pressure of 1060 GPa, and remains the most favorable solid structure of W at higher *P*. Hence, in the case of W, the s–s transition boundary is the bcc–dhcp one and its starting point is (*P*, *T*)=(1060, 0).

**Figure 6.** The *T* = 0 free energies of different solid structures of W (listed in the legend) from ab initio calculations using VASP (Vienna Ab initio Simulation Package): The free energy of bcc-W is taken to be identically zero.

#### *3.2. Equations of State of bcc-W and dhcp-W*

Next, we calculate the cold equations of state of both bcc-W and dhcp-W. For dhcp-W, we also determine the density dependence of the *c*/*a* ratio, where *a* and *c* are the lattice constants of the dhcp unit cell.

Our ab initio results on the cold equation of state (EOS) of bcc-W are described by the third-order Birch–Murnaghan (BM3) form:

$$P(\rho) = \frac{3}{2} B\_0 \left( \left( \frac{\rho}{\rho\_0} \right)^{7/3} - \left( \frac{\rho}{\rho\_0} \right)^{5/3} \right) \left[ 1 + \frac{3}{4} (B\_0' - 4) \left( \left( \frac{\rho}{\rho\_0} \right)^{2/3} - 1 \right) \right],\tag{4}$$

where *B*<sup>0</sup> and *B* <sup>0</sup> are the values of the bulk modulus and its pressure derivative at the reference point *ρ* = *ρ*0. In our case of bcc-W,

$$
\rho\_0 = 19.3 \text{ g/cm}^3, \quad B\_0 = 314.1 \text{ GPa}, \quad B\_0' = 4.07. \tag{5}
$$

Since the *P* = 0 values of the density of W at *T* = 0 and 300 K differ by ∼0.3% (19.31 vs. 19.26 g/cc) and *T* = 300 K introduces a negligibly small thermal pressure correction, the *T* = 0 and *T* = 300 K isotherms can be described by the same values of *B*<sup>0</sup> and *B* <sup>0</sup>. Consequently, we can compare room-*T* isotherm data of References [44,50–53] to our *T* = 0 isotherm as determined with VASP. A comparison is shown in Figure 7. It is seen that our EOS is virtually identical to the experimental isotherm of

Reference [53], for which *B*<sup>0</sup> = 317.5 GPa and *B* <sup>0</sup> = 4.05 are very similar to ours, as well as to the theoretical calculations of Reference [51] in the so-called mean-field potential (MFP) approach.

For dhcp-W, the *P* dependence of the *c*/*a* ratio for the lattice constants of the dhcp unit cell is accurately described by

$$\frac{c}{a} = 3.2660 - 0.4912 \left( \frac{a}{a\_0} \right)^3 + 0.7484 \left( \frac{a}{a\_0} \right)^5,\tag{6}$$

where *a*<sup>0</sup> = 2.8525 Å coresponds to the ideal (*c*/*a* = 4 <sup>√</sup>2/3 <sup>≈</sup> 3.2660) unit cell of the same reference density *ρ*<sup>0</sup> as that for the corresponding cold EOS. Since the unit cell volume is *ca*<sup>2</sup> <sup>√</sup>3/8 (*a*3/ <sup>√</sup><sup>2</sup> for the ideal structure), the above relation can be translated into the *P* dependence of *c*/*a* via the corresponding cold EOS of dhcp-W that we obtained. It is described by the BM3 form, Equation (4), with

$$
\rho\_0 = 18.6 \text{ g/cm}^3, \quad B\_0 = 290.8 \text{ GPa}, \quad B\_0' = 3.95. \tag{7}
$$

Above 1000 GPa, the dhcp-W structure is virtually ideal. At the transition pressure of 1060 GPa, the two density values predicted by these EOSs are (in g/cm3) 40.38 for bcc-W and and 40.53 for dhcp-W; hence, the bcc–dhcp transition corresponds to a small volume change of ∼0.4%.

**Figure 7.** The *T* = 0 equation of state of bcc-W: our own ab initio calculations using VASP vs. experimental data (Dubrovinsky et al. [44]; Ruoff et al. [50]; Dewaele et al. [52]; and Mashimo et al. [53]) and other theoretical calculations (Wang et al. [51]).

#### **4. Melting Curves of Different Solid Structures of W**

We now discuss the calculation of the melting curves of different solid structures of W. For this calculation, we used the Z method, which is described in detail in References [9,36,54]. We calculated the melting curves of bcc-W as well as a number of other solid structures of W which have been mentioned in the literature in connection with transition metals: all the close-packed structures with different layer stacking (fcc, hcp, dhcp, thcp, and 9R), open structures (simple cubic, A15, and hex-*ω*), and different orthorhombic structures (Pnma, Pbca, Pbcm, Cmca, Cmcm, etc.). Again, in all the non-cubic cases, we relaxed the structure to determine its unit cell parameters; those unit cells were used for the construction of the corresponding supercells. We used systems of 400–500 atoms in each case.

#### *4.1. Ab Initio Melting Curve of bcc-W*

For the calculation of the melting curve of bcc-W (with the Z method), we used a 432-atom (6 <sup>×</sup> <sup>6</sup> <sup>×</sup> 6) supercell with a single <sup>Γ</sup>-point. Full energy convergence (to <sup>&</sup>lt; ∼ 1 meV/atom) was verified by performing short runs with 2 × 2 × 2 and 3 × 3 × 3 *k*-point meshes and by comparing their output with that of the run with a single Γ-point. For each of the five points listed in Table 1, we performed ten *NVE* runs (five values of *V* correspond to five values of densities in Table 1) of 10,000–20,000 time steps of 1 fs each, with an increment of the initial *T* of 250 K for the first datapoint and 625 K for the remaining four. Since the error in *T* is half of the increment [36], the *T* errors of our five values of *Tm* are within 5% each. The *P* errors are negligibly small: less than 1 GPa for the first point and 1–2 GPa for the remaining four. Hence, our melting results on bcc-W are very accurate.


**Table 1.** The five ab initio melting points of bcc-W, (*Pm*, *Tm* ± Δ*Tm*), obtained from the *Z* method implemented with VASP.

The best fit to the five bcc-W melting points gives the melting curve of bcc-W in the Simon–Glatzel form (*Tm* in K, *P* in GPa):

$$T\_m(P) = 3695 \left( 1 + \frac{P}{41.8} \right)^{0.50} \text{.} \tag{8}$$

Its initial slope, *dTm*(*P*)/*dP* = 44.2 K/GPa at *P* = 0, is in excellent agreement with 44 K/GPa from isobaric expansion measurements [26]. Both the five bcc-W melting points and the melting curve (Equation (8)) are shown in Figure 8 and compared to the experimental results in Reference [27,30], the calculated Hugoniot of W [55], and other theoretical calculations [6,56–58].

**Figure 8.** The melting curve of bcc-W: QMD simulations using VASP vs. other theoretical calculations (Xi an Cai [56]; SNL, 2008 [57]; and Liu et al., 2012 [6] and 2017 [58]), the low-pressure melting data (experiment in Reference [30]), and the experimental shock melting datapoint (Hugoniot melting [27]). The calculated Hugoniot is shown as a thin black curve.

#### *4.2. Ab Initio Melting Curve of dhcp-W*

The five melting points of dhcp-W that we obtained are listed in Table 2. The dhcp-W melting simulations were carried out exactly the same way as the bcc ones, namely, 10 runs per point, 10,000–20,000 time steps of 1 fs each per run, and an increment of the initial *T* of 300 K for the first datapoint and 625 K for the remaining four. Just as for bcc-W, a 432-atom (6 × 6 × 3) supercell was used to simulate each of the five dhcp-W melting points. Also, as for bcc-W, the *T* errors of our five values of *Tm* are within 5% each.


**Table 2.** The five ab initio melting points of dhcp-W, (*Pm*, *Tm* ± Δ*Tm*), obtained from the *Z* method implemented with VASP: The lattice constant values correspond to the ideal dhcp structure.

The best fit to the five dhcp-W points gives the melting curve of dhcp-W in the Simon–Glatzel form (*Tm* in K, *P* in GPa):

$$T\_m(P) = 2640 \left( 1 + \frac{P}{23.0} \right)^{0.51} \,\text{.}\tag{9}$$

Figures 9 and 10 demonstrate the time evolution of *T* and *P*, respectively, in the *Z*-method runs of the bcc-W melting point (*P* in GPa, *T* in K) (*P*, *T*)=(947, 18, 070). Figures 11 and 12 demonstrate the same for the dhcp-W melting point (970, 18, 000). These two points are chosen as examples and are shown in Figure 13 as open blue and green circles, respectively.

All the other solid structures that we considered melt below bcc. As an example, one of the fcc-W melting points, namely, (*P*, *T*) = (942, 16,920), along with a short segment of the corresponding melting curve, are shown in Figure 13. The corresponding time evolution of *T* and *P* in the *Z*-method runs are shown in Figures 14 and 15, respectively.

The melting curves of bcc-W and dhcp-W cross each other at (*P* in GPa, *T* in K) (*P*, *T*) = (1675, 23,680), which is the bcc–dhcp-liquid triple point. The choice of *<sup>T</sup>*(*P*) = *<sup>a</sup>*(<sup>1060</sup> − *<sup>P</sup>*)*b*, 0 < *<sup>b</sup>* < 1 as a functional form for the bcc–dhcp phase boundary leads to

(iii) the bcc–dhcp solid–solid phase transition boundary:

$$T(P) = 73.2 \left(P - 1060\right)^{0.90} \text{\AA} \tag{10}$$

which crosses the triple point and lies within the bounds imposed by the solidification simulations using the inverse Z method. We discuss these inverse Z simulations in the following section.

**Figure 9.** The same as in Figure 1 for bcc-W.

Although the rigorous derivation of the thermal equations of state of bcc-W and dhcp-W goes beyond the scope of this work, we note that the finite-*T* counterparts of the above two EOSs can be written approximately as *<sup>P</sup>*(*ρ*, *<sup>T</sup>*) = *<sup>P</sup>*(*ρ*) + *<sup>α</sup> <sup>T</sup>*, where *<sup>α</sup>*bcc = 7.3 × <sup>10</sup>−<sup>3</sup> and *<sup>α</sup>*dhcp = 6.8 × <sup>10</sup>−3. The resulting "approximate" thermal EOSs turn out to be quite accurate. For example, for the five bcc-W melting points in Table 1, the corresponding thermal EOS gives pressures of −17.8, 90.3, 258, 543, and 948, which are basically identical to those in the third column of Table 1. For the five dhcp-W

melting points in Table 2, the five *P* values are 7.9, 150, 509, 970, and 1476, in excellent agreement with those in the third column of Table 2.

**Figure 10.** The same as in Figure 2 for bcc-W.

**Figure 11.** The same as in Figure 1 for dhcp-W.

**Figure 12.** The same as in Figure 2 for dhcp-W.

**Figure 13.** Ab initio phase diagram of tungsten.

**Figure 14.** The same as in Figure 1 for fcc-W.

**Figure 15.** The same as in Figure 2 for fcc-W.

#### **5. Inverse Z Solidification Simulations of Liquid W**

To constrain the location of the bcc–dhcp solid–solid phase boundary on the *P*–*T* plane between the points (*P* = 1060, *T* = 0) and (*P* = 1675, *T* = 23,680), we carried out two sets of independent inverse Z runs [9] to solidify liquid W and to confirm that liquid W solidifies into bcc on one side of this boundary and into dhcp (or any other solid structure) on the other side, such that the location of this phase boundary may be constrained. Since the solidification kinetics is approximately governed by the factor exp{Δ*F*/*T*∗}, where Δ*F* ≡ *Fl* − *Fs* is the liquid–solid free-energy difference at the solidification temperature *T*∗, in the case of several energetically competitive solid phases, the most thermodynamically favorable solid phase has the largest Δ*F* and is therefore the fastest to solidify. Hence, the inverse Z method yields the most stable solid phase at a given (*P*, *T*).

For our inverse Z simulations, we used a computational cell of 512 atoms prepared by melting a 8 × 8 × 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.0, 1.935, 1.915, 1.895, and 1.870 Å; the dimensions of bcc unit cells having the same volume as the sc ones are 2.520, 2.438, 2.413, 2.388, and 2.356 Å, respectively, which corresponds to the bcc-W pressures of ∼870, 1225, 1355, 1505, and 1715 GPa and to slightly lower pressures for dhcp-W.

We carried out *NVT* simulations using the Nosé–Hoover thermostat with a timestep of 1 fs, with the initial *T* increment of 2500 K. Complete solidification typically was required from 15 to 25 ps or 15,000–25,000 timesteps. The inverse Z runs indicate that liquid W only solidifies into bcc at ∼900 GPa in the whole temperature range from 0 to essentially the corresponding *Tm*. However, at ∼1200–1400 GPa, it solidifies into bcc above the transition boundary in Figure 13, 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 seven structures (seven green bullets in Figure 13 in the 1200–1400 GPa range) 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 ∼1200 GPa above the transition boundary are shown in Figure 16 and of those solidified below the transition boundary are shown in Figure 17. The 7500 K state virtually lies on the boundary. We tentatively assign it to bcc because it definitely has features of bcc (RDF peaks at *R* ∼ 55, 65, and 85, etc.). At the same time, it certainly has some features that are both uncharacteristic of bcc (e.g., the disappearance of the bcc peak at *R* ∼ 95) and characteristic of dhcp (peaks at *R* ∼ 90 and 120, small peaks at *R* ∼ 100 and 130, etc.). Most likely, this 7500 K state is bcc with some admixture of dhcp.

A few more comments are in order. The 17,500 K state at ∼1250 GPa did not solidify, most likely for the reason of not being supercooled enough to intiate the solidification process [9]. Indeed, 17,500 K constitutes ∼0.8 of the corresponding *Tm* of ∼22,000 K (20% of supercooling) while, e.g., for the other set of points at ∼900 GPa, the highest solidification *T* of 12,500 K constitutes ∼0.75 of the corresponding *Tm* of ∼17,000 K (25% of supercooling), which apparently allows for the solidification process to go through in this case.

The phase diagram of W to 2500 GPa is shown in Figure 13. It includes the two bcc and dhcp melting curves, the bcc–dhcp solid–solid phase transition boundary, and the results of the solidification of liquid W into final states of either solid bcc or solid dhcp using the inverse Z method.

**Figure 16.** Radial distribution functions (RDFs) of the final states of the solidification of liquid W at ∼1200 GPa at higher temperatures.

Both the Mo and W phase diagram figures, Figures 5 and 13, do not show the previous experimental DAC melting curves of Errandonea et al. [13,16] and are both flat, with the initial slope of ∼7–8 K/GPa. We did not include those in the two figures because we do not consider them to be relevant. The very recent experimental study by Hrubiak et al. [23] 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. [13]. Hence, most likely, Errandonea's curve is an intermediate DAC transition boundary, while our curve as well as Alfe's [41], which is basically identical to ours, is the actual melting curve of Mo.

**Figure 17.** RDFs of the final states of the solidification of liquid W at ∼1200 GPa at lower temperatures.

A behavior similar to that of Mo (a texture change) was recently observed in the high-*PT* melting experiments on V [25]. Hence, this phenomenon of a texture (microstructure) change can be common to a number of transition metals, including W. It is therefore natural to assume that, in the case of W, too, the corresponding Errandonea's melting curve is an intermediate DAC transition boundary while our curve is the actual melting curve of W.

#### **6. Concluding Remarks**

The phase diagrams of Mo and W, Group 6B partners in the periodic table, shown in Figures 5 and 13, respectively, are topologically equivalent. In both cases, the ambient bcc solid structure transforms into dhcp on increasing *P*. Only these two solid phases are confirmed as being present on the two phase diagrams, via ab initio QMD simulations using the Z methodology. Even the critical aspect ratio of the two phase diagrams is identical. Indeed, the ratio of the pressure values for the bcc–dhcp-liquid triple point and the *T* = 0 bcc–dhcp transition (it determines the curvatures and the locus of the bcc–dhcp transition boundary on the *P*–*T* plnae) is the same in both cases: 1045/660 = 1.583 for Mo and 1675/1060 = 1.580 for W. The two transition boundaries, Equations (3) and (10), are described by the same exponent 0.9 which further confirms the topological equivalence of the two phase diagrams.

It is interesting to note another case of the phase diagrams of Nb and Ta, group 5B partners in the periodic table, which appear to be topologically similar (look-alike) but not equivalent. In both cases, a bcc-Pnma solid–solid phase transition occurs [9,59], but in contrast to dhcp in Mo and W, the orthorhombic Pnma phase exists at high *T* only. In Ta, the bcc-Pnma-liquid triple point is at a pressure which is an order of magnitude higher than that in Nb [9,59]. Also, the bcc-Pnma transition *T* increases with *P* in Nb but decreases with *T* in Ta.

Inverse Z solidification simulations confirm the phase diagrams of Mo and W shown in Figures 5 and 13, respectively, and validate the corresponding bcc–dhcp solid–solid phase boundaries. Additional support for the existence of the high-*P* dhcp structure in both Mo and W comes from Reference [60]. The electronic wave functions of Mo and W and of their group 5B partner Cr are

such that, with decreasing volume (increasing *P*), they become more spherically symmetric, which favors a higher atomic coordination number than 8 as that for bcc and, correspondingly, a more symmetric structure than bcc. It is well known that the highest possible coordination number is 12, and it corresponds to all of the hexagonal polytypes, namely, hcp (stacking sequence *AB*), fcc (*ABC*), dhcp (*ABAC*), triple-hcp (*ABCACB*), etc. It is therefore natural to expect that, at high *P*, both Mo and W and perhaps Cr as well will have one of the hexagonal polytypes as their equilibrium solid structure. Our study demonstrates that this polytype is in fact dhcp. It is interesting to note that, in Cr, as theoretical calculations show [61], hcp becomes energetically more favorable than bcc with decreasing volume (increasing *P*). Since different hexagonal polytypes are usually close to each other energetically, one can therefore expect them to be the high-*P* equilibrium solid structure of Cr. Thus, it is quite possible that the phase digram of Cr is topologically similar, if not equivalent, to those of both Mo and W. In contrast to both Mo and W, the study of Reference [62] shows that the electronic wave functions of V, Nb, and Ta are such that, with decreasing volume, they become less spherically symmetric, which favors a lower atomic coordination number than bcc's 8 and, correspondingly, a less symmetric structure than bcc, e.g., an orthorhombic one. This appears to be in agreement with the fact that, in both Nb and Ta, the high-*P* equilibrium solid structure is orthorhombic Pnma.

Fast recrystallization observed in HP-HT experiments on transition metals may imply either a microstructural transition (discussed above) at which the sample texture changes but its crystal structure remains the same or a true solid–solid phase transition. The former was recently proven to be the case in Mo [23], and it is very likely the case in V [25] as well as in W. The latter is apparently the case in Nb [59] and likely in Re as well [10]. In either case, the corresponding fast recrystallization lines have been misinterpreted as flat melting curves in the former laser-heating DAC experiments.

**Author Contributions:** All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The work was done under the auspices of the US DOE/NNSA. The QMD simulations were performed on the LANL clusters Pinto and Badger as parts of the Institutional Computing projects w17\_molybdenum, w18\_meltshear, and w19\_thermoelast.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
