Dynamical Avenues for Mercury's Origin. I. The Lone Survivor of a Primordial Generation of Short-period Protoplanets

Both the isotopic record derived from meteorites (Budde et al. 2016; Kruijer et al. 2017; Brasser & Mojzsis 2020) and the predictions of modern disk models (e.g., Morbidelli & Nesvorny 2012; Levison et al. 2015; Johansen et al. 2015; Dra̧żkowska et al. 2016; Lambrechts et al. 2019; Bitsch et al. 2019) rather strongly disfavor super-Earth formation having been achieved in the inner parts of the solar nebula. This staunchly contradicts the original model of Volk & Gladman (2015). Moreover, it is difficult to reconcile the high inferred occurrence frequencies of super-Earths (∼50% or higher around Sun-like stars; Chiang & Laughlin 2013; Zhu et al. 2018) with a scenario where such a chain of planets in the solar system is cataclysmically ground to dust. Indeed, the orbital period ratios of such multiplanet systems seem to imply that their architectures were reshaped by delayed epochs of giant impacts and instability (Pu & Wu 2015; Izidoro et al. 2017, 2019). Thus, it is unclear why, unlike the vast majority of known systems, a hypothetical population of super-Earths in the solar system would be catastrophically demolished. Indeed, recent work in Poon et al. (2020) using a collisional fragmentation algorithm concluded that the formative episodes of giant impacts and instabilities within chains of hot super-Earths typically do not yield collisional velocities in the cataclysmically destructive regime (e.g., Kokubo & Genda 2010; Genda et al. 2012; Leinhardt & Stewart 2012; Stewart & Leinhardt 2012; Gabriel et al. 2020). For these reasons, we intentionally deviate from the original proposition by Volk & Gladman (2015) and remind the reader that our work does not invoke super-Earth formation in the solar system or exotic dynamical processes (e.g., planetesimal shepherding) to generate our proposed initial conditions (see additional discussion in Raymond et al. 2016; Deienno et al. 2020; Lenz et al. 2020; Bean et al. 2021).

Figure 1 summarizes the systematic inability of dynamical models to generate Mercury–Venus systems. For simplicity, we define a "Mercury–Venus" system as any numerical simulation finishing with exactly two planets inside of 0.85 au with M_M < M_V. The black points illustrate a set of 100 simulations from Clement et al. (2018) assuming the so-called classic initial conditions (∼5 M_⊕ of terrestrial-planet-forming material distributed uniformly between 0.5 and 4.0 au, e.g., Chambers & Wetherill 1998; Raymond et al. 2009). Here, Mercury forms directly within the massive disk in the same manner as Earth and Venus. The gray points depict a batch of 100 simulations from Clement et al. (2019b) assuming the same classic initial conditions and also including a prescription (Chambers 2013) accounting for the effects of collisional fragmentation. In successful realizations, Mercury forms small as the result of a series of imperfect accretion events. The gold points represent a suite of 16 simulations reported in Chambers (2001) where the Mercury-forming region is modeled with a linearly increasing surface density of embryos and planetesimals between 0.3–0.7 au (hence the disk is "peaked"). Successful outcomes occur when Mercury forms directly from within this mass-depleted inner disk component. The tan points plot 40 Mercury–Venus analog systems formed in 89 simulations in Lykawka & Ito (2019) considering similar inner disk components. The work of Izidoro et al. (2015), comprising 60 numerical simulations where the slope of the disk surface density profile is varied, is denoted by blue points. Successful Mercury analogs are produced when an embryo is scattered out of the disk. The light green points illustrate a sample of 71 Mercury–Venus systems formed in 118 total simulations of terrestrial-planet formation in the Grand Tack scenario (Walsh et al. 2011) from Jacobson & Morbidelli (2014). As in the depleted asteroid belt framework of Izidoro et al. (2015), successful Mercury analogs are typically embryos that scatter and dynamically decouple from the disk. The green points plot a set of 125 simulations in Clement et al. (2019b), including a fragmentation prescription where the terrestrial planets form out of a narrow annulus of material (∼2 M_⊕) between 0.7 and 1.0 au as envisioned in Hansen (2009). Successful Mercury analogs are occasionally generated when an embryo is scattered out of the annulus. Conversely, the magenta points depict a subset of Mercury–Venus systems generated in 800 simulations from Clement et al. (2018) where the giant-planet instability (Tsiganis et al. 2005; Nesvorný & Morbidelli 2012) occurs within the first 10 Myr of the process terrestrial-planet formation. Appropriate Mercury analogs are produced when the giant planets' resonant perturbations liberate an embryo from the planet-forming disk and strand it on a Mercury-like orbit. Finally, the pink points represent a set of 600 simulations from Clement et al. (2019b) modeling the same early instability scenario and also incorporating the effects of collisional fragmentation.

