**Fe Melting Transition: Electrical Resistivity, Thermal Conductivity, and Heat Flow at the Inner Core Boundaries of Mercury and Ganymede**

#### **Innocent C. Ezenwa** † **and Richard A. Secco \***

Department of Earth Sciences, University of Western Ontario London, London, ON N6A5B7, Canada

**\*** Correspondence: secco@uwo.ca

† Now at Institute for Planetary Materials, Okayama University, Misasa, Tottori 682-0193, Japan.

Received: 10 June 2019; Accepted: 12 July 2019; Published: 15 July 2019

**Abstract:** The electrical resistivity and thermal conductivity behavior of Fe at core conditions are important for understanding planetary interior thermal evolution as well as characterizing the generation and sustainability of planetary dynamos. We discuss the electrical resistivity and thermal conductivity of Fe, Co, and Ni at the solid–liquid melting transition using experimental data from previous studies at 1 atm and at high pressures. With increasing pressure, the increasing difference in the change in resistivity of these metals on melting is interpreted as due to decreasing paramagnon-induced electronic scattering contribution to the total electronic scattering. At the melting transition of Fe, we show that the difference in the value of the thermal conductivity on the solid and liquid sides increases with increasing pressure. At a pure Fe inner core boundary of Mercury and Ganymede at ~5 GPa and ~9 GPa, respectively, our analyses suggest that the thermal conductivity of the solid inner core of small terrestrial planetary bodies should be higher than that of the liquid outer core. We found that the thermal conductivity difference on the solid and liquid sides of Mercury's inner core boundary is ~2 W(mK)−1. This translates into an excess of total adiabatic heat flow of ~0.01–0.02 TW on the inner core side, depending on the relative size of inner and outer core. For a pure Fe Ganymede inner core, the difference in thermal conductivity is ~7 W(mK)<sup>−</sup>1, corresponding to an excess of total adiabatic heat flow of ~0.02 TW on the inner core side of the boundary. The mismatch in conducted heat across the solid and liquid sides of the inner core boundary in both planetary bodies appears to be insignificant in terms of generating thermal convection in their outer cores to power an internal dynamo suggesting that chemical composition is important.

**Keywords:** melting transition; Fe; electrical resistivity; thermal conductivity; high pressure; heat flow; thermal and chemical convection

#### **1. Introduction**

The processes of magnetic field generation and sustainability in planetary bodies depend on the composition and thermal state of their cores. Among the rocky planetary bodies with an active dynamo, Mercury has the weakest internally generated magnetic field, with a surface field strength of ~0.3 μT or ~1% compared with the Earth's field. Though a possible remnant crustal magnetization has been suggested [1], a self-sustained dynamo in Mercury's Fe core is the most plausible source of its global magnetic field [2]. A recent study suggests that a double-diffusive convective regime operates, where both thermal and compositional convection drive the system [3]. Earth-based radar measurements of subtle deviations from the mean resonant spin rate of Mercury demonstrate that Mercury's mantle is decoupled from its liquid or partially molten core [4–7]. This supports earlier assertions that Mercury has a molten outer core [8,9]. Recent geodetic constraints on the interior of Mercury from the MESSENGER spacecraft are consistent with a high degree of internal differentiation and a solid inner core with a radius of 0.4–0.7 times the outer core radius [10]. While the role of a solid inner core and its contribution to chemical composition convection in a liquid outer core was recognized long ago for Mercury [9], the possibility of Mercury's weak surface magnetic field resulting from dynamo action in a thin shell geometry has been shown more recently [11]. For the case of Ganymede, although remnant magnetization cannot be completely ruled out as the source of its magnetic field, magneto-convection in its core like that in the Earth has been suggested [12]. However, dynamo action in Ganymede differs from that of the Earth due to the presence of the strong nearby Jovian magnetic field. Thus, the magneto-hydrodynamic equation is variant under the transformation of *B* → *–B* with a directional preference for the self-generated field which could lead to a non-reversing magnetic field unlike the reversing nature of the geomagnetic field [13].

Convection in a terrestrial planetary core can arise from heat transport in excess of the conducted heat (i.e., by thermal convection) or from exsolution and precipitation of core components (i.e., by chemical convection) such as Fe at the inner core boundary, ICB [14], or SiO2 [15], or MgO [16] at the core-mantle boundary as suggested for Earth. Recent studies have both challenged [17] and supported [18] the MgO precipitation model. There is continuing debate about the relative contributions of thermal vs. chemical convection throughout the thermal and chemical evolution of terrestrial-like planetary cores [19–27]. In a purely thermally driven core, as expected in the early stages of core evolution where a solid inner core and chemical convection are absent, thermal convection is the only source of energy to power the dynamo. Thus, knowing the relative contribution of thermal conduction and thermal convection to thermal transport in the core is essential to understanding the source of energy of a core-generated magnetic field, inner core age, and thermal evolution of the core.

The contribution of conductive heat flow in the core requires the thermal conductivity of core material to be known. Thermal conductivity for metals can be approximated using the Wiedemann-Franz relation if values of electrical resistivity of Fe at high pressure *(P)* and temperature *(T)* conditions are known. This approach is often adopted [19–21] over direct measurement of thermal conductivity due to the enormous challenges in maintaining a well-controlled *T*-gradient in a small sample at very high *T* and *P* conditions [22].

Much recent attention in attempts to determine core conductive heat flow is focused on Earth. The electrical resistivity of the Earth's core was estimated to be 350–450 μΩcm from analysis of low *P* static and high *P* dynamic shock compression data [23,24], leading to calculated values of core thermal conductivity of 30–50 W(mK)−<sup>1</sup> that are generally consistent with the only experimental measurements of thermal conductivity made on Fe at core *P,T* in the diamond anvil cell [25]. However, theoretical [26,27] and experimental investigations [19,20] have suggested a much lower core resistivity (and thermal conductivity values greater than 90 W(mK)<sup>−</sup>1) for the outer core because of the effect of resistivity saturation at high *T*.

Theoretical investigation by Wagle and Steinle-Neumann [28] used a thermodynamic model and the Ziman approximation to determine the resistivity of solid and liquid Fe up to core *P* and *T* conditions. They found a decreasing resistivity change (ρ*liq* − ρ*sol*) on melting with increasing *P*. From their experimental resistivity data on *hcp* Fe at high *P* and room *T* in the diamond anvil cell (DAC) Gomi et al. [19] asserted that Fe resistivity at core conditions is close to saturation and therefore the resistivity change on melting should be negligible. From their DAC measurements, Ohta et al. [20] reported, ~20% change in Fe resistivity on melting from the *fcc* Fe phase at 51 GPa. However, lower pressure measurements in the multi-anvil press of the *T*-dependent electrical resistivity of Co up to 5 GPa [29], Ni up to 9 GPa [30], and Fe up to 12 GPa [21] demonstrated an increasing change of resistivity on melting with increasing *P*. This lower pressure regime is relevant for thermal transport at the ICB in the small planetary bodies Mercury and Ganymede.

From these multi-anvil studies, the resistivity of liquid Co and Ni along their respective *P*-dependent melting boundaries remained invariant while Fe showed a decreasing trend of resistivity below the δ-γ-liquid triple point at ~5 GPa but then remains constant above the triple point *P*. Although

experimental investigation of electrical resistivity of the α, γ, ε phases of Fe at combined static high *P* and *T* conditions have been made [20,21,31–33], its behavior through the melting transition is still contentious, hence, a detailed discussion is needed.

Generally, for the *3d* ferromagnetic metals Fe, Co, and Ni, the weak interaction of *d* electrons gives rise to an ordered magnetic state characterized by different numbers of electrons with up and down spins. Since the electronic state of a metal can be probed through the investigation of its electrical resistivity, and since electronic state and magnetism in a metal are interwoven [34], electrical resistivity can also provide information about the magnetic state of these metals. We discuss qualitatively our observation of the effect of decreasing magnon-induced electron scattering with increasing *P* on the *T*-dependent electrical resistivity of these metals at the solid–liquid transition. In addition, we discuss the possible implications of this behavior on the thermal conductivity and heat flow at the ICB of Mercury and Ganymede.

#### **2. Electronic Scattering in Ferromagnetic Metals**

