Terrestrial planets are differentiated into a silicate mantle and a metallic core. The core-forming liquid is required to segregate from the mantle over a planetary lengthscale. Many core formation scenarios have been proposed based on geophysical and geochemical constraints, as summarized by Rubie and Jacobson.1 In the stage of planetary accretion, a magma ocean forms at the surface because of the large amount of potential energy released by accretion. Segregation of metal and silicate liquids easily occurs in the magma ocean simply because of the difference between their densities. Beneath the magma ocean, a segregation mechanism of liquid metal from solid silicate plays a key role in core formation. If the interfacial energy between liquid metal and solid silicate is sufficiently small, the core-forming melt can percolate through the solid silicate mantle. Although the percolation mechanism is unlikely at least up to the top of the lower mantle,2,3 it has been reported that Fe–Ni–S liquid can percolate in silicate perovskite if the pressure is greater than 50 GPa.4,5 This suggests that percolative core formation may occur in the deep mantle. However, if the pressure at the base of the magma ocean is less than 50 GPa, segregated liquid metal accumulates at the top of the solid mantle. In this case, the large metal cumulates descend through the solid mantle because of gravitational [Rayleigh–Taylor (RT)] instability or negative diapirism,1,6,7 which is downward migration of large metal blobs (diapirs) through a lighter solid mantle.

The RT instability is a gravitational instability that occurs when a heavy fluid overlies a light fluid under gravity. The growth of the RT instability is generally classified into four stages.8 In the early stage (stage 1), perturbations, i.e., the onset of instability, occur at the interface of the heavy and light fluid layers. The amplitude l of a sinusoidal perturbation grows exponentially with time following the classical theory:9l=l0expγt,(1)where l₀, γ, and t are the initial amplitude, growth rate, and time, respectively. The growth rate γ can be expressed asγ=Akg,(2)where k, g, and A respectively denote the perturbation wavenumber (k = 2π/λ), gravitational acceleration, and Atwood number. The Atwood number is defined asA=ρH−ρLρH+ρL,(3)where ρ_H and ρ_L are the densities of the heavy and light fluids, respectively. The growth of a single mode of l follows the classical theory [Eq. (1)] up to l of the order of 0.1–0.4λ, where λ is the wavelength of the perturbation. Then, in stage 2, the perturbation amplitude begins to grow nonlinearly until l becomes approximately λ.10 The nonlinear growth phase is characterized by mode mixing of the perturbation and high-harmonic generation. The heavy fluid begins to form downwelling spikes (diapirs) between the upwelling diapirs of the light fluid. In stage 3, the interacting up- and downwelling diapirs develop to form stable shapes with a stable wavelength. The up- and downwelling diapirs then move in a steady state. In the final stage (stage 4), the spikes and diapirs are broken up, causing turbulent or chaotic mixing.

A core formation mechanism related to the RT instability has been reported based on the results of numerical simulations6,7,11–13 and experiments using analog materials.14 These simulations and experiments have mainly considered the evolution and segregation of a large metal droplet as a consequence of RT instability (stages 3 and 4). Several studies have considered RT instability as a possible scenario for core formation and have indicated that the segregation timescale and characteristic length vary considerably depending on the viscosity model adopted for the mantle. However, experiments using actual core-forming melts, i.e., liquid Fe alloys, have not been performed at high pressures to study the RT instability, and the onset of the RT instability and its initial behavior (stages 1 and 2) for these materials have not been well constrained. In the study reported here, we perform in situ observations of the RT instability of liquid Fe and Fe–Si alloys at high pressure using a high-power laser, and we measure the growth rate of the RT instability.

