Formation of super-Mercuries via giant impacts

ABSTRACT

During the final stage of planetary formation, different formation pathways of planetary embryos could significantly influence the observed variations in planetary densities. Of the approximately 5000 exoplanets identified to date, a notable subset exhibits core fractions reminiscent of Mercury, potentially a consequence of high-velocity giant impacts. In order to better understand the influence of such collisions on planetary formation and compositional evolution, we conducted an extensive set of smoothed particle hydrodynamics giant impact simulations between two-layered rocky bodies. These simulations spanned a broad range of impact velocities from 1 to 11 times the mutual escape velocity. We derived novel scaling laws that estimate the mass and core mass fraction of the largest post-collision remnants. Our findings indicate that the extent of core vaporization markedly influences mantle stripping efficiency at low impact angles. We delineate the distinct roles played by two mechanisms – kinetic momentum transfer and vaporization-induced ejection – in mantle stripping. Our research suggests that collisional outcomes for multilayered planets are more complex than those for undifferentiated planetesimal impacts. Thus, a single universal law may not encompass all collision processes. We found a significant decrease in the mantle stripping efficiency as the impact angle increases. To form a 5 M_⊕ super-Mercury at 45°, an impact velocity over 200 km s⁻¹ is required. This poses a challenge to the formation of super-Mercuries through a single giant impact, implying that their formation would favour either relatively low-angle single impacts or multiple collisions.

hydrodynamics, methods: numerical, planets and satellites: composition, planets and satellites: formation, planets and satellites: terrestrial planets

1 INTRODUCTION

Several of rocky exoplanets discovered to date, including Kepler-107 c (Bonomo et al. 2019), GJ 367 b (Lam et al. 2021), Kepler-406 b (Marcy et al. 2014), K2-38 b (Toledo-Padrón et al. 2020), and K2-229 b (Santerne et al. 2018), have very high bulk densities, suggesting that they are iron-rich and may have a composition similar to that of Mercury. These planets are often referred to as ‘super-Mercuries’ due to their measured mass being at least one order of magnitude larger than Mercury’s. The formation mechanisms of Mercury and super-Mercuries are still debated and could be the result of mantle photoevaporation (Cameron 1985; Fegley & Cameron 1987), the separation of silicate and iron particles in the initial inner disc (Weidenschilling 1978; Wurm, Trieloff & Rauer 2013; Johansen & Dorn 2022), or giant impacts (Benz, Slattery & Cameron 1988; Benz et al. 2007; Chau et al. 2018; Bonomo et al. 2019; Reinhardt et al. 2022). According to Ito & Ikoma (2021), photoevaporation of mantle material is not significant enough to greatly affect the bulk composition of a 1 M_⊕ rocky planet. Although photophoretic forced separation could explain the formation of Mercury, it cannot accurately model Kepler-107 c, which is more than twice as dense (about 12.6 g cm⁻³) as the innermost Kepler-107 b (about 5.3 g cm⁻³). Instead, the dissimilar densities are consistent with a giant impact event that resulted in striping the silicate mantle off of the precursor planet resulting in a dense Kepler-107 c (Bonomo et al. 2019).

Giant impacts (violent collisions between planet-sized bodies) are extremely energetic events that could strip away a large fraction of mantle and generate a post-collision body with a high core fraction. Benz et al. (1988, 2007), Asphaug & Reufer (2014), and Chau et al. (2018) previously examined several impact conditions that could have contributed to the formation of Mercury. Mercury’s current core fraction could have been achieved by a single giant impact, in which the proto-Mercury was hit by a smaller projectile (Benz et al. 2007), by a massive one like Venus or Earth after which Mercury was the left remnant from the impactor (Asphaug & Reufer 2014), or by a sequence of less energetic giant impacts (Chau et al. 2018). The newly discovered super-Mercury exoplanets exhibit a wide range of mass and core fractions, and little is known about their formation history. There is currently no reliable way to link an observed planet with the specific conditions of the impact that led to its formation.

Scaling laws provide a means to predict collision outcomes based on the impact parameters, and can thus reveal which kind of collisions can produce high-density planets. Stewart & Leinhardt (2009) derived a universal law for the post-collision mass of the largest remnant for planetesimal–planetesimal collisions up to radii of 50 km. Based on about 60 smoothed particle hydrodynamics (SPH) simulations of giant impacts with targets ranging from 1 to 10 M_⊕, Marcus et al. (2009) found that the universal law was also in agreement with head-on impacts in the super-Earth regime. Additionally, they derived a scaling law to predict the iron mass fraction of the post-collision remnant for head-on impacts. Following this work, Leinhardt & Stewart (2012) made one of the most comprehensive and widely used scaling laws that can predict the outcomes of planetesimal impacts with different impact angles and impactor-to-target mass ratios. However, the data for collisions between planet-sized bodies were adopted from several previous studies that had different target composition set-ups. Therefore, Leinhardt & Stewart (2012) did not provide scaling laws to predict the iron mass fraction of post-collision remnants. In contrast, Carter et al. (2018) simulated impacts between differentiated targets of sizes that fall between the size regimes explored by Marcus et al. (2009) and Leinhardt & Stewart (2012). Carter et al. (2018) showed that the masses of the largest remnant at both head-on and oblique impacts show excellent agreement with the updated universal law from Leinhardt & Stewart (2012). None the less, Carter et al. (2018) reported that the core mass fraction generally reaches a plateau at relatively high impact energies and deviates from the prediction laws proposed by Marcus et al. (2009). However, very few hit-and-run simulations were performed in Carter et al. (2018). Gabriel et al. (2020) also developed empirical relationships for the accretion and erosion of gravity-dominated bodies of various compositions during various impact conditions, but the largest target mass was around 1 M_⊕. Studies by Marcus et al. (2009), Leinhardt & Stewart (2012), Carter et al. (2018), and Gabriel et al. (2020) have mainly focused on low-velocity collisions and lack simulation data for extreme impact conditions where super-Mercuries could be formed. Additionally, all simulations except for those by Marcus et al. (2009) are designed to study smaller bodies than Earth, and focus more on the collision outcome of planetesimals.

As a replacement for the scaling laws proposed by Marcus et al. (2009), Reinhardt et al. (2022) proposed new scaling laws based on super-Earth head-on impact simulations with the new M-ANEOS equation of state (Stewart et al. 2019; Stewart 2020) at high enough impact velocities to form very dense planets. They showed that in the critical disruption regime, above a certain normalized impact energy, the mass of the largest remnant will rapidly decrease and the iron mass fraction will increase rapidly, both following a power law. These new scaling laws deviate significantly from those proposed by Marcus et al. (2009) and Carter et al. (2018) at high impact energies.

In addition to the scaling laws derived from planetesimal impacts and two-layered super-Earth impacts, Denman et al. (2020) and Denman, Leinhardt & Carter (2022) conducted a comprehensive study on the loss of atmosphere during both head-on and oblique impacts between super-Earth-sized planets with thick atmospheres. They derived scaling laws to predict the total mass and remaining atmosphere mass of the largest remnant. Although the target planets used in their study had three layers and massive atmospheres, at high impact energies, most of the atmosphere is expelled, and the scaling law provides a useful reference for super-Earth impacts without atmospheres.

In this paper, we aim to bridge the gap between previous studies and provide a more systematic survey of giant impacts that could potentially form super-Mercuries. We present the results of over 1200 SPH giant impact simulations, with target masses ranging from 0.06 to 20 M_⊕ and impact parameters from 0 (head-on) to 0.7 (45°). We propose new scaling laws for both head-on and oblique impacts that can be used to constrain the collision conditions of super-Mercuries.

This paper is organized as follows. In Section 2, we describe the methods used to create planetary bodies and their respective giant impact SPH simulations. We explain the tools and methods used to analyse the results of simulations. In Section 3, we present the results of head-on and oblique impacts and the derived disruption criteria. In Section 4, we discuss the applicability of our scaling laws to predict the outcome of giant impacts and the two important mechanisms involved in the mantle stripping process. In Section 5, we summarize the results presented in the paper.

2 METHODS

We ran around 1200 SPH simulations with various target masses, impact velocities, and impact angles using swift (v 0.9.0) (Schaller et al. 2016, 2023; Kegerreis et al. 2019). We focus only on equal-mass giant impacts as these are the most efficient collision type to form dense planets in terms of total impact energy for both head-on impacts and oblique impacts (Leinhardt & Stewart 2012; Denman et al. 2020). Thus, we set the impactor-to-target mass ratios (γ = M_imp/M_targ) to 1 for all simulations in this work. Fig. 1 shows the parameter space spanned in our simulations. We simulated targets spanning the range of rocky planet masses: 0.06, 0.12, 0.35, 0.6, 1.0, 2.0, 3.0, 5.0, 5.8, 7.0, 8.5, 11.9, 15.0, and |$20.0\, \mathrm{M_{\oplus }}$|⁠, mostly consisting of 2–5 × 10⁵ SPH particles.

Figure 1.

The parameter space of impact simulations in this work. The colour of each point represents the maximum impact velocity normalized by the mutual escape velocity of a specific target mass and impact parameter.

2.1 Initial conditions