For unfilled *d* band transition metals, *s–d* scattering dominates over normal *s–s* electron scattering as *T* increases due to the high density of *d*-band states. This is generally understood in Mott's *s–d* scattering model [35]. For diamagnetic metals at 1 atm, those with filled *d*-bands (e.g., Cu, Ag, and Au), the combined results from many studies show that their *T*-dependent resistivity in the solid state follows a linear dependence on *T* [36]. Similarly, for some paramagnetic metals at 1 atm (e.g., Pt, Pd, etc), their resistivity follows a near-linear dependence on *T* [37,38]. However, for the ferromagnetic metals, resistivity follows a *T2*-dependence below the Curie point and *T*-dependence above the Curie point [39–43]. With increasing *T*, the increasing phonon and spin-disorder induced scattering (magnon-induced scattering) of the highly mobile *s* conduction electrons into unfilled *d*-band states leads to decreased mobility of *s* electrons and higher resistivity. Below the Curie *T*, electron scattering is caused by a combination of phonon- and magnon-induced scattering, as well as a contribution from the asymmetry of the Fermi surface (Mott, 1964). Above the Curie *T* in the paramagnetic state, paramagnon-induced scattering tends toward a constant value while the phonon-induced scattering continues to increase with increasing *T* and therefore controls the *T*-dependent resistivity trend. Even if only qualitatively known, the relative contribution of the different scattering mechanisms is important for our study.

Probing band structure effects through resistivity investigation of the ferromagnetic metals under *P* and *T* conditions may provide an understanding of the complex electron scattering mechanisms which can occur due to topological features of the Fermi surfaces, Fermi level position, and energy gap between the spin sub-bands (δ*Eex*). Experimental studies mapping the Fermi surfaces of Fe, Co, and Ni have been accomplished primarily through the use of de Hass-van Alphen (dHvA) oscillatory effects [44] along with magnetoresistance investigations that have confirmed the existence of a complicated open orbit topology of the Fermi surfaces of these metals [45–47]. In *3d* ferromagnetic metals, magnetism is largely caused by electrons in the high density of states *3d* bands at the Fermi level. Angle-resolved photoemission studies demonstrated that the decrease in δ*Eex* above the Curie *T* for Fe, Co, and Ni is due to the energy of the spin-down sub-band shifting ~2–3 times faster than the spin-up sub-band [47–49]. Interestingly, values of δ*Eex* for Fe, Co, and Ni, and the population of the *3d*-band at ambient conditions correlate with the magnitude of the abrupt change in the electrical resistivity on melting as shown in Figure 1. Ni has the highest number of *3d* electrons (least number of unoccupied *3d* states), lowest value of δ*Eex*, and it has the greatest change in resistivity on melting as shown in Figure 1. Fe has the least number of *3d* electrons (highest number of unoccupied *3d* states), highest value of δ*Eex*, and it has the smallest change in resistivity [38] on melting. This implies that Fe, having the highest number of unoccupied *3d* states with the greatest contribution of *s–d* electron scattering induced by phonons and magnons, should show a smaller change in the resistivity on melting. The small change in the resistivity on melting can thus be explained by the extensive pre-melting scattering relative to the additional scattering arising from atomic structural change on

melting. Conversely, just prior to melting, Ni has the least contribution of scattering from phonon- and magnon-induced *s–d* electron scattering and therefore shows a larger jump in resistivity arising from the relatively larger scattering contribution on melting due to the effect of atomic structural change.

