Earth and Planetary Science Letters Assessing the inner core nucleation paradox with atomic-scale simulations

We investigate the conditions required to freeze liquid iron and iron alloys near the centre of Earth’s core. It is usually assumed that inner core growth begins once the ambient core temperature falls below the melting temperature of the iron alloy at Earth’s centre; however, additional (under)cooling is required to overcome the energy barrier associated with creating a solid–liquid interface. Predictions based on Classical Nucleation Theory (CNT) have estimated a required undercooling of ∼ 1000 K, which cannot be reconciled with predicted core cooling rates of ∼ 100 K Gyr − 1 . This apparent contradiction has been called the ‘inner core nucleation paradox’. Here we address three major uncertainties in the application of CNT to inner core nucleation using atomic-scale simulations. First, we simulate freezing in Fe and Fe–O liquids at core conditions to self-consistently constrain all parameters required by the CNT equations. Second, we test the basic validity of CNT by directly calculating the waiting time to observe freezing events in Fe and Fe–O liquids. Third, we investigate the inﬂuence of wave-like forcings applied to the atomic simulations, which have been suggested as a means to signiﬁcantly reduce the energy barrier. Our results are consistent with CNT in the computationally accessible parameter regime, though error estimates on the waiting time can reach 50% of the measurement at the largest undercooling temperatures. Using CNT to extrapolate to inner core conditions yields estimated undercooling of 730 ± 20 K for the pure iron system and 675 ± 35 K for the Fe–O system. Forcings corresponding to large pressure variations of O ( 10 ) GPa reduce these values by ∼ 100 K. While our undercooling estimates are signiﬁcantly lower than previous estimates they are not low enough to resolve the inner core nucleation paradox.