Given the trends depicted in Figure 1, we begin our study with a simple question that will, unfortunately, remain unanswered in our manuscript: what was the structure of the terrestrial disk in the Mercury-forming region? Perhaps the disk conditions that gave rise to Mercury are still completely unknown (e.g., Bean et al. 2021), and the planets' origin will remain a mystery until future exploration of the surface shines a light on the nature of the Mercury-forming planetesimals, and detailed disk observations and models conclusively reveal the elusive initial conditions that generate Mercury-like planets. It is also possible that the standard setups employed by contemporary formation studies are qualitatively correct (for example, Clement et al. 2019b; Deienno et al. 2019; Lykawka & Ito 2019; Walsh & Levison 2019), and the apparent inability of such models to reproduce Mercury (Figure 1) is the result of a flaw in methodology or approach (i.e., insufficient resolution, inappropriate collisional treatments, etc.). Regardless of which proposition (or both) is correct, we argue that the scope of the Mercury problem permits some margin for experimentation and exploration outside of the conventional bulwark of solar system evolutionary narratives.

In this paper, we investigate a scenario where a quasi-stable chain of four to six large, well-spaced embryos emerge in the region interior to the terrestrial-forming disk. While these initial conditions might seem largely contrived, we believe that they are not particularly fictitious given the diversity of disk structures generated in hydrodynamical and N-body models. For instance, Dra̧żkowska et al. (2016) developed a semianalytical disk model to follow the evolution of dust and pebbles in the Sun's gaseous nebula and found that planetesimals form and pile up in a narrow annulus of width ∼0.3–3.0 au beginning around a ∼ 0.3–1.0 au. Thus, it is important to refrain from grounding any investigation in a narrow assumption that planetesimal growth and, by extension, planet formation transpired identically throughout the terrestrial-forming disk. This is particularly relevant given the tendency of disk models to find planetesimal formation to be highly localized and sensitive to variations in initial conditions (e.g., Levison et al. 2015; Dra̧żkowska et al. 2016; Yang et al. 2017; Dra̧żkowska & Dullemond 2018; Lenz et al. 2019, 2020). Additionally, the conversion of the planetesimals that do form into large embryos (similar to those studied in this paper; Kokubo & Ida 1996, 2000; Chambers 2006) occurs with similarly variable efficiencies at different locations in the disk (Morishima et al. 2008; Walsh & Levison 2019; Wallace & Quinn 2019). In particular, high-resolution N-body simulations of this epoch in Clement et al. (2020a) found the innermost regions of the terrestrial-forming disk to be highly processed after the runaway growth phase and virtually devoid of planetesimals near the Mercury-forming region (bins of ≃0.45–0.6 au in that work). Moreover, it seems plausible to envision a population of embryos traversing the disk under the influence of Type I migration becoming stranded in the Mercury-forming region during the disk's photoevaporation phase. Indeed, new work in Woo et al. (2021) considering a variety of conceivable disk parameters and giant-planet configurations reveal similar chains of dynamically isolated embryos (although their dynamical separations are not as broad in terms of mutual Hill radii as presumed in our simulations; see Figures 2 and 12 in that work).

Zoom In Zoom Out Reset image size

We select a modified version of the Mercury6 hybrid integrator (Chambers 1999) described in detail in Chambers (2013) for our numerical simulations. This augmented package includes algorithms designed to incorporate the various fragmentation regimes mapped in Leinhardt & Stewart (2012) and Stewart & Leinhardt (2012). In short, when an erosive collision is detected, the code calculates the mass of the largest remaining remnant and divides the remaining material into a number of equal-mass fragments that are ejected at v ≃ 1.05v_esc in uniformly spaced directions in the collisional plane. While this selection of ejection velocity is not precisely physically motivated, it represents a reasonable compromise between the premature reaccretion of displaced material and the creation of excessively dynamically excited fragments (see Clement et al. 2019a for a recent study that tested different ejection velocities). The total number of fragments is determined by the user-established minimum fragment mass (MFM), which we set to 10% of Mercury's modern mass (0.0055 M_⊕) in our simulations. Establishing an MFM mass is necessary to avoid injecting an unreasonable number of particles in a simulation such that the calculation is no longer computationally feasible.

As the parameter space of possible primordial configurations of short-period protoplanets is extensive, we limit our current study to systems of equal-mass objects in the vicinity of Mercury's modern orbit. In four independent batches of integrations (Table 1), we investigate systems of M_emb = 0.3, 0.1, 0.05, and 0.025 M_⊕ embryos, distributed approximately uniformly between 0.25 and 0.6 au (we impose a random variation of δ_a ≲ 0.02 au on each semimajor axis to generate different initial conditions). Eccentricities and inclinations for the embryos are selected at random from nearly circular, coplanar Rayleigh distributions (σ_e = 0.01, σ_i = 01), and the remaining angular orbital elements are assigned randomly by sampling from uniform sets of angles. While our initial conditions are clearly somewhat contrived in the absence of high-resolution dynamical investigations of runaway growth at small radial distances (e.g., Section 2.1), our higher values of M_emb (0.1 and 0.3 M_⊕) are loosely based off recent high-N studies of planetesimal accretion in the gas-disk phase (Walsh & Levison 2019; Clement et al. 2020a; Woo et al. 2021). In contrast, our lower presumed values of M_emb (0.05 and 0.025 M_⊕) are selected to boost the probability of forming Mercury analogs with the correct mass (i.e., they are not particularly physically motivated beyond the fact that they are similar to the actual mass of Mercury). It is also important to note that, in all cases, our embryos begin the simulation separated from one another by ∼30–50 mutual Hill radii (R_h ). Thus, our initial conditions are selected such that the system of short-period planets should possess a relatively high degree of Hill stability (e.g., Chambers et al. 1996; Kokubo & Ida 1998; Petit et al. 2018) in the absence of perturbations from the other planets.