The planetary profiles and initial conditions were generated with woma (Kegerreis et al. 2019; Ruiz-Bonilla et al. 2020) initializing each layer with an isentropic temperature profile. The initial planets were differentiated with 30 per cent iron core and 70 per cent forsterite mantle. We used the iron (Stewart 2020) and the forsterite (Stewart et al. 2020) M-ANEOS equations of state (EoS) and regenerated the EoS tables with higher maximum densities (100 |$\mathrm{g\, cm}^{-3}$| for forsterite and 200 |$\mathrm{g\, cm}^{-3}$| for iron) using Stewart et al. (2019) and Stewart (2020). We chose entropies of the mantle to ensure that the mantles are hot but entirely below the melt curve, which varied between 2800 and 3027 J K⁻¹ kg⁻¹. The core entropies were set to give iron close to the solidus with temperatures larger than mantle temperature at the core–mantle boundary, which were between 1700 and 1800 J K⁻¹ kg⁻¹.

Before the impact simulations, each body was equilibrated in a cooling simulation for 20 h of simulation time in isolation to reach a stable hydrostatic state. During the first 10 h, the entropy of the core and mantle were fixed to the desired values at each time-step (see Carter et al. 2018) in order to produce planets with isentropic layers. Over the next 10 h, planets evolved without this additional damping, towards a hydrostatic profile. After the equilibration process, particles have a root-mean-squared velocity that is less than 1 per cent of the planet’s escape velocity.

We ran simulations for a range of impact parameters and target masses; additional head-on impact simulations were run for target masses of 2, 5, 5.8, 8.5, and 15 M_⊕. For oblique impact simulations, the impact parameter (b = sin θ, where θ is the angle between the centres of the bodies and the velocity vector at the time of contact; see Leinhardt & Stewart 2012) ranges from 0.1 to 0.7. The impact velocities range from 1 to 11 times the mutual escape velocity (V_esc), which is defined by

$$\begin{eqnarray} V_\mathrm{esc}=\sqrt{\frac{2G(M_\mathrm{targ}+M_\mathrm{imp})}{R_\mathrm{targ}+R_\mathrm{imp}}}, \end{eqnarray}$$

(1)

where G is the gravitational constant, and M_targ, R_targ and M_imp, R_imp are the mass and radius of the target and the impactor, respectively. For equal-mass oblique impacts with larger impact parameters, the required impact speed to generate a relatively dense planet (e.g. 50 per cent iron mass fraction) could be more than seven times the mutual escape velocity of the two-body system. In this work, higher impact velocities were chosen for larger impact parameters. Head-on impacts have impact velocities up to 3V_esc, while impacts with b = 0.7 can have impact velocities as high as 11V_esc.

Bodies were initially separated by a distance such that contact occurs 1 h after the start of the simulation to allow tidal deformation of the SPH planets (Kegerreis et al. 2020). For head-on impact simulations, the simulation time was approximately 15–20 h. We observed that the properties of post-collision remnants became relatively stable after 6–10 h of simulation time at the particle resolution we used. For oblique collisions, the simulations typically ran for 30 h to distinguish between accretion impacts and erosive hit-and-run impacts. The simulations were run in cubic boxes with sides ranging from 1000 to 5000 R_⊕, depending on impact velocities. Any particles leaving the box were removed from the simulation.

swift uses a parameter h_max (see swift documentation¹) setting the maximal allowed smoothing length of SPH particles. Particles in the simulations can only have smoothing length less than h_max, which will in turn set a density floor (minimum density of SPH particles). A maximum smoothing length or a minimum SPH density has a great influence on the post-collision disc (Hull et al. 2024) as relatively few particles are left in the disc region. Most of the simulations in this study (N = 2 × 10⁵) used an h_max = 0.2 R_⊕ resulting in a density floor of ∼0.1 g cm⁻³. To understand the influence of maximum smoothing length on post-collision remnants, we tested a subset of simulations with a larger h_max ranging from 0.5 to 5 R_⊕. The mass and iron fraction of the largest post-collision remnant agrees within 5 per cent for all the investigated h_max values. The larger the h_max the longer a simulation will take to run. Since this study focuses on the large-scale properties of post-collision remnants (rather than the properties of the disc), our smaller value for h_max provides an acceptable balance between computing cost and accuracy.

Finally, we also investigated the dependence of numerical resolution on the mass and iron mass fractions of the largest post-collision remnants. The resolution tests focused on the lower and higher mass range as these are most susceptible to phase changes and are thus the most sensitive to numerical resolution. We tested target masses below 1 M_⊕ and above 15 M_⊕ with higher resolution (10⁶ particles) and did not find significant resolution dependence for the properties of post-collision remnants, which is consistent with Meier, Reinhardt & Stadel (2021).

2.2 Remnant properties

We used the same method (Dou 2023) as in Marcus et al. (2009) and Carter et al. (2018) to determine particles bound to a post-collision remnant. Particles’ potential and kinetic energies were first computed in relation to the seed particle that was closest to the potential minimum. Then, for the bound particles, the centre of mass position and velocity were computed and used as the seed for the next iteration. This process was then repeated for the remaining unbound particles until the calculation converged. We do not take account of the re-accretion of ejecta unbound to the largest post-collision body but bound to the star, as in the previous similar studies (Marcus et al. 2009; Leinhardt & Stewart 2012; Carter et al. 2018; Reinhardt et al. 2022). The mass of the resulting iron core was found from the mass of iron particles in the largest remnant. The iron mass fraction of a remnant is the ratio between its iron mass and remnant mass.

We calculate the radii of post-collision planets corresponding to the radius of a super-Earth, using interior models with an iron mass fraction and mass determined as described in Section 2.2. We assume that all bound particles (debris) will re-accumulate and that the planets will eventually become differentiated again after a long cooling process. Benz et al. (2007) studied the long-term evolution of debris after potential Mercury-forming impacts using analytical calculations and found that around 40 per cent of particles would be re-accreted by Mercury after several Myr. Therefore, our simple assumption of no re-accumulation of unbound debris should provide a lower limit for the radii of super-Earths. The evolution of ejected debris will vary for different impact scenarios and planetary systems, and the inclusion of the effect of re-accumulation is beyond the scope of this work. For each post-collision remnant mass and iron mass fraction, we use the planet interior model code magrathea (Huang, Rice & Steffen 2022) to calculate the final planet’s radius. We use magrathea’s default core and mantle equations of state as described in Huang et al. (2022) (note that these are different from the equations of state used for the SPH simulations). We use isentropic temperature gradients and set the surface temperature to 1000 K; we find that the final radii vary little with the surface temperature.

2.3 Catastrophic disruption criteria

Scaling laws are typically functions of the mass or iron mass fraction of the largest remnant and the normalized impact energies. When the impact energies from different target masses are normalized by the catastrophic disruption criteria, |$Q^{*}_\mathrm{RD}$| (the impact energy needed to disperse half of the system’s total mass), they can be effectively summarized using similar equations. We follow the same variable definition and annotation as in Leinhardt & Stewart (2012). The specific impact energy Q_R is defined by

$$\begin{eqnarray} Q_\mathrm{R}=0.5\mu \frac{V^2_{i}}{M_\mathrm{tot}}, \end{eqnarray}$$

(2)

where M_tot = M_targ + M_imp is the total mass of the system, μ = M_targM_imp/M_tot is the reduced mass, and V_i is the impact velocity. Therefore, M_lr/M_tot = 0.5 (where M_lr is the mass of the largest post-collision remnant) for a collision with specific reduced impact energy equal to |$Q^{*}_\mathrm{{RD}}$|⁠. For head-on impacts, to determine the critical specific impact energy for catastrophic disruption we use linear interpolation between the two data points that surround the specific impact energy where about half of the total colliding mass remains bound.

During head-on impacts, most materials from the target and impactor will be compressed and shocked, in the resulting collisions. At impact energies meeting the criteria for catastrophic disruption, half of the masses are lost due to the high pressure or shock propagation experienced during the impacts. However, in oblique impacts, especially equal-mass impacts, half of the system’s mass can be lost during hit-and-run impacts, where the impactor deposits little energy into the target. A significant amount of mass can be lost due to the misalignment between the target and the impactor. As a result, if traditional catastrophic disruption criteria are still used to describe oblique impact systems, the criteria may not accurately reflect the energy behaviour of such systems. Therefore, for oblique impacts, we define a new disruption criterion |$Q^{^{\prime }*}_\mathrm{TD}$|⁠, which is the specific impact energy required to disperse half of the target’s mass. For each group of target mass and impact parameter, the target catastrophic disruption criterion |$Q^{^{\prime }*}_\mathrm{TD}$| was estimated by fitting the equation

$$\begin{eqnarray} \frac{M_\mathrm{lr}}{M_\mathrm{targ}} = c_1Q_\mathrm{R}^{c_2}+c_3, \end{eqnarray}$$

(3)

with c₁, c₂, and c₃ as free parameters and using data points of M_lr/M_targ between 0.8 and 0.2, and evaluating Q_R when M_lr/M_targ is 0.5. Fig. A1 shows an example of fitting process when M_targ is 1 M_⊕.

2.4 Hit-and-run velocity

