Phase stability and thermodynamic properties of FeS polymorphs Journal of Physics and Chemistry of Solids

We present a Hubbard-corrected density functional theory investigation of FeS polymorphs based on the quasi-harmonic theory of lattice vibrations. We show that the ﬁ rst temperature transition of troilite FeS cannot involve the MnP-type FeS phase, as sometimes reported in the literature, while the sequence of polymorphs in the pressure range 0 – 100 GPa at room temperature supports the experimental observations. Although with some differences in the critical pressures, our ab initio phase diagram is in line with those derived from X-ray diffraction studies. The thermodynamic properties of troilite FeS are in good agreement with those measured, which lends support to the accuracy of our predictions for the other FeS phases that are less accessible experimentally.


Introduction
FeS polymorphs are of significant relevance to condensed matter physics and planetary science. Troilite (FeS I) is an antiferromagnetic insulator with a peculiar semiconductor-metal-semiconductor transition induced by pressure [1][2][3] and a remarkably high magnetoelectric response suitable for electromagnetic devices [4]. FeS I is also the stoichiometric end-member of a pyrrhotite group of iron sulfides with promising applications in catalysis [5].
FeS phases are thought to form the cores of Earth and Mars, which is suggested by the presence of FeS I in many meteorites [6]. This possibility has sparked interest in investigating the numerous polymorphs forming at different thermodynamic conditions, especially with the goal of modelling the cores of the terrestrial planets [7][8][9][10][11][12][13][14][15].
At ambient conditions, FeS I has a hexagonal NiAs-type superstructure (P62c; a ¼ ffiffiffi 3 p A, c ¼ 2C), where A and C are the unit cell parameters of the NiAs-type, and it transforms to an MnP-type structure (FeS II) with an orthorhombic cell (Pnma; a ¼ C, b ¼ A and c ¼ ffiffiffi 3 p A) by room temperature compression at 3.4 GPa [16]. A further increase in pressure results in a second transition at 6.4 GPa to FeS III, which has a non-magnetic monoclinic structure (P2 1 =a) [17,18]. Upon heating at low pressures, FeS I transforms to FeS IV [19,20] and FeS V [6,17] in this order, having respectively hexagonal (P6 3 mc; a ¼ 2A, c ¼ C) and simple NiAs-type superstructures (P6 3 =mmc; a ¼ A, c ¼ C). At pressures above 36 GPa and room temperature, a non-magnetic MnP-type phase (FeS VI) becomes stable over FeS III [14,21]. First principles results have predicted a phase transition from FeS VI to a phase FeS VII with an orthorhombic cell (Pmmn) [22], but this polymorph has never been observed experimentally. Sata et al. have instead reported the existence at 180 GPa of a CsCl-type phase (FeS VIII), which is not related to the NiAs-type structure and has a cubic unit cell (Pm3m) [23]. In all the FeS I-VIII polymorphs, Fe atoms present octahedral coordination, as illustrated in Fig. 1.
While previous first-principles studies by Martin et al. [24] and Ono et al. [22] have provided fundamental insights into the stability of the FeS polymorphs at 0 K, in this work, we go beyond these earlier studies to introduce the effect of temperature. The calculation of the Gibbs free energies by means of the quasi-harmonic approximation allows us to derive the ab initio phase diagram as well as the thermodynamic properties of the FeS polymorphs, knowledge of which is key to advance our understanding of the cores of the terrestrial planets.

Thermodynamics within the quasi-harmonic approximation
We have used the quasi-harmonic approximation to derive the thermodynamic properties of the FeS phases at finite temperatures and pressures, as detailed in a comprehensive recent review [25]. In this approach, the Helmholtz free energy F of a vibrating crystal is separated as the electronic energy E 0 , obtained by DFT, and the vibrational free energy F vib of the ions: If we denote with ω ks the sth oscillation mode of the wavevector k in the harmonic approximation, F vib can be expressed as a sum of all the vibrational modes in the crystal: Given FðV; TÞ, we can derive the pressure p of the system through the equation of state: Equations (1)-(3) imply a dependence of the pressure on the vibrational modes. Since these decrease when the volume of a crystal is expanded, it is possible to take into partial account some anharmonic effects by calculating the phonon spectrum at different volumes (in a fully harmonic crystal the frequencies would not change with the volume). The pressure can then be used in conjunction with the Helmholtz free energy to derive the Gibbs free energy G: Finally, thermodynamic properties like the heat capacity at constant volume C v or the entropy S can be derived from the phonon spectra as a function of temperature:

