Two-step nucleation of the Earth's inner core.

The Earth's inner core started forming when molten iron cooled below the melting point. However, the nucleation mechanism, which is a necessary step of crystallization, has not been well understood. Recent studies have found that it requires an unrealistic degree of undercooling to nucleate the stable, hexagonal, close-packed (hcp) phase of iron that is unlikely to be reached under core conditions and age. This contradiction is referred to as the inner core nucleation paradox. Using a persistent embryo method and molecular dynamics simulations, we demonstrate that the metastable, body-centered, cubic (bcc) phase of iron has a much higher nucleation rate than does the hcp phase under inner core conditions. Thus, the bcc nucleation is likely to be the first step of inner core formation, instead of direct nucleation of the hcp phase. This mechanism reduces the required undercooling of iron nucleation, which provides a key factor in solving the inner core nucleation paradox. The two-step nucleation scenario of the inner core also opens an avenue for understanding the structure and anisotropy of the present inner core.

The core plays a key role in the Earth’s evolution. The present core contains two major parts, a solid inner core and a liquid outer core. Iron dominates both parts with a small amount of light elements (1). The solid core is generally believed to be hexagonal, close-packed (hcp) iron, while the possible existence of body-centered, cubic (bcc) iron has also been suggested (2–5). The growth of the solid inner core is believed to be the major driving force of the present geodynamo, providing the main power source for convection in the liquid core (6, 7). Despite its importance, the initial formation of the solid core, which directly relates to its thermal evolution and Earth’s history, is far from being completely understood (8–12). Most of Earth’s thermal history models assume that the inner core started to crystallize when molten iron cooled right below its melting temperature at the Earth’s center (7). However, in practice, nucleation does not happen at the melting point but requires some undercooling because of the formation of a solid–liquid interface (SLI) that accompanies it. While the bulk solid phase is thermodynamically favored, the SLI costs energy. These two factors lead to a nucleation barrier , which is described in classical nucleation theory (CNT) (13) as

where is the nucleus size, is the free energy difference between the bulk solid and liquid, is SLI free energy, and is the SLI area. The liquid must be cooled sufficiently below the melting temperature to overcome the free-energy barrier during thermal fluctuations. After considering this mechanism, it was found that a very large undercooling of ∼1,000 K is required for the nucleation of hcp iron in the Earth’s core (14). However, considering the slow cooling rate of ∼100 K/Gyr throughout the core history (15), it is impossible to reach such a large degree of undercooling inside the Earth within the inner core’s age. This “inner core nucleation paradox,” recently described by Huguet et al. (14), strongly challenges the current understanding of the inner core formation process. While Huguet et al.’s argument relies on a few estimations of thermodynamic quantities, Davies et al. also confirmed the paradox with atomic-scale simulations (16). Even considering the effect of light elements on the nucleation process, it still requires 675 K undercooling to nucleate hcp iron, nearly impossible to reach in the Earth core (16).

CNT was proposed more than a century ago, and its formalism is the most widely used to describe nucleation phenomena nowadays. The simplest scenario in CNT assumes a single-nucleation pathway where only the nucleus of the thermodynamically stable phase forms and grows toward the bulk phase. This was the situation considered in refs. 14 and 15), in which the authors assumed that the melt in the Earth’s core crystallized directly into the hcp phase. Recent studies have shown that nucleation can be a multistep process that includes multiple intermediate stages and phases (17–19). While the CNT concept of nucleus formation is still valid under these situations, phase competition must be considered (18, 19). Therefore, instead of the single-pathway scenario, we can consider a complex process in which nucleation is facilitated by forming an intermediate phase with a high-nucleation rate. For example, it has been observed that the bcc phase can nucleate before the face-centered, cubic (fcc) or hcp phases in a few alloys in which the fcc/hcp phase is the most stable one (20–24). Could the bcc phase also facilitate hcp iron nucleation and relate to the inner core nucleation paradox? Making a quantitative prediction on such complex nucleation processes is a challenging problem. In addition to the extreme conditions in the core, nucleation involves microscopic-length scales that are extremely hard to probe in real time, even with state-of-the-art measurements (25). Hence, it requires computer simulations, particularly large-scale molecular dynamics (MD), to reproduce the temporal evolution of the liquid into the crystal (26). Unfortunately, nucleation under Earth’s core conditions is a rare event that occurs on the geological time scale, far beyond the reach of conventional MD simulations. Besides, large-scale MD simulations require semiempirical potentials to describe atomic interactions, and the outcome may depend heavily on the potential’s quality (27). In this work, we assess the inner core nucleation process with the account of competition between bcc and hcp phases during the nucleation process using the persistent embryo method (PEM) (28) to overcome the significant time limitation in conventional MD simulation of nucleation.