For equal-mass oblique impacts, especially at low impact parameters, we find that it is unlikely for a ‘clean’ hit-and-run impact (where target masses remain unchanged) to occur, i.e. the mass of the largest remnant is equal to the target mass (M_lr = M_targ and M_lr/M_tot = 0.5). As impact speeds increase from the mutual escape velocity, equal-mass oblique impacts undergo a sudden transition from accretion-dominated impacts to erosive hit-and-run or graze-and-merge. The transition region is quite sensitive to the initial impact scenario and makes it challenging to determine the exact impact speed at which a ‘clean’ hit-and-run occurs. We therefore check the collision outcome more closely for impacts located in this transition region by running the SPH simulations for a longer time, up to 200 h. We then ran simple N-body simulations using the centre of masses and velocities of bound remnants with a duration of up to 10 yr as an additional way to determine the outcome of an impact. All the N-body simulations in this work were run using the rebound  N-body code (Rein & Liu 2012) and integrated using the hybrid symplectic mercurius integrator (Rein et al. 2019). We take the median speed of the two data points between which the ratio of the largest remnant mass to total mass (M_lr/M_tot) drops from above 0.5 to below 0.5 as the hit-and-run velocity. We used a finer speed grid around the speed where the hit-and-run transition begins, and the difference between the speeds we use to calculate median speeds is no larger than 1 km s⁻¹.

3 RESULTS

In this study, we simulate equal-mass collisions between target planets ranging from 0.06 to 20 M_⊕, encompassing a range of impact angle from head-on to oblique impacts. Our focus is on high-energy erosive impacts, where more than half of the total material in the system is ejected for head-on impacts, and more than half of the target mass is ejected for oblique impacts. Overall, our simulations reveal that the mantle stripping process is primarily influenced by two mechanisms: vaporization-induced ejection and kinetic momentum transfer. Vaporization plays a crucial role in head-on and nearly head-on impacts, while kinetic transfer significantly affects oblique impacts at relatively large impact angles.

We introduce new scaling laws for head-on impacts that accurately capture the subtle changes in mass and iron mass fraction of the largest remnant, thereby demonstrating that collisions between multilayered planets should be modelled separately at different impact energy regimes to better represent the distinct shock experiences encountered by each layer. Regarding oblique impacts, we employ a new criterion for target catastrophic disruption to normalize impact energies and present new scaling laws that can predict the outcomes of oblique impacts at different angles. Additionally, we present target mass-dependent equations that predict the onset of equal-mass hit-and-run impacts. Based on the new scaling laws for both head-on and oblique impacts, we derive the radii of post-collision remnants at different impact velocities for both types of impacts and illustrate the new mantle stripping curves alongside potential super-Mercuries in a mass–radius diagram.

3.1 Head-on impact simulations

Fig. 2 shows the mass and iron mass fraction of the largest remnant against the specific impact energy relative to the catastrophic disruption criteria for head-on impact simulations. Each symbol represents an impact simulation with different colours for different target masses. We find that for head-on collisions, the impact results separate into two regimes at a target mass of ∼2 M_⊕: Figs 2(a) and (c) show low target mass impacts (M_targ < 2 M_⊕) and Figs 2(b) and (d) show high target mass impacts (M_targ ≥ 2 M_⊕).

Figure 2.

Normalized mass (panels a and b) and iron mass fraction (panels c and d) of the largest post-collision remnant, plotted against the normalized impact energy Q_R/|$Q^{*}_\mathrm{RD}$| for all head-on collisions. Panels (a) and (c) display the results of impacts with target masses below 2 M_⊕, while panels (b) and (d) show the impact results with target masses above 2 M_⊕. The dashed lines in panels (a) and (b) represent the universal law from Leinhardt & Stewart (2012), where there is a break in the slope at Q_R/|$Q^{*}_\mathrm{RD}$| equal to 1.8, above which is the supercatastrophic disruption power law. The dashed curves in panels (c) and (d) show the scaled down version of the fit from Marcus et al. (2009) updated by Carter et al. (2018). The solid lines represent the scaling laws fitted from this work. Unfilled symbols show simulations with M_lr/M_tot less than 0.1.

3.1.1 Sharp change of mass and iron mass fraction

M_lr/M_tot shows several sharp changes in gradient: in Fig. 2(a) at around Q_R/|$Q^{*}_\mathrm{RD}$| ≃ 0.55 and 1.15, and in Fig. 2(b) at around 0.625. In Figs 2(c) and (d), the increasing trend of iron mass fraction also changes at the normalized impact energies where the first sharp change in remnant mass occurs.

In Fig. 2(c), the iron mass fraction of low target mass impacts starts to increase more slowly above Q_R/|$Q^{*}_\mathrm{RD}$| ≃ 0.55, and the maximum iron mass fractions of low target mass impacts are unlikely to exceed 50 per cent, which suggests a low mantle stripping efficiency for impacts in the low target mass regime. We note that Carter et al. (2018) found that the iron mass fraction in their differentiated planetesimal impacts reached a plateau at a value of around 0.4, which is similar to iron mass fraction trend shown here in Fig. 2(c). However, their break was around Q_R/|$Q^{^{\prime }*}_\mathrm{RD}$| = 1.04 (Carter et al. 2018 use |$Q^{^{\prime }*}_\mathrm{RD}$| to denote the catastrophic disruption threshold in their study), which is larger than the break 0.55 found in this work. The difference could be related to the different choice of initial core fraction of the target planets: the initial core fraction in Carter et al. (2018) is 22 per cent, while our targets have a core fraction of 30 per cent. According to Leinhardt & Stewart (2012), Denman et al. (2020), and Reinhardt et al. (2022), head-on collisions should be the most efficient way to strip off mantle materials in terms of impact energy, but this might not be the case for low target mass impacts. Therefore, it is unlikely that Mercury and Mercury-like exoplanets that are small and dense with a core fraction larger than 50 per cent were formed by a single roughly equal-mass head-on impact.

At impact energies higher than the sharp change of Q_R/|$Q^{*}_\mathrm{RD}$| at 0.625 (Fig. 2d), the iron mass fraction for high target mass impacts initially shows a slowly increasing trend similar to the low target mass impacts for a similar energy (Fig. 2c). However, at energies greater than Q_R/|$Q^{*}_\mathrm{RD}$| around 0.8–0.9 the iron mass fraction increases rapidly with increasing energy and can reach 80 per cent or greater.

3.1.2 Vaporization-enhanced mantle stripping

The difference between the iron mass fraction of low and high target mass impacts at high impact energies is likely related to the vaporization of core materials. Fig. 3 shows the vaporized fraction of all the core materials at the final step of each impact simulation for the low and high target masses. The vapour fractions are calculated using the lever rule, and the entropies and pressures on the iron vapour dome are taken from Stewart (2020). When normalized impact energy (Q_R/|$Q^{*}_\mathrm{RD}$|⁠) is small, the core vapour fraction of both low and high target mass impacts barely exceeds 10 per cent. As target mass increases, when Q_R/|$Q^{*}_\mathrm{RD}$| is larger than around 0.8–0.9, more than 10 per cent core materials in high target mass impacts can be vaporized, as shown in the bottom panel of Fig. 3.

Figure 3.

Vapour fraction of the core material right after the end of each impact simulation. Top panel shows the vapour fraction of all the impacts with target mass below 2 M_⊕, and bottom panel shows vapour fraction for impacts with target mass above 2 M_⊕. Different colours represent different target masses in each panel. Note that the scales of y-axes of the two panels are different.

Therefore, we propose that the mantle stripping process of differentiated planets’ impacts is dominated by two processes: kinetic momentum transfer and vaporization-induced ejection. At relatively low impact speeds, starting from V_i = V_esc, mantle materials begin to be expelled while only a small amount of core materials are kicked fast enough to overcome the gravitational binding energy of the system, since they are deeper and denser in the centre of planets compared to mantle materials. With increasing impact energy, some core materials begin to be expelled, as shown by the red dots and lines in Fig. 4. Core materials begin to be lost at around Q_R/|$Q^{*}_\mathrm{RD}$| = 0.5 for M_targ = 0.6 M_⊕ (representing low target mass impacts) and Q_R/|$Q^{*}_\mathrm{RD}$| = 0.6 for M_targ = 7.0 M_⊕ (representing high target mass impacts), which is consistent with the normalized energy where the first sharp change in remnant mass occurs (Fig. 2). The iron mass fractions start to increase more slowly when core materials begin to be lost, since partial of the impact energy is deposited into core materials. If the core vaporization fractions remain very low, the increasing trend in iron mass fraction will appear as shown in Fig. 2(c). In these simulations, we find that if the core vaporization fraction exceeds |$\sim\!\! 10{{\ \rm per\, cent}}$|⁠, the vaporized core material will enhance the mantle stripping efficiency and make the iron mass fractions increase again, as shown in Fig. 2(d).

Figure 4.

The mass-loss fraction of core and mantle materials, plotted against the normalized impact energy for head-on collisions of target mass 0.6 M_⊕ (top panel) and 7.0 M_⊕ (bottom panel).

3.1.3 Fragmentation of the largest remnant