Computational details
All geometry optimisations were performed with VASP 5.3 [26,27] using the PBE functional [28]. It has been reported that a DFT þ U approach with a moderate on-site Coulomb interaction reproduces the experimental equilibrium volume of FeS I, while values of U larger than 1 eV lead to structural distortions [29]. In addition, Ricci et al. [4] have shown that the PBE functional with a Hubbard correction of 1 eV on the Fe-d orbitals eliminates the imaginary frequencies from the phonon spectrum predicted by functionals based on the density gradient [15]. In view of these findings, we have applied the same U parameter of 1 eV in all the calculations. We mention that this value is relatively close to that of 1.6 eV resulting from minimising the differences between the calculated and the experimental band gaps, unit cell volumes and bulk moduli of three different FeS x compounds (including FeS I and FeS II) [30].
We performed a first volume scan including eight data points spanning a typical range of ±10% from the equilibrium volume, during which we relaxed the shape and the atoms of the unit cells while keeping the volumes fixed. The initial volume range was expanded with additional points to ensure the inclusion of all the thermodynamic conditions considered in the study [31]. With the optimised lattice vectors, we performed a second scan to obtain the ground state energies by relaxing the internal atomic positions within each space group symmetry. At this point, we obtained the unit cell volumes per formula unit and the bulk moduli at zero pressure, V 0 and B 0 , respectively, by fitting the energy-volume data of all phases to a third order Birch-Murnaghan equation of state [32]. Next, in order to derive the vibrational free energies (equation (2)), we calculated the phonon spectra at each volume. Finally, we also fitted the Helmholtz free energy-volume data (equation (1)) of the phases that were identified as dynamically stable to a third order Birch-Murnaghan equation of state.
We have employed the projector augmented wave method to model the core-electron interaction [33], treating explicitly the 4s, 3d and 3p electrons of Fe, and the 3s and 3p of S. The spin configurations of FeS I and FeS II were taken from a previous work [24]. Non-polarised spin calculations were performed for FeS III, FeS VI, FeS VII and FeS VIII, in agreement with the loss of magnetism reported in the literature [22,34]. The spin configuration of FeS I, antiferromagnetic along the c-axis, was employed for FeS IV and FeS V [8]. All the optimisations were performed with a plane wave cutoff of 600 eV and stopped when the forces acting on the ions were less than 10 À3 eV/Å. After finding that a 4 Â 4 Â 2 Monkhorst-Pack grid [35] ensures that absolute energies of FeS I are converged to better than 1 meV per formula unit, we scaled the grids inversely with their unit cells dimensions as follows: Following the finite-displacement method implemented in the PHO-NOPY package [36,37], we have used its VASP interface to calculate the phonon frequencies. Test calculations with a 1 Â 1 Â 1 and a larger 2 Â 2 Â 1 supercell of FeS I gave absolute changes in the vibrational free energy of 6 meV per formula unit at the highest temperature of 600 K considered in this study. Based on this, force constants were calculated with the following supercells:

Structural properties
We report in Table 1 the volumes per formula unit and the bulk moduli of the FeS polymorphs at p ¼ 0 GPa and T ¼ 0 K. Our calculated V 0 are within 6% of the experimental observations, which supports the choice of the U parameter to model all FeS polymorphs considered (see subsection 2.2). Consistent with all the previous literature [8,22,24,38], the agreement is less satisfactory for the experimental bulk moduli at zero pressure. Estimating B 0 by fitting the Helmholtz free energy, rather than the static DFT energy, could not explain the differences reported. For example, at the experimental room temperature of references [16] and [39], the calculated bulk moduli of FeS II and FeS III reduce to 65.8 and 144.1 GPa, from the 0 K values of 70.8 and 149.3 GPa, respectively. However, the accuracy of the experimental data is very dependent on the stability field of the phase investigated, and an interval of only a few GPa can result in large uncertainties in the experimental functions used to fit the compression behaviour [22]. Accordingly, Nelmes et al. did not consider that the pressure range of their data was large enough to provide any meaningful estimates of B 0 in FeS I-III [18]. Finally, we note that the calculated bulk modulus of FeS VI (158.4 GPa) is remarkably close to the experimental value at room temperature (156 GPa), reflecting the fact that volume-pressure data could be collected for a window of several tens of GPa in the work of reference [21].