The compositions of the sample used in this study were Fe, Fe₈₂Si₁₈ (Fe–10 wt. % Si), and Fe₇₇Si₃₃ (Fe–20 wt. % Si). The Fe sample was a thin foil with a thickness of 6.2 ± 0.4 or 11.6 ± 0.6 µm (99.85% purity, Goodfellow). The Fe–Si samples were synthesized from a mixture of Fe (99.98% purity, Strem Chemicals, Inc.) and Si (99.9% purity, Wako Pure Chemical Industries, Ltd.) with an arc furnace in argon gas. The synthesized Fe–Si samples were sliced and polished to a thickness of 7.5–16 µm. A sine-wave-like perturbation was imposed on the sample surface as an initial perturbation of the RT instability using micro-laser machining (L.P.S. Works Co., Ltd.). With this imposition of an initial perturbation, the growth of perturbations with specific wavelengths can readily be detected. This technique has been used for observations of RT instability in inertial confinement fusion research using laser shocks.15–19 The shapes and sizes of the initial perturbations were measured precisely using a 3D laser scanning confocal microscope (VK-X Series, Keyence Corp.). A typical microscope image of an initial perturbation on a sample surface is shown in Fig. 1(a), the horizontal line profiles of the perturbation at the vertical center of the sample are shown in Fig. 1(b), and the amplitudes of each mode of the perturbation are shown in Fig. 1(c). The initial wavelength λ₀ and amplitude l₀ of the perturbation in the fundamental mode are 62–80 and 0.9–2.3 µm, respectively, and the details of the sample conditions are listed in Table I. The values of l₀ in Table I correspond to l₀ in the range 0.01λ₀ to 0.04λ₀.

Some of the foil samples were coated with forsterite (Mg₂SiO₄) to study the effect of the presence of a silicate layer on the RT instability. Forsterite was coated onto the metal samples using a radio-frequency sputtering machine (SVC-700RF I, Sanyu Electron Co., Ltd.). In the coating process, argon ions bombarded a charged polycrystalline forsterite target plate in the sample chamber, and the target forsterite was ejected and coated the substrates (sample foils). The coating was performed for 4 h, and the thickness of the coated layer was ∼0.6 µm. After forsterite coating, the wavelength and amplitude of the perturbation on the sample were measured using the laser microscope.

Laser-shock experiments were performed using the high-power GEKKO-XII (GXII) laser, at the Institute of Laser Engineering of Osaka University.20 The experimental setup is shown in Fig. 2. The drive laser directly irradiated the sample, causing an intense shock to propagate in the sample. As a result, the sample was accelerated in a strong effective gravitational field. The direction of this effective gravity is shown in Fig. 2. The total energy of drive laser was about 1 kJ. The drive laser was stacked by delaying three coherent laser beams of 527 nm wavelength to create a nearly flat top with a relatively long pulse duration [3.2 ns full width at half maximum (FWHM)]. The diameter of the laser beam was ∼250 µm. A gravity measurement using the same drive laser condition as this study was performed based on the sample trajectory from x-ray shadowgraphs obtained by Sakaiya et al.21 The typical gravitational acceleration in this study was 1.5 ± 0.6 × 10¹³ m/s², which corresponds to a pressure of ∼1000 ± 400 GPa at the target thickness of 10 µm. The error in the experimental pressure was derived from the error in the gravity estimation. After shock wave propagation, the sample could be approximated as being in a state of constant acceleration (i.e., g = constant).21 Since the constantly accelerated sample was then under constant pressure, the state of compression and the temperature of the sample could also be regarded as constant. Therefore, there was no density heterogeneity in the sample during this period. This period corresponded to the timing of in situ x-ray observation of the RT instability in this study (t = 0.8–2.5 ns). The timing of the RT instability is discussed in Sec. III.

To evaluate the validity of the pressure estimation and the timing of the RT instability from the experiments, we also performed a 1D hydrodynamic simulation (ILESTA-1D)22 for the sample conditions used in this study and compared the simulation results with the measured ones. The sample trajectory from the 1D simulation for a sample thickness of 10 µm was consistent with that from the gravity measurement. In the 1D simulation, the estimated pressure in the RT instability period (1–2 ns) was 800 GPa, which is close to the experimentally obtained pressure of 1000 ± 400 GPa. In terms of the sample state, it has been reported from ramp compression experiments (e.g., those by Remington et al.23) that melting of Fe occurs above 400 GPa. The pressure condition in this study of 1000 ± 400 GPa corresponds to the Fe alloy sample being in a liquid state. Furthermore, since the 1D simulation was based on a hydrodynamic code,22 the agreement between its results and those from experiments suggests that the sample behavior can be treated using a fluid dynamical approach.