Head-on impacts with low target masses will enter a fragmented regime when |$Q_\mathrm{R}/Q^{*}_\mathrm{RD}$| is larger than approximately 1.15. In this regime, the post-collision largest remnant will split into several smaller pieces, as shown in the bottom panel of Fig. 5. The fragmentation causes a rapid decrease in the largest remnant mass, as shown in Fig. 2(a), around the second sharp change (⁠|$Q_\mathrm{R}/Q^{*}_\mathrm{RD}$| = 1.15), which was also observed by Reinhardt et al. (2022). In Fig. 5, each colour represents a bound remnant. In the top panel, for a head-on impact at velocity 30 km s⁻¹ (Q_R/|$Q^{*}_\mathrm{RD}$| = 1.15), there is only one bound remnant. However, the bottom panel shows that the bound remnant splits into six smaller remnants at early times for the same mass impact with an impact speed of around 31 km s⁻¹ (Q_R/|$Q^{*}_\mathrm{RD}$| = 1.23).

Figure 5.

Evolution of bound particles in the post-collision remnants after two equal-mass 1.0 M_⊕ planets collided head-on at 30 km s⁻¹ (top panel) and 31 km s⁻¹ (bottom panel). Different colours represent different post-collision remnants. The green particles represent the largest remnant and the cyan particles are bound particles of the second largest remnant. In the bottom panel, from 10 to 20 h, the largest remnant changes from the top to the lower left due to the merging of two remnants.

We adopted several different approaches to verify this surprising fragmentation phenomenon. First, we tracked the position of the remnants for 10 yr using simple N-body simulations (for details on code, see Section 2.4), taking the centre of mass and centre of velocity of each remnant from the last snapshot (simulation time of 60 h) in the bottom panel and ignoring the effect of the central star. N-body simulations suggest that these remnants will move far away from each other and will not re-accrete to merge together again. Then, to check whether the fragmentation is caused by the h_max or density floor set in the swift simulations, we tested the same simulations with a larger h_max (up to 50 R_⊕). We found that the fragmentation did not change, which suggests that it is not due to h_max being too small or the minimum resolved mass density being too large. As an additional comparison, using the same initial condition file and equation of state tables, we also performed equivalent impact simulations with the gadget-2 SPH code (Ćuk & Stewart 2012) that does not have a density floor and found that similar fragmentation occurred (although the mass of each remnant varied slightly – within 10 per cent – from that found in the swift simulations). Finally, we ran a friend-of-friend search (Creasey 2018) with a relatively large linking length (e.g. 5 R_⊕) and found that all particles were loosely bound together. We then applied our bound particle searching method to these particles to check whether the fragmentation could be an artefact of the searching algorithm. Despite this, we still found the same splitting scenario with multiple remnants. Since all these tests show the fragmentation phenomenon, we infer that the fragmentation is likely due to the equation of state used and might be related to the phase transition from liquid to vapour.

It is noteworthy that the fragmentation of the largest remnant does not only occur within the small target mass regime. Fragmentation occurs when M_lr/M_tot is around 0.2 for target masses of 2.0 and 3.0 M_⊕. For target masses higher than 7.0 M_⊕, fragmentation occurs when M_lr/M_tot is less than 0.1, which is defined as the supercatastrophic disruption region (Leinhardt & Stewart 2012). Simulation data points of target masses 1.0 M_⊕ (red dots) and 1.58 M_⊕ (pink dots) in Figs 2(a) and (c) suggest that, with increasing target mass, the normalized energy required to trigger fragmentation also increases, and the mass and iron mass fraction of low target mass impacts will become increasingly similar to those of high target mass impacts. The fragmentation of the largest remnant in the small target mass regime might be crucial for the smaller planetesimal collision stage, where target masses are small but impact velocities are fairly large, which makes head-on or nearly head-on impacts more likely to happen. Previous N-body simulations have tended to use the universal law derived by Stewart & Leinhardt (2009) and Leinhardt & Stewart (2012) to predict the masses of remnants after collision. This law defines that fragmentation or supercatastrophic impacts only occur when M_lr/M_tot is less than 0.1. However, based on our new simulations, fragmentation could occur very early when M_lr/M_tot is as large as 0.4. We plan to conduct a more detailed study of the fragmentation in our following papers.

3.1.4 Scaling laws

We fit piecewise-defined scaling laws for the mass and iron mass fraction of the largest post-collision remnant separately for low and high target mass head-on impacts based on the locations of the first and second sharp changes. For low target masses,

$$\begin{eqnarray} \frac{M_\mathrm{lr}}{M_\mathrm{tot}}= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 - 0.91 \Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}\Bigl)^{2.27} & \text{if } \frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}} \le 0.55, \\\\ 0.75 - 0.50\Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}-0.50\Bigl) & \text{else}, \\ \end{array}\right. \end{eqnarray}$$

(4)

$$\begin{eqnarray} \frac{M_\mathrm{Fe}}{M_\mathrm{lr}}= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0.3 + 0.39\Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}\Bigl)^{2.55} & \text{if} \frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}} \le 0.55, \\\\ 0.38 + 0.12 \Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}-0.5\Bigl) & \text{else}, \\ \end{array}\right. \end{eqnarray}$$

(5)

and for high target masses,

$$\begin{eqnarray} \frac{M_\mathrm{lr}}{M_\mathrm{tot}}= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1 - 0.86 \Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}\Bigl)^{2.77} & \text{if } \frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}} \le 0.625, \\\\ 0.75 - 0.86 \Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}-0.625\Bigl)^{1.24} & \text{else}, \\ \end{array}\right. \end{eqnarray}$$

(6)

$$\begin{eqnarray} \frac{M_\mathrm{Fe}}{M_\mathrm{lr}}= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0.3 + 0.25\Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}\Bigl)^{2.6} & \text{if} \frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}} \le 0.625, \\\\ 0.37 + 0.66 \Bigl (\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{RD}}-0.625\Bigl)^{2} & \text{else}. \\ \end{array}\right. \end{eqnarray}$$

(7)

These scaling relations are shown using solid lines in Fig. 2.

For low target mass impacts, the first sharp change is around Q_R/|$Q^{*}_\mathrm{RD}$| equaling 0.55, and we only use simulation data where Q_R/|$Q^{*}_\mathrm{RD}$| is less than 1.15. Beyond this energy, impacts tend to become fragmented and M_lr/M_tot shows a sudden drop from 0.4 to 0.1. For high target mass impacts, the sharp change occurs around Q_R/|$Q^{*}_\mathrm{RD}$| equaling 0.625, and we only use simulations where M_lr/M_tot > 0.1, as impacts tend to be supercatastrophic below this mass and the largest remnant might not be resolved due to low number of particles. The boundary between low and high target mass scaling laws is not strict. Although fragmentation can still happen when the target mass is above 2.0 M_⊕, the fitted scaling laws can predict the mass of the largest remnant even after splitting when M_targ is above 2.0 M_⊕. Therefore, we recommend using high target mass scaling laws for M_targ equal to or above 2.0 M_⊕ and low target mass scaling laws for target mass below 2.0 M_⊕.

3.2 Oblique impact simulations

Fig. 6 presents the mass and iron mass fraction of the largest post-collision remnants versus normalized specific impact energies for oblique impacts. Note that for oblique impacts, we define target catastrophic disruption criteria, |$Q^{^{\prime }*}_\mathrm{TD}$|⁠, as the specific impact energy required to dispel half the mass of the target. Therefore, the mass of the largest remnant is now represented by M_lr/M_targ. We select and plot simulation data with remnant masses M_lr/M_targ greater than 0.1 with filled symbols, below which the number of particles in the post-collision remnants may not be sufficient for convergence. Figs 6(a) and (c) show results with target mass less than 2.0 M_⊕, and Figs 6(b) and (d) show results for high target mass above 2.0 M_⊕. Oblique impacts do not show a substantial difference in the mass of the largest remnant for different target mass regimes. However, we notice that high target mass impacts tend to have slightly lower remnant masses compared to low target mass impacts at higher impact energies with Q_R/|$Q^{^{\prime }*}_\mathrm{TD}$| larger than 1.0.

Figure 6.

Normalized mass (panels a and b) and iron mass fraction (panels c and d) of the largest post-collision remnant, plotted against the normalized impact energy for all oblique collisions with low (left panels a and c) and high (right panels b and d) target masses. Unfilled symbols show the simulations with M_lr/M_targ less than 0.1.

In particular, at high impact energies with low impact parameters (b = 0.2, blue symbols, and b = 0.3, green symbols), the iron mass fractions of low target mass impacts (Fig. 6c) appear to be lower compared to those of higher target mass impacts, especially M_targ = 0.06 and 0.12 M_⊕ (squares and circles). Among all impact parameters tested, the mass and iron mass fraction of b = 0.2 impacts follow a different pattern compared to others. Also, at b = 0.2 and b = 0.3, an impact tends to transition abruptly from accretion to erosive hit-and-run without going through a ‘clean’ hit-and-run impact, and thus M_lr/M_targ starts below 1.0 at these impact parameters. On the other hand, at b = 0.7, impacts with target mass larger than 3.0 M_⊕ at higher energies tend to have much lower iron mass fraction (M_Fe/M_lr around 60 per cent) compared to low target mass high-energy impacts (M_Fe/M_lr around 90 per cent), as shown in Fig. 6(d).