Stability at finite temperatures
The transition sequence of the FeS polymorphs at increasing temperature under ambient pressure has long been debated. According to Kings et al., FeS I transforms at 420 K to orthorhombic FeS II [16]. Similar Table 1 Volumes per formula unit and bulk moduli of the FeS polymorphs at zero pressure. All data from this and other DFT works correspond to T ¼ 0 K.  results have been obtained by Kruse et al., who concluded that the so-called α transition occurs in an interval of at least 40 K and is complete at 413 K [40]. Other works have instead suggested that hexagonal FeS IV is the phase forming at the α transition [6,19]. There is a general agreement, however, that further increases in temperature result in a transition to hexagonal FeS V [6,7,13,14]. Fig. 2 shows the Gibbs free energy per formula unit of FeS I and FeS II at p ¼ 0 GPa. The curve of FeS I is below that of FeS II by approximately 50 meV, which suggests that the formation of FeS II is not allowed by thermodynamics. We did not calculate the Gibbs free energies of FeS IV and FeS V, due to the strong presence of imaginary frequencies (up to 10%) in their phonon spectra at all the volumes investigated. Approaches such as the renormalisation or the neglect of the imaginary frequencies can be adopted when these constitute a negligible fraction of the quasiharmonic phonon spectrum [41,42]. However, the application of these schemes to systems like FeS IV and FeS V, where the fraction of unstable modes is very large, may not be reliable, as at high temperatures the vibrational free energy is very sensitive to the lowest phonon frequencies (equation (2)). However, our findings are consistent with FeS IV being the result of the α transition and FeS V occurring at even higher temperatures. In this scenario, it is reasonable to expect the phonon branches of FeS IV and FeS V to be stabilised by anharmonic contributions that a quasi-harmonic treatment cannot reproduce. Fig. 3 shows the relative Gibbs free energy per formula unit of the FeS polymorphs as a function of pressure at 300 K. FeS I is the most stable phase up to a pressure of 7.0 GPa, after which FeS II becomes more favoured. A second transition takes place at about 15.0 GPa, when FeS II transforms to monoclinic FeS III. The two transitions are accompanied by volume reductions of 3 and 11% respectively. Experimentally, pressures of approximately 3 and 7 GPa have been found at room temperature for the two transitions above [16,39]. According to our results, FeS III undergoes a subsequent transformation to FeS VI at 61.2 GPa, which is characterised by an approximate 1% volume reduction. As in a previous DFT study [22], the order of appearance of the polymorphs as well as the successive reductions in volume agree with previous experimental works [6,16,18,21,39]. However, DFT þ U improves the theoretical accuracy both by predicting FeS I and FeS II to be stable at positive pressures and by shifting the pressure of the FeS III to FeS VI transition from 90 to 61.2 GPa, i.e. closer to the experimental value of 36 GPa [14]. We confirm, moreover, the prediction of a transition from FeS VI to FeS VII at 145.9 GPa [22], which has never been observed experimentally, as FeS VIII has been reported to form instead [23]. Unfortunately, the presence of unstable modes in the spectrum of FeS VIII, in addition to those reported above for FeS IV and FeS V, did not make it possible to include it in this analysis (similar to FeS IV and FeS V). Further studies incorporating anharmonic contributions, for example within the self-consistent phonon theory of anharmonic lattice dynamics [43], are necessary to shed light  on this controversial aspect [44].

Stability under pressure
In order to include FeS IV, FeS V and FeS VIII in the comparison of the relative stabilities under pressure, we have calculated the enthalpies of FeS I-VIII at 0 K neglecting the zero point energies. Fig. 4 shows that despite the fact that the transition pressures shift significantly to 5.5, 12.5, 76.7 and 144.7 GPa, the picture emerging from the analysis above has not altered. In particular, we find that it is not possible for FeS I to transform to FeS IV by only increasing pressure, as the two curves are separated by a significant 70 meV difference, which seems to contradict X-ray diffraction experiments [7], while indirectly supporting the presence of a MnP-type FeS II intermediate [6]. Finally, we note that there is no pressure range of stability for either FeS V, which is in line with all previous literature reports [6,39,45], or for FeS VIII, which agrees with previous calculations [22], but it is at odds with experiment [23]. However, we stress that these findings should be treated with a degree of caution, as thermal effects have been neglected in Fig. 4.

Phase diagram
By evaluating the Gibbs free energy of the dynamically stable phases at different temperatures and pressures, we have derived the phase diagram for the FeS system. We show the results in Fig. 5, where they are completed with the boundaries involving FeS IV and FeS V, which we have speculated on the basis of the results by Fei et al. [6]. We have included data only up to a pressure of 100 GPa, as we could not investigate the appearance at higher pressures of FeS VIII (see subsection 3.3) [23]. During the transition from FeS I to FeS II, contrary to the current understanding in the literature [16], the critical temperature increases with pressure. Similarly, we find a positive slope for the transition from FeS II to FeS III. In contrast, the phase boundary between FeS III and FeS VI has a negative pressure dependence, which supports X-ray diffraction measurements [14].
Changes in entropy ΔS and volume ΔV across the transitions can be related to the gradient of the phase boundaries by the inverted Clausius-Clapeyron equation: Given the monotonic decrease in volume with pressure, our results imply that the entropy decreases across the transitions from FeS I to FeS II and from FeS II to FeS III, whereas it increases across the transition from FeS III to FeS VI, thus playing a role in the destabilisation of monoclinic FeS III at high temperatures.