Introduction
The formation of the solid inner core was a defining moment in Earth's history. Continual freezing of the liquid core releases latent heat and lighter elements, the latter providing the gravitational energy that is the main power source maintaining the present geomagnetic field (Braginsky, 1963;Gubbins, 1977;Nimmo, 2015). Prior to inner core formation, higher cooling rates are needed to explain the existence of the geomagnetic field for the last 3.5 Ga (Tarduno et al., 2010). If the inner core is 400-700 Myrs old, as suggested by high core thermal conductivity values (Pozzo et al., 2012;de Koker et al., 2012;Gomi et al., 2013), then conventional thermal history models predict early core temperatures far above current estimates of the lower mantle solidus (Nimmo, 2015;Davies et al., 2015;Labrosse, 2015), implying widespread melt-* Corresponding author.
E-mail address: c.davies@leeds.ac.uk (C.J. Davies). ing in early times. However, some studies have suggested that this rapid cooling scenario is unsustainable and argue that novel crystallization mechanisms powered the early magnetic field (O'Rourke and Stevenson, 2016;Badro et al., 2016;Hirose et al., 2017), slowing core cooling and allowing an inner core age of >1 Ga, though the efficiency of these processes has been questioned (Du et al., 2017). Clearly, the preferred scenario depends critically on the age of the inner core.
All current core thermal history models assume that the inner core nucleates when the melting curve for the liquid iron alloy falls below the core (adiabatic) temperature at Earth's centre (e.g. Nimmo, 2015;Davies, 2015;Labrosse, 2015). However, freezing of solid from liquid requires the creation of a solid-liquid interface, with an associated free energy barrier. The size of the barrier is determined by the competition between the (negative) volumetric and (positive) interfacial free energy; the top of the barrier corresponds to the critical radius, beyond which crystals will continue to grow (e.g. Christian, 2002). The time required for crystals to reach the critical radius is infinite if the system is at the melt-ing temperature and therefore some undercooling is needed for the system to reach equilibrium with solid as the stable phase below the melting temperature. The effect of undercooling is to delay the onset of inner core nucleation compared to predictions from existing models; larger undercooling leads to greater delay. Determining the amount of undercooling is therefore crucial for predicting the inner core age, which in turn impacts estimates of the power available to the dynamo and the thermal evolution of the core.
Consequences of undercooling at the inner core boundary (ICB) have long been discussed, mainly in connection with the possible existence of a mushy zone ahead of the ICB (e.g. Loper and Roberts, 1977;Fearn et al., 1981). Shimizu et al. (2005) used classical nucleation theory (CNT) to investigate the formation of crystals in a slurry layer above the ICB. They found that the probability of forming a slurry by homogeneous nucleation (nucleation in the absence of solid surfaces) would be essentially zero and suggested that nucleation must begin heterogeneously (i.e. on a pre-existing surface). Huguet et al. (2018) used CNT to investigate the origin of the inner core and found that a critical undercooling δT c ∼1000 K is needed to bring the waiting time for generating a stable nucleus below O (1) Gyr. Since the core is cooling at O (100) K Gyr −1 (Davies, 2015), such a large δT c is clearly impossible. Moreover, upon reaching an undercooling of ∼1000 K most of the core would rapidly freeze since a large region would be below the melting point. The existence of the inner core coupled with the apparently impossibility of homogeneously nucleating solid near Earth's centre led Huguet et al. (2018) to label the problem the 'inner core paradox'. Huguet et al. (2018) considered several limitations of CNT, but still found that homogeneous nucleation cannot explain inner core nucleation. They considered various catalysts for heterogeneous nucleation, all with caveats; the least implausible scenario was considered to be delivery of solid to Earth's centre from the mantle. The difficulty here is that the solid must be dense enough to sink to Earth's centre, which requires that it is composed of pure iron or iron alloyed with an element that is heavier than the main light elements in the core, thought to be silicon, sulphur and oxygen (Alfè et al., 2002b;Hirose et al., 2013). All viable alloys are thought to be highly soluble in liquid iron and so the solid would need to be large enough to avoid complete dissolution in either the lower mantle or core during its descent. Given the significant uncertainties with the existence and dynamics of this process it cannot be yet be considered a robust explanation for the inner core paradox.
In this paper we address three major uncertainties involved in the application of CNT to inner core nucleation using molecular dynamics simulations. First, we calculate the quantities needed to evaluate CNT for pure iron and iron alloys at core conditions. Huguet et al. (2018) were forced to use highly uncertain estimates for the prefactor I 0 and interfacial energy γ (defined precisely below) required by the CNT equations since relevant estimates for iron alloys were not available. We focus on iron-oxygen alloys since O does not enter the solid at core conditions (Alfè et al., 2002b) and is therefore the element primarily responsible for the seismically-observed density difference between the solid and liquid core and the 600-700 K depression in the core melting temperature compared to the melting temperature of pure iron (Davies, 2015).
Second, we test the basic validity of CNT at core conditions using two sets of simulations for each composition. In the first set we model freezing from a pure liquid and calculate the waiting time required to observe freezing events for a given undercooling without recourse to CNT. In the second set we model freezing from a seed, which is much less computationally demanding but requires CNT to relate the undercooling to the waiting time.
Third, we investigate the effect of external wave-like forcings applied to our simulations. Previous studies have found that such forcings can significantly reduce the nucleation barrier (see Huguet et al., 2018, for a discussion). A wide variety of waves exist in the liquid core, ranging from slow waves that arise from magnetic and rotational effects to seismic and electromagnetic waves (Moffatt, 1978). We consider an arbitrary forcing with prescribed amplitude that is not directly related to a particular wave-like motion in the core. The goal at this stage is simply to understand whether a wave-like forcing can reduce the undercooling for core materials at core conditions. This paper is organised as follows. In Section 2 we present relevant results from CNT and describe the atomic potentials used in our calculations of the Fe and Fe-O liquid systems. In Section 3 we present results for nucleation from a pure liquid and nucleation in the presence of a solid seed in both Fe and Fe-O systems and compare the results to CNT. We discuss the application of these results to Earth in Section 4 and provide conclusions in Section 5.

Classical nucleation theory
At equilibrium a thermodynamic system is found in the state of minimum Gibbs free energy g. In particular, equilibrium between a solid and a liquid is obtained at the melting temperature T m , where g sl = g s − g l = 0, with g s and g l the free energies of solid and liquid, respectively. If a liquid is cooled below T m it will eventually solidify, but homogeneous solid nucleation is an activated process, because the formation of a solid embryo inside the liquid introduces a solid-liquid interface, with an associated positive free interfacial energy γ . The process is usually described by classical nucleation theory (CNT) (Christian, 2002). If the embryo is spherical with radius r then the Gibbs free energy change as it forms is: G = 4 3 πr 3 g sl + 4π r 2 γ . (1) The number of embryos I(r) formed in the liquid per unit volume per unit time is proportional to exp (− G/k B T ), where k B is the Boltzmann constant and T is the temperature.
Above T m both terms in equation (1) are positive and so G increases with r. The distribution of embryos is therefore confined to small radii, as the probability of forming an embryo of radius r is a decreasing function of r. Below T m the first term in equation (1) is negative, and therefore G increases at first, reaches a maximum, and then decreases monotonically. If we define r c to be the critical embryo radius corresponding to the maximum of G, then I(r c ) is minimum and an embryo of that size will be unstable, having equal probability of either remelting or growing. If the embryo starts to grow beyond r c , the probability of growing further increases exponentially, and the whole system solidifies.
The condition to determine r c is obtained by setting d G/dr = 0, which gives r c = −2γ /g sl . Near T m we can write: where h f is the enthalpy of fusion, δT = T − T m is the undercooling and h c (δT ) is a correction that accounts for the non-linear behaviour of the solid-liquid free energy difference as δT increases.
In standard CNT h c = 1, but we found that for δT ∼ 1000 K h c can be as large as 1.06 and is therefore not completely negligible. By analysing our free energy data for pure iron (Alfè et al., 2002c) we found that a correction of the type h c (δT ) = 1 − 7.046 × 10 −5 × δT reduces the discrepancy between the left-and right-hand sides of equation (2) to less than 1% up to δT = −2000 K. Note that the correction is positive, meaning that the same value of G is obtained for somewhat lower δT .
The critical radius can thus be re-written as and is therefore a function of the interfacial energy γ and the critical undercooling temperature δT c : the closer the temperature to T m the larger the critical radius. By substituting equations (2) and (3) into equation (1) we obtain the value of the height of the free energy barrier to nucleation: The number of embryos of critical size r c per unit volume per unit time is where I 0 is a prefactor which is related to the number of atoms per unit volume and an attempt frequency of bringing atoms together to form a solid embryo.
Since they are at the top of the free energy barrier, on average half of these embryos will re-melt, and the other half will grow and lead to a freezing event. 1 The average waiting time to observe a freezing event, τ v (units of s m 3 ), is therefore This discussion suggests two possible ways to study homogeneous freezing via molecular dynamics simulation: freezing from a pure liquid and freezing with a seed. In each case we consider both pure iron and iron-oxygen alloys. In the case of freezing a pure liquid, described in Section 3.1, we simulate the liquid at different temperatures below the melting temperature until it freezes and measure the time for these freezing events to occur. This method provides a direct test of CNT using equation (6). In the case of freezing with a seed, described in Section 3.2, we fix the seed size and vary temperature until the seed melts. This method does not directly test CNT via equation (6); instead it tests the linear relation given by equation (3), which can be used to infer γ . In the following two subsections we describe the technical details of the pure iron and the iron oxygen potentials, including the prediction of their melting properties.

Pure iron
In the pure iron system we use an embedded atom model (EAM) (Daw and Baskes, 1983;Sutton and Chen, 1990) developed in Alfè et al. (2002a), in which the total energy has the form where the atomic energies E i consist of two parts: 1 In reality the probability of re-melting or growing into a larger solid are not exactly equal, as the free energy profile in equation (1) is not symmetric with respect to r c , but deviations will only marginally affect the pre-factor and can be ignored.
Here E rep i is a purely repulsive function of the interatomic distances r ij and F (ρ i ) is an embedding term accounting for the metallic bonding, with the densities ρ i = j =i (a/r ij ) m . The values of the parameters are: n = 5.93, = 0.1662 eV, C = 16.55, a = 3.4714 Å and m = 4.788 (Alfè et al., 2002a). We performed simulations with two system sizes, including 7776 and 40960 atoms. The volume was 7.04285 Å 3 /atom, giving pressures p close to 323 GPa at the temperatures of interest. We prepared the systems by melting and equilibrating them close to the melting temperature. These initial configurations were then given random velocities drawn from a Maxwellian distribution. All simulations were performed in the microcanonical ensemble (constant number of atoms N, constant volume V and constant energy E), using the Verlet algorithm with a time step of 1 fs. We extended each simulation by up to 2 ns, which resulted in a maximum accumulation of the drift in the constant of motion of no more than 10 K.

Iron-oxygen mixtures
To simulate iron oxygen mixtures we have constructed a more general EAM potential of the form where i Fe runs over the iron atoms, i O over the oxygen atoms and i FeO over the whole list of atoms. The atomic energies are and with densities The bonding term F O is not essential, but we keep it for generality.
To obtain the values of the parameters we used a 30 ps trajectory in a system containing 127 iron and 30 oxygen atoms at a temperature of 6350 K and a pressure of 240 GPa using density functional theory (DFT). The simulation was performed with the vasp code (Kresse and Furthmüller, 1996) using the technical parameters in Davies et al. (2018). The accuracy of the DFT description of iron at Earth's core conditions has been demonstrated in previous work, which showed good predictive power for its high pressure static, dynamic and thermodynamic properties (Alfè et al., 2000(Alfè et al., , 2001(Alfè et al., , 2002aGillan et al., 2006). We extracted N l = 6000 configurations (one every 5 fs) and calculated DFT energies U (l) and pressures p(l) for each configuration l. With these we constructed the quantities δU (l) = E tot (l) − U (l) and δ p(l) = p EAM (l) − p(l), where p EAM (l) is the pressure computed with the EAM model. The parameters of the EAM were then obtained by minimising the quantities . 1. Example of solid-liquid coexistence in an iron-oxygen mixture at a pressure of 323 GPa. Initially, the cell is built with solid and liquid iron in coexistence at the pure iron melting temperature of 6215 K. The simulation is then stopped, 1000 iron atoms (out of a total of 40,000) in the liquid are transformed to oxygens (blue dots), and the simulation is resumed with new random velocities (top panel). In this particular example the velocities are drawn from a Maxwellian distribution with a temperature of 5000 K. Since the total energy of the system has been reduced, the pure iron part of the cell starts to freeze (middle), and the release of latent heat increases the temperature. The temperature is above the melting temperature of the mixture, and this starts to melt (left side of the cell), while the pure iron side of the cell keeps freezing until the temperature reaches the melting temperature of the mixture (bottom). The snapshots are taken at t = 0, 0.1, 0.3 ns, respectively. (For interpretation of the colours in the figure, the reader is referred to the web version of this article.) and is about 50% larger than the value found for the EAM for pure iron, but still small.

Iron-oxygen melting point depression
To study the waiting time as a function of undercooling in the Fe-O system we first computed T m as a function of oxygen concentration. We used the coexistence approach, in which solid and liquid are simulated side by side in a large simulation cell containing 40,000 atoms in the NVE ensemble. Starting with the pure iron system to establish solid-liquid coexistence, we then transmuted 1000 iron atoms into oxygen in a portion of the liquid neighbouring the solid on one side (see Fig. 1(a)) and removed some kinetic energy from the system by restarting the simulation with lower velocities. This process causes the temperature to fall below the melting temperature of the pure system and the iron solidliquid interface to move into the liquid as the solid phase grows ( Fig. 1(b)), while oxygen also starts to diffuse in the liquid. Release of latent heat causes the system to heat up: solid continues to grow on the pure iron side of the simulation, which remains below its melting temperature, while solid iron neighbouring the liquid-oxygen mixture starts to melt due to the reduced melting temperature (Fig. 1(c)). Once oxygen is well mixed throughout the liquid and the temperature has settled to the melting temperature of the mixture the system reaches equilibrium, with the amount of solid and liquid remaining roughly constant. Note that almost no oxygen penetrates the solid, in agreement with the predictions of Alfè et al. (2002c).
Changing the amount of energy given to the system alters the fraction of solid and liquid, and therefore the concentration of oxygen in the liquid, c O . Simulations were repeated with initial temperatures of 6200, 5000, 4000, 3000, 2600 and 2200 K, for which the values of c O were 3.1%, 4.8%, 7.4%, 12.5%, 14.3% and 19.6% respectively. The resulting melting temperatures are plotted in Fig. 2, corrected by adding a term (p −323) ×dT m /dp K to account for the changing pressure, where, dT m /dp = 9 K GPa −1 is the slope of the melting curve at these conditions (Alfè et al., 2002d). The melting temperature of the mixture with respect to the melting temperature of pure iron, 6215 K at 323 GPa, is described well by the linear approximation T m = 6215 × (1 − c O /S ls ) (equation (12) in Alfè et al., 2002c) up to liquid oxygen concentrations c O 0.075, where S ls = 0.82 k B /atom is the entropy of melting of pure iron.

Freezing the pure liquid
In this approach we simulate the liquid at different temperatures T below T m until it freezes and measure the time for these freezing events to occur. We can therefore directly relate the rate of freezing to the undercooling without requiring CNT. We assume that a freezing event is a stochastic process with a distribution of waiting times τ w . At a given T we typically performed a set of M = 128 simulations, though for the highest temperatures M = 2300 was used. Within each set all simulations had the same kinetic energy, but statistically independent initial velocities. An example of 4 simulations using the 7776-atom system is shown in Fig. 3. After a short transient (not visible in the graph), the temperature settled around an initial value T 4190 K. All simulations froze at different times and settled around different finite temperature values because of different distributions of defects and stacking faults present in the solids, which change the average potential energy of the system. The effect of these defects and stacking faults also affects the pressure, though to a lesser extent.
Parameter values at the reference pressure of 323 GPa are presented in Table 1. To define the undercooling temperature δT c = T − T m requires T m at the pressure of the simulation, which can vary by a few GPa, so we always correct the melting temperature  for each simulation using T m = 6215 +(p −323) ×dT m /dp K. When the system freezes the pressure drops, and the release of latent heat causes a temperature increase as can be seen clearly in Fig. 3. The rapid temperature increase allows us to accurately identify the onset of freezing in each simulation, from which we obtain realisations, τ j , of τ w . The waiting time multiplied by volume defined in equation (6) can be obtained as τ v = τ w × V , where V is the volume of the simulation cell. Similarly to the case of homogeneous melting (Alfè et al., 2011), the distribution of waiting times is well approximated by an exponential form 1/τ 0 exp(−τ w /τ 0 ) with average waiting time τ 0 (Fig. 3). This form is expected for a random process with a constant probability per unit time 1/τ 0 of occurring, given that it has not already occurred.
For each undercooling δT c we obtained the average waiting time as follows: where τ j is either the time taken by simulation number j to freeze, or the total time of the simulation if no freezing event is observed, and N frozen is the number of simulations that ended up freezing. The error bar is σ 0 = τ 0 / √ N frozen . The approach described above was applied to pure Fe and Fe-O liquid. In both cases high (negative) values of δT c yielded short waiting times and all simulations froze in the total time allowed.
As δT c is reduced the waiting time grows very quickly, and some simulations ran their whole course without freezing. At the lowest δT c = −1720 K we observed only two freezing events in 2300 separate simulations, for a total simulation time of 1.8 μs, on the 40960-atom system. Fig. 4 shows waiting times as a function of undercooling for pure Fe and Fe-O systems (data are provided in Supplementary  Table 1). For pure Fe the systems with 40960 and 7776 atoms show only minor differences and so to study Fe-O we only used the 40960-atom system. It is clear that the behaviour of the Fe and Fe-O systems is very similar. Fig. 4 also shows the fit to equation (6) with γ and I 0 treated as free parameters. The fitted parameters and values of T m and h f used in the fits are listed in Table 1. The values of I 0 are somewhat lower than typical values of 10 −42 s −1 m −3 (Christian, 2002), while the values of γ are slightly lower than those obtained by Zhang et al. (2015). Nevertheless, CNT provides a reasonable fit to the data.

Freezing with a seed
A seed of radius r c has roughly 50% probability of growing into a solid at the corresponding critical undercooling δT c defined by equation (3). A series of simulations performed at different temperatures can therefore be used to identify δT c . Computationally, this approach is much less demanding than the one described in Section 3.1 because the seed will either melt or grow into a full solid on a time scale of O (10) ps. However, it does not yield a direct relation between nucleation rate and undercooling temperature and therefore does not provide results that can be extrapolated to core conditions. Instead, this method provides a test of equation (3), which predicts a linear relationship between r c and δT −1 c . If this relationship is borne out by the data then equation (3) can be used to compute γ .
We performed simulations using cells containing 40960 atoms, and seeds with radii equal to 10 (575 atoms), 15 (1965 atoms) and 20 Å (4679 atoms). At each temperature we performed at least 10 statistically independent simulations, and obtained the fraction f that ended up freezing. Since it represents uncorrelated events, f is a random variable following a Poisson distribution, which has mean value 0.5 at δT c , and approaches zero and one above and below δT c , respectively.
For the pure iron system f changes from 0 to 1 roughly linearly (Fig. 5a), which yields the results for δT c shown in Fig. 5b (data are provided in Supplementary Table 2). While our data are not obviously in conflict with the linear relationship predicted by equation (3), the deviation of the fit from the data is quite large; the standard deviation based on the sum of squares of errors is 2.2 Å, or 10-20% of the measurement. The small number of points precludes a robust error estimate, but also means that this test does not provide compelling support for equation (3).
If we assume that equation (3) holds then we obtain γ = 1.13-1.26 J m −2 , for the three size embryos. There is some variability, which arises because the product r c δT c is not a monotonic function, as can be seen from Fig. 5b; again, this behaviour likely reflects the low number of data points. Taken together, the range of values is consistent with the value of 1.2 J m −2 reported by Zhang et al. (2015) and used by Huguet et al. (2018).
To study freezing of Fe-O we used the 40960-atom system containing 36960 iron atoms and 4000 oxygen atoms (concentration near 10%) randomly distributed and inserted a pure iron solid seed of radius r c = 10 Å. The melting temperature of this system at 323 GPa is 5600 K (Fig. 2). We found that the system always melted for T > 4150 K and the seed grew into a solid for T < 4100 K, giving δT c −1475 K, which is very close to the value found for pure iron for a similar sized seed (Fig. 5b).

Freezing with an added perturbation
To explore the effect of perturbations on the nucleation process we repeated the investigations of Section 3.1 by adding an imposed sinusoidal force. We used the 7776-atom cell and as the simulations progressed we subjected each atom at position 0 ≤ (x, y, z) ≤ 1 to an additional planar force of the form A sin(2π z) with A = 1.0 eV/Å. These simulations were performed using the canonical ensemble (NVT) and a Nosé thermostat (Nosé, 1984) to fix the temperature. With these arbitrary values, chosen purely to illustrate the potential effect of a perturbation on the nucleation rate, the amplitude of the additional force is less than 10% of the typical forces acting on each atom due to thermal motion. A full analysis of the behaviour of the system as function of A is beyond the scope of this study, but we did perform a limited number of simulations with A = 0.75 and A = 0.5 eV/Å.  Table 3 show that the presence of the perturbation has a dramatic effect on the nucleation rate. For A = 1.0 similar nucleation rates are obtained at almost half the undercooling temperature compared to the unforced system. Increasing A by 0.25 eV/Å decreases δT c by around 300 K at similar nucleation rates, a significant effect.
The effect of the perturbation on the nucleation rate is shown clearly by considering the atomic density. The change in density across the box is large in all simulations with a perturbation, about 10% of the mean value for the example shown in Figs. 6b and c. Freezing occurs where the density is highest (in the middle of the box in Fig. 6c), which corresponds to a region of high pressure and hence high melting temperature. Since the temperature in the simulation is fixed, regions of high atomic density correspond to regions of anomalously high δT . Nucleation occurs preferentially in these regions because the energy barrier is lowest: the increase in T m that would enhance G is more than compensated by the increase in δT c which lowers G.
The calculated δT c for each simulation (as shown in the figures) is based on the average pressure (see Section 3.1) and therefore underestimates the actual δT c required for nucleation. The perturbation does not act to reduce the overall energy barrier to nucleation; instead it provides an additional mechanism for locally increasing the undercooling to the required level. This result may depend on the thermostat used in our simulations to control the temperature. Without the thermostat, perturbations of the kind considered here could generate an increase in temperature, which would reduce the predicted undercooling.

Discussion
Computational limitations mean that accessible waiting times are limited to values of τ v ∼ 10 −31 s m 3 . The value of τ v appro- priate for the inner core, assuming a waiting time of 1 billion years which is probably near the upper limit of plausible estimates, is approximately τ ic v ∼ V ic τ ic = 7.6 × 10 18 × 3.155 × 10 16 = 2.3 × 10 35 s m 3 . Evidently the simulated waiting times are over 60 orders of magnitude lower than appropriate for the inner core; with no hope of achieving these extreme values in the near future using our current approach, a means of extrapolating to core conditions is required. It is therefore important to consider whether classical nucleation theory (CNT) is suitable for this task.
Our results lend some support to the use of CNT for estimating nucleation barriers at core conditions. The calculations with a solid seed of radius r c in Section 3.2 are compatible with the CNT relation between r c and δT c , but the limited size of the dataset and significant (though also uncertain) error of 10-20% on the bestfitting line means that we cannot rule out substantial deviations from CNT. The pure liquid calculations in Section 3.1 also show that the relation between τ v and δT c predicted by CNT can be fit by our data, though our estimate of I 0 ∼ 10 −46 s −1 m −3 (Table 1) is quite different from the CNT prediction of I 0 ∼ 10 −42 s −1 m −3 .
The pure liquid calculations only confirm that the data can be fit by CNT; they do not require that the CNT relation G c ∼ (δT c ) −2 holds. Tests have revealed that empirical relations of the form G c ∼ δT c and G c ∼ (δT c ) −1 can fit our data just as well as the CNT relation. The relation G c ∼ δT c could be realised if the growing solid has a surface area to volume ratio that is a positive power of r, rather like a sponge. However, in this situation the energy barrier grows with the crystal size, which is unphysical. The relation G c ∼ (δT c ) −1 does not seem to have a physical interpretation and so it does not seem sensible to use it for extrapolation. We therefore recourse to CNT to extrapolate to core conditions, noting that the reservations above still stand. Fig. 7 shows extrapolations of our data to Earth's core conditions using CNT. In the absence of perturbations, nucleation of the solid requires undercooling of 730 ± 20 K for the pure iron system and δT c = 675 ± 35 K for the Fe-O mixture. Interestingly, δT c /T m ≈ 0.12 in both cases, demonstrating the similar behaviour between Fe-O alloy and pure metal systems. On the one hand, this is a significant reduction from the value of δT c ∼ 1000 K predicted by Huguet et al. (2018) since the extra 300 K requires an additional 3 billion years of core cooling based on a cooling rate of 100 K Gyr −1 (Davies, 2015). However, our values are still much too large to be achieved by core cooling alone since the inner core is thought to be less than 1 billion years old (Nimmo, 2015;Davies et al., 2015).
Extrapolating the simulations that were subjected to wave-like perturbations shows a drastic reduction in the undercooling, with δT c ≈280 K for the most extreme case (Fig. 7). However, in the simulations the density fluctuations induced by the perturbation are ∼10%, which is comparable to the density difference across the liquid outer core. The density fluctuations associated with convection are thought to be about 9 orders to magnitude smaller than the mean density (Stevenson, 1987), while even the density changes caused by seismic waves are insignificant compared to this value. The problem is that the perturbations we have considered do not significantly reduce the surface energy γ , which is the key factor in equation (1); instead, nucleation is induced locally in regions of anomalously high density. Perturbations of this magnitude seem unlikely to exist in Earth's core. Moreover, such perturbations would likely be accompanied by heating, an effect that is removed by the thermostat in our simulations, thus reducing the predicted undercooling. The extrapolations based on these simulations shown in Fig. 7 are therefore unrealistic and do not resolve the large required undercooling. Perturbations driving pressure changes of several GPa reduce the undercooling by ∼100 K.
As expected, our simulations show that the nucleation rate with a solid seed is significantly faster than from a pure liquid, effectively removing the barrier if the undercooling δT is higher than the critical undercooling δT c defined in Eq. (3). Extrapolating the results in Supplementary Table 2 shows that a seed of O (10 3 ) Å is needed to reduce the undercooling to a few tens of Kelvin. The problem remains to explain how such a seed could have survived inside the core. If this seed was introduced to the core from the mantle, the analysis of Huguet et al. (2018) shows that it would rapidly dissolve before reaching the centre of the planet.

Conclusions
We have investigated homogeneous nucleation of iron and iron alloys at core conditions using atomic-scale simulations. We have directly calculated the time to observe freezing events as a function of undercooling below the melting temperature in simulations of a pure liquid and established the relationship between the critical radius for nucleation and undercooling in a suite of simulations containing a solid seed. Our simulations are limited to conditions that are far-removed from those in Earth's core and so we have considered classical nucleation theory (CNT) as a means to bridge the gap. In both cases we find reasonable agreement with CNT, but also some important discrepancies. Future work should look to further test CNT at lower undercooling temperatures than we have been able to achieve here.
Assuming CNT we find undercooling temperatures of δT c = 730 ± 20 K for the pure iron system and δT c = 675 ± 35 K for the Fe-O system at inner core boundary conditions (Table 1). These values are about 300 K lower than those obtained by Huguet et al. (2018), but are still much too large to be achieved in the Earth's core, which cools at only ∼100 K Gyr −1 and so our calculations confirm the existence of the 'inner core paradox' as described by Huguet et al. (2018).
Of the ∼600-700 K required for homogeneous nucleation of the inner core, at most 200 K could plausibly come from the long-term cooling of the core, leaving another mechanism to facilitate the remainder. Our calculations have shown that wave-like perturbations represent another mechanism, independent from cooling, for providing undercooling. A plausible future direction for resolving the inner core paradox is to consider the role of different kinds of perturbations in the nucleation process. It may be that further reductions in the required homogeneous supercooling can be achieved via mechanisms we have not considered.