Results

Melting Curve of hcp and bcc Phases.

The nucleation rate of hcp iron was previously estimated (14) based on the driving force and the SLI free energy obtained in ref. 29, with the semiempirical potential developed by Ackland et al. (30). However, this potential was developed to simulate iron at ambient conditions such that no high-pressure data were used in the potential development (30). In the present study, we developed a potential, explicitly considering its application at Earth’s core conditions. contains details of the present potential development. One of the vital target properties in the potential development is the latent heat , because along with the melting temperature it defines the driving force for solidification. Fig. 1 shows excellent agreement between the latent heat calculated using the developed potential and ab initio MD (AIMD) for both the hcp and bcc phases. Besides, elastic properties and liquid structure predicted with the developed potential also agree well with those calculated with AIMD (), making this potential suitable for simulations of the iron crystallization process under inner core conditions.

Fig. 1.

Melting curves and thermodynamic properties of hcp and bcc iron. (A) The latent heats of hcp and bcc iron at 323 GPa from AIMD and CMD with the developed semiempirical potential. (B) Comparison between the current melting curve of hcp iron and several others obtained by previous simulations and experiments (16, 31–38). (C) Melting curves of hcp and bcc iron from CMD with the developed semiempirical potential. The inset shows a projected snapshot of the bcc solid–liquid coexistence simulation with 22,500 atoms in the simulation cell. (D) Change in bulk free energy upon solidification (nucleation-driving force) at 323 GPa.

In Fig. 1, we compare the current melting curve of the hcp phase with previous measurements and simulation results (16, 31–38). This is especially important because melting curves from previous classical MD (CMD) simulations deviate considerably from each other (31, 32), which points out the importance of employing a thoroughly developed, semiempirical potential. The melting temperatures in the present study were determined at several pressures using the solid–liquid coexistence approach (39). The current melting curve agrees well with the experimental curve obtained using fast X-ray diffraction in the laser-heated diamond anvil cell in the pressure range between 130 and 200 GPa (37). It also reasonably agrees with the recent estimation of melting boundary from shock compression measurements in the higher-pressure range from 250 to 360 GPa (38). Compared to previous simulations, our current melting curve provides the closest agreement to the recent, high-pressure experiments for the hcp phases. It also validates the calculations of latent heat, which is directly related to the slope of the melting curve according to the Clausius–Clapeyron equation.

The melting curves of hcp and bcc phases are compared in Fig. 1. The hcp phase has higher melting temperatures in the pressure range from 130 to 330 GPa. Thus, the hcp phase is thermodynamically stable, while the bcc phase is metastable in the range from the core–mantle boundary to the inner core boundary. Interestingly, the melting points of hcp and bcc phases are predicted to cross at 350 GPa, which suggests a similar free energy of bcc and hcp phases at pressures near the inner core center.

In the present study, we chose to conduct nucleation simulations at 323 GPa, the pressure at the inner core boundary. We used the Gibbs–Helmholtz equation and MD simulations to calculate the free-energy difference between the bulk solid and liquid phases, as described in ref. 40. Fig. 1 shows as a function of undercooling with respect to the hcp melting temperature, , for both the hcp and bcc phases. The absolute value of for the hcp phase is always larger than that for the bcc phase.

Nucleation of hcp and bcc Phases.