Thermodynamic properties
The quasi-harmonic approximation is based on the dynamic stability of the phase under investigation. As discussed in subsections 3.2 and 3.3, because of the significant presence of imaginary frequencies in their spectra, we could not calculate the vibrational free energies of FeS IV, FeS V and FeS VIII, which precluded us from deriving their thermodynamic properties. Given that these polymorphs do not present a stability window for the enthalpy at 0 K (Fig. 4), this section could be considered as an investigation of the thermodynamic properties of the FeS polymorphs existing at low temperatures (we recall that thus far FeS VII has only been predicted theoretically). In particular, we studied the volume expansion behaviour, heat capacity and entropy of FeS I-III, FeS VI and FeS VII at representative values of pressure falling within the stability ranges resulting from the Gibbs free energy plots of Fig. 3 (Table 2). Fig. 6 shows the dependence on temperature of the relative volume expansion V rel of the FeS polymorphs at different pressures with respect to their reference volume V 0 at 300 K, i.e: For FeS I, the agreement with experiment is very good in the range between 300 and 400 K. The discrepancy at higher temperatures is in line with the occurrence of the α transition to a different polymorph, which our results in subsection 3.2 suggest should be FeS IV. At a temperature of 400 K, we have obtained for FeS I a value of the thermal expansion coefficient, defined as   of 4.9 Â 10 À5 K À1 , which is significantly different from that of 14.1 Â 10 À5 K À1 measured in the range 373-573 K [46]. Finally, we note that the thermal expansion of the FeS phases reduces with pressure ( Table 2). The upper panel of Fig. 7 shows the molar heat capacity at constant pressure as a function of temperature, which we have estimated by means of the thermodynamic relation [48]: where BðTÞ is the bulk modulus at temperature T. The predicted heat capacity of FeS I at 0 GPa agrees very well with the experimental value up to 200 K, after which some differences between the two start to be visible. In particular, our data form the typical sigmoidal curve, while the experimental points mark an upward trend. At 350 K, our underestimation reaches its maximum value, with the theoretical prediction of 47.9 J mol À1 K À1 against the experimental value of 57.2 J mol À1 K À1 . We stress that the interatomic force constants calculated within the finitedisplacement method can be very sensitive to the size of the supercell adopted [49]. It is therefore important to ensure that the observed discrepancy is not caused by such an effect. To this end, in the inset of Fig. 7, we show that the heat capacity of FeS I at 0 GPa is converged well with respect to the supercell size, as absolute differences between 1 Â 1 Â 1 and 2 Â 2 Â 1 supercells, containing respectively 24 and 96 atoms, are less than 1 J mol À1 K À1 . Thus, the deviation of the theoretical heat capacity from experiment parallels the behaviour of the volume expansion with temperature and is in line with the attribution of the experimental upward trend to the approach of the α transition [47]. The lower panel of Fig. 7 shows the dependence on temperature of the molar entropy. For FeS I, the agreement between calculated and measured results is excellent up to the highest experimental value of 68.9 J mol À1 K À1 at 350 K, with our value of 67.6 J mol À1 K À1 differing by less than 2%. The perfect match between the two curves suggests that the differences encountered in the heat capacities should be ascribed to the increasing fraction of FeS IV with temperature before the α transition, rather than to inaccuracies in our methodology or neglected anharmonic terms. Most importantly, it gives us significant confidence about the thermodynamic predictions for the high pressure phases (Table 2), for which anharmonic effects are expected to become even smaller [50].

Conclusions
In this work, we have presented an investigation of the FeS polymorphs based on Hubbard-corrected density functional theory within the quasi-harmonic approximation. We have conclusively shown that, at ambient pressure, FeS I cannot transform to MnP-type FeS II as reported previously [16,40], which indirectly supports a transition towards FeS IV [6]. Under pressures up to 100 GPa, although some discrepancies between calculated and measured transition pressures exist, the order of appearance of the polymorphs as well as the discontinuities in the volumes are in full agreement with the literature. Our study suggests that some caution should be taken when neglecting temperature effects, as they can affect non negligibly the transition pressures. We found that the hexagonal FeS IV and FeS V, as well as the cubic FeS VIII structures are dynamically unstable at 0 K, which prevented us from assessing their stability in the phase diagram. However, this does not contradict experiments, as these polymorphs occur in nature at high temperatures, where anharmonic contributions may not be negligible anymore. Future investigations based on a level of theory beyond the quasi-harmonic approximation may help to elucidate this aspect. Finally, we stress that the good agreement between theory and experiment for the dependence on temperature of volume expansion, heat capacity and entropy at zero pressure of FeS I lends support to our predictions for the other, less experimentally accessible, FeS phases which are speculated to have a role in geophysics.