Table 1. Summary of Initial Conditions for Our Various Simulations

Set	Nemb	Memb (M⊕)	MVenus (M⊕)	Nsim	${N}_{\mathrm{instb}}/{N}_{\mathrm{sim}}$ (%)
1	3	0.3	0.6	200	31
2	4	0.1	0.6	200	88
3	5	0.05	0.815	200	100
4	6	0.025	0.815	200	100

Note. The columns are as follows: (1) the simulation set, (2) the total number of embryos placed interior to Venus' orbit, (3) the masses of the individual embryos, (4) the initial mass of Venus, (5) the total number of simulations, and (6) the percentage of all simulations that experience an instability.

Download table as:  ASCIITypeset image

In all of our simulations, we include Venus, Earth, Mars, Jupiter, and Saturn in their modern orbits. As Venus typically accretes at least one of the embryos in our simulations (Section 3.4), we reduce its initial mass to M_Venus = 0.6 M_⊕ in our two sets, testing the highest values of M_emb in order to boost the probability that Venus finishes with the correct mass. Each of our systems is integrated for 100 Myr utilizing a 1.5 day time step and algorithms designed to account for the effects of relativity.

The fraction of simulations that experience instability resulting in the loss of one or more protoplanets is given in Table 1. For simplicity, we define the "instability time" as the first collisional encounter involving an embryo. Instabilities typically transpire within a few tens of Myr in simulations evaluating lower values of M_emb (Figure 2) and occur at lower frequencies and after longer delays in our integrations investigating M_emb = 0.3 and 0.1 M_⊕. As our various configurations of protoplanet masses begin our integrations with roughly equivalent mutual Hill spacings (∼30–50 R_h ), the lengthy instability delays in our M_emb = 0.3 M_⊕ batch are somewhat conspicuous (e.g., Chambers et al. 1996). We find that this higher degree of stability in our more massive embryo configurations is a result of the larger protoplanets being initialized on more widely separated orbits. We expound further on this discrepancy in Section 3.1.

We evaluate the compositional evolution of the objects in each simulation using the methodology of Chambers (2013) for a first-order approximation of the final Mercury analogs' Fe/Si ratios. We begin by assuming that each inner planet's initial composition is approximately chondritic and fully differentiated such that 30% of the total mass is in an iron-rich core and the silicate mantle makes up the remaining 70%. When an erosive collision occurs between two bodies in the simulation, fragments are first generated from the mantle of the projectile, then its core, then the mantle of the target body, and finally, from the target's (largest remnant) core material.

An alternative methodology for assembling the final planets' bulk compositions would be to assign each fragment with an equal mixture of the cumulative ejected material. One could imagine this as a crude approximation of the fragments themselves having reaccreted from a sea of vaporized material. However, we experimented with such an implementation and determined that it did not significantly alter the final core mass fraction (CMF) distributions in the majority of our systems. This is because collisional fragments are rarely involved in prolonged, subsequent episodes of CMF-altering collisions. While such fragments do occasionally experience repeated follow-on hit-and-run events with the larger remnant particle, such interactions almost exclusively result in the eventual accretion of the fragment with no net transfer of mantle or core material outside of the interacting system. We elaborate further on whether or not our fragmentation algorithm is an accurate physical representation of CMF alteration during erosive collisions in Section 3.3.

The pitfall of our chosen CMF-tracking implementation is the obvious degeneracy that arises when both the mantle and core of an object are eroded. In this case, the various fragments are essentially randomly assigned CMFs of either 0.0 or 1.0, with one fragment typically possessing a mixture of the two to ensure mass convergence. If one of these fragments survives the simulation and becomes a Mercury analog, the simulation's success or failure in terms of Mercury's mean density would be determined by what material was randomly assigned to the individual fragment that went on to become Mercury. While we find that this situation rarely arises as the majority of fragments are reaccreted by either the remnant particle or another embryo in the system (and seldom become Mercury analogs themselves), we still focus the majority of our analysis in the subsequent sections on our various simulation sets' collective statistical properties. However, we also discuss individual systems that are excellent solar system analogs as examples of how our proposed scenario might have transpired.

