Quenching of bcc-Fe from high to room temperature at high-pressure conditions: a molecular dynamics simulation

The new high-temperature (T), high-pressure (P), body-centered cubic (bcc) phase of iron has probably already been synthesized in recent diamond anvil cell (DAC) experiments (Mikhaylushkin et al 2007 Phys. Rev. Lett. 99 165505). These DAC experiments on iron revealed that the high-PT phase on quenching transforms into a mixture of close-packed phases. Our molecular dynamics simulation and structural analysis allow us to provide a probable interpretation of the experiments. We show that quenching of the high-PT bcc phase simulated with the embedded-atom model also leads to the formation of the mixture of close-packed phases. Therefore, the assumption of the stability of the high-PT bcc iron phase is consistent with experimental observation.


2
The problem of establishing the iron phase diagram has long been central to many fields of metallurgy and materials science. The reason for this is obvious because of the practical importance of this material. For the last 50 years the 'extreme' high-PT part of the iron phase diagram has attracted significant attention. Apart from military applications, the importance of the high-PT diagram is related to the Earth's core. It is established that the Earth's inner core (IC), a solid spherical body about 2400 km in diameter, consists mainly of iron. It has long been thought that at pressures (P) 330-360 GPa and temperatures (T) above at least 5000 K, i.e. conditions corresponding to the IC, iron is stable in the hexagonal close-packed (hcp) phase [2]. However, the experimental work by Brown and McQueen [3] suggested that iron might transform to another, different from the hcp, phase at the conditions of IC. That discovery [3] was not confirmed by later shock-wave experiments [4]; however, the difference in compositions of iron samples could be the reason for that (see [5] for discussion). It was suggested, on the basis of non-empirical [6] molecular dynamics (MD) simulations, that under high pressure iron transforms on heating to the body-centered cubic (bcc) phase [7]. A number of IC properties can be explained within the paradigm of the bcc iron phase stability under the IC PT conditions, such as the low shear modulus [8,9] and elastic anisotropy [10]. However, computations performed by Vocadlo et al [11] have not confirmed the bcc stability.
Except for the measurements of velocities of rarefaction waves behind the shock-wave in iron [3,4], no other experiments have either confirmed or ruled out the stability of another, different from the hcp, iron phase at the IC conditions. Recently, the hcp-bcc transition was observed in iron slightly alloyed with nickel at pressures around 200 GPa on temperature increase. This experiment was performed in a diamond anvil cell (DAC) [12]. The temperatures of the transition are in agreement with earlier prediction [7], provided that the correction for alloying is made [12]. However, the confirmation of the bcc stability in pure iron at the conditions in static high pressure devices is still lacking. Recently, another DAC experiment, this time on pure iron, was reported [1]. The structural measurements were only performed on quenched samples. The samples were pressurized, heated, and then the heating was switched off. After that, x-ray spectra were taken from cold samples under pressure. Depending on the temperature of heating, the x-ray diffraction (XRD) patterns of quenched samples showed the presence of either only hcp phase or both hcp and face-centered cubic (fcc) phases (if heated above some PT boundary and then quenched). It was suggested that the appearance of the fcc peaks is related to the stabilization of the high-T fcc phase proposed on the basis of firstprinciples calculations [1].
We note, however, that such a straightforward relating of phases in both in situ and quenched samples might not always be correct. For example, Zr, stable in the hcp phase at ambient conditions, transforms into the bcc phase on heating. The Zr bcc phase, when quenched from high to room temperature, transforms into a martensitic structure with hcp and fcc fragments [13], very much like what was observed in the experiment on iron [1].
Therefore, the appearance of the fcc phase in a quenched sample is not necessarily evidence for the stability of the fcc phase at high temperature. Moreover, we know from other experiments on iron [14] that when iron is quenched from the stability field of the fcc phase, the XRD pattern of the quenched sample shows lines belonging to hcp and double hcp phases. Therefore, it is not necessary that the fcc peaks in XRD patterns of quenched samples are due to the stabilization of the fcc phase at high temperature. On the other hand, robust and reproducible changes in the quenched samples are indeed evidence that a structure change takes place on heating. When we placed the experimental points on the phase diagram of Fe, as was proposed in 2003 [7], we noticed the following. If the temperature of heating was above the hcp-bcc boundary [7] (figure 1) then the quenched samples consisted of hcp + fcc mixture. Otherwise, only the hcp phase was present. This suggested that the bcc phase could have been synthesized in the DAC experiments. We decided to simulate the experimental process of quenching, assuming that the high-T phase is bcc.
The simulation of quenching the non-magnetic high-PT bcc phase of Fe was performed by the MD method. Fe atoms interacted via the embedded-atom model (EAM) potential [6] parameterized solely on the basis of first principles calculations, without any experimental data. This EAM has been proved to perform very well. For example, it allows one to reproduce the structure of liquid iron [14]- [16] as well as its viscosity [7,17]. It also allows one to explain a number of features of IC [9,10]. The melting temperature of Fe at the pressure of the outer core-inner core boundary (330 GPa), according to the EAM, is consistent with the temperature computed from first principles [18] (later revised to a somewhat lower temperature) within the mutual error bars. The EAM melting temperature is in agreement with the shock-wave data [3] at a pressure of 243 GPa (figure 1). The most recent melting data by Ma et al [19] was predicted almost exactly back in 2000 [6]. Therefore, we are confident in this method.
Starting from the bcc structure, we equilibrated it at the experimental PT of 165 GPa and 3700 K [1]. The computational cell consisted of 1,024,000 atoms and was obtained by 80 × 80 × 80 multiplications of the bcc unit cell that contained two atoms. That is, the initial structure represented a single bcc crystal without grain boundaries. The 3D periodic boundary conditions were applied. The details of simulations are the same as reported earlier [6]. The  [20,21]. Blue (108 072 atoms)-fcc up to fourth nearest neighbour (nn), dark blue (141 100 atoms)-fcc up to first nn, red (403 850 atoms)-hcp up to fourth nn, dark red (272 642 atoms)-hcp up to first nn, grey (34 196 atoms)-12 coordinated atoms that are not bcc, fcc, or hcp, dark grey (64 140 atoms)-8 coordinated atoms that are not bcc, fcc, or hcp. Hcp and fcc regions delineated by 111 planes, suggesting that, indeed, during the quenching a displacive transition has occurred from bcc to the hcp+fcc phase. Such a phase transition is normally associated with the development of a soft phonon mode. timestep in our simulations was 1 fs. After the equilibration in the NPT ensemble, we saved the last configuration. The structure was clearly the bcc one. Then, starting from the saved configuration and preserving the volume, we ran 10 000 timesteps scaling temperature to 300 K. After that, the scaling was switched off and the system was simulated for another 500 000 timesteps. As a result of quenching, the pressure dropped to 128 GPa. The atoms showed practically no diffusion during the last 100 000 timesteps. Therefore, the last configuration can be considered as representative of the quench process.
This last configuration was analyzed by applying the medium range order (MRO) analysis [20,21]. The atoms were treated one by one and assigned a certain structure according to their local crystal structure. The simulated sample (see figure 2) represents a granular network, with the grain boundary network being defined by the grey and dark grey atoms which do not have an assigned crystallography. Inspection of the sample shows that grains which are predominantly hcp (red and dark red atoms) and fcc (blue and dark blue atoms) are both present. This is in full agreement with the structure analysis performed in the experiment on quenched samples [1]. We also computed x-ray spectrum for the quenched sample and compared it to the measured one ( figure 3). The spectra are similar; however, there are differences in the intensity of peaks as well as in the lattice constant. We note that the resulting pressure after quench was not measured and that the pressure could in principle be different between the simulation and experiment. Also, grains in the experiment are most likely larger than in the simulation. The  [22,23] (red) and experimental (blue). Note that in calculations the atomic form factor is assumed to be independent of scattering angle and set equal to unity. The calculated pattern is shifted to the left, indicating that the hcp + fcc lattice constant is about 3% larger than that seen in experiment. When this correction is taken into account, all simulated peaks are present in the experimental data, although with significant differences in both peak heights and widths. Such differences must be expected given different sizes of experimental and simulated samples. Some of the computed peaks are, therefore, appearing as shoulders.
ratio of fcc and hcp phases is clear from the numbers of atoms in the different coordinations (figure 2). The atoms in the hcp phase are the most abundant. The typical linear size of the grains is about 40 Å.
In principle, the results of an MD simulation of quenching might be affected by the size of a sample and a quenching rate. We did quenching simulation for considerably smaller samples [7] and did not notice any impact on structure. As for the quenching rate, one might expect different phase ratio if quenching is done slower, but the general picture will likely remain the same because of the bcc dynamical instability.
Why is the quenched structure a mixture of close-packed phases for the samples quenched from above the hcp-bcc boundary? The bcc phase is dynamically unstable at low temperature. It becomes dynamically (and eventually thermodynamically) stable above a certain critical temperature [7] (figure 1). As soon as the samples are quenched from the PT field of bcc-Fe stability, they become dynamically unstable. The atoms fall from the positions corresponding to the state with high potential energy in the bcc structure into the positions corresponding to the states with low potential energy in the hcp and fcc structures. Because the quench is a very fast process, the atoms fall from the positions corresponding to the high energy states to those corresponding to the low energy states without probing the entropy of the states. This is why the quenched structure is a metastable mixture of hcp and fcc atoms, contrary to the stable hcp structure. The speed of structural changes is clear from figure 4. The fragments of the fcc and The microstructure of the sample after first 500 steps of quenching. Green, red and blue correspond to the bcc, hcp, fcc structures, correspondingly; radial distribution functions computed at different times of quenching as indicated in the legend. After the first 500 timesteps the first peak is still split, indicating the initial bcc structure is still present. After 10 000 timesteps the splitting disappears; however, the magnitude of the RDF between the first and second peaks is still substantial, indicating a significant amount of disorder. Yet, no bcc phase is present. This means that after 10 000 timesteps, that is, after 10 ps, the collapse of the unstable bcc structure is essentially completed. This time (10 ps) is much shorter than the kinetics of the hcp-fcc transition. This is why both phases are present in the quenched sample. The RDF averaged over the following 500 000 timesteps is the same as the one computed for the last configuration (the corresponding curves are on top of each other). This indicates that the structure no longer changes. 7 hcp phases appears almost immediately within the first 500 MD steps ( figure 4(b)). One can see that just after the first 10 ps the structure has no characteristic features of the bcc structure. Namely, the splitted first peak of the radial distribution function (RDF) disappears and the uniform sharp first peak is formed (see the RDF for 10 ps in figure 4(a)). This is clear evidence for the close-packed structure. The structure becomes less disordered in the next 0.5 ns (note the tail of RDF at 10 ps and 0.5 ns, figure 4(a)). Also, the atoms move from the range between the first and second peaks. The collapse of the bcc structure is much faster than the kinetics of the solid-solid transitions in iron.
Our major finding is that the quenched samples in both experiment and simulation consist of a mixture of hcp + fcc, provided that the conditions before quenching are above the hcp-bcc boundary (figure 1). We conclude, therefore, that the assumption of the bcc stability allows one to reproduce the experimental observations. Considering very close proximity of the experimental points to the hcp-bcc phase boundary [7], such an assumption seems very probable. We hope that future experiments will fully clarify this issue.