**Figure 1.** Data of 1 atm of Fe, Co, and Ni. (**a**) Resistivity discontinuity [38] on melting (note the differences in resistivity scale for Fe, Co, and Ni, whose melting T's are 1809 K, 1768 K, and 1728 K, respectively); (**b**) resistivity discontinuity on melting and number of occupied *d-*electron band, and; (**c**) sub-band energy gap and magnetic moments [34,50].

Focusing on the magnitude of resistivity of Fe, Co, and Ni in the solid state just prior to melting, the *T*-dependent resistivity of paramagnetic Fe above the Curie *T* is similar in trend to the *T*-dependent resistivity of paramagnetic Pd as phonon-induced scattering dominates in both cases, as shown in Figure 2. On the other hand, an x-ray magnetic circular dichroism study [51] showed that the net magnetic moment of Fe decreases with increasing *P* and vanishes at ~18 GPa at ambient *T* while both Ni and Co remain ferromagnetic to well over 100 GPa. The increasing population of *d-*band electrons due to *s–d* hybridization with increasing *P* [52–54] will lead to termination of magnetism. It is expected that the relative change in the positions of *s* and *d* bands in Fe, Co, and Ni with increasing *P* control the rate of *d-*band population and loss/retention of magnetism. Theoretical investigation demonstrated that the non-spin state of Fe is the most energetically favored electronic state at high *P* [51]. Through *P*-induced reduction of magnetism and tendency toward spin disorder saturation above the Curie *T*, these two effects combine to reduce or eliminate the contribution of paramagnon-induced electron scattering in the *T*-dependent resistivity region of ferromagnetic metals at high *P*.

**Figure 2.** *T*-dependent electrical resistivity of Fe at different *P* compared with Pd at 1 atm.

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

#### *3.1. Electrical Resistivity and Thermal Conductivity at the Melting Transition*

As shown in Figure 3, recent experimental investigations of the *T*-dependence of resistivity of Co [29], Ni [30], and Fe [21] at high *P* demonstrate that the effect of *P* on resistivity is greater in the high *T* region (*T*-dependent resistivity) above the Curie *T* than in the low *T* region (*T2*-dependent resistivity) below the Curie *T*. This suggests that magnon-induced scattering is less sensitive to *P* than is scattering caused by simple phonon scattering or phonon scattering that results in *s*-electrons being scattered into *d*-states. This appears intuitively expected as phonon scattering or phonon-induced *s–d* scattering arise from atomic vibration whereas magnon-induced scattering is operative at the electronic level.

**Figure 3.** Temperature dependence of resistivity of solid and liquid Fe, Co, and Ni at 1 atm and at various high pressures.

*Crystals* **2019**, *9*, 359

The *P*-dependence of liquid resistivity of Co and Ni along the melting boundary appears constant up to 5 GPa and 9 GPa, respectively, with values of resistivity on melting - ρ*liq* comparable to their corresponding values at 1 atm. The resistivity of Fe on melting decreases up to 5 GPa as it melts from the *bcc* phase but then resistivity on melting remains constant up to 12 GPa as it melts from the *fcc* phase. With a constant value of resistivity on the melting boundary, ρ*liq*, and a decreasing value of solid resistivity just before melting (ρ*sol*) with increasing *P*, ρ*liq* − ρ*sol* increases with increasing *P* up to the maximum pressures investigated in these studies as shown in Figure 4. Although, these data for Fe show an increasing ρ*liq* − ρ*sol* with increasing *P* in this low *P* range, theoretical calculation [28] up to core *P* and *T* show that the ρ*liq* − ρ*sol* for Fe melting from the *hcp* phase decreases with increasing *P*. Further experimental work is needed at higher *P* to assess the trend of ρ*liq* − ρ*sol* for Fe shown here within the context of a larger pressure range.

**Figure 4.** The difference in electrical resistivity value of solid and liquid Fe, Co, and Ni at the melting transition with increasing *P*. The least squares fits of ρ*liq* − ρ*sol* vs. *P* for Fe, Co, and Ni are, respectively, (11.33 ± 3.19) + 0.74 *P*, (12.99 ± 0.0016) + (1.35 ± 0.6) *P*, and (31.84 ± 0.61) + (1.21 ± 0.1) *P*.

Theoretical calculations demonstrate that *d*-resonance scattering dominates the electrical resistivity of unfilled *d*-band liquid metals [55–58]. Experimental study using a flash heating technique in the DAC showed the electrical resistivity of Pt along its high *P* melting boundary is constant [59]. The constancy in the liquid resistivity on the melting boundary may be understood based on the expectation that increasing *P* brings the Fermi level closer to the *d*-resonance site, hence, decreasing conduction electron mobility and increasing resistivity. However, increasing *P* also decreases phonon amplitudes and thus phonon-induced scattering which decreases resistivity. The combined antagonistic effects of *P* on these scattering mechanisms on melting could compensate each other in such a way that resistivity remains constant along the melting boundary, especially in closed packed structures [29,60]. Upon loss of, or reduction in, paramagnon-induced electron scattering at high *P* and *T* conditions, one might infer that the *T*-dependence of resistivity of ferromagnetic metals Fe, Co, and Ni in the solid state could eventually, at high enough *P*, mimic that of paramagnetic metals such as Pt and Pd, at 1 atm [36,37] or perhaps Cu, Ag, Au, and Zn [61–64], where there is a constant paramagnon-induced scattering contribution to its *T*-dependent resistivity.

#### *3.2. Heat Flow at the Inner Core Boundaries of Mercury and Ganymede*

The electronic thermal conductivity, *ke*, at planetary inner core conditions can be estimated using the Wiedemann-Franz law - *ke* = *LT* ρ where *L* is the Lorenz number. The total thermal conductivity of metals is dominated by electronic thermal conductivity [65] and one can reasonably assume they are similar in value. Mercury is thought to have a solid inner and liquid outer core with the *P* and *T* conditions at the ICB of ~5 GPa [6] and (1800–2000) K [66], respectively. Parameter values are provided in Table 1. For a pure Fe core in Mercury, using measured resistivity and melting *T* data of Fe at 5 GPa by Silber et al. [21] and the Sommerfeld value (Lo = 2.445 <sup>×</sup> 10−<sup>8</sup> WΩ/K2) of the Lorenz number [67], we compute a value of *ke* of 39 W(mK)−<sup>1</sup> for the solid just before melting and 37 W(mK)−<sup>1</sup> for the liquid side. The errors on these calculated values are mainly due to the errors on the experimentally derived values of ρ*sol* and ρ*liq* which are ~5% [21] and the *T* at the ICB; however, we used the same value of 1880 K to calculate both values. The difference in the calculated *ke* values suggest ~5% difference in thermal conductivity across the ICB of a pure Fe core in Mercury. While the choice of Lorenz number may also be debated, a single value of *L* is appropriate for calculating *ke* on both sides of the ICB which are at a single set of *P,T* conditions. An *L* value different than the one used here will not change the relative values of *ke* across the ICB. For Ganymede with *P*ICB of ~9 GPa [68] and using measured Fe resistivity and melting *T* data by Silber et al. [21], we calculate *ke* on the solid side of the ICB in a pure Fe core to be 46 W(mK)−<sup>1</sup> and on the liquid side to be 39 W(mK)<sup>−</sup>1, a difference of more than 7%. For Mercury and Ganymede, this analysis suggests that their thermal conductivity on the solid side of their ICB is likely to be higher than on the liquid side of their ICB, but only marginally when errors are considered. This difference is likely to be higher in Ganymede with *P*ICB of ~9 GPa compared with Mercury with *P*ICB of ~5 GPa.

The heat flow (*Q*) along the adiabatic *T* gradient in a liquid outer core can be expressed as:

$$Q\_{cond} = k \frac{(dT)}{dz} \Big|\_{adiabatic} = k \frac{\arg T}{C\_p} \tag{1}$$

where - *dT dz adiabatic* is the adiabatic *<sup>T</sup>* gradient and <sup>α</sup>, *<sup>g</sup>*, and *Cp* are thermal expansion, gravitational acceleration, and heat capacity at constant *P*, respectively. Heat flow transported away from the inner core that exceeds the conducted heat flow in the liquid outer core is transported by thermal convection, which in turn is available for driving a dynamo. Here, we concentrate on the heat transport across the solid and the liquid sides of the ICB of Mercury and Ganymede. At the solid side of Mercury ICB of ~5 GPa, we use a melting *T* for Fe at ~5 GPa of 1880 K [21] and we adopt an average value of 8.9 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>K</sup>−<sup>1</sup> for <sup>α</sup> from the range of values (6.4–11.4) <sup>×</sup> 10−<sup>5</sup> K−<sup>1</sup> estimated at the top of Mercury's core by Secco [67], a value of 4.0 ms−<sup>2</sup> for *g* [69], a value of 39 W(mK)−<sup>1</sup> for *ke*, and for *Cp* a value of

835 J/Kg/K which is assumed independent of *P* and *T* [70], we calculate a value of 31 mWm−<sup>2</sup> for the heat flow conducted down the adiabat on the solid side of Mercury ICB. To calculate the total adiabatic heat flow on the solid side of the ICB, we use a total core radius of 2004 km [69] along with the recently obtained estimates of inner core radius of 0.4–0.7 times the outer core radius [10]. These values yield a total adiabatic heat flow of 0.25–0.77 TW. On the liquid side of Mercury ICB, using the calculated *ke* of 37 W(mK)−<sup>1</sup> and keeping other quantities constant, we calculate a value of 30 mWm−<sup>2</sup> for the conducted heat flow and a range of total adiabatic heat flow of 0.24–0.75 TW. This analysis suggests that for a pure Fe core in Mercury, the difference in the heat conducted along the adiabat across the ICB is small and in the range of 0.01–0.02 TW and likely too small to generate significant thermal convection in the liquid outer core.


**Table 1.** Parameter values for Mercury and Ganymede used in heat flow calculations \*.

\* Values in table without references are calculated in this study.

We calculate the heat flow on the solid side of a pure Fe core in Ganymede where *P*ICB is taken as ~9 GPa and *T*ICB ~2200 K [71]. We determine *ke* on the solid and liquid side of Ganymede's ICB to be 46 W(mK)−<sup>1</sup> and 39 W(mK)<sup>−</sup>1, respectively. The size of Ganymede inner core, *r*, is not well determined, however, its core size has been estimated to lie between 650–900 km [73] and we assume an ICB radius of 650 km in our calculations. We estimate gravity *g(r)* by 4π*G*ρ*cr* where, *G* is the gravitational constant, ρ*<sup>c</sup>* is core density ~8000 kg/m2 [12] to be 4.36 m/s2. From the research of Jeanloz [72], we determine α*Fe* at 9 GPa to be 4.8 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>K</sup><sup>−</sup>1. The melting *<sup>T</sup>* of Fe at ~9 GPa is ~1990 K [21]. Using these parameters in equation 1, we estimate the heat flow on the solid and liquid sides of Ganymede's ICB to be 23 mW/m<sup>2</sup> and 19 mW/m2, respectively. For an inner core radius of ~650 km, this yields a total adiabatic heat flow of ~0.12 TW and ~0.10 TW conducted on the solid and liquid side of Ganymede ICB, respectively. This analysis shows that the larger thermal conductivity difference on the solid and liquid sides of Ganymede's ICB of ~7 W(mK)−<sup>1</sup> compared to Mercury only causes a difference of ~0.02 TW in the heat flow conducted along its adiabat, which is similar to the value for Mercury.

#### **4. Conclusions**

The *T* variation of the electrical resistivity of solid and liquid Fe, Co, and Ni through the melting transition at high *P* was discussed using experimentally measured data from previous studies. These findings were examined on the basis of reduction of magnon-induced electron scattering (quadratic dependence on *T*) at high *P* and *T*. The scattering of *s*-electrons to *d*-states in Fe, Co, and Ni above their Curie *T* can be related to the increasing phonon-induced scattering to empty *d*-states (linear dependence on *T*) and the diminishing relative effect of constant magnon-induced scattering. Relative increases of resistivity on melting in these three metals are self-consistently interpreted within this model. The *ke* of solid and liquid at the onset of melting was calculated using the Wiedemann-Franz law with the Sommerfeld value of Lorenz number. These analyses suggest that the thermal conductivity of the solid inner core of small terrestrial planetary bodies could be higher than the liquid outer core. Analysis of the thermal conductivity difference on the solid and liquid side of a pure Fe Mercury and

Ganymede inner core were performed. We found that the thermal conductivity difference on the solid and liquid sides of Mercury's ICB at ~5 GPa is ~2 W(mK)<sup>−</sup>1, which translates into a difference in total adiabatic heat flow of ~0.01–0.02 TW, depending on the size of the inner core relative to the outer core. For a pure Fe Ganymede inner core at ~9 GPa, the difference in thermal conductivity is ~7 W(mK)<sup>−</sup>1, corresponding to difference in total adiabatic heat flow of ~0.02 TW across its ICB. The cores of both planetary bodies appear to have a difference in conducted heat across their ICB that is insignificant in terms of generating thermal convection to power an internal dynamo suggesting that chemical composition is important.

**Author Contributions:** I.C.E. conceptualized the work and carried out heat flow analyses; I.C.E. and R.A.S. contributed equally to the interpretations and writing of the manuscript.

**Funding:** RAS acknowledges the Natural Sciences and Engineering Research Council of Canada (NSERC) (grant number 2018-05021) for research funding.

**Acknowledgments:** The authors thank Wenjun Yong for comments on the manuscript. ICE thanks Takashi Yoshino for the opportunity to continue research on transport properties applied to planetary bodies at the Institute for Planetary Materials, Okayama University, Misasa, Japan.

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

#### **References**


© 2019 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* **Stepwise Homogeneous Melting of Benzene Phase I at High Pressure**

**Ravi Mahesta <sup>1</sup> and Kenji Mochizuki 1,2,\***


Received: 26 April 2019; Accepted: 27 May 2019; Published: 28 May 2019

**Abstract:** We investigate, using molecular dynamics simulations, the spontaneous homogeneous melting of benzene phase I under a high pressure of 1.0 GPa. We find an apparent stepwise transition via a metastable crystal phase, unlike the direct melting observed at ambient pressure. The transition to the metastable phase is achieved by rotational motions, without the diffusion of the center of mass of benzene. The metastable crystal completely occupies the whole space and maintains its structure for at least several picoseconds, so that the phase seems to have a local free energy minimum. The unit cell is found to be unique—no such crystalline structure has been reported so far. Furthermore, we discuss the influence of pressure control on the melting behavior.

**Keywords:** benzene phase I; homogeneous melting; Ostwald's step rule; molecular dynamics simulation; high pressure; metastable phase

#### **1. Introduction**

Benzene is a renowned small and simple molecule. Regardless, its polymorphism [1–3] and local liquid structure [4–6] have been extensively explored until now. Five crystalline phases of pure benzene have been reported so far experimentally [7–14], and computational studies predicted many other potential crystalline structures [1,15–18]. Although the phase diagram shows the most stable phase at a given condition, there can be rich intermediate phases on the transition pathway from the initial to the finally prevailed phase [19–23]. Ostwald argued that a phase transition can proceed via an intermediate metastable phase due to the reduction in the surface energy of nucleation (Ostwald's step rule) [24] in contrast to the direct nucleation described in the classical nucleation theory [25]. Indeed, a metastable phase or structure has been shown to play a key role in the melting of colloids [26], copper and aluminum [27], and ice [28]. Further, under high pressures, the melting pathway can be much more complicated, because of the competition between potential energy, entropy, and packing effect (negligible at ambient pressure).

The existence of intermediate metastable phase in the melting of benzene phase I is still controversial, although the phase transition dynamics of benzene has been investigated [29–31]. Tohji et al. claimed the existence of a premelting stage at 10 K below melting point using powder X-ray diffraction [30]. Furthermore, they predicted that a plastic crystal transiently appears in the melting. This premelting stage was later disproven by the investigation over the complete temperature range of 4 to 280 K [31]. Very recently, we performed the molecular dynamics (MD) simulations for the homogeneous melting of benzene phase I crystal near the limit of superheating and statistically demonstrated that there is no intermediate transient state between the crystal and liquid phases [32]. While these studies have been conducted under ambient pressure, the melting dynamics under high pressure have not been explored.

In this study, we perform MD simulations of the homogeneous melting of benzene phase I at a high pressure of 1.0 GPa and observe the stepwise transition via a metastable crystal phase. We show that the formation of the metastable phase is achieved by rotational motions, and the unit cell structure differs from any other crystalline phases reported so far.

#### **2. Methods**

#### *2.1. Force Field*

Benzene was modelled with full atomistic detail using CHARM22 [33,34]. This force field reproduces the crystal structure of benzene phase I [29,35] and the physical properties of liquid benzene, such as density and self-diffusion coefficient [36]. The intermolecular nonbonded interactions were described by Lennard–Jones plus Coulomb potential. The intermolecular interactions were truncated at 1.20 nm. The Lennard–Jones parameters for cross-interactions were obtained using Lorentz–Berthelot mixing rules and the long-range Coulombic interactions were evaluated using particle-mesh Ewald algorithm [37].

#### *2.2. MD Simulations*

MD simulations were performed using GROMACS 2019 package [38], in which the equations of motion are integrated using leapfrog algorithm with a time step of 2 fs. The temperature *T* for equilibration was controlled using a Berendsen thermostat [39] with a damping constant of 1.0 ps, while a Nose-Hoover thermostat [40,41] was used for the production runs. The pressure *p* was isotropically controlled by a Berendsen barostat [39] for both equilibration and production runs, with a damping constant of 2.0 ps. The pressure was set to 1.0 GPa for all the simulations. The periodic boundary condition was applied in all three directions.

To evaluate the influence of pressure control on the melting behavior, we performed some additional simulations using anisotropic pressure control, where the box dimensions and the box angles can change.

#### *2.3. Crystals*

The structure of benzene phase I was obtained from Cambridge Structural Database (refcode: BENZEN) [42]. Phase I is the most stable crystalline phase up to 1.2 GPa [9,14]. The unit cell consists of four benzene molecules with space group *Pbca.* The lattice parameters are *a* = 0.7460 nm, *b* = 0.9666 nm, and *c* = 0.7034 nm at 270 K [11]. Lattice axes *a*, *b*, and *c* corresponded to *x*, *y*, and *z* axes, respectively, in this study.

For the spontaneous homogeneous melting, we replicated the unit cell to manufacture an approximately cubic system with a box size of 5.21 nm × 4.78 nm × 4.84 nm. The number of molecules, *N*, was 980. To equilibrate the crystal structure, energy minimization was first carried out using the steepest descent algorithm. Subsequently, an isochoric isothermal (*NVT*) MD simulation was performed at 270 K for 500 ps. Thereafter, an isobaric isothermal (*NpT*) MD simulation for 5 ns at 270 K and 0.1 MPa was performed, followed by a 5 ns *NpT* MD simulation at 270 K and 1.0 GPa. Furthermore, an *NpT* MD simulation at 500 K and 1.0 GPa was performed for 5 ns. The final configuration obtained was used as the initial structure for the following heating simulations.

To determine the equilibrium melting temperature we used the two-phase method [43–46]. First, we developed two identical boxes of benzene phase I crystal containing 448 molecules with a box size of 2.68 nm × 6.45 nm × 2.86 nm. We randomly erased 48 molecules from one of these crystal structures, and we melted it by performing a short *NVT* MD simulation at 800 K. Then, we joined these crystal and liquid boxes. The box dimension was 2.68 nm × 13.5 nm × 2.86 nm. To equilibrate the solid-liquid coexistence structure, energy minimization was first carried out using the steepest descent algorithm. Subsequently, an *NVT* MD simulation was performed at 200 K for 20 ps. Thereafter, the system was gradually heated up by performing each 20 ps *NpT* MD simulations at two temperatures of 200 and

400 K and at 1.0 GPa, and the equilibrated configuration was used for the production runs at 1.0 GPa. The configuration equilibrated at 200 K and 1.0 GPa was gradually decompressed from 1.0 to 0.0 GPa in increments of 0.2 GPa at 200 K by performing each 20 ps *NpT* MD simulations. The systems were further equilibrated at 400 K for 0.8 GPa, at 300 K for 0.6 and 0.4 GPa, and at 200 K for 0.2 and 0.0 GPa by performing each 20 ps *NpT* MD simulation, then the obtained configurations were used for the production runs at each pressure.

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

We initially evaluated the equilibrium melting temperature *T*<sup>m</sup> of benzene phase I at 1.0 GPa in our model using the two-phase method [43–46]. The initial configuration for production runs was the coexistence between liquid and phase I crystal, as shown in Figure 1A. We performed each 10 ns *NpT* MD simulation at 430K, 435K, and 440 K and at 1.0 GPa. Figure 1B depicts the time evolution of the volume for these runs. At *T*m, the phase I crystal and liquid exhibit the same stability in free energy, and the total volume does not change. On the other hand, the volume increases at *T* > *T*<sup>m</sup> due to melting, while it decreases at *T* < *T*m. The results show that *T*<sup>m</sup> is 435 ± 5 K at 1.0 GPa for our computational model, which is approximately 20 K lower than the experimental estimation of 457 K [9]. We also computed *T*m at 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, and 1.2 GPa in the same manner. Figure 1C shows that the melting curve for our computational model qualitatively reproduces the experimental result [9,47].

**Figure 1.** Two-phase method for *T*m: (**A**) initial configuration consisting of liquid and phase I structures, (**B**) time evolution of the total volumes in *NpT* MD simulations at 1.0 GPa and at three different temperatures (430, 435 and 440 K), (**C**) the resulting melting curve (*T*m) for the phase I (blue filled circles) with the experimental data by Bridgman (open triangles) [47] and by Akella et. al (open squares) [9]. The panel C also presents the homogeneous melting temperature (red sharp), and the stable (filled diamonds) and unstable (open diamonds) conditions for the metastable phase, estimated in this study.

Benzene phase I, equilibrated at 500 K, was heated up in the increments of 10 K at 1.0 GPa, and an *NpT* MD simulation was performed for 10 ns at each *T*. We first observed the melting of crystal at 590 K. This result indicates that the limit of superheating achieved by our heating rate (1 K/ns) is 590 K, which is considerably higher than the equilibrium melting temperature of 435 K for our computational model. The large superheating for homogeneous melting has also been reported for other materials both experimentally [26,48,49] and computationally [28,50–55], because of the absence of an apparent trigger, such as surfaces and impurities, to initiate the disordering. To observe the melting under moderate conditions, we focused on the melting at 584 K, heated up from 580 K. At 584 K, 6 of 95 independent trajectories melted within 10 ns when we gave different momenta to the same configuration.

Figure 2A,B present the time evolution of potential energy (PE) and density for a typical melting trajectory at 584 K. We find a plateau between 490 and 510 ps (shown by green shade), where both the PE and density remain constant. The preservation of these properties over a time range implies the existence of a metastable phase. Figure 3A–C depicts the molecular structures at each stage of 100, 490, and 575 ps, respectively. In contrast to the disordered liquid structure in Figure 3C, layered structures along the *y* axis are depicted in Figure 3A,B. However, their molecular orientations are clearly different. The same structure as that in Figure 3B was also observed in all the other melting trajectories at 584 K, although the structure sometimes did not fill the whole space of the simulation box and its lifetime was in the order of several 10 picoseconds.

**Figure 2.** Time dependence of (**A**) PE, (**B**) density, and (**C**) MSD for the melting of benzene phase I at 584 K. The time range when the metastable phase is formed is shaded in green.

**Figure 3.** Molecular structures of (**A**) phase I (100 ps), (**B**) metastable phase (490 ps), and (**C**) liquid phase (575 ps). The two distinct layers, M1 and M2, and the *y* axis are shown in (B).

Although the metastable structure exists for a very short time in the melting pathway at 584 K, the same structure might be (meta) stable at different *pT* conditions. To explore such a *pT* region, we performed *NpT* MD simulations for each 10 ns at pressures ranging from 0.8 to 1.2 GPa and at temperatures ranging from 450 to 550 K. At a given *pT* condition, we gave 10 different momenta to the structure obtained at 490 ps in the melting trajectory (Figure 3B). The metastable structure was found to be preserved for 10 ns in the all trajectories at several thermodynamic conditions: 450 K and 0.8, 0.9, 1.0, 1.1, and 1.2 GPa; 500 K and 0.9, 1.0, 1.1, and 1.2 GPa; 550 K and 1.1 and 1.2 GPa. At the other conditions, the metastable structure melted into liquid. These stable and unstable *T-p* conditions for the metastable structure are plotted in Figure 1C. These results imply that the metastable structure can be a "phase" rather than a transiently appearing fragile structure. Its free energy is likely to be between phase I and liquid benzene, in accordance with Ostwald's step rule [23].

In contrast to the short lifetime of the metastable phase, the waiting times for the phase transition from phase I to the metastable phase range from hundreds of picoseconds to several nanoseconds at 584 K. The significant difference in their lifetimes indicates that the free energy barrier from the metastable phase to liquid is much lower than that from phase I to the metastable phase at 584 K.

To characterize the self-diffusion over melting, we computed the mean square displacement (MSD) of center of mass (COM) of benzene from the initial configuration at 584 K (*t* = 0 ps), using the expression of |→ *<sup>r</sup>* (*t*) <sup>−</sup> <sup>→</sup> *<sup>r</sup>* (0)|<sup>2</sup> (Figure 2C). The time-dependent MSD shows a small hump at around 490 ps, arising from the re-organization of molecules to form the metastable phase. The difference in the MSD of COM is negligibly small (0.04 nm2). The MSD does not change between 490 and 510 ps, where the PE and density are also preserved. After 510 ps, the MSD gradually increases, during which a liquid nucleus forms and erodes the metastable phase. Finally, the MSD increases rapidly after 530 ps, at which point the liquid structure covers approximately half the space of the simulation box. Figure 4 shows the radial distribution functions (RDFs) between COMs of benzene molecules in the three phases. The results show that the metastable phase maintains a long-range order structure in comparison with the liquid phase. Although the metastable phase exhibits broader peaks than phase I, the peaks in the two phases are located at a similar radial distance. This indicates that there is only a minor change in the position of COM between phase I and the metastable phases. Hence, the analyses of MSD and RDF demonstrate that the transition from phase I to the metastable phase arises from the rotational motion, and COM in the metastable phase exhibits an ordered structure.

**Figure 4.** RDFs, g(r), of the COM of benzene molecule for phase I, metastable phase, and liquid phase, computed in the time range 100–150 ps, 490–510 ps, and 9000–9050 ps, respectively.

We next characterize the molecular orientation in the metastable phase. Figure 3B demonstrates that benzene molecules in every other layer of the metastable phase (we name them M1 and M2 layers) tend to have the same alignment. We project the probability density of the normal vector on a sphere using HEALPix algorithm [56–58], in which the order parameter (θ, ϕ) is used to define the vector orientation (Figure 5A). In phase I, benzene molecules are most likely to orient at (90◦, 40◦) and (90◦, 140◦), as shown in Figure 5B. It is to be noted that the molecular orientations of (90◦, 40◦) and (90◦, 140◦) are identical to (90◦, −140◦) and (90◦, −40◦), respectively. The molecular orientation in phase I is slightly different from that obtained experimentally [11] (see the green circles in Figure 5B); the difference probably arises from the limitations of the force field or the effect of high pressure. In the metastable phase, there are distinct orientations in both M1 and M2 layers (Figure 5C,D). More specifically, the molecules in the M1 layer are most likely to orient at (90◦, 0◦) and (90◦, 180◦) (or (90◦, −180◦)), while those in the M2 layer are most likely to orient at (90◦, 90◦) and (90◦, −90◦). These preferential orientations are obviously different from those for phase I. These results indicate that all molecules change the orientations in the transition from phase I to the metastable phase. In the previous study, we have shown that at ambient pressure the flipping motion played an important role in the formation of the critical nucleus [32]. Hence, the rotational motion triggers the melting transition of phase I at both the pressures. Further, it is primarily the high pressure that induces the formation of the metastable phase in the melting dynamics, because there is no intermediate metastable phase at ambient pressure.

**Figure 5.** Probability density of the normal benzene vector (**A**) whose orientation is defined by angles θ and ϕ. The normalized distributions for (**B**) phase I (100–120 ps), (**C**) M1, and (**D**) M2 layers in the metastable phase (490–510 ps) are depicted using HEALPix, such that each pixel represents an equal area on the spherical surface. The experimentally obtained orientation for phase I [11] is shown in green circles.

Furthermore, we find that the probability densities in the metastable phase are more broadly distributed than phase I. The loose orientation results in the increase of conformational entropy, which may compensate for the increase in PE (Figure 2A) and stabilize the metastable phase in term of free energy. Thus, based on the preferential orientation in each layer and the negligible diffusivity of COM, we conclude that the metastable phase is a crystal.

Using the same method, we evaluated the probability density for the two-fold axis of the benzene molecule, which is in parallel to the molecular plane and determined its most probable orientation. The average box size was calculated between 490 and 510 ps. Based on the repeating structural pattern, we determined that the unit cell is tetragonal and contains two molecules with a dimension of 0.51 nm × 0.51 nm × 0.96 nm (Figure 6). The space group is *P*42*nm*. Interestingly, as far as we know,

the obtained unit cell structure is different from any other crystalline structures of benzene reported so far [1,7–18]

**Figure 6.** Schematic representation of the metastable phase and its unit cell (dashed black lines). Every other layer along the *y* axis has the same color code. The *x*, *y*, and *z* for the simulation box and the *a*, *b*, and *c* parameters for the unit cell are shown.

It is worth investigating the influence of pressure control on the melting behavior. In the above simulations, we used the isotropic pressure control, where the dimension of box was isotropically changed. For comparison, we performed some *NpT* MD simulations using anisotropic pressure control, which allows the deformation of the simulation box. First, we performed 10 independent *NpT* MD simulations at 500 K for each 10 ns under the anisotropic pressure control, initiated from the metastable structure obtained under the isotropic pressure control. In all the trajectories, the metastable structure transformed into a different crystal form—phase III with line defect [16]. The formed phase III structure was preserved for 10 ns. Note that on the *pT* phase diagram, the phase III locates at very high pressures over 4 GPa and does not have the phase boundary with phase I [9,14]. Our results imply that the phase III can be a metastable phase between phase I and liquid. Second, we heated up the phase I crystal under the anisotropic pressure control until it melted at 580K. However, we did not observe a stepwise homogeneous melting via phase III at the limit of superheating. The phase III may appear on the melting pathway at lower temperatures. These comparisons indicate that the metastable structure that we found in this study is likely to appear under the isotropic pressure control, which can be realized by the recently developed experimental techniques [59–61].

#### **4. Summary**

We carried out the MD simulations to investigate the homogeneous melting of benzene phase I under a high pressure of 1.0 GPa. We found a stepwise change in the PE, density, and MSD in the melting trajectories. The metastable phase preserves the layered structure as phase I but the molecular orientations are obviously different from those of phase I. The formation of the metastable phase is attributable to the flipping motion of benzene rather than the diffusion of COM. Although the flipping motion has been previously shown to lead the formation of the critical nucleus, no intermediate metastable phase was observed at ambient pressure. Furthermore, the probability density of the normal vector demonstrates the preferential orientation in the metastable phase. Their orientations are more relaxed than those in phase I, indicating the gain of conformational entropy compensates the increase of potential energy. We further determined the unit cell and found that such a crystalline structure has not been reported so far.

**Author Contributions:** Conceptualization, K.M.; molecular simulations, R.M.; writing, K.M. and R.M.; funding acquisition, K.M.

**Funding:** We acknowledge the support from the Japan Society for the Promotion of Science (KAKENHI 18K19060) and Inoue Foundation for Science (Inoue Science Research Award).

**Acknowledgments:** We thank M. Matsumoto for determining the unit cell. A part of calculations was performed at the Research Center for Computational Science in Okazaki, Japan.

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

#### **References**


© 2019 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* **Pressure E**ff**ects on the Optical Properties of NdVO4**

#### **Enrico Bandiello 1,\* , Josu Sánchez-Martín 1, Daniel Errandonea <sup>1</sup> and Marco Bettinelli <sup>2</sup>**


**\*** Correspondence: enrico.bandiello@uv.es

Received: 25 April 2019; Accepted: 4 May 2019; Published: 6 May 2019

**Abstract:** We report on optical spectroscopic measurements in pure NdVO4 crystals at pressures up to 12 GPa. The influence of pressure on the fundamental absorption band gap and Nd3<sup>+</sup> absorption bands has been correlated with structural changes in the crystal. The experiments indicate that a phase transition takes place between 4.7 and 5.4 GPa. We have also determined the pressure dependence of the band-gap and discussed the behavior of the Nd3<sup>+</sup> absorption lines under compression. Important changes in the optical properties of NdVO4 occur at the phase transition, which, according to Raman measurements, corresponds to a zircon to monazite phase change. In particular, in these conditions a collapse of the band gap occurs, changing the color of the crystal. The changes are not reversible. The results are analyzed in comparison with those deriving from previous studies on NdVO4 and related vanadates.

**Keywords:** vanadate; zircon; high pressure; band gap; phase transition; optical absorption

#### **1. Introduction**

Trivalent metal orthovanadates have been constantly studied during the last decade because of their optical properties, which makes them appropriate for many different technological applications. Such applications go from photocatalytic water purification and hydrogen production to scintillators, phosphors, solid-state lasers, magnetoelectric, and medical applications [1–6]. Neodymium vanadate (NdVO4) is one of the most promising orthovanadates for the development of technological applications [1–6], for which an accurate knowledge of NdVO4 optical properties is a fundamental requisite. Several studies have been carried out on the optical properties of NdVO4 [6–11], mainly focused on the study of the band-gap energy and the optical absorptions associated to internal transitions of Nd3<sup>+</sup> ions from 4I9/<sup>2</sup> to 4G5/2, 4F7/2, 4F5/2, and 4F3/<sup>2</sup> [12]. Discrepancies on the value of the fundamental band gap (*E*gap) can be found in the literature [6–11], with values of the *E*gap ranging from 2.95 to 3.72 eV. Consequently, new studies are needed to accurately determine this parameter, which is crucial for the understanding of the properties of this wide gap semiconductor.

A second focus of research on NdVO4 has been the influence of pressure on its properties and crystal structure. Numerous studies have been carried out, reporting the existence of a phase transition around 5–6 GPa [13–16]. Interestingly, it has been found that the crystal structure of the high-pressure (HP) phase strongly depends on the conditions of the experiment (namely, hydrostatic or non-hydrostatic conditions). On the other hand, a second phase transition has been found to occur at 18–20 GPa, according to Raman and X-ray diffraction (XRD) experiments [13,14]. In contrast, optical-absorption experiments located this transition at 12 GPa. These results also suggest that additional studies of the pressure effects on the optical properties of NdVO4 are required. This could give not only information about the pressure dependence of the band gap but also on the structural stability. In particular, the study of the behavior under compression of the internal absorption transitions

of Nd3<sup>+</sup> ions, which strongly depend on their local environment, can be useful to characterize the HP phase transitions [17,18].

Here we report a systematic study of the optical properties on NdVO4 in the ultraviolet (UV)–visible (VIS) range. The studies have been carried on single crystals growth by the flux growth method. The maximum pressure achieved is 12 GPa. Our results confirm that NdVO4 is a wide-gap semiconductor with *E*gap = 3.72(2) eV. The pressure dependence of *E*gap is reported and discussed, as well as the pressure dependence of the absorption bands associated to internal transition of Nd3<sup>+</sup>. We demonstrate the existence of only one phase transition in the pressure range covered by the experiment. This transition is non-reversible and causes a band-gap collapse, along with important changes in the absorption bands associated to Nd3<sup>+</sup>, as shown and discussed in the following.

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

NdVO4 single crystals, up to 10 mm long, were obtained by the flux growth method [19–21]. A 1:1.923 mol. mixture of V2O5 and PbO (99%, Carlo Erba, Val-De-Reuil, France, and 99.9%, Sigma Aldrich, St. Louis, MO, USA, respectively) plus borax (Na2B4O7·10H2O, > 99%, JT Baker, Waltham, MA, USA) was used as the flux. Borax was added as a flux modifier, in order to increase the size of the crystals [22]. Nd2O3 (0.7394 g, purity 99.99%, Sigma Aldrich, St. Louis, MO, USA) was added as the crystal precursor. All the reagents were in form of fine powders. The flux composition was 21.3082 g of PbO, 9.1235 g of V2O5, and 1.25 g of Na2B4O7·10H2O. Platinum crucibles were filled with the flux and the precursor, sealed with a platinum lid, and put in a programmable oven. To dehydrate the mixture, the temperature was ramped up at 0.6 ◦C/hour and kept at 250 ◦C for 2 hours. Afterwards, the temperature was kept fixed for 15 hours at 1300 ◦C, after a gradual increment at a rate of 105 ◦C/hour. During this step, the flux melts and acts as a solvent for Nd2O3. The temperature was then slowly decreased until 950 ◦C (−1.8 ◦C/hour), in order to promote the formation of the crystals by precipitation and spontaneous nucleation. The crucible was then removed from the oven, rapidly reversed (to facilitate the recovery of the crystals) and allowed to cool down to room temperature. Subsequently, the platinum lid was removed and the crucible was immersed in hot nitric acid (1.5 M), which was continuously stirred and renewed multiple times until the complete dissolution of the flux. Finally, the clean crystals were washed with deionized water and recovered using a paper filter.

In order to characterize the growth of NdVO4 crystals we performed energy-dispersive X-ray spectroscopy (EDXS) in a transmission electron microscopy (TEM, FEI Company, Hillsboro, OR, USA) operated at 200 kV (TECNAI G2 F20 S-TWIN, Waltham, MA, USA), powder XRD measurements using a Rigaku D/Max diffractometer (Tokyo, Japan) with Cu-K<sup>α</sup> radiation (λ = 1.5406 Å), and Raman measurements in backscattering geometry with a Jobin Yvon THR1000 single spectrometer (Kyoto, Japan) equipped with an edge filter, a 10X objective (Mitutoyo, Kawasaky, Japan) with a numerical aperture of 0.28, and a thermoelectric-cooled multichannel charge-coupled device (CCD) detector, using a 632.8 nm HeNe laser (10 mW power).

Optical-absorption experiments were accomplished in the UV–VIS range using the optical set-up previously described by Segura et al. [23], which is equipped with Cassegrain objectives (6X magnification and 0.60 numerical aperture, Edmund Optics, Mainz Germany) and an UV–VIS spectrometer (USB4000-UV–VIS - Ocean Optics, Duiven, Netherlands). Optical absorption spectra were obtained from the transmittance of the sample, measured using the sample-in sample-out method [24]. For these measurements, platelets of size 60 × 60 μm and approximately 5 μm thick were cleaved along the {110} plane. The experiments were carried out both at ambient and at high pressure. The HP experiments were performed using a diamond anvil cell (DAC, Chervin-type, Sorbonne University, Paris, France) equipped with IIA-type diamonds with a culet size of 480 μm. The crystals were loaded in a 200 μm hole of a steel gasket, pre-indented to a thickness of 50 μm. As the pressure medium, we employed a 16:3:1 methanol–ethanol–water mixture, which remains quasi-hydrostatic up to 10.5 GPa [25]. The NdVO4 crystals were carefully located at the center of the gasket hole to avoid bridging [26]. Small ruby chips were loaded next to the sample for pressure determination [27] with

an accuracy of 0.05 GPa. The full width half maximum of R1 and R2 lines of ruby fluorescence and the splitting between them confirm the quasi-hydrostatic conditions of the experiments [28].

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

#### *3.1. Ambient-Pressure Characterization*

The obtained single crystals are transparent with a strong violet hue. The growth is preferentially along the *c*-axis, with size of the order of 5 × 1 × 0.2 mm (see inset in Figure 1). EDXS shows the presence of only Nd, O, and V, indicating that the grown crystals are free from impurities. According to EDXS analysis, the ratio of Nd to V is found to be ~1, confirming the formation of stoichiometric NdVO4. Results from powder XRD measurements are shown in Figure 1 together with the results of the structural refinement and residuals. The phase purity of NdVO4 has been confirmed by the Rietveld refinement, along with the expected zircon-type crystal structure (space group I41/amd) [29]. A total number of 31 reflections have been observed and all of them have been used in the refinement. Nd and V are located at high-symmetry atomic positions 4a (0, 3/4, 1/8) and 4b (0, 1/4, 3/8), respectively. This facilitates the determination of the positions of oxygen atoms, which were established to be at 16 h (0, 0.4292(5), 0.2055(5)). The unit-cell parameters are *a* = 7.3311(8) Å and c = 6.4359(7) Å. According to these results, the zircon structure is composed of regular VO4 tetrahedral units with a V-O distance of 1.7076(6) Å and NdO8 dodecahedra with two different Nd-O distances, 2.4082(9) and 2.5001(9) Å. These results are in full agreement with the literature [4,29]. The small residuals in Figure 1 and the R-factors of the refinement (RP = 5.8% and RWP = 7.11%) fully support the assignment of the zircon-type structure.

**Figure 1.** Powder XRD pattern measured in NdVO4. The dots are the experiment. The solid lines the refinement and residual. Ticks indicate the positions of Bragg peaks. The most intense peaks are labeled with their indexes. A picture of the grown crystals is show in the inset, showing their size and transparency. For reference, the squares are 1 × 1 mm.

Results from the Raman experiment at ambient pressure are shown in Figure 2. Group theory predicts 12 Raman active modes for the zircon-type structure, whose irreducible representation is 2A1g + 4B1g + B2g + 5Eg [30]. All the expected Raman active modes could be identified in our experiment. Notice that the two broadest peaks (those at ~240 and 378 cm<sup>−</sup>1) are doublets consistent of modes very close in frequency (see Table 1). No other modes than those representative of zircon-type NdVO4 have been observed. This supports that only the stable polymorph of NdVO4 is present in our crystals, a conclusion that agrees with the results of XRD experiments. The wavenumbers of the different modes are summarized in Table 1. They agree with those reported by Santos et al. [30], Nguyen et al. [31], and Panchal et al. [13]. The mode assignment was done according to these previous works. ω ω ω

**Figure 2.** Raman spectrum of low-pressure zircon NdVO4 at ambient pressure (bottom) and high-pressure NdVO4 at 6.9 GPa (top). The modes of the zircon phase are identified.


**Table 1.** Frequencies (ω) in cm−<sup>1</sup> of the Raman-active phonon modes at room conditions in NdVO4.

Results from optical-absorption experiments are shown in Figure 3. The absorption spectrum has an abrupt absorption edge at low wavelengths plus eight absorption bands, named as A1 to A8. These bands resemble those of the transitions between the levels of Nd3<sup>+</sup>, while the strong low-wavelength absorption corresponds to the fundamental band gap of NdVO4, with the typical Urbach-tail shape of this family of compounds [32,33]. This can be clearly seen in the inset of the figure, where the absorption is plotted as a function of the photon energy (*E*). In particular, it closely resembles the absorption edges of the direct band gap of zircon-type YVO4, LuVO4, and YbVO4 [11]. By assuming that the absorption coefficient (α) of NdVO4 obeys an Urbach-tail energy dependence, α(*E*) = α0*e*−(*Egap*−*E*)/*Eu*, we have been able to explain the observed strong low-wavelength (high-energy) absorption edge. In the equation, *E*<sup>u</sup> is the Urbach energy, which is related to the steepness of the absorption edge. In the inset of Figure 3 we show the goodness of this dependence fit for the high-energy part of the absorption spectrum (red line). The best fit has been obtained with *E*gap = 3.72(2) eV and *E*u = 0.081(5) eV. The value of the band-gap agrees with the results reported by Panchal et al. [11] and the value of *E*u is comparable with the values reported for related oxides [34].

**Figure 3.** Optical-absorption of zircon-type NdVO4. The inset shows the high-energy region as a function of energy including the fir to the step fundamental absorption described in the text (red line).

The rest of the absorption features observed in our experiment (A1 to A8 in Figure 3) can be correlated to transitions between 4f3 levels of Nd3<sup>+</sup>. All the absorption bands are composed by several overlapping peaks, due to the expected broadening at room temperature. They are centered at 359, 439, 474, 528, 593, 688, 753, and 817 nm. They likely correspond to transitions from 4I9/<sup>2</sup> to 4G11/2, 4G9/2, 4G7/2, 4F9/2, 4G5/2, 4F7/2, 4F5/2, and 4F3/2, respectively [35]. The features at long wavelengths are the strongest and have been previously reported [6–10]. Among the others, only the absorption band centered at 359 nm has been previously reported, even if mistakenly assigned to the fundamental band-gap of NdVO4 [8]. This leads to the underestimation of the band gap, with values from ranging 2.95 to 3.5 eV [6–10]. The band gap of NdVO4 measured in our experiment is larger (3.72 eV), compatible with the band gap of rare-earth vanadates, and its value is determined by the classical field splitting of the VO4 <sup>3</sup><sup>−</sup> ion [26], the top of the valence band and bottom of the conduction band being dominated by V 3d and O 2p states. A possible cause for the aforementioned underestimation of the band-gap can be the fact that the previous studies were carried out on nanoparticles and nanowires and not on single

crystals. The diffusion of light through nanomaterials used in previous experiments, along with the poorer crystallinity, would lead to a larger Urbach energy, making the fundamental absorption edge less steep and causing it to merge with the highest energy absorption band of Nd3+, therefore leading to an underestimation of *E*gap.

#### *3.2. High-Pressure Studies*

We will now present the results from HP experiments. In Figure 4, we show the high-energy portion of the absorption spectrum, which we will use first to discuss the effect of pressure on the band gap. Afterwards, we will focus on the influence of pressure on the absorptions related to Nd3<sup>+</sup>. These results will be correlated with the information known on the HP behavior of the crystal structure and Raman modes of NdVO4 [13–16]. The first phenomenon we observe is that the absorption edge gradually blue-shifts from ambient pressure to 4.7 GPa. In particular, we found that in this range of pressure *E*gap increases from 3.72(2) eV to 3.79(2) eV. The shape of the absorption edge does not change with pressure, which indicates that the band gap remains a direct one. This result is in good agreement with the report by Panchal et al. [11]. In Figure 5, we present the pressure dependence of *E*gap and compare it with the literature [11]. Both experiments show a linear increase of *E*gap, which is a consequence of the increase of the repulsion of bonding and anti-bonding states as the V-O distances of the VO4 tetrahedron are reduced upon compression [13].

**Figure 4.** Optical-absorption NdVO4 at different pressures (indicated in the figure). The spectra assigned to the low(high)-pressure phase are represented with solid (dashed) lines.

**Figure 5.** Pressure dependence of the band-gap energy of NdVO4. Circles (diamonds) correspond to the zircon (monazite) phase. Squares correspond to the post-monazite phase. Results from the present study are shown with solid symbols and results from the literature [11] are shown with empty symbols.

At 5.4 GPa we found an abrupt decrease of *E*gap. This can be seen in Figure 4 by the sudden red-shift of the absorption edge when moving from 4.7 to 5.4 GPa. This corresponds to a decrease of the band gap from 3.79(2) eV to 3.17(2) eV. This fact is indicative of the occurrence of a phase transition. From 5.4 up to 12 GPa, the absorption edge gradually red-shifts, with no indication of the existence of a second phase transition. Upon decompression, the transition is found to be irreversible in our experiment, being *E*gap = 3.26(2) eV at ambient pressure in the recovered sample. According to the literature, the observed transition can be either from zircon to scheelite or monazite [13–16]. In order to clarify this issue, we have carried out a Raman experiment at 6.9 GPa. The results are shown in Figure 2. For the scheelite phase, according to group theory, 13 Raman-active modes are expected (3Ag + 5Bg + 5Eg) [36]. For the monazite phase, 36 Raman-active modes are expected (18Ag + 18Bg) [14]. The difference in number of modes is very noticeable in the high-frequency region, where scheelite has a very characteristic feature with only three modes. In our experiment, the Raman spectrum has manifestly more than thirteen peaks. In addition, it very much resembles the Raman spectrum associated to monazite-type NdVO4 by Panchal et al. [14]. Therefore, the changes observed in the absorption spectrum are likely to be triggered by the occurrence of the zircon–monazite transition. This is a first-order reconstructive transition that involves a volume collapse and an increase in the coordination number of Nd from 8 to 9. In the HP monazite phase, Nd is surrounded by nine oxygens forming NdO9 pentagonal interpenetrating tetrahedral units, while in zircon Nd has eight neighboring oxygens forming NdO8 dodecahedral units.

From the absorption spectra shown in Figure 4 we obtained the pressure dependence of *E*gap, using the same procedure to determine it at ambient pressure. The results are shown in Figure 5 and compared with those reported by Panchal et al. [11]. For the zircon-type phase we obtained a very similar linear pressure dependence. In our case d*E*gap/d*P* = 11.2(9) meV/GPa. In the figure, the abrupt collapse of the band gap from 4.7 to 5.4 GPa can be seen, which is associated to the zircon–monazite

transition. This transition involves a collapse of the volume, a lowering of the crystal symmetry, and a distortion of the VO4 tetrahedron [16]. These changes of the crystalline structure are reflected in the electronic structure of NdVO4, producing the above described collapse of *E*gap. Notice that the variation of the electronic band gap is much more rapid and sharp than the changes observed in the Raman and XRD experiments, suggesting that optical experiments at the intrinsic absorption edge are more accurate for the determination of the transition pressures. This is because the electronic properties are much more sensitive than the structural and vibrational properties to the influence of pressure. A small change in the crystal structure can induce large changes in the band structure and consequently in the electronic properties. In our case, the band-gap collapse at the transition is 0.1 eV larger than in previous experiments [11].

After the phase transition, the band gap follows a linear behavior but red-shifts under compression with a rate d*E*gap/d*P* = 19.5(9) meV/GPa. The HP monazite phase remains stable up to 12 GPa, the highest pressure covered by our experiments. Notice that we did not find evidence of the second transition reported by Panchal et al. [11]. However, we observed the formation of cracks in the crystal at 12 GPa, which allowed the diffusion of light and precluded the performance of accurate experiments at higher pressure. We consider the formation of cracks a possible precursor effect of the second transition reported by Panchal et al. [11] at 11. 4 GPa [37], which could be triggered at a lower pressure than in our experiment because of the use of a thicker crystal, causing the bridging of the crystal between diamonds [23]. Interestingly, we found that after releasing pressure, our samples remain in the HP phase. This allows us to determine the band gap of the metastable monazite-type polymorph at ambient pressure, with *E*gap = 3.26(2) eV. This value is probably more suitable for many photocatalytic applications than that of zircon-type NdVO4. The red-shift of the band gap in the HP phase is consistent with the results obtained for monazite-type LaVO4 [18]. The explanation for the closing of the gap with pressure in monazite (in contrast to the opening of the gap in zircon) comes from the contribution of Nd 4f and 5d orbitals to the bottom of the conduction band in monazite NdVO4. These states are much more localized than O 2p and V 3d states. As a consequence, the bottom of the conduction band does not shift with pressure, but the top of the valence band moves toward higher energies, causing the observed band gap decrease.

We will comment now on the absorption bands related to Nd3+. In Figure 4, it can be seen that at the phase transition the absorption band around 359 nm (A1 in Figure 3) assigned to 4I9/<sup>2</sup> to 4G11/<sup>2</sup> disappears because the band gap of the HP phase is at a lower energy (longer wavelength). In the low-pressure zircon phase, the effect of pressure on the A1–A8 absorption bands is qualitatively the same for all of them. We will therefore focus on bands A4 and A5 ( 4I9/<sup>2</sup> to 4F9/<sup>2</sup> and 4I9/<sup>2</sup> to 4F5/<sup>2</sup> transitions), two of the better resolved bands, to illustrate the pressure behavior of internal Nd3<sup>+</sup> transitions. Their pressure dependence is shown in Figure 6. The absorption bands consist of several multiplets caused by the Stark splitting. In the low-pressure phase (pressure ≤ 4.7 GPa), there are no qualitative changes as pressure increases. Basically, the absorption bands show the typical red shift which can be attributed to an increase of the covalence in NdVO4 [37]. The position of each of the identified peaks as function of pressure is presented in Figure 7. The behavior is qualitatively similar to related oxides doped with Nd3<sup>+</sup> or other lanthanides [17,18,38,39]. In particular, our results suggest that the increase of the crystal field interaction with pressure in NdVO4 is qualitatively similar to that of isomorphic Nd-doped YVO4 [38].

**Figure 6.** Absorbance relative to the 4I9/<sup>2</sup> to 4G11/<sup>2</sup> and 4I9/<sup>2</sup> to 4F9/<sup>2</sup> transitions measured at room temperature at different pressures. The spectrum collected after pressure release is marked with (r). The spectra are plotted as function of wavenumber (top) and wavelength (bottom). They have been vertically shifted to facilitate comparison.

It is noteworthy that, for pressures below the phase transition, some band splitting and merging is present in the 1.5–3.8 GPa range. This may be caused by the well-established distortion of the NdO8 tetrahedron under high-pressure conditions [16]. For pressures higher than 4.7 GPa, the absorption bands are totally different, confirming the phase transition. The optical spectroscopy associated to Nd3<sup>+</sup> is thus quite sensitive to pressure, proving also to be a very efficient tool for the detection of phase transitions. Again, the changes induced by pressure in the absorption spectra are not reversible upon decompression (see Figure 4). The changes of the absorption bands are consistent with the coordination change of Nd3<sup>+</sup> at the zircon-monazite transition. In the HP monazite phase Nd3<sup>+</sup> is coordinated by nine oxygen atoms forming a distorted polyhedron. As in the low-pressure phase, in the HP phase the absorption bands also mainly red shift under compression. However, the pressure shift is slightly reduced, which is caused by the decrease of the compressibility due to the density increase associated to the transition [16]. To conclude, we would like to add that in both phases we observed the splitting of some of the multiplets. We believe this is the results of enhancement of exchange interaction in Nd3<sup>+</sup>-Nd3<sup>+</sup> pairs as pressure increases [18,40].

**Figure 7.** Peak positions observed under pressure increase. Lines are a guide for the eye.

#### **4. Conclusions**

In this study, single crystals of NdVO4 where synthetized and characterized by XRD, UV–VIS absorption, and Raman spectroscopy at ambient conditions and at high pressure (up to 12 GPa). In particular, optical spectroscopy allowed a precise determination of the band-gap energy of NdVO4 (3.72(2) eV), which is found to be larger than that reported in previous studies. We compared the pressure evolution of the fundamental band-gap and of the absorption bands associated with transitions of Nd3<sup>+</sup> energy levels. The shape of the fundamental absorption edge does not change with pressure and a phase transition around 5 GPa is evidenced by a sudden collapse of the band gap (Δ*E*<sup>g</sup> = −0.69 eV). These HP experiments are thus an accurate tool for establishing the band gap. On the other hand, the absorption bands of Nd3<sup>+</sup> show a more gradual evolution but the phase transition is equally put in evidence by their reduced number in the high-pressure phase. The band-gap collapse

near 5 GPa can be correlated to the well-known zircon–monazite transition of NdVO4. The transition is irreversible, as the HP phase could be recovered at ambient pressure, with a band gap of 3.26(2) eV. The obtained results have been compared with previous studies. We hope that these conclusions will contribute to deepening the understanding of the optical properties of NdVO4 and its behavior in high-pressure conditions.

**Author Contributions:** All authors contributed equally to this work; being all involved in experiments, analysis, interpretation of results, and writing of the manuscript.

**Funding:** This work was supported by the Spanish Ministry of Science, Innovation and Universities, the Spanish Research Agency, and the European Fund for Regional Development under grant MAT2016-75586-C4-1-P and by Generalitat Valenciana under grant Prometeo/2018/123 (EFIMAT). E. Bandiello thanks the Generalitat Valenciana for his postdoctoral contract (ValI+D, APOSTD2017).

**Acknowledgments:** The authors thank Erica Viviani (Department of Biotechnology, University of Verona) for the technical assistance in the synthesis of NdVO4 crystals.

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

#### **References**


© 2019 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/).