Simultaneously with the drive laser, we used a backlight laser to irradiate the Ti foil to produce x-rays for radiographic measurements (Fig. 2). The backlight laser beam was generated by stacking two delaying lasers with a laser pulse of 2.2 ns FWHM, and its total energy was about 470 J. The behavior of the perturbations on the Fe alloys was observed based on temporal variations in the transmitted x-ray intensity from the samples.15,24 The drive laser and the x-rays emitted by the backlight laser came from the same direction with respect to the sample surface, as shown in Fig. 2 (face-on radiography). The intensity of the x-rays transmitted through the sample was measured using an x-ray streak camera (XSC) (Hamamatsu C7700-31). The spatial and temporal resolutions of the XSC were ∼14 µm and 130 ps, respectively. The energy of the backlight x-rays was about 4.7 keV. Since spatial profiles of backlight x-ray intensity at each time were obtained [as shown in Figs. 3(c) and 3(d)], the spatial and temporal intensity profiles of the backlight x-rays are taken into account in the image analysis of the x-ray shadowgraphs of the sample. Details of the image resolution and detection of perturbation amplitudes are given in the supplementary material.

Images obtained from face-on radiography are shown in Figs. 3(a) and 3(b). These images reveal the temporal variations of the sample surface perturbations. Horizontal (spatial) line-scan profiles from the x-ray radiography are shown in Figs. 3(c) and 3(d). The perturbations of the sample surfaces can be observed as cyclic variations in transmitted x-ray intensity. The fundamental modes of the perturbations in Figs. 3(c) and 3(d) indicate wavelengths of 62 µm for Fe and 80 µm for Fe–10 wt. % Si, which correspond to the initial perturbation wavelengths. Thus, the temporal variation in amplitude of this fundamental mode shows the behavior of the RT instability. We analyzed the temporal variation in perturbation amplitude using the image analysis procedure given by Sakaiya et al.21 The perturbation amplitude of the fundamental mode increases with time in the time range of 0.8–2.5 ns in Fe–10 wt. % Si, whereas that of the third mode remains almost constant.

To understand the variation in perturbation amplitude from the initial state, we consider the growth factor G. This is defined as the measured surface perturbation normalized by the initial value and can be expressed as follows:25G≡ρlρ0l0∼F(I/Ii)μMρ0l0,(4)where l, ρ, μ, and M are respectively the amplitude of the fundamental mode, the density, the mass absorption coefficient, and the modulation transfer function of the diagnostic system used. Subscript 0 indicates the initial state at t = 0, and F(I/I_i) denotes the Fourier-analyzed ratio I/I_i, where I_i and I are the intensities of incident (backlight) and transmitted x-rays, respectively. Since I/I_i is measured by x-ray radiography, G can be obtained using Eq. (4) together with the information about the sample initial state given in Table I. μ of the samples can be obtained using the reported μ of the end-member components (Fe and Si) at 4.75 keV.26M is estimated following the procedure of Azechi et al.24 The extraction procedure for M has also been reported by Smalyuk et al.27 The temporal variation of the growth factor G thus obtained is shown in Fig. 4. The temporal variation of the growth factor of the perturbation on a laser-shocked sample is generally divided into three regimes:21 lateral mass flow, RT instability, and decrease in the drive laser. The details of these regimes are explained in the supplementary material. G increases exponentially with time, corresponding to the regime associated with RT instability and suggesting that the perturbation growth follows the classical theory expressed by Eq. (1). From Eqs. (1) and (4), G can be written in terms of a growth rate γ asG=ρl0exp(γt)ρ0l0=ρρ0exp(γt).(5)G in the RT instability regime is then fitted with this exponential function to obtain the growth rate γ of the RT instability. The values obtained for γ are in the range 0.30–0.34 ns⁻¹ for the Fe–Si samples and are summarized in Table S1 in the supplementary material.

To check the timing of the RT instability regime, we performed 1D hydrodynamic simulations22 of samples with a density of 6.8 g/cm³ and thicknesses of 5 and 10 µm. In the simulation results, both the pressure and acceleration become almost constant after 1 ns up to at least 2 ns. Before 1 ns, the pressure and acceleration change significantly with time, which corresponds to a period of shock-wave propagation before sample acceleration. As described in the supplementary material, the stage of shock-wave propagation (the lateral mass flow regime) is characterized by propagation of a rippled shock. After the lateral mass flow regime, acceleration of the sample occurs, and the RT instability regime starts. Thus, the period of 1–2 ns corresponds to the RT instability regime in the simulation, which is consistent with the observed RT instability regime (0.8–2.5 ns).