The region of the modern solar system interior to ≃0.65 au is quite complex dynamically (outside of 0.65 au, in the vicinity of Earth and Venus, NEA (near-Earth asteroid) dynamics are primarily governed by Kozai cycles; Michel & Thomas 1996). Overlapping mean motion resonance (MMR) and secular resonance in the region results in large swaths of parameter space where our hypothetical protoplanets' orbits evolve chaotically (e.g., Figure 3, see also Michel & Froeschlé 1997; Laskar & Gastineau 2009; Batygin et al. 2015). Most significant among these various unstable zones are objects with semimajor axes near the intersection of the 3:1 MMR with Venus and the powerful ν₅ secular resonance at ∼0.35 au (leftmost blue regime in Figure 3) which are highly unstable and can easily be excited in eccentricity to the point of loss or collision with another protoplanet in relatively short timescales. Orbits around a = 0.46 and a = 0.55 au are also unstable in our simulations (in that objects attain Venus-crossing orbits within 10 Myr)—a consequence of overlaps between ν₂, ν₁₂, and the 2:1 MMR with Venus, and the ν₁₆ ν₆, ν₃, ν₄, and 3:2 MMR, respectively.

Zoom In Zoom Out Reset image size

Figure 4 demonstrates how perturbations from Jupiter's g₅ mode can easily excite objects with masses akin to those of our embryos on to crossing orbits (see also Batygin et al. 2015). The plot depicts a 10 Myr integration of the solar system (utilizing the Mercury6 Hybrid integrator and a 1.5 day time step; Chambers 1999), where Mercury's semimajor axis is shifted to the unstable region in Figure 3 at a = 0.34 au. It is important to recognize that the secular architecture of this system (or those of any of the other integrations presented in this manuscript) is not necessarily akin to that of the solar system. For instance, Mercury's nodal precession rate induces a significant perturbation on Venus' inclination evolution in the modern solar system (e.g., Laskar 1990) and thus the location and precise dynamics of the ν₁₂ resonance is clearly different in our simulations than in the solar system by virtue of moving Mercury or replacing it with a system of short-period protoplanets. However, it is clear from Figures 3 and 4 that the influence of g₅ in the region interior to Mercury's modern orbit is quite significant in terms of its ability to excite the orbits of Moon–Mars-mass objects. It is also important to note that this analysis only considers the giant planets on their modern orbits. If our hypothetical epoch of short-period planet destabilization transpired in conjunction with the giant-planet instability (Tsiganis et al. 2005; Nesvorný & Morbidelli 2012; Deienno et al. 2017; discussed further in Sections 3.2 and 4), it is plausible that the inward sweeping of ν₅ toward its ultimate modern location might serve to destabilize our chains of embryos (e.g., Brasser et al. 2009; Roig et al. 2016).

Zoom In Zoom Out Reset image size

While it is beyond the scope of our present manuscript to fully scrutinize each individual instability and determine the precise destabilization mechanism, it is obvious that instabilities are fairly ubiquitous in our simulations (particularly for the lower-mass systems of embryos; Table 1). Moreover, from a simple analysis of the encounter histories of our systems of protoplanets, it is readily apparent that the highly perturbed region of orbital parameter space in the vicinity of the ν₅ secular resonance plays a key role in driving our configurations toward an instability. Indeed, ≳45% of the first collisions in our simulations involve the embryo that originates closest to the 3:1 resonance with Venus. Another ∼30% of our systems' instabilities unfold after the outermost embryo's eccentricity is gradually excited until the object inhabits a crossing orbit with Venus (at least partially as a consequence of being in close proximity to the chaotic region in Figure 3 around the 3:2 MMR at 0.55 au). Additionally, we used a set of preliminary runs to determine that placing the central embryo in our M_emb = 0.3 M_⊕ sets (the most resilient to destabilization; Table 1) at a semimajor axis roughly corresponding to the 3:1 MMR with Venus boosts the total fraction of systems experiencing an instability by ∼50%.

The highly chaotic orbital evolution of protoplanets in the radial regime interior to Mercury's modern semimajor axis (e.g., Figure 4) is indeed somewhat conspicuous and remarkably consistent with the hypothesis that the solar system once possessed additional planet-mass bodies interior to Venus (Volk & Gladman 2015). The modern Mercury–Venus period ratio is peculiar in that it is quite larger than that of the other neighboring pairs of terrestrial planets and is seldom matched in planet formation models (Clement et al. 2019a). Furthermore, Mercury is positioned at the edge of one of the more regularly behaved islands in a/i space (Figure 3) between Venus and the ν₅ resonance. It may be impossible to ascertain whether or not additional planets formed inside of Venus' orbit. However, it is easy to imagine a simplified version of our scenario, where an additional planet once existed between Mercury and Venus in the vicinity of the 3:2 MMR, ν₁₆, ν₃ overlap, and was destabilized and lost via a collision with Venus. Further additional planets interior to Mercury might also have existed and been eventually lost via dynamical interactions with the powerful ν₅ resonance as discussed above. Such hypothetical short-period planets might present an intriguing explanation for Mercury's peculiar composition if they underwent a series of erosive impacts with proto-Mercury en route to being lost. To ascertain the likelihood of this hypothetical scenario, we turn to our suites of full dynamical simulations (Table 1).