We fit scaling laws for erosive hit-and-run oblique impacts, as these are the impacts that can form dense planets. Impacts with velocities lower than the hit-and-run velocity can be approximated as perfect-merging events. Therefore, only simulation data with M_lr/M_targ less than 1.0 (beginning of erosive hit-and-run) and larger than 0.1 are selected to fit the scaling laws. The general equations for mass and iron mass fraction for oblique impacts are as follows:

$$\begin{eqnarray} \frac{M_\mathrm{lr}}{M_\mathrm{targ}} =1.4\left(\alpha _{M,\,b} \right)^{\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{TD}}}+(0.5-1.4\alpha _{M,\,b}), \end{eqnarray}$$

(8)

$$\begin{eqnarray} \frac{M_\mathrm{Fe}}{M_\mathrm{lr}} =0.3+\alpha _{\mathrm{ Fe},\,b}\left(\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{TD}}\right)^{\beta _{\mathrm{ Fe},\,b}}. \end{eqnarray}$$

(9)

Except for impacts at b = 0.2, all impacts follow a similar trend for both mass and iron mass fraction. Thus, we first fit general scaling laws by combining all the simulation data from b = 0.3 to b = 0.7, as shown in Fig. 7, giving α_{M, b} = 0.576, α_{Fe, b} = 0.238, and β_{Fe, b} = 1.356 with a coefficient of determination r² of 0.973 and 0.976 for mass and iron mass fraction, respectively. In order to balance the error, especially for iron mass fraction, we exclude the iron mass fraction data points with Q_R/|$Q^{^{\prime }*}_\mathrm{TD} \gt 1.3$| for target mass of 0.06 and 0.12 M_⊕ at b = 0.3, and target mass larger than 3 M_⊕ at b = 0.7, which are shown as unfilled symbols in Fig. 7.

Figure 7.

General fitted scaling laws of mass (top panel) and iron mass fraction (bottom panel) of the largest remnant for oblique impacts. Open symbols are excluded when fitting the scaling laws. Data with b = 0.2 are not plotted.

When M_lr/M_targ is equal to 1.0, Q_R/|$Q^{^{\prime }*}_\mathrm{TD}$| is not near zero, especially at low impact parameters. M_lr/M_targ starts to change from greater than 1.0 to less than 1.0 when impacts transition from accretion to erosive hit-and-run. When M_lr/M_targ is close to 1.0, the value of Q_R/|$Q^{^{\prime }*}_\mathrm{TD}$| is determined by the hit-and-run velocity. For equal-mass oblique impacts, lower impact parameters tend to result in larger hit-and-run velocities (discussed in Section 3.4), and thus larger Q_R/|$Q^{^{\prime }*}_\mathrm{TD}$| values. As a result, the data points for mass and iron mass fraction for each impact parameter slightly diverge from each other when M_lr/M_targ is above ∼0.7. To more accurately model mass and iron mass fraction, we found that the parameters α_{M, b}, α_{Fe, b}, and β_{Fe, b} in equations (8) and (9) can be described by the following equations at different impact parameters:

$$\begin{eqnarray} \alpha _{M,\,b} = 1.602\log _{10}(b)-0.848b+1.478, \end{eqnarray}$$

(10)

$$\begin{eqnarray} \alpha _{\mathrm{ Fe},\,b} = 0.085b^2+0.218, \end{eqnarray}$$

(11)

$$\begin{eqnarray} \beta _{\mathrm{ Fe},\,b} = 2.5\times (0.1)^b+b, \end{eqnarray}$$

(12)

where b is the impact parameter between 0.3 and 0.7. Since impacts at b = 0.2 are in the transition regime, and have relatively different trends of mass and iron mass fraction, we fit scaling laws at b = 0.2 individually. At b = 0.2, the mass follows

$$\begin{eqnarray} \frac{M_\mathrm{lr}}{M_\mathrm{targ}} =2.8\times \left(0.23 \right)^{\frac{Q_\mathrm{R}}{Q^{*}_\mathrm{TD}}}+ (0.5-2.8\times 0.23), \end{eqnarray}$$

(13)

and the iron mass fraction follows the same form as equation (9) with α_{Fe, 0.2} = 0.199, and β_{Fe, 0.2} = 1.904. Iron mass fraction data points for target masses below 0.6 M_⊕ are excluded from fitting at b = 0.2.

3.3 Catastrophic disruption criteria

Leinhardt & Stewart (2012) reported that, for equal-mass head-on impacts, the catastrophic disruption criteria or the principal disruption curve follows

$$\begin{eqnarray} Q^{*}_{\mathrm{RD},\,\gamma =1} = c^{*}\frac{4}{5}\pi \rho _{1}GR^2_\mathrm{C1}, \end{eqnarray}$$

(14)

where c* is a measure of the catastrophic disruption threshold in units of the gravitational binding energy, ρ₁ = 1000 kg m⁻³, and R_C1 = |$(\frac{3 M_\mathrm{{tot}}}{4\pi \rho _1})^{\frac{1}{3}}$| is the radius a spherical body would have if it had the total system mass and a density of ρ₁. Leinhardt & Stewart (2012) found c* to be 1.9 ± 0.3 (for planetary-mass bodies), while Denman et al. (2020) fitted c* for planets with atmospheres to be around 2.52. Based on all of our head-on impact data, we found a c* value of 2.926. The top panel of Fig. 8 shows the relationship between |$Q^{*}_\mathrm{RD,\,\gamma =1}$| and R_C1, along with the fitting lines from our work and previous work.

Figure 8.

Relationship between R_C1 and catastrophic disruption criteria |$Q^{*}_\mathrm{RD}$| for head-on impacts (top panel) and target catastrophic disruption criteria |$Q^{^{\prime }*}_\mathrm{TD}$| for oblique impacts (bottom panel). The unfilled symbol at b = 0.2 and M_targ = 0.35 M_⊕ was excluded from fitting.

In Denman et al. (2020), the target planets have thick and massive atmosphere layers, which can account for up to 30 per cent of the planet’s total mass. For planets with the same mass, including an atmosphere increases the planet’s radius and decreases the total gravitational binding energy. This means that less impact energy is required for catastrophic disruption. This explains why their value of c* is slightly lower than ours.

The difference between our results and those of Leinhardt & Stewart (2012) may also be related to the different internal structures of the target planets used in the simulations. Leinhardt & Stewart (2012) fitted their c* by combining previous hydrodynamic simulation results, including collisions between differentiated iron-rock bodies, differentiated ice-rock bodies, and pure rock bodies. The lower c* value, compared to the fitted c* from iron-rock bodies of this study, is a result of combining a range of different density planets, as ice-rock and pure rock bodies of the same masses have lower gravitational binding energies. In addition, the choice of equations of state and number of fitting points used may also affect the final result. Our simulation data can be directly used for fitting since they are all equal-mass impacts. However, in Leinhardt & Stewart (2012), a correction related to the mass ratio, γ, needs to be applied first to the raw data to convert the results to the equal-mass equivalent, which necessarily introduces some errors. The new M-ANEOS equation of state tables used in this work are better at modelling the phase transition from liquid to vapour state, and therefore will have a more accurate representation around the catastrophic disruption criteria.

For oblique impacts with the same mass, we found that the newly defined target catastrophic disruption criteria |$Q^{^{\prime }*}_\mathrm{TD}$|⁠, which represent the specific impact energies required to expel half of the target’s mass, follow a similar trend as equation (C3) for head-on impacts but with larger constant parameters. |$Q^{^{\prime }*}_\mathrm{TD}$| at different impact parameter can be modelled as

$$\begin{eqnarray} Q^{^{\prime }*}_\mathrm{TD} = \left(0.175 \left(\frac{1}{1-b} \right)^{3.3}+0.635\right)c^{*} \frac{4}{5}\pi \rho _{1}GR^2_\mathrm{C1}, \end{eqnarray}$$

(15)

where c* is the same as fitted for head-on impacts. The bottom panel in Fig. 8 shows target catastrophic disruption criteria for oblique impacts and the corresponding fitting lines.

3.4 Hit-and-run velocities

Fig. 9 shows the normalized velocities where hit-and-run impacts begin to happen for different impact parameters for our equal-mass oblique impacts. Hit-and-run velocities become more target mass-dependent with decreasing impact parameter. Previously published hit-and-run criteria, as shown in the top panel of Fig. 9, fit the general trend of hit-and-run velocities relatively well, but cannot account for the variation with target mass.

Figure 9.

The relationship between normalized hit-and-run velocities and impact parameter. V_esc is the mutual escape velocity. Colours represent the mass of the target in Earth masses. Top panel: The dashed line represents the hit-and-run velocity criteria of Kokubo & Genda (2010) at γ = 1; the dash–dotted line represents the criteria from Denman et al. (2022); and the dotted line represents the fitting results from Gabriel et al. (2020). Our fitted criteria at M_targ = 1 M_⊕ are represented by the solid black line. Bottom panel: Solid coloured lines show the target mass-dependent hit-and-run criteria described by equation (C5).