Temporal variations in the growth factor G of Fe, Fe–10 wt. % Si, and Fe–20 wt. % Si are shown in Figs. 4(a) and 4(b). The RT instability regime corresponds to the period shown by the filled symbols. The variations in G earlier than the RT regime are caused by the lateral mass flow associated with rippled shock propagation, as discussed in Sec. III A and by Sakaiya et al.21 After the RT regime, G saturates and gradually decreases as a result of reduction in drive laser intensity. In the RT instability regime, the G values of both Fe–10 wt. % Si and Fe–20 wt. % Si clearly increase exponentially with time [Fig. 4(a)]. This suggests that the RT instability of the Fe–Si alloys grows in the observed time ranges following classical theory [Eq. (1)]. By contrast, the G values of Fe do not show a clear trend with time. The G values of the two Fe samples starting with different initial perturbation wavelengths (62 and 80 µm) show similar behavior (no clear trend with time), as shown in Fig. 4(b). With regard to the sample thickness, we measured both thin (6.2 µm) and thick (11.6 µm) Fe samples, and the time window of 1–2 ns overlaps the RT instability regime estimated from the 1D simulation for both sample thicknesses. The absence of any growth of RT instability for both sample thicknesses [Fig. 4(b)] implies that the growth rate for Fe is too small to be detected in the measured time range. These results indicate that dissolution of Si into Fe increases the growth rate γ, i.e., leads to fast growth of RT instability. The effect of Si content on γ is shown in Fig. 5. Although the error is relatively large, γ tends to increase with Si content. This increase in γ might be caused by the reduction in sample density due to the addition of Si. The density of the sample decreases with increasing Si content, and this density reduction may cause a difference in the gravitational effect. However, since the gravitational acceleration g is given by g = P/ρd (where P and d are respectively pressure and sample thickness), g will increase by 10% when ρ is reduced by 10% through addition of Si at the same pressure. However, this increase in g leads to an increase in γ that is less than the error in γ, and hence the extent of the increase in γ cannot be explained solely by the difference in density. A difference in diffusion rate between Fe and Si might also influence the behavior of γ. However, atomic mixing has been reported to occur in the later stage of nonlinear growth of RT instability28 and so the difference in diffusion rate is not relevant in the present case (the early stage of RT instability). Differences in other transport properties, such as viscosity, between Fe and Fe–Si might also influence γ. Furthermore, the effect of interfacial tension of the Fe alloy on RT behavior has also been considered using the dimensionless Bond number B, which is the ratio between buoyancy force and interfacial tension. Calculations of B for the present samples and experimental conditions (details of which are given in the supplementary material) gives values that are large enough to suggest that the gravitational force is dominant and the interfacial tension is negligibly small.

The effect of the forsterite layer on the temporal variations of the growth factor G of Fe–20 wt. % Si is shown in Fig. 4(c). It can be seen that the values of G of Fe–20 wt. % Si with the forsterite layer do not show any obvious trend with time, which is different from the behavior of the sample without the forsterite layer (γ = 0.34). Therefore, the existence of a forsterite layer suppresses the growth of RT instability. The RT instability of the Fe sample does not grow in the observed time range [Fig. 4(b)]. The Fe + forsterite sample also does not show any growth of the RT instability [Fig. 4(d)].