Figure 5 plots the final orbits of Mercury analogs in systems that finished with the correct number of planets (referred to throughout the subsequent text as one-Mercury systems) for each of our different simulation batches. It is clear that, regardless of initial conditions, Mercury's modern eccentricity and inclination are rather typical outcomes of our modeled scenario. The cumulative distributions of eccentricities and inclinations plotted in Figures 7 and 8 demonstrate how our simulations do tend to slightly underexcite Mercury's orbit, though it is clear that our scenario yields a wide spectrum of outcomes. It is still unclear whether Mercury's modern dynamical excitation is largely a consequence of the giant-planet instability or the terrestrial-planet formation process itself (e.g., Kaib & Chambers 2016; Roig et al. 2016; Clement et al. 2019b). Our work presumes that the so-called Nice Model instability (Tsiganis et al. 2005; Nesvorný & Morbidelli 2012; Deienno et al. 2017) occurred prior to the loss of our hypothetical population of short-period planets (consistent with an early instability as argued for by multiple recent studies; Clement et al. 2018; Morbidelli et al. 2018; Nesvorný et al. 2018; Brasser et al. 2020). While we leave the full study of our scenario's timing within the various proposed timelines for the terrestrial planets' early evolution to future work (discussed further in Section 4; see recent reviews in Izidoro & Raymond 2018; Raymond et al. 2018), we note that the Nice Model instability possibly played a role in sculpting Mercury's modern orbit. Thus, while many of our final systems are remarkable solar system analogs, our scenario is by no means a comprehensive model for the early evolution of the inner solar system.

Zoom In Zoom Out Reset image size

We compare the final mass and orbital period ratios of our Mercury–Venus analog systems to those of previous dynamical studies (Figure 1) in Figure 6. As in Figure 5, it is clear that the final semimajor axes of Mercury analogs in one-Mercury systems decrease with decreasing initial protoplanet masses. In the most extreme case (M_emb = 0.025 M_⊕), Mercury's modern semimajor axis lies outside of the range of values generated in our simulations. In contrast, our higher-mass configurations of embryos (Table 1) are plagued by a low rate of instabilities and preponderance of overly massive Mercury analogs (e.g., top-left panel of Figure 5). Thus, it is apparent that our configurations assuming more moderate initial embryo masses (∼0.05–0.10 M_⊕) are the most viable in terms of their ability to simultaneously match both Mercury's mass and orbit. The magnitude of the eccentricity excitation induced on our protoplanets via secular resonant perturbations (Section 3.1) from the other planets varies inversely with the mass of the primordial embryo (e.g., Murray & Dermott 1999). This is also evidenced by the fact that Mercury analogs produced in systems of less-massive protoplanets possess hotter e and i distributions (e.g., Figure 7). Thus, our M_emb = 0.025 M_⊕ embryos attain higher eccentricities during the instability and are therefore more easily accreted by the other terrestrial planets. Indeed, Venus and Earth accrete approximately three times as many embryos in our simulations investigating lower embryo masses than the higher-mass runs. Thus, embryos in the lower-mass cases that begin interior to the ν₅ resonance tend to merge and survive on an orbit that is closer to the Sun than in the actual solar system. The embryos beginning these simulations exterior to ν₅ are largely lost via collisions with Earth or Venus (discussed further in Section 3.4). The initial semimajor axes of the innermost one to two embryos also play a minor role in producing the trend of final Mercury analogs surviving closer to the Sun in systems with less-massive protoplanets. The inner planets in these systems originate further from ν₅. Thus, they are more likely to begin the integration with orbits that inhabit the more stable regions depicted in Figure 3. However, on closer inspection, we find that the innermost embryos never survive an instability without experiencing at least one collision with either a fragment or another embryo, thus we conclude that this is a minor effect.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast to our less-massive systems, Mercury analogs produced from our M_emb = 0.3 M_⊕ runs inhabit a wider range of final semimajor axes. Often in these cases either all three embryos survive or the protoplanets combine with one another to form a single Mercury analog that is far too massive. Finally, our M_emb = 0.05–0.10 M_⊕ sets occasionally finish with small planets exterior to ν₅ (a ≳ 0.35) like Mercury as the protoplanets are able to avoid collisions with Venus by attaining smaller eccentricities. Moreover, the final analog planet mass in these simulations is closer to that of Mercury by virtue of the system beginning with less total mass in the region.

An example evolution of a successful simulation from our M_emb = 0.05 M_⊕ batch is plotted in Figure 9. The instability ensues at t = 60 Myr when the second, third, and fourth embryos (in order of increasing distance from the Sun) undergo a series of six hit-and-run collisions. Through this process, the most distant embryo is excited in eccentricity until it ultimately experiences a hit-and-run collision with Venus at t = 63 Myr that produces two fragments. All three remnant objects from this interaction are re-accreted by Venus within the subsequent ∼1.5 Myr, boosting the planets' mass by 0.05 M_⊕. Three more fragments are generated in a hit-and-run impact between the second and third embryo at t = 66 Myr. Swiftly following the injection of these additional particles, the system transforms rapidly through a series of ∼30 collisional interactions that produce an additional 12 fragments. During this chaotic period of evolution, the embryo that goes on to become Mercury (second embryo) experiences a series of several imperfect accretion events that cumulatively boost its CMF. In particular, a fragmenting collision between the ultimate Mercury analog and the systems' fourth embryo (the mantle of which had already been significantly eroded in the t = 66 Myr series of impacts) at t = 68 Myr produces five new, mantle-only fragments. Through this series of interactions, the CMF of the eventual Mercury analog is permanently altered to ∼0.60. The mantle-rich fragments ultimately merge with the innermost embryo, which collides with the Sun at t = 78 Myr. It is also apparent from Figure 9 that the orbits of Earth and Venus remain mostly quiescent during the instability. We integrated the system for an additional 1 Gyr to confirm its dynamical stability.