The dependence on target mass at low impact parameter is related to the core radius ratio, R_core/R_planet. The core radius ratio of each target planet decreases with increasing mass, for example, R_core/R_planet is 0.52 for a target planet with a mass of 0.06 M_⊕, while it is 0.48 for a target with a mass of 20 M_⊕. In an impact where the target and impactor are the same size, the impact parameter at which the core material from both planets no longer intersect along the impact direction is the core radius ratio (b = R_core/R_planet, if the deformation of the planets when they approach each other is not accounted for, Fig. B1). Below this limit, the projected length of the impactor’s core overlapping the target’s core is

$$\begin{eqnarray} l_\mathrm{core} = 2R_\mathrm{planet} \left(\frac{R_\mathrm{core}}{R_\mathrm{planet}}-b \right). \end{eqnarray}$$

(16)

Therefore, at the same impact parameter, lower mass targets have slightly larger core overlap and thus require higher impact velocities to transition into hit-and-run.

We fit the hit-and-run velocity as a function of target mass, as shown in the bottom panel of Fig. 9. We do not include the data for b = 0.1, as the impact angle is too small and the tidal deformation of the planets makes impacts behave very similarly to head-on impacts. Our target mass-dependent hit-and-run criterion is

$$\begin{eqnarray} \frac{V_\mathrm{HnR}}{V_\mathrm{esc}} = \Gamma (1-b)^{3.6}+1.2, \end{eqnarray}$$

(17)

where Γ = −0.3 log₁₀M_{targ, ⊕} + 2.18, and M_{targ, ⊕} is the mass of target planet in Earth masses. Since all our simulations have impactor-to-target mass ratio γ = 1, the hit-and-run criterion does not account for the effect of varying γ.

4 DISCUSSION AND CONCLUSIONS

4.1 Prediction of impact conditions

Prediction of impact conditions required to achieve a particular mass and core fraction can be made by combining the mass and iron mass fraction scaling laws derived in the previous sections. As we are motivated to model the formation of dense planets, for the head-on scaling laws, we only use the high target mass scaling law equations (6) and (7) when Q_R/|$Q^{*}_\mathrm{RD}$| is larger than 0.625 because low target mass head-on impacts can at most form a remnant with 50 per cent iron mass fraction. Applying our new c* = 2.926 and combining equations (6), (7), and (C3), Figs 10(a) and (c) show the prediction of head-on impact velocities and iron mass fractions derived from both our scaling laws and the combination of those from Leinhardt & Stewart (2012) and Marcus et al. (2010). For the same remnant mass and iron mass fraction, we require higher impact velocities compared to Leinhardt & Stewart (2012) as our c* are larger. At different impact velocities, our iron mass fraction predictions converge at around 0.9. From the prediction of this work, the impact velocities necessary to produce a super-Mercury with mass above 5 M_⊕ and iron mass fraction larger than 60 per cent are greater than 60 km s⁻¹, which suggests an impact occurring close to the star.

Figure 10.

Prediction of impact conditions for head-on impacts (left panels a and c) and oblique impacts (right panels b and d). (a) and (b): The impact velocity versus the mass of the largest remnant for iron mass fraction of 60, 70, and 80 per cent. (c) and (d): The iron mass fraction versus the mass of the largest remnant for different impact velocities. Dashed lines show the prediction derived by combination of Leinhardt & Stewart (2012) and Marcus et al. (2010).

For oblique impacts, combining equations (8), (9), and (C4) in a similar way, we derive the prediction of impact conditions for b = 0.3, 0.5, and 0.7 shown in Figs 10(b) and (d). For the same remnant mass and iron mass fraction, larger impact parameter requires higher impact velocities. At b = 0.7, a remnant with mass above 5 M_⊕ and 60 per cent iron mass fraction would need an impact velocity greater than 200 km s⁻¹. As shown in Fig. 10(d), for a remnant mass of 5 M_⊕ and a velocity of 80 km s⁻¹, an impact at b = 0.7 could only produce an iron mass fraction less than 40 per cent; while at b = 0.3, the remnant could have an iron mass fraction of around 65 per cent. With increasing impact parameter, the mantle stripping efficiency drops significantly.

The scaling laws can also be used to predict the pre-impact target masses required. For an exoplanet with a given mass and core fraction, Fig. 11 provides a good reference for potential target masses and impact velocities that could form the planet via a giant impact. Fig. 12 illustrates an alternative prediction of the mass of target bodies based on the iron mass fraction at different impact velocities. Comparing it to oblique impacts with b equal to 0.5, head-on impacts require a higher target mass but lower impact velocities for a remnant with a specific mass and iron mass fraction.

Figure 11.

The relationship between the mass of the largest remnant and the inferred target mass at different impact velocities. Here, we present predictions for impact parameters of b = 0.0, b = 0.3, and b = 0.4.

Figure 12.

The relationship between the mass of the largest remnant and the inferred target mass at different iron mass fractions. To avoid overlap of lines, we present only predictions for impact parameters of b = 0.0 and b = 0.5.

Using the method discussed in Section 2.2, we derived the radius² of post-collision remnants at different impact velocities for both head-on and oblique impacts, as shown in Fig. 13. These radii do not correspond to the radii of the simulated largest remnants immediately after the collision, but rather to those of a super-Earth with an assumed surface temperature of 1000 K and a core fraction Z_Fe after a long cooling period. Compared to the maximum stripping line [originally from Marcus et al. (2010), and later corrected by Reinhardt et al. (2022)], our maximum stripping line predicts slightly smaller radii for planets. Note that several super-Mercuries, such as Kepler-406 b and Kepler-107 c, fall just below the 80 km s⁻¹ head-on stripping line, suggesting that they could require a higher impact speed than 80 km s⁻¹ to form via impacts.

Figure 13.

Mass–radius diagram for Mercury and super-Mercuries. The red solid line shows the Reinhardt et al. (2022) corrected maximum stripping line generated using the method described in Section 2.2 for head-on impacts at 80 km s⁻¹ and the black solid line shows our new maximum stripping line. The dashed lines show the maximum stripping lines of oblique impacts. Blue solid line represents radii of planets made by 100 per cent core generated using the same method described in Section 2.2. The data points show several possible super-Mercuries, GJ 367 b, Kepler-406 b, Kepler-107 c, K2-38 b, and K2-229 b.

4.2 Comparison with Denman et al. (2020)

We note that Denman et al. (2020) found a similar trend of M_lr/M_tot, as shown in their fig. 4. Denman et al. (2020) found that, above an impact energy Q_piv at higher impact energies Q_piv/|$Q^{*}_\mathrm{RD}$|⁠, M_lr/M_tot will have a steeper decreasing trend, and the loss efficiency of the atmosphere will decrease, which is due to the onset of mantle and core erosion. Their Q_piv/|$Q^{*}_\mathrm{RD}$| is close to our normalized impact energies where the first sharp change happens. The difference is that their Q_piv/|$Q^{*}_\mathrm{RD}$| is target mass-dependent, and the break in M_lr/M_tot slope occurs at slightly different Q_R/|$Q^{*}_\mathrm{RD}$| for different target masses. We infer that their changing of Q_R/|$Q^{*}_\mathrm{RD}$| is due to the different amount of atmosphere (ranging from 8 to 33 per cent) in their three different target planets. For targets with low atmosphere mass ratio, the impact energy required to expel enough atmosphere for the core and mantle to start being lost is smaller compared to targets with higher atmosphere mass ratio. Combined results from their work and this work show that a universal law having a single trend may not be enough to describe all the collision processes that happen within multilayered planets well. The universal law from Stewart & Leinhardt (2009) and Leinhardt & Stewart (2012) works well for undifferentiated planetesimals, but when dealing with larger planets with two, three, or even four layers with oceans and atmospheres, due to the different shock impedances of different materials, the mass-loss shows different trends. In Fig. A2, we compare scaling laws from this work to scaling laws derived from Marcus et al. (2009), Denman et al. (2020), and Reinhardt et al. (2022).

4.3 Polar head-on impact

Fig. 12 illustrates that, according to the current scaling laws, a target mass of at least 10 M_⊕ is typically required to generate a relatively dense remnant with mass greater than 5 M_⊕ and core fraction greater than 60 per cent. To form a significantly denser planet, the required target mass would need to exceed 20 M_⊕. However, Crida & Bitsch (2017) suggest that planets with a mass above a critical value of 20 M_⊕ must accrete gas and grow large enough to create a gap in the disc, in order to avoid being consumed by the star during fast Type-I migration. Given that planets with greater mass are more likely to undergo runaway accretion, it would be preferable to have targets with smaller mass as candidates for forming dense planets.

For high target mass impacts, head-on impacts still would be the most efficient way to strip off large amounts of mantle. To form a massive super-Mercury (e.g. 5.0 M_⊕ with 80 per cent core fraction) via a giant impact requires either an extremely energetic head-on impact resulting in catastrophic disruption of a target or some new impact configurations. If pre-impact spin is included, the effective gravitational potential is reduced and in extreme spinning cases (Lock & Stewart 2017; Lock et al. 2018) the target is deformed making the core more exposed on the poles. To generate super-Mercuries with relatively low target masses and impact velocities compared to common head-on impacts, we propose an alternative impact configuration involving a fast-spinning target to increase the post-collision remnant iron mass fraction.