According to the CNT (13), the nucleation barrier, , is the key quantity to determine the nucleation rate. can be computed (28) aswhere is the critical nucleus size. To obtain , we employ the PEM (28), which invokes the central CNT concept that homogeneous nucleation happens via the formation of a critical nucleus in the undercooled liquid. Fig. 2 shows a typical result of the PEM-MD simulation. The plateau on the N(t) curve indicates the appearance of the critical nucleus (see technical details in Materials and Methods). The critical nucleus sizes of both the hcp and bcc phases at several moderate undercooling temperatures are shown in Fig. 2. The hcp phase shows a systematically larger critical nucleus size than the bcc phase. Fig. 3 shows the free-energy barriers of both bcc and hcp phases computed using Eq. . The hcp phase has a larger nucleation barrier than the bcc phase at all undercooling temperatures considered here, although the hcp phase has a larger bulk-driving force, (Fig. 1). To explain this, we compare the SLI free energies of these phases obtained from the PEM simulation. Fig. 3 shows that this quantity is much larger for the hcp phase. This is consistent with a previous study which suggests that bcc metals show lower interface free energy than hcp metals at ambient conditions (41). Because the scales with (see Materials and Methods), the difference in can significantly change the ratio of the nucleation barriers. The temperature dependences of obtained for both hcp and bcc phases are almost linear as shown in Fig. 3, which is similar to the ones found for Ni and Al in ref. 42. Therefore, the SLI free energy can be linearly extrapolated to smaller undercooling (higher temperatures), in which the critical nucleus size is too large to be simulated directly. In Fig. 3, we further compare values for the hcp phase with previous results. Davies et al. estimated based on a different, semiempirical potential and freezing simulations (16). While the temperature dependence was not considered in that work, the value of is very similar to our results. In contrast, the determined in ref. 29 highly deviates from our results and those in ref. 16. This can be attributed to the fact that the empirical potential used in ref. 29 was not designed to simulate iron under Earth’s core conditions.

Fig. 2.

PEM-MD simulation and critical nucleus size. (A) The nucleus size versus time from a PEM-MD trajectory at (i.e., ). The blue dashed line shows the size of the persistent embryo, N0, and the green dashed line indicates the threshold for spring removal, Nsc (see Materials and Methods). The inset enlarges the plateaus at the critical size. The red line shows the plateaus to determine the critical nucleus size, N*. (B) The critical nucleus size as a function of undercooling temperature for hcp and bcc phases.

Fig. 3.

Temperature dependence of the free-energy barrier and SLI free energy. (A) Free energy barrier as a function of undercooling temperature for hcp and bcc at 323 GPa. (B) SLI free energy of hcp and bcc phases at 323 GPa. The dashed dot line is from ref. 14. The dashed line with the CI (green band) is from ref. 16.

The nucleation rate, , can be calculated as , where is the Boltzmann constant, and is a kinetic prefactor. The kinetic prefactor can be determined from MD simulations, based on the fluctuations of the nucleus size around the critical value (). Thus, we obtained all essential parameters to compare the bcc and hcp nucleation rates using the PEM in MD simulations. Fig. 4 shows that the bcc phase has a much higher nucleation rate than the hcp phase in a broad undercooling regime. For example, at an undercooling of 660 K, the nucleation rate of the bcc phase is 31 orders of magnitudes higher than that of the hcp phase. With such a vast difference, the bcc phase should nucleate much quicker than the hcp phase under the Earth’s core conditions.

Fig. 4.

Nucleation rate and waiting time at 323 GPa. (A) Nucleation rate as a function of undercooling for the hcp and bcc phases. (B) Waiting time as a function of undercooling. The dotted lines are the extrapolation with CNT from current PEM results. The dashed line is from ref. 16. The thick black line indicates the range of nucleation waiting time in the Earth’s core.