Zoom In Zoom Out Reset image size

We now turn our attention to our Mercury analogs' collisional histories. As our fragmentation algorithm (Chambers 2013) and corresponding CMF-tracking methodology (Section 2) are based on the collisional regimes mapped in Leinhardt & Stewart (2012), objects that emerge from our simulations with highly processed CMFs axiomatically experience one or more energetic, mantle-stripping collisional encounters (Benz et al. 1988; Chau et al. 2018; Jackson et al. 2018). Nevertheless, it is worthwhile to discuss the collisional environment within our computations in the context of hydrodynamical studies specifically devised to transform proto-Mercury's bulk density. As discussed in the introduction, one possibility is that Mercury was the projectile striking a much larger body (either proto-Venus, proto-Earth, or an object that no longer resides in the solar system; Asphaug & Reufer 2014). We analyzed the dynamical problems with this scenario in detail in Clement et al. (2019a). A second plausible circumstance is one where proto-Mercury (near its current orbit) is struck by a smaller body. Benz et al. (2007) analyzed this possibility and found that the collisional velocities required to erode the mantle of a 0.12 M_⊕ (2.25 times Mercury's modern mass) version of Mercury are typically high (≳6v_esc; diamond in Figure 10). However, the authors obtained one moderately successful case from a nearly head-on impact at ∼4v_esc. More recently, Chau et al. (2018) thoroughly investigated the applicable collisional parameter space in the classic Benz et al. (2007) scenario and favored a lower-velocity (∼4v_esc) collision and less-extreme impact geometry (impact parameter, b ≃ 0.4). Additionally, the authors explored the possibility that Mercury's mantle was eroded in a series of similar collisions and tabulated the number of impacts required for given combinations of ${M}_{\mathrm{proj}}/{M}_{\mathrm{targ}}$, b, and v_imp/v_esc. Figure 10 plots the relative frequency of impact velocities and ${M}_{\mathrm{proj}}/{M}_{\mathrm{targ}}$ for all collisions on protoplanets in our various simulation batches, compared with a range of successful scenarios from Benz et al. (2007) and Chau et al. (2018). While the typical target particle masses in our simulations investigating M_emb = 0.3 and 0.025 M_⊕ are different than what is invoked in the classic single-erosive impact model (a fully differentiated body with ∼2.25 times Mercury's modern mass), those in our more dynamically successful batches testing M_emb = 0.1 and 0.055 M_⊕ are quite similar. Interestingly, the most common combinations of ${M}_{\mathrm{proj}}/{M}_{\mathrm{targ}}$ and v_imp/v_esc in these sets (blue regions in Figure 10) necessitate a head-on collision (black square) to adequately reshape Mercury's mantle. However, a multiple-impact scenario with a less-extreme orientation (upside-down triangle) is certainly plausible in our model. In such a scenario, Mercury would experience multiple erosive encounters with other protoplanets and collisional debris in the region.

Zoom In Zoom Out Reset image size

To evaluate our scenario's success in terms of its ability to collisionally strip Mercury's mantle, we plot the final CMFs of all surviving Mercury analogs in our simulations in Figure 11. It is clear that, with the exception of our M_emb = 0.3 M_⊕ set where instabilities occur quite infrequently, our systems' populations of Mercury analogs are highly processed and altered in CMF via the collisional interactions that occur during the instability. Indeed, only ∼20% of surviving Mercury analogs in our three sets testing lower values of M_emb retain their original CMFs. Thus, the one-Mercury systems depicted in Figure 5 are quite diverse in terms of Mercury's CMF. However, we remind the reader that these results are partially a consequence of our code's fragmentation algorithm (Section 2). In particular, when a fragmenting collision occurs between two lower-mass embryos, the resulting collisional fragments have masses closer to those of the embryos themselves as we utilize the same MFM (0.0055 M_⊕) in all of our simulations. Follow-on collisions in these systems can still be erosive and thus still CMF-altering over a wider range of collisional velocities (e.g., Genda et al. 2012; Leinhardt & Stewart 2012; Gabriel et al. 2020).

Zoom In Zoom Out Reset image size

At this point, it is reasonable to question whether our fragmentation model provides an accurate representation of the problem we seek to investigate. Hydrodynamical models tend to find that the majority of the disrupted material in giant impacts similar to those studied in this manuscript should be ejected in the form of dust. Indeed, observed debris disks (e.g., Lisse et al. 2009; Weinberger et al. 2011) interpreted to be potential tracers of the giant-impact process of terrestrial-planet formation (Genda et al. 2015) are largely consistent with this assumption. Thus, while a prescription for tracking fragmentation events within our N-body code (Leinhardt & Stewart 2012; Chambers 2013) is undoubtedly necessary for modeling mantle loss on our Mercury analogs, it is at best unclear whether tracking the ejecta as massive (∼lunar-mass), dynamically interacting bodies is physical. However, we note that, in all four of our simulation batches, ≳90% of produced fragments are subsequently removed via a perfectly accretionary collision without ever undergoing subsequent fragmentation events. In this manner, the vast majority of our high-CMF Mercury analogs acquire their large Fe/Si ratios through a small number of fragmenting collisions that tend to occur with low impact parameters (similar to the head-on scenario of Benz et al. 2007; see Figure 12). The mantle-rich fragments produced in these impacts dynamically interact with the other objects in the simulation until they are eventually accreted by another planet. If this ejected material were in the form of dust rather than massive bodies, it is reasonable to argue that it, too, might be carried away from an eventual Mercury analog via PR drag (Melis et al. 2012) or interactions with the solar wind (Spalding & Adams 2020) rather than being reaccreted (e.g., Gladman & Coffey 2009). If this were the case, our simulations should systematically overestimate both the masses and underestimate the CMFs of Mercury analogs (as disrupted material would be reaccreted less frequently). Thus, we contend that our CMF results are not an artifact of our computation methodology as repeated collisions between fragments do not appreciably contribute to the final CMFs of our Mercury analogs.

Zoom In Zoom Out Reset image size

While our simulations are quite successful at generating a diverse range of Mercury analog Fe/Si compositions (Figures 5 and 11), as discussed in Section 3.2, the surviving analogs themselves tend to be positioned on the sunward side of the ν₅ resonance. However, our integrations begin with the giant planets on their modern orbits (i.e., after the giant-planet instability; Tsiganis et al. 2005; Nesvorný & Morbidelli 2012). Thus, if our formation scenario occurred either prior to, or in conjunction with, the Nice Model instability (e.g., Clement et al. 2018), the ν₅ resonance would be located farther from the Sun by virtue of Jupiter's increased rate of perihelia precession (induced by stronger mutual interactions with Saturn; Clement et al. 2020b). This might boost the likelihood of Mercury's survival near its modern semimajor axis (0.387 au) exterior to ν₅, provided the diminutive planet could withstand the dynamical perturbations involved with ν₅ migrating through its orbit during the instability (Roig et al. 2016). In this same manner, if our scenario predated the Moon-forming impact (e.g., Touboul et al. 2007; Kleine et al. 2009; Canup 2012), the resonant perturbations in the a < 0.6 au region induced by the terrestrial planets would obviously be different from those modeled in our simulations by virtue of the system possessing additional large terrestrial bodies. Thus, a complete validation of our hypothetical evolutionary sequence must study its viability throughout the various stages of solar system evolution and particularly focus on the consequences of the giant-planet instability.

A significant consequence of our scenario is a delayed accretion event on Venus. While not ubiquitous, these delayed mass additions occur on timescales of some 20–30 Myr (Figure 13) after the onset of the instability (Figure 2). However, as our simulation batches investigating lower values of M_emb litter the inner solar system with significantly more debris in the form of fragments generated in repeated imperfect collisions (Section 3.3), the ultimate merger with Venus can occur ≳60–80 Myr after the onset of the instability. In general, the total mass of additional material accreted by Venus (and, to a lesser extent, Earth) is proportional to the total initial mass of the system's protoplanets. Venus gains an average of 0.13 M_⊕ of material after the instability in our batch testing a configuration of three, M_emb = 0.3 M_⊕ objects. Conversely, the planet is only impacted by an average of 0.043, 0.022, and 0.017 M_⊕ worth of embryos and debris fragments in our systems testing protoplanet masses of 0.1, 0.05, and 0.025 M_⊕, respectively. In addition to these delayed accretion events, Venus also occasionally experiences a small number of hit-and-run collisions (∼35%–42% of all collisions on Venus, depending on the simulation batch) en route to eventually accreting the impacting material. While Venus' ultimate structure would be rather unaffected by this additional delivery of mass (presuming a chondritic composition) if our scenario occurred early in the solar system's history (in conjunction with terrestrial-planet formation; Clement et al. 2018), the fact that Venus does not possess an internally generated magnetic field or natural satellite is potentially at odds with primordial stratification in its core having been disrupted by a late giant impact as was the case on Earth (Jacobson et al. 2017). On closer inspection, it is clear that the vast majority (≳90%) of the material impacting Venus in our simulations begins the simulation outside of Mercury's modern orbit. Thus, it might be easy to prevent late impacts on Venus by altering the range of initial semimajor axes for our protoplanets (i.e., not placing any embryos at a ≳ 0.5 au). However, such a mass configuration is quite ad hoc, and it is unclear how such a wide gap in the surface density profile of our protoplanets would have developed. Thus, we favor a uniform distribution of embryos as our modeled evolutions provide a compelling explanation for the modern Mercury–Venus period ratio.

Zoom In Zoom Out Reset image size

While Venus' mass is somewhat altered in our scenario, the effect on the other terrestrial planets' orbits is rather minimal, though not negligible. Figure 14 plots the fractional change in angular momentum deficit (AMD) for Earth and Venus in our various simulation sets:

$$\begin{eqnarray}&&\mathrm{AMD}=\displaystyle \frac{{\displaystyle \sum }_{i}{M}_{i}\sqrt{{a}_{i}}[1-\sqrt{(1-{e}_{i}^{2})}\cos {i}_{i}]}{{\displaystyle \sum }_{i}{M}_{i}\sqrt{{a}_{i}}}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

The majority of the terrestrial-planet AMDs in our simulations testing higher values of M_emb are altered (almost exclusively increased) by over a factor of 2. In some ∼5% of cases, AMD_EV is increased by a whole order of magnitude. Thus, our results indicate that larger initial values of M_emb are potentially incompatible with the modern, dynamically cold orbits of Earth and Venus. While a small fraction of systems in all four of our simulation batches do experience only minor changes in AMD_EV, a majority of the configurations considering smaller values of M_emb undergo only mild alterations in the orbits of Earth and Venus. Thus, our results indicate that a version of the Volk & Gladman (2015) scenario incorporating smaller protoplanets is perhaps more promising given the improved likelihood of systems experiencing smaller net changes in AMD_EV, and the systematically lower masses of Mercury analogs in systems finishing with the correct number of planets (Figure 5).

It is worthwhile to investigate our systems of protoplanet stability during the chaotic giant-impact epoch of terrestrial-planet formation and similarly understand whether it is reasonable to disentangle our proposed scenario from Earth's and Venus' formation as assumed in Section 3.2. To accomplish this, we perform an additional suite of 50 simulations where we place our systems of protoplanets interior to a terrestrial-forming disk of embryos and planetesimals. For simplicity, we distribute the planetesimals and embryos in a narrow annulus between 0.75 and 1.05 au (e.g., Hansen 2009; Walsh & Levison 2016; Lykawka 2020). While the process of terrestrial-planet formation is rather complex and a subject of ongoing research, we select these initial conditions to maximize the possibility of forming Venus-, Earth-, and Mars-like planets at the proper radial locations and minimize the computational cost of the calculation (for a recent review, see Raymond et al. 2018). Each system considers 2.0 M_⊕ of planet-forming material, distributed such that half the mass is in 12 equal-mass embryos and the remaining half is in 200 equal-mass planetesimals. Orbits for these particles are selected using the methodology described in Clement et al. (2019b). In 25 simulations, we study the evolution of our 4 × 0.1 M_⊕ protoplanet systems, and in an additional 25 runs, we follow the dynamics of our 5 × 0.05 M_⊕ protoplanet initial conditions (Table 1). Each system is integrated for 100 Myr using the same Mercury6 integrator and time step employed throughout our manuscript (Chambers 1999). However, these additional simulations do not incorporate a fragmentation algorithm.

In general, our additional simulations evince a broad spectrum of outcomes. In fact, the final states of four of our additional runs are qualitatively similar to those attained by our embryos integrated within the modern solar system (Section 3.2). Thus, the outer one to two protoplanets incorporate into Venus, and the inner two to three merge in a series of high-energy collisions to form a dynamically isolated Mercury analog. More typically, however, the outer one to two embryos merge with the forming Venus, and the innermost protoplanets survive the 100 Myr integration on relatively unaltered orbits. Indeed, the median number of surviving protoplanets among our 50 terrestrial-planet formation simulations is 3. While we observe a few cases where the outermost protoplanet does not merge with Venus and instead combines with one of the interior protoplanets, the overwhelming number of systems largely retain their initial protoplanet structures, with the exception of objects that merge with Venus. Indeed, one proto-Mercury system of four 0.1 M_⊕ embryos survives the entire process of terrestrial-planet formation intact (aside from each protoplanet accreting material scattered out of the terrestrial-forming disk). The evolution of this successful simulation is plotted in Figure 15.

Zoom In Zoom Out Reset image size

It is far beyond the scope of this paper to comprehensively study the survivability of our Mercury-forming protoplanets within the various proposed terrestrial-formation (e.g., Walsh et al. 2011; Bromley & Kenyon 2017; Raymond & Izidoro 2017; Clement et al. 2018) and giant-planet migration sequences (e.g., Nesvorný & Morbidelli 2012; Pierens et al. 2014; Raymond & Izidoro 2017; de Sousa et al. 2020; Clement et al. 2021). Nevertheless, we present these initial simulations as a proof of concept in order to demonstrate that the dynamical evolution of the embryos in our proposed scenario (particularly the innermost ones that tend to form the final Mercury analog) is not strongly coupled to the process of terrestrial-planet formation.