An erosive equal-mass polar impact on to a fast-spinning protoplanet could potentially create a planet with a higher core fraction and larger remnant mass with lower impact velocity compared to common head-on giant impacts in the equatorial plane. A fast-spinning target with a spin period close to its instability rotational limit will deform greatly from a sphere to an oblate spheroid with mantle material concentrated at the equator. The core material is then closer to the surface in the polar regions. The core from an impactor can more easily penetrate through the thin mantle near the polar region of the target efficiently mixing with the target core. The impact shock generated from the polar impact will tend to kick more mantle material off than an impact on to a non-spinning target.

For example, run 2 in Table 1 (a pole-on impact with a fast-spinning target, P = 1.9 h, mass = 11.9 M_⊕, and an equal-mass non-spinning projectile with an impact velocity of 76 km s⁻¹) generates a post-collision body with a mass of 5.15 M_⊕ and core fraction of 78.3 per cent. The pole-on impact requires less total impact mass (23.8 M_⊕) and lower impact speed (76 km s⁻¹) in contrast with the equatorial head-on impact that requires 32.4 M_⊕ and an impact speed of 83.8 km s⁻¹ (run 1 in Table 1).

Table 1.

Potential giant impact formation of an example super-Mercury with mass 5.0 M_⊕ and 80 per cent core fraction. M_targ mass of initial pre-impact target in the unit of Earth mass; γ mass ratio between target and projectile; V_imp impact speed in km s⁻¹; b impact parameter (zero indicates head-on equatorial impact if hit direction is −X, pole-on impact if hit direction is −Z); P spin period of a pre-impact target in hours; M_lr mass of the largest remnant after 200 h simulation time; and Z_Fe iron element mass fraction in the largest remnant after a collision. Impacts are set up with the centre of mass of two impactors located at the origin with zero net velocity. For equatorial head-on impacts, impactors are moving towards each other along the X-axis with zero Y and Z direction velocities, while for polar impacts, impactors are moving towards each other along Z-axis with zero X and Y direction velocities. Hit direction −X and −Z means that the projectile is colliding with the target with the velocity pointing to the negative side of the X-axis and Z-axis, respectively.

Run	Mtarget/M⊕	γ	Vimp (km s−1)	b	P (h)	Mlr/M⊕	ZFe (per cent)	Hit direction
1	16.21	1.00	83.8	0	0	5.10	78.9	−X
2	11.90	1.00	76.0	0	1.9	5.15	78.3	−Z

Through stochastic accretion and migration stages, embryos and protoplanets can have a significant tilt with respect to their orbit. Also, rapidly rotating bodies are a common outcome of off-axis collisions, and average angular velocity after accretion or several giant impacts could be as high as the critical angular velocity (Kokubo & Genda 2010; Nakajima & Stevenson 2015; Lock & Stewart 2017). Adopting a realistic accretion model, N-body simulations (Kokubo & Genda 2010) show that spin obliquity could range from 0 to 180° following an isotropic distribution. Hence, the pole-on impact would not need a projectile embryo to have an orbit perpendicular to the orbital plane of the target. More detailed study of polar impacts in a wider parameter space will be carried out in future work.

4.4 Vaporization versus kinetic kicking

Transiting from head-on to oblique impacts, at low impact parameters, there will still be a significant amount of intersection between the target and impactor cores along their initial velocity vector. Therefore, core vaporization could still contribute to the mantle stripping process at relatively low impact angles. We examine core vaporization in Fig. 14 for impacts at different impact parameters. At b = 0.2 and b = 0.3, impacts with target mass M_targ = 0.06 and 0.1 M_⊕ have nearly zero core vaporization and therefore only kinetic kicking contributes to the mantle stripping process, leading to relatively lower mantle stripping efficiency compared to higher target mass impacts that have contributions from both kinetic kicking and core vaporization. At higher impact parameters (b = 0.4 and b = 0.5), low target mass impacts still have a low core vapour fraction. However, their iron mass fraction can be as high as high target mass impacts. This suggests that kinetic momentum transfer dominates the mantle stripping process at high impact angles.

Figure 14.

Fraction of vaporized core material of different target masses at various impact parameters.

The iron mass fraction of the very small target mass impacts (M_targ = 0.06 and M_targ = 0.1) at b = 0.2 and b = 0.3 (blue and green square and circle symbols in Fig. 6c, respectively) is roughly 5–10 per cent lower than the higher target mass impacts. However, these low target mass impacts still generally have higher iron mass fraction than the corresponding head-on impacts that have maximum iron mass fraction around 50 per cent (see Fig. 2c). This also shows the limited role of core vaporization in oblique impacts. Impacts at b = 0.2 represent the transition from head-on impacts to normal oblique impacts. Below this impact parameter, both kinetic momentum transfer and vaporization-induced ejection can greatly affect the result of impacts. As the impact parameter increases, the intersection between the cores becomes smaller, and mantle stripping is gradually dominated by the kinetic ejection of mantle materials. For head-on impacts, the difference between mantle loss and core loss is around 30 per cent after the initial loss of the core, as shown in Fig. 4. However, due to the misalignment nature of oblique impacts, the difference could be as high as 70 per cent. As a result, unlike vaporization-dominated ejection of mantle for head-on impacts, mantle loss for oblique impacts is mainly dominated by kinetic kicking. Therefore, even low target mass impacts can form dense post-collision remnants.

At high impact parameters (b = 0.5 and b = 0.7) and high impact energy (Q_R/|$Q^{^{\prime }*}_\mathrm{TD}\gt 1.3$|⁠), the iron mass fraction of impacts with a target mass larger than 7 M_⊕ is surprisingly lower compared to smaller target mass impacts (Fig. 6d). Due to the extremely high impact energies, the core and mantle materials in these impacts with high target masses and high impact parameters are mostly in a supercritical state after the collision. The difference in iron mass fraction at these high impact parameters may be due to the different phase transition processes experienced by the core and mantle materials during normal and extremely high energy impacts. For these high target mass impacts, a greater amount of kinetic energy is dissipated due to the phase transition, which might in turn lead to lower mantle stripping efficiency. Furthermore, since these impacts occurred at extremely high impact velocities (V_i > 200 km s⁻¹) and energies, the EoS tables we used may not be sufficient to model the state of the material under these extreme conditions. We plan to investigate this further in our future work.

4.5 Caveats

Although the revised M-ANEOS equations of state (Stewart et al. 2019; Stewart 2020) can better model phase transitions, we extended these EoS tables to higher maximum densities (100 |$\mathrm{g\, cm}^{-3}$| for forsterite and 200 |$\mathrm{g\, cm}^{-3}$| for iron), which exceed the phase space that the tables were originally designed to deal with into a region with no experimental data to confirm material properties. This will mainly affect the high target mass impacts at high impact energies. The mass and iron mass fraction for high target mass impacts shown in Figs 2(b) and (d) do not show significantly different behaviour at high impact energies compared to relatively low target mass impacts, which suggests that any error introduced by extending the EoS tables is negligible.

All simulations in this work used the same target and impactor for each impact. Therefore, the scaling laws derived from these simulations apply conservatively only to equal-mass impacts. To better understand the impact conditions that can form dense planets, we tested oblique impacts with impact parameters up to b = 0.7, so all the fitted scaling laws are valid below b = 0.7. For impact parameters above b = 0.7, the probability of having energetic erosive hit-and-run events is very low. Additionally, the impact velocities required to form a post-collision remnant with an iron mass fraction larger than 60 per cent could easily exceed 10 times the mutual escape velocities.

As we expect lower mass projectiles to cause less erosion for a given impact energy, equal-mass impacts provide a useful ‘lower limit’ and our study can provide useful reference to constrain the formation conditions of super-Mercuries. Although using only equal-mass impacts narrows the parameter space, the data are consistent and allow for more accurate fitting of scaling laws.

All impact simulations in this work only represent the collision outcome immediately after giant impacts, within about 30 h. This does not account for cooling or the effect of re-accretion of the ejecta after a long period of time. Benz et al. (2007) reported that up to approximately 40 per cent of ejected materials will accumulate back to Mercury after several million years. However, the amount of re-accumulation depends on the configuration of an impact and the specific structure of the planetary system. In addition, Spalding & Adams (2020) showed that the primordial solar wind could provide sufficient drag on ejected debris to remove them from Mercury-crossing trajectories before re-accreting back to the planet’s surface. This suggests that the giant impact hypothesis remains a viable pathway towards the formation of dense planets, particularly those close to young stars. Therefore, the results derived from this work give an upper limit of the core fraction and a lower limit of the mass of super-Mercuries formed by giant impacts.

4.6 Dilemma for the formation of super-Mercuries

So far, more than 5000 exoplanets have been discovered but less than 10 planets are confirmed to be candidate super-Mercuries. Head-on super-Earth impacts theoretically are the most efficient mantle stripping method; however, compared to oblique impacts, they are very rare. Thus, impact-generated super-Mercuries would most likely be formed by one or a sequence of oblique giant impacts. Oblique impacts are more likely to happen around 45° or b = 0.7, with probability decreasing with decrease of impact angle. However, at high impact parameter as b = 0.7, the mantle stripping efficiency is extremely low, as shown in Figs 6 and 10. At b = 0.7, the required impact velocities to form a super-Mercury with 60 per cent core even above 1 M_⊕ are unrealistically high. Balancing the low probability of impacts at low impact parameters and low mantle stripping efficiency at high impact parameters, we conclude that impact-generated super-Mercuries would more likely occur with impacts between b = 0.3 and b = 0.5 (17–30°).