Using the obtained nucleation rate, we are now able to estimate the nucleation waiting time. Because the critical nucleus only has half a chance to grow at the top of the nucleation barrier, the waiting time in a fixed volume can be expressed as (16). As shown in Fig. 4, the present waiting time for hcp nucleation is smaller than the previous estimation by Davies et al. (16). This is mainly because the semiempirical potential used in ref. 16 shows a different melting point and bulk free-energy difference than the present values (). The reachable in the Earth’s core can be estimated as follows. The nucleation incubation time is approximated as 1 billion y, probably the upper limit of plausible inner core age (16). The volumes of the inner core and the entire core are ∼ to ∼, respectively. Therefore, of the Earth’s core should be in the range between to , which is indicated by the solid black line in Fig. 4. The intersections of this line with the hcp and bcc’s in Fig. 4 provide the required undercoolings of 470 K for bcc nucleation and 610 K for hcp nucleation. Therefore, the bcc phase significantly reduces the required undercooling of inner core nucleation by 140 K, which corresponds to ∼1.4 billion y of cooling based on a cooling rate of 100 K/Gyr (8).

Discussion

As illustrated in Fig. 5, we have shown that the bcc nucleation is likely to be the first step of iron crystallization under core conditions, which effectively decreases the large nucleation barrier of the hcp phase. We note that this two-step nucleation process can also be observed from the brute force MD simulation by cooling the iron melts with ultrahigh cooling rates (). While this cooling simulation is far from the core condition, it provides a qualitative validation of the PEM results. As the two-step nucleation process reduces the required undercooling of inner core nucleation, it can be a key to solving the core nucleation paradox. While the necessary undercooling of the bcc phase (470 K) is still more significant than the inner core’s maxima undercooling of ∼200 K (14), a few factors remain to be accounted for. First, the core contains ∼10 weight percent light elements. Davies et al. (16) have shown that oxygen can reduce the required undercooling by ∼55 K for the hcp nucleation. It was also found that the light elements can further stabilize the bcc phase with respect to the hcp phase (4, 43, 44). Therefore, their presence in the melt should cause a similar, if not larger, reduction of the required undercooling for the bcc nucleation. Moreover, stresses associated with pressure variations can also reduce the undercooling by ∼100 K (16). Considering all these effects, it should be possible to solve the nucleation paradox with the future development of semiempirical or machine-learning potentials for iron and light elements and more sophisticated PEM nucleation simulations of a multicomponent melt. Our simulations consider homogenous nucleation only. While heterogeneous nucleation might decrease further the required undercooling, it is challenging to propose plausible scenarios to explain a preexisting stable substrate in the early deep core. This issue has been thoroughly discussed along with the nucleation paradox proposal in Huguet et al. (14). Besides, for any proposal of a preexisting substrate in the melt, the competition between bcc and hcp phases should still be considered to address the heterogeneous nucleation of the solid core. Then, our free-energy results for bcc and hcp phases will also be relevant because the free-energy barrier for heterogeneous nucleation, , includes the term , where is the homogenous nucleation barrier, and is the catalytic factor that depends on the contact angle, size, and substrate shape (13). Recent high-pressure and high-temperature experiments have reported the formation of bcc iron on the MgO surface (45), which indicates that the bcc should be the first nucleated phase in a heterogeneous nucleation process.

Fig. 5.

Schematics of the two-step nucleation process. The insert shows a bcc nucleus spontaneously formed during the brute force MD simulations described in .

The two-step nucleation with the intermediate bcc phase can impact the present inner core structure. While the metastable bcc phase should eventually transform into the stable hcp phase at 323 GPa, the melting curves in Fig. 1 suggest that the bcc phase could stabilize over the hcp phase with increasing pressure when approaching the core center. Therefore, the initially formed bcc phase in this region may remain at the core center nowadays. This coincides with hypothesis of a different innermost core structure proposed to explain the anomalous anisotropy of the inner core (46–48). Therefore, the two-step nucleation scenario for inner core formation opens a path to understanding Earth’s deepest interior.

Materials and Methods

AIMD Simulations.