Since the Atwood number A, representing the density difference between the two phases, decreases from 1 for the sample without the forsterite layer to 0.35–0.41 for that with such a layer [see Eq. (3)], the RT instability of the sample with the forsterite layer grows more slowly than that of the sample without the layer. This suggests that a longer measurement time is necessary to observe the perturbation growth of the sample with the forsterite layer. Alternatively, the Atwood number of the laser-shocked sample might correspond to the density difference between the Fe alloy and the ablation plasma. If the density varies smoothly between Fe-alloy and ablation plasma and minimum density gradient scale L is finite, then the expression for γ is different depending on the wavelength magnitude compared with L (see, e.g., Betti et al.29). If the product of the perturbation wavenumber k and L is much smaller than 1 (kL ≪ 1, i.e., the long-wavelength mode), then this mode is not affected by the finite L and grows according to the classical growth rate [Eq. (2)]. However, if kL ≫ 1 (the short-wavelength mode), then this mode is localized inside the smooth interface and grows at a rateγ=Akg1+AkL.(6)Thus, if kL ≫ 1 and ablation of the forsterite layer increases L compared with the sample without such a layer, then γ could decrease depending on the effect of L. We have estimated the value of kL using λ = 62–80 µm [the initial wavelength (= 2π/k) in this study] and L = 1.5 µm, which is the density gradient scale reported by Betti et al.29 The calculated values of kL are in the range 0.12–0.15 and thus approximately one order of magnitude smaller than 1. Thus, L can be considered to have no influence in this study unless it is an order of magnitude larger than 1.5 µm. This suggests that the growth rate γ in this study can be approximated by Eq. (2). If the values of g and k in Eq. (2) are the same between the samples with and without a forsterite layer, then the difference in γ between these samples can be explained only by the difference in the Atwood number A. Furthermore, in the case of a laser-shocked sample with an ablation stability effect, γ can be expressed as the classical growth rate [Eq. (2)] plus an ablation stability term βkV_a:γ=Akg−βkVa,(7)where β and V_a are a constant and the ablation velocity, respectively (see, e.g., Betti et al.29). According to Eq. (7), even if the ablation stability term is taken into account, γ can vary only by a factor of order unity, not by an order of magnitude. Hence, we have not taken the ablation stability term for γ into account when estimating the core formation timescale.

The growth rate γ obtained in this study is applicable to estimation of the timescale of onset and initial growth of RT instability (stage 1) in planetary core formation. When l reaches 1/10 of λ, it saturates. The timescale t₁ of stage 1 corresponds to the grow-up time for realistic values of λ in the range 0.1–100 km (i.e., l in the range 0.01–10 km) at the bottom of the magma ocean in the Earth’s interior.7,13 After this, the growth becomes nonlinear10 (stage 2). In stage 2, the descending spike (downward perturbation) has been reported to reach a terminal velocity approximated by Stokes’s law8 or by a growth constant, or by the free-fall velocity in the case of large Atwood number (A > 0.6).28 Thus, the timescale t₂ of stage 2 corresponds to the falling time of the descending spike according to the terminal velocity or Newtonian free fall in the solid mantle (500–3000 km in depth) beneath the magma ocean. Based on the γ values obtained in this study with a simplified scaling (see the supplementary material), the estimated t₁ and t₂ are proportional to the respective characteristic lengths (l for t₁ and the descending length for t₂). Therefore, compared with the entire core formation timescale (30–40 × 10⁶ years),30 the initial growth of RT instability is likely to occur on quite a short timescale (more than two orders of magnitude faster than the entire formation timescale).

In situ observations of the RT instability of liquid Fe and Fe–Si alloys were performed using a high-power laser-shock technique combined with x-ray radiography. The perturbation amplitude increases exponentially with time for the Fe–Si samples. By contrast, the growth of perturbations on the Fe sample is too small to be detected in the measured time range. The obtained growth rate of RT instability ranges from 0.30 to 0.34 ns⁻¹ and increases on the addition of Si. It is also found that the growth of RT instability is suppressed by coating with a forsterite layer. Based on the present results with simplified scaling, the timescale of the initial evolution of RT instability (stage 1) in a planetary interior is much faster than the timescale of the later growth of RT instability in the nonlinear regime. The observations and results of this study could be used as constraints for future numerical or analog models of planetary core formation driven by RT instability.

Acknowledgment. The authors thank R. Hosogi, T. Ueda, and T. Fujikawa for their technical assistance and are also grateful to the staff at the GEKKO-XII facility for their technical support. The authors acknowledge four anonymous reviewers for their constructive comments. This study was performed under a joint research project of the Institute of Laser Engineering, Osaka University and was partly supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Science, and Sport and Technology (MEXT), Japan to H.T. (Grant Nos. 26610141, 26247089, 15H05828, and 20H02008).