It is also worth noting that previous N-body simulations normally have an inner semimajor axis cut-off between 0.1 and 0.5 au, while several potential super-Mercuries have even smaller orbits around star; e.g. Kepler-107 c has orbit radius 0.06 au and GJ 367 b has orbit around 0.007 au. The conditions of impacts that occur at such short orbital distances need further study in order to determine the likelihood of super-Mercury-forming collisions.

5 CONCLUSIONS

In this work, we present over 1000 giant impact SPH simulations, ranging from head-on to oblique impacts. All simulations feature equal-mass impacts, with target mass spanning a wide range from 0.06 to 20 M_⊕. We find that the vaporization of core materials during collisions can greatly affect the outcome of giant impacts. We propose that the stripping of mantle materials is mainly controlled by two mechanisms: (1) core vaporization-enhanced stripping and (2) kinetic momentum transfer. For target planets with a mass less than 1 M_⊕, head-on impacts fail to form remnants with an iron mass fraction larger than 50 per cent, due to the lack of core vaporization. For small target mass impacts, we find that the largest remnant tends to fall into a new fragmentation regime, splitting into several smaller pieces at a small Q_R/|$Q^{*}_\mathrm{RD}$| value of around 1.15, which corresponds to much smaller impact energies compared to those predicted by previous studies. We plan to further investigate the reasons causing fragmentation in our follow-up paper. For target masses above 2 M_⊕, the mass and iron mass fraction show a clear trend, and we fit a new scaling law for these impacts.

Oblique impacts do not show significant differences between the low and high target mass regimes. We fit mass and iron mass fraction scaling laws for oblique impacts and propose the new target mass-dependent hit-and-run velocity criteria. Based on our head-on data, we fit a more accurate free parameter c*, 2.926, used to define the catastrophic disruption criteria for two-layered super-Earth impacts. We also propose a further angle correction for target catastrophic disruption criteria. Combining the catastrophic disruption criteria and the new scaling laws, one can easily predict the outcome for both head-on and oblique impacts. Based on our new scaling laws, we describe the necessary impact conditions required to form planets with various mass and iron mass fraction. Combining these with interior structure models for super-Earths, we update the ‘maximum mantle stripping line’ at different impact velocities and angles. Additionally, we provide a prescription on how to use these scaling laws in Appendix  C.

We demonstrate that for oblique impacts, the efficiency of mantle stripping decreases significantly as the impact angle increases. At the most common impact angle of 45°, an exceedingly high impact velocity is required to create a relatively dense planet. This presents a challenge in the formation of super-Mercuries through a single giant impact. However, the majority of currently observed dense planets have extremely short orbital periods and consequently high orbital velocities. This mitigates the restrictions on the impact conditions and makes the scenario of impact formation plausible.

These simulations provide important insights into the formation and evolution of dense planets. Future studies could build on our work by exploring a wider range of parameters and conditions, including different impactor-to-target mass ratio, the effect of spin, different chemical composition, different stellar type, and long-term evolution of the post-collision remnants.

Our work provides a useful framework for interpreting observational data of exoplanets and understanding the formation history of planetary systems. With more planets discovered in the future, we hope that we can better understand the formation of our own Solar system and the prevalence and diversity of dense rocky exoplanets throughout the Galaxy.

ACKNOWLEDGEMENTS

We thank the anonymous reviewer for the constructive comments that have improved the quality of this manuscript. JD acknowledges funding support from the Chinese Scholarship Council (no. 202008610218). PJC and ZML acknowledge financial support from the Science and Technology Facilities Council (grant number ST/V000454/1). The giant impact simulations were carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol (http://www.bristol.ac.uk/acrc/) and Isambard 2 UK National Tier-2 HPC Service (http://gw4.ac.uk/isambard/) operated by GW4 and the UK Met Office, and funded by the Engineering and Physical Sciences Research Council (EP/T022078/1).

DATA AVAILABILITY

Full simulation output is available from the authors on reasonable request.

Footnotes

1

https://swift.strw.leidenuniv.nl/docs/index.html

2

The line data can be found at https://github.com/JingyaoDOU/Forming_S_Mercuries_2024.

APPENDIX A: FITTING RESULTS

Below, we show the fitted scaling laws for head-on impacts and oblique impacts.

Figure A1.

Demonstration of the fitting process for equation (3) to determine the target catastrophic disruption criterion |$Q^{^{\prime }*}_\mathrm{TD}$| for the target with mass 1 M_⊕ at various impact parameters. The star symbols represent the fitted Q_R when M_lr/M_targ is 0.5 at different impact parameters.

Figure A2.

Comparison of scaling laws of mass (top panel) and iron mass fraction (bottom panel) of the largest post-collision remnant. In the top panel, the pink dashed line shows the universal law first derived in Stewart & Leinhardt (2009). The blue dashed line shows the fitting result from Denman et al. (2020) at target mass around 10.5 M_⊕. The green dashed line shows the scaling laws from Reinhardt et al. (2022). In the bottom panel, the pink dashed line shows the iron mass fraction scaling law from Marcus et al. (2009) with the initial iron mass fraction shifted to 30 per cent.

APPENDIX B: OTHER SUPPLEMENTARY PLOTS

Figure B1.

Schematic of the collision geometry. θ is the impact angle and the impact parameter is b = sin θ.

APPENDIX C: PRESCRIPTION TO PREDICT MASS AND IRON MASS FRACTION

Below, we describe the prescription that can be used to predict the mass and iron mass fraction of the largest remnant based on the scaling laws we derived in this work. The prescription should only apply to equal-mass single giant impacts with a target mass above 0.01 M_⊕. The prescription for an impact at impact velocity v_imp and impact parameter b is as follows:

• (i) For a target and an impactor planet with masses and radii M_targ, M_imp, R_targ, and R_imp, first calculate the mutual escape velocity of the system:

  $$\begin{eqnarray} V_\mathrm{esc}=\sqrt{\frac{2G(M_\mathrm{targ}+M_\mathrm{imp})}{R_\mathrm{targ}+R_\mathrm{imp}}}. \end{eqnarray}$$

  (C1)
• (ii) Next, calculate the specific impact energy:

  $$\begin{eqnarray} Q_\mathrm{R}=0.5\mu \frac{V^2_\mathrm{i}}{M_\mathrm{tot}}, \end{eqnarray}$$

  (C2)

  where μ = M_targM_imp/M_tot.
• (iii.a) Head-on impact (b = 0): Calculate catastrophic disruption criteria for head-on impacts:

  $$\begin{eqnarray} Q^{*}_{\mathrm{RD},\,\gamma =1} = c^{*}\frac{4}{5}\pi \rho _{1}GR^2_\mathrm{C1}, \end{eqnarray}$$

  (C3)

  where c* equals 2.926, ρ₁ = 1000 kg m⁻³, and R_C1 = |$(\frac{3 M_\mathrm{{tot}}}{4\pi \rho _1})^{\frac{1}{3}}$|⁠.

• (1) If a target mass is less than 2 M_⊕, use equations (4) and (5) to predict the mass and iron mass fraction, respectively. If Q_R/|$Q^{*}_\mathrm{RD}$| is larger than 1.2, label the impact as fragmentation, indicating that no main remnant is left in the system.

• (2) If a target mass is larger than 2 M_⊕, use equations (6) and (7) to predict the mass and iron mass fraction, respectively.

• (iii.b) Oblique impact (b ≥ 0.2): Calculate the target catastrophic disruption criteria:

  $$\begin{eqnarray} Q^{^{\prime }*}_\mathrm{TD} = \left(0.175 \left(\frac{1}{1-b} \right)^{3.3}+0.635\right)c^{*} \frac{4}{5}\pi \rho _{1}GR^2_\mathrm{C1}, \end{eqnarray}$$

  (C4)

  where b is impact parameter and c* equals 2.926.

1. Find the impact velocity where hit-and-run starts to happen:

   $$\begin{eqnarray} V_\mathrm{HnR} = \left(\Gamma (1-b)^{3.6}+1.2\right)V_\mathrm{esc}, \end{eqnarray}$$

   (C5)

   where Γ = −0.3 log₁₀ M_{targ, ⊕} + 2.18, and M_{targ, ⊕} is the mass of target planet in Earth masses.

2. If V_i is less than V_HnR, then an impact should be treated as perfect merging.

3. Else, if V_i is larger than V_HnR and impact parameter b ≥ 0.3, equations (8) and (9) should be used to predict the mass and iron mass fraction. The parameters α_{M, b}, α_{Fe, b}, and β_{Fe, b} in equations (8) and (9) can be calculated by equations (10), (11), and (12).

4. If b is between 0.2 and 0.3, equations (13) and (9) (with α_{Fe, 0.2} = 0.199, and β_{Fe, 0.2} = 1.904) should be used to approximate the mass and iron mass fraction of the largest remnant.

5. If b is between 0.0 and 0.2, then an impact should be approximated as head-on impact and using prescription described above for head-on impact.