AIMD simulations were performed to obtain the input data for the semiempirical potential development. The Vienna ab initio simulation package (49) was employed for the density functional theory (DFT) calculations. The projected augmented wave method was used to describe the electron–ion interaction and the generalized gradient approximation in the Perdew–Burke-Ernzerhof form was employed for the exchange–correlation energy functional. The Mernin functional was used to equilibrate electrons and ions at the same temperature in AIMD simulations (50). This functional includes electronic entropic effects on the DFT energy. The Γ point was used to sample the Brillouin zone. The AIMD simulations were performed for the constant number of atoms, volume, and temperature ensemble. The Nosé–Hoover thermostat was employed to control the temperature. A time step of 2.0 fs was used to integrate Newton’s equations of motion. Supercells with 288, 250, and 256 atoms were used to simulate hcp, bcc, and liquid models, respectively. To fit the potential for a sizeable pressure-temperature (P-T) range, three P-T conditions near the melting curves were investigated: 140 GPa at 4,000 K, 250 GPa at 5,500 K, and 350 GPa at 6,000 K. No phase transitions were observed during the simulation of either hcp, bcc, or liquid model under these conditions during the AIMD simulations. To save the equilibration time, the AIMD simulations were run iteratively with the potential development, such that configurations equilibrated with the semiempirical potential were used as initial configurations for AIMD. We monitored the energy and pressure as a function of time to determine when the models were equilibrated enough to start collecting data. Data collection took place in the last 10 ps of each AIMD run. The potential was fitted to reproduce AIMD results using the Mermin functional, which includes electronic entropic effects.

CNT Formulae.

Based on the CNT (13), the competition between the bulk and interface energies leads to the nucleation barrier. From Eq. , one can derive the nucleation barrier aswhere is the crystal density, and is a shape factor so that the interface area in Eq. can be written as . This equation can be reduced to Eq. by introducing the critical nucleus size

The SLI free energy, , is obtained from Eq. . All quantities to determine can be obtained from PEM-MD simulations.

Persistent Embryo Method.

The PEM (28) is employed to measure the critical nucleus size with the CMD simulations. It utilizes the main CNT concept that homogeneous nucleation happens via the formation of the critical nucleus in the undercooled liquid. During the simulation, a small crystal embryo containing atoms (should be much smaller than the critical nucleus) is constrained by spring forces to prevent melting (28). These forces are only applied to the original embryo atoms. The spring constant of the harmonic potential decreases with increasing nucleus size as if and otherwise. Here, is a subcritical threshold. No spring forces are applied to the embryo if . This strategy ensures that the system is unbiased at the critical point, such that a reliable critical nucleus can be obtained. If the nucleus melts below , the harmonic potential is gradually enforced, preventing the complete melting of the embryo. When the nucleus reaches the critical size, it has an equal chance to melt or to further grow, causing fluctuations around . As a result, the curve tends to display a plateau during the critical fluctuations, giving a unique signal to detect the appearance of the critical nucleus, as shown in Fig. 2. The critical nucleus size is directly measured by averaging the nucleus size at the plateau (28). We repeated the PEM-MD simulation to collect at least four plateaus to obtain sufficient statistics and the CI of the critical nucleus size . During the nucleation process, the bcc and hcp nuclei are stable, and no transformations to other solid phases were observed. The classical PEM-MD simulations were performed with the graphic processing unit–accelerated large-scale atomic/molecular massively parallel Simulator code (51). The interatomic interaction was modeled using the semiempirical potential developed in this work based on the embedded atom method (52). The Nosé–Hoover thermostat/barostat was used in our NPT (constant number of atoms, pressure, and temperature) simulations. The damping time in the Nosé–Hoover thermostat was set as which is frequent enough for heat dissipation during the crystallization (42). The time step in the simulation was 1.0 fs. In some cases, in which the nucleation was very fast, the time step was reduced to 0.8 fs. The simulation cell contained 31,250 atoms, which are at least 15 times larger than the critical nucleus size. The simulation cell is sufficiently large to avoid the finite size effect (see ). The semiempirical potential of iron developed in this work is publicly accessible at National Institute of Standards and Technology (NIST) potentials repository (https://www.ctcms.nist.gov/potentials/).