Dependencies of Mantle Shock Heating in Pairwise Accretion

Collisions between large (⪆300 m), gravity-dominated bodies ("giant impacts") are a ubiquitous feature of nearly all planet formation models; however, the timing, number, and conditions of these impacts can vary considerably. "Classical" accretion (Weidenschilling 1977; Wetherill 1980) involves pairwise accretion of material from submeter-sized bodies to the current planets observed today. Planetesimals mutually accrete in a process called "hierarchical growth" (Tanga et al. 2004) until the establishment of large planetary "embryos." Stirring by large bodies then causes mutual collisions between small bodies to be catastrophic (shutting off their growth); simultaneously, the larger protoplanets sweep up the remaining small bodies in a process termed "oligarchic growth" (Kokubo & Ida 1998; Chambers 2006). Although this model produces some inconsistencies with the current solar system configuration, e.g., overly massive Mars analogs (Chambers 2001; Raymond et al. 2009), it remains an important benchmark (Chambers 2013; Quintana et al. 2016). The "grand-tack" model (Walsh et al. 2011, 2012) features the inward migration of gas giants around 2–5 Ma, which importantly depletes the Mars orbital zone of mass. The terrestrial zone still undergoes pairwise accretion, albeit under an excited state that produces higher impact velocities (O'Brien et al. 2014). The "low-mass asteroid belt" model (Raymond & Izidoro 2017) features a depleted Mars region, potentially provided by earlier, nebular processes, but again, planets in the terrestrial zone would undergo accretion via pairwise interactions. In most planet formation models, however, the formation of Jupiter remains a conundrum; the mutual accretion of kilometer-sized bodies does not form a large enough core (∼10 Earth masses) early enough to attract the waning gas disk to form a sufficiently large envelope (Lambrechts & Johansen 2012). Thus, a new "pebble accretion" paradigm (Lambrechts & Johansen 2012) was developed, whereby planetary embryos form rapidly (within the first 1–10 Ma) through accretion of small pebbles with assistance from gas drag, providing the seeds for gas giants. The overall size distribution of embryos, however, is sensitive to model parameters (Chambers 2016) and thus so is the nature of giant impacts that occur in this model. Considering that giant impacts are often an unavoidable feature of solar system formation models, we can expect these collisions play important roles in exoplanetary systems as well, which include a more exotic array of inferred formation conditions.

The rich diversity of planet formation models have predictable effects on the nature and number of giant impacts potentially with measurable consequences. In "self-stirred" dynamical systems, i.e., where the accreting protoplanets are the dominant gravitational perturbers, collision velocities are just above the mutual escape velocity (freefall speed) of the protoplanets (Wetherill 1980). Thus, impact velocity is proportional to the mass of the protoplanets, causing them to collide at higher absolute velocities as they grow. A natural consequence is that the escape velocity at sufficient scales can exceed the sound speed of mantle silicates, causing unavoidable shocks in even the gentlest collisions (Asphaug et al. 2015). Thus, the planetesimal/embryo size distribution and degree of dynamical stirring can modulate the energy budget and partitioning, which can involve large exchanges between gravitational potential, kinetic, and internal energy (Okeefe & Ahrens 1977; Carter et al. 2020). For example, migrating gas giants would naturally shift velocities to higher multiples of their mutual escape velocity, reducing accretion efficiency through nonmerging ("hit-and-run") collisions that poorly couple impact energy to the target (Asphaug et al. 2006); however, the stirred velocities would also induce shock heating in smaller bodies. A recent pebble accretion model for terrestrial planet formation (Johansen et al. 2021) may mostly avoid an extensive period of pairwise accretion altogether, an important aspect that we address later.

Formation models have a strong lever arm on the nature of giant impacts; however, the role of shocks and their thermodynamic consequences in various giant impacts is relatively understudied. From experimental and theoretical work in classical cratering context, shocks are understood to induce devolatilization (Tyburczy et al. 1986, and other works) and modern experimental approaches continue to produce new insights (Kraus et al. 2013). In light of modern advances in equations of state (e.g., Melosh 2007; Kraus et al. 2015), classical scaling laws for vapor and melt production in the cratering regime (e.g., Okeefe & Ahrens 1977; Pierazzo et al. 1997) continue to undergo revision (Kraus et al. 2015). These revisions also help resolve inaccuracies from the use of simplified equation-of-state models (e.g., Tillotson) in impact shock physics codes, which has been widely employed in giant impact literature despite the lack of explicit treatment of phase boundaries and the nonphysical behavior during cold decompression (Stewart et al. 2020). However, despite equation-of-state advances, the scaling of thermodynamic processes in giant impacts has yet to be established. Using up-to-date thermodynamically consistent equations of state and new shock physics experimental data, we systematically explore the thermal consequences of giant impacts at various scales and velocities in a suite of full 3D hydrocode simulations. This allows us to understand how global-scale thermodynamic processes in colliding protoplanets may play a role in compositional diversity and volatile inventory observed in solar and exosolar systems, as well as help identify consequences of various planet formation paradigms. We also briefly explore numerical effects that are important to consider when making predictions using a common giant impact modeling method.

To determine if colliding materials surpass thermodynamic thresholds for melting and vaporization in giant impacts, we use the "peak entropy" method (Zeldovich & Raizer 1966; Ahrens 1972) as a postprocessing step for our impact simulations described below; however, we largely use these thresholds as metrics for the extent of thermodynamic processing. This method allows for the tracking of melting and vaporization conditions by inspection of the evolution of entropy through shock loading and unloading. The approach is ideal in scenarios where shear strength can be neglected (Quintana et al. 2015), which is the case at giant impact scales where gravitational stresses readily overwhelm material coherence (Asphaug et al. 2006). We used the experimentally derived entropy thresholds from Davies et al. (2020) for complete melting, incipient vaporization, and 50% (by mass) ^{1 , 2} vaporization for forsterite (Mg₂SiO₄) starting at 1200 K: S_CM = 3474(±197) J K⁻¹ kg⁻¹, S_IV = 4270(±279) J K⁻¹ kg⁻¹, and S_50V = 7616 J K⁻¹ kg⁻¹, respectively, for release at P = 10⁵ Pa. We also use the entropy values of S_CM = 3474(±197) J K⁻¹ kg⁻¹, S_IV = 3474(±197) J K⁻¹ kg⁻¹, and S_50V = 6635 J K⁻¹ kg⁻¹, respectively, for release at P = 5.2 Pa. ³ We note that our planets are initialized at ∼1500 K and feature a temperature gradient, which may contribute to reaching entropy thresholds at lower absolute impact velocities than predicted in Davies et al. (2020), as discussed later.

To model collisions, we use the smoothed particle hydrodynamics (SPH) method, which is a Lagrangian code that discretizes modeled planets into mass nodes with certain properties numerically computed using a smoothing function (kernel). We use the standard three-dimensional cubic spline kernel (Monaghan & Lattanzio 1985), a traditional form of artificial viscosity to resolve shocks (Monaghan 1992), and a resolution of 10⁵ nodes in each body. We use the SPHLATCH code (Reufer 2011; Emsenhuber et al. 2018) to perform our SPH simulations. We arrange nodes according to the HEALPix algorithm (Gorski et al. 2005), which was developed to generate discretized spherical maps for the analysis of the cosmic microwave background, but has since been extended to giant impact applications (Reinhardt & Stadel 2017). SPH bodies were simulated in free space for several hours of simulation time to achieve hydrostatic equilibrium (Reufer 2011), before running the collision. We primarily explore head-on geometries of two equal-mass bodies to ensure the entire impactor intersects the target during the collision, mitigating the issue of inefficient coupling of impact energy (Canup 2008) and simplifying the analysis, but we include some off-axis collisions to gain insight into more probable geometries. The masses of our colliding bodies span 10⁻⁵–1 M_⊕, in log-normal increments. Mantle silicates and the cores are initialized at 70 wt% and 30 wt% abundance, respectively, broadly reflecting the bulk silicate-to-iron ratios of differentiated bodies composed of primitive chondritic materials; this is a common assumption in giant impact literature as it represents planets that are not fractionated from primordial composition. We utilize the ANEOS/M-ANEOS equation of state (Melosh 2007) for iron (core materials; Thompson & Lauson 1974) and updated parameters for forsterite (mantle materials; Davies et al. 2020; Stewart et al. 2020).

We simulate collisions between two end-member planet formation scenarios: (1) a self-stirred orbital dynamical system and (2) a system excited by a large gravitational perturber (e.g., a migrating gas giant). In the former, the expected impact velocity is just above the mutual escape velocity (Wetherill 1980):

$$\begin{eqnarray}&&{v}_{\mathrm{esc},{\rm{m}}}=\sqrt{\displaystyle \frac{2{ \mathcal G }({M}_{\mathrm{tar}}+{M}_{\mathrm{imp}})}{{R}_{\mathrm{tar}}+{R}_{\mathrm{imp}}}},\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal G }$ is the gravitational constant, M is the body mass, and R is its radius. In orbital dynamical systems, impact velocities lower than the mutual escape velocity are impossible (without gas drag) and velocities appreciably greater than ∼2–3 v_esc are exceedingly improbable without retrograde orbits and/or strong gravitational stirring from much larger bodies (Wetherill 1980; Ida & Makino 1992; Agnor et al. 1999), the latter of which can shift velocity distributions to several times the mutual escape velocity (O'Brien et al. 2014). Thus, we predominantly simulate impact velocities at ∼1.1 and 3 times the mutual escape velocity: the former is near the median impact velocity in late-stage planet formation simulations (Chambers 2013; Gabriel et al. 2020) and the latter is just above the median impact velocity of Grand-Tack planet formation simulations (∼2.1 v_imp/v_esc; O'Brien et al. 2014). This allows us to understand near end-member, destructive outcomes unique to formation paradigms that involve strong dynamical excitement. We also simulate collisions at 2.36 v_esc as they have the same absolute velocity as a 1.1 v_esc collision between bodies 10 times more massive. We also perform collisions at 7.5 v_esc since they have equal energy to 1.1 v_esc collisions between bodies 10 times more massive. Additionally, we simulate a few collisions at 5 v_esc to fill in our analysis.

Given our focus on thermodynamic transitions, we chose impact scales (total mass) that span escape velocities from subsonic (a few hundred meters per second in our smallest bodies) to supersonic (up to 11 km s⁻¹ escape velocities in our largest bodies). ⁴ Since any material competency (e.g., tensile strength) would be readily overcome in simulations of bodies with R ≫ ∼100–300 m (Benz & Asphaug 1999) and our smallest bodies (M = 10⁻⁵ M_⊕) have radii much greater than this transition (D ∼ 300 km), we assume a fluid rheology. However, Jutzi & Asphaug (2015) showed that a variety of material strength models allow for the transition to purely gravity-dominated collisions to extend to larger sizes (D ∼ 100 km), and this is just below the size of our smallest bodies. Real collisions would also involve plastic work that may not reflect changes in post-impact dynamics studies (Jutzi & Asphaug 2015), but may influence the heating budget. One might intuit from this that our results may underestimate the amount of impact heating, but for numerical reasons described later, this is not strictly the case.

We also note the importance of using the ANEOS/M-ANEOS equation of state as the Tillotson Hugoniot for forsterite is considerably inaccurate at pressures beyond ∼200 GPa (Stewart et al. 2020). For context, peak pressures of ∼ TPa and several tens of GPa are readily reached in even our gentlest collisions (v_imp = 1.1 v_esc) between Earth-mass and Mars-mass bodies, respectively. Furthermore, Tillotson unloading paths do not cross a phase boundary, resulting in far lower expansion velocities (affecting dynamics) and result in inaccurate densities (affecting internal energy and sound speeds; Stewart et al. 2020), which can be especially consequential during gravitational resettling events after the initial collision. We also note that along with most impact hydrocode studies, our study assumes materials are in instantaneous thermodynamic equilibrium and does not account for mixed phase, chemical evolution/reactions, and radiative transfer, common sources of inaccuracy.

In Figure 1 we show simulations with similar impact velocities (∼12 km s⁻¹) but at two different scales (roughly Mars and Earth scales). In the center, we show mass fraction of silicates that reach 3474 J K⁻¹ kg⁻¹ (entropy at complete melting and incipient vaporization assuming release pressure of 5.2 Pa, and complete melting assuming release at 10⁵ Pa). In both simulations, the initial compression generates in a precipitous increase in thermodynamic processing within the first hour; a process that predominantly scales with peak pressure and thus absolute velocity. The post-impact phase, however, involves additional heating through intense flows and reaccreting debris. Weakly bound ejected material returns at reimpacting velocities that scale with the mass of the remaining remnant, which can exceed the sound speed—especially in the impact-heated (potentially supercritical) remnant.

In Figure 2 we show the results for two entropy thresholds after 8 hr. ⁵ For each criterion the fraction of the mantle that reaches each threshold shows log-linear trends as a function of impact energy (right) and a diffuse log-linear trend as a function of impact velocity (left). The sublinear data near 100% mass fraction are expected, as the amount of material available for processing decreases considerably and only stray mass nodes are unprocessed. Values of 1% represent only ∼1400 mantle nodes and data near this region may be subject to numerical effects. Notably, almost none of the thresholds for melting or vaporization are reached in a meaningful fraction of planets smaller than 10⁻² M_⊕, which is to be expected as the escape velocity of these bodies is subsonic. Our results are reported in Table 1 in Appendix A, and we discuss the influence of numerical effects in Appendix B, which likely play a role in our results.

Zoom In Zoom Out Reset image size

For the log-linear trends in Figure 2, we provide the fit values for 3474 J K⁻¹ kg⁻¹ and 4270 J K⁻¹ kg⁻¹, respectively:

$$\begin{eqnarray}\begin{array}{rcl}{\mathrm{log}}_{10}\left({M}_{\mathrm{CM}},{M}_{\mathrm{IV}}\right) & = & \left(0.63\pm 0.06\right){\mathrm{log}}_{10}\left({{Mv}}_{\mathrm{imp}}^{2}\right)\\ & & +\left(1.43\pm 0.10\right)\\ {\mathrm{log}}_{10}\left({M}_{\mathrm{IV}}\right) & = & \left(0.64\pm 0.06\right){\mathrm{log}}_{10}\left({{Mv}}_{\mathrm{imp}}^{2}\right)\\ & & +\left(0.87\pm 0.08\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where M is in units of Earth masses and M_CM,IV are in percent. Note that the complete melting thresholds at both release pressures are identical to the incipient vaporization threshold at the low release pressure.

To understand the origin (provenance) of thermodynamically processed materials, we track the initial locations in the pre-impact bodies that reached various entropy thresholds (Figure 3). The region of greatest thermodynamic processing is expectedly established near the impact point. Two notable regions of strong thermodynamic processing arise: (1) a shell-like region along the surface of the planet that decreases in size with distance from the impact point and (2) a roughly cylindrical region reaching from the surface toward the core. Both regions are most readily observed in the M = 10⁻¹ M_⊕, 1.1 v_esc collision. Additional zones of processed material at the antipode are also seen in the cross-sectional view, as also observed in Nakajima et al. (2021). Higher-energy collisions feature more extensive melt and vapor production in each of these zones, and the zones themselves grow larger and more diffuse due to the increasing strength of shocks and post-impact effects. These symmetries break down in off-axis collisions (not shown), with more processing initially occurring along the velocity vector direction and post-impact processing occurring more globally due to the rotation of the remnant and spiraling in of accreting materials.

Using the algorithm of Emsenhuber & Asphaug (2019) to identify gravitationally bound clumps, we find the amount of escaping debris that reaches incipient vaporization in our head-on Mars-scale collisions is 77% in merging (1.1 v_esc) and 95% in disruptive (3.0 v_esc) collisions for 5.2 Pa release pressure, and 56% and 67%, respectively, at 10⁵ Pa release pressure. In collisions at the most probable geometry (45°), a dominant fraction of the debris is still melt and vapor. The exact amount varies with impact outcome, with 73%, 46%, and 49% in graze-and-merge (1.1 v_esc), hit-and-run (1.7 v_esc), and disruptive (3 v_esc) scenarios at 5.2 Pa, respectively, and 63%, 28%, and 16% at the higher release pressure, respectively. In head-on Earth-scale collisions, we find ∼99% and ∼100% of the debris reach incipient vaporization (release at 5.2 Pa) in 1.1 and 3 v_esc scenarios, respectively, and 97% and 100% at the higher release pressure, respectively. At these scales at 1.1 and 3 v_esc, respectively, we also predict 59% and 75% (low release pressure), 14% and 65% (high release pressure) to reach the 50% vaporization threshold in the 1.1 and 3.0 V_esc collisions. Across all cases at the Mars scale, the fraction of escaping debris that reaches the 50% vaporization threshold is no more than 23% and is often less than 10%. The role of numerical resolution in these results is discussed in Appendix B.

Zoom In Zoom Out Reset image size

We find that the absolute scale (total mass) of the collision influences the thermodynamic consequences of collisions between sufficiently large bodies as predicted by the SPH method, primarily due to post-impact effects. In Figure 1, we show that after the initial compression and excavation, giant impacts can feature a complex series of reimpacts and resettling events when impact velocity is near the escape velocity. These events can induce shock heating, depending on the gravitational potential of the remnant, and overprint the initial compression stage. To test the strictness of energy scaling, we performed simulations at 7.5 v_esc (highly disruptive collisions), which are greater in absolute velocity, but equal in energy to collisions at 1.1 v_esc between bodies 10 times more massive. In the M = 10⁻³ M_⊕, 7.5v_esc collision, impact heating exceeds that of its equal energy, yet lower scale counterpart at M = 10⁻² M_⊕. This deviation from energy scaling at very high multiples of v_esc can be expected. Here the collision is highly disruptive, experiencing none of the post-impact heating; however, the shocks are sufficiently disruptive to exceed the contribution of post-impact heating experienced at equal-energy collisions at larger scales. Under most planet formation scenarios, however, impact velocities are near the mutual escape velocity, where post-impact dynamics are highly consequential.

Importantly, results in Figure 2 are not only subject to numerical effects (Appendix B), but they are simply metrics for the extent of thermodynamic processing; they cannot strictly be used to determine the actual mass fraction—release and initial temperature pressure is a critical factor (Davies et al. 2020). In reality, bound silicates settle at relatively high pressure in the (sometimes supercritical) remaining remnant, or as a diffuse transient envelope; only a small fraction of the total colliding mass produced to escaping debris in 1.1 v_esc scenarios (into presumed near-vacuum). Shock-inducing reimpacts could also act to further erode the transient atmosphere (Ahrens 1993), depending on the dynamical environment and size distribution of the released debris, the latter of which is governed in part by the scale of the collision. This effect could provide a natural explanation for the considerable diversity in exoplanet densities (and inferred envelope sizes) at the super-Earth scale (Inamdar & Schlichting 2016), depending on the nature (scale, velocity, and geometry) of the disruptive collision and the efficiency in which reimpacting debris erodes the nascent, likely volatile-rich atmosphere.

We can expect that at certain scales, debris unloading into vacuum or nebula pressures is almost guaranteed to be predominantly melt and/or vapor. At roughly Mars-scale collisions, we found self-stirred velocities at common (off-axis) geometry produce a greater proportion of molten and vaporized debris than the dynamically stirred scenarios. Stirred velocities tend to produce hit-and-run, with debris much less thermodynamically processed due to poor impact energy coupling. The ultimate fate of debris in our simulations, however, is subject to processes that we do not resolve (condensation, collisional grinding, chemistry, radiative transfer, etc.) and warrants further study. Vaporized giant impact debris would be less likely to contribute a meaningful fraction of silicates to the main belt and thus would not exacerbate the "missing mantle" paradox, a long-standing problem regarding the apparent shortage of mantle-derived asteroids (Chapman 1986) and the relatively small population of asteroids thought to be collisional remnants (DeMeo et al. 2019). "Pebble accretion" (Lambrechts & Johansen 2012), for example, may rapidly establish large bodies not strictly through pairwise accretion potentially avoiding subsonic impacts between growing bodies almost entirely, especially due to the warmer (and thus lower sound speed) conditions of early bodies. It is thus likely that classical pairwise accretion and pebble accretion, and the vapor-rich debris we infer that they produce, are more consistent with the lack of mantle material in the main belt (Chapman 1986). Exoplanetary systems that quickly form into tightly packed resonant dynamical configurations, such as Kepler-223, Kepler-60, and Trappist-1 (Leleu et al. 2021), may avoid these vaporizing giant impacts, helping to explain their volatile-rich inventory inferred from density estimates. In contrast, the high-energy mantle stripping events that may explain high-density (large core-mass fraction) planets, for example the large-scale ∼3.3 v_esc collision posited for the Kepler-107 system (Bonomo et al. 2019), would be nearly entirely vapor. This would reduce the tendency for ejected material to be preserved dynamically and reaccrete, retaining the large inferred core-mass fraction of the remnant. Finally, a natural consequence of the transition to vapor- and melt-dominated debris is likely in the interpretation of giant impact candidates in the debris disk catalog (Lisse et al. 2009); observed amorphous silicates are predominantly formed in sufficiently large impacts that generates quickly quenched melt and vapor. It thus becomes clear the timing and scale/size distribution of colliding bodies and their dynamical environment is critical to interpreting debris disk observations, understanding exoplanet diversity, and evaluating various planet formation models.

It is important to note that the provenance of giant impact heating is not entirely global, especially when considering common geometries, which is generally consistent with Nakajima et al. (2021). High-speed hit-and-run collisions preferentially heat one hemisphere, leaving only graze-and-merge and highly disruptive collisions as candidates for a global heating provenance. This is especially meaningful as impact angle distributions are centered about θ_imp = 45°, according to $P({\theta }_{\mathrm{imp}})=\sin (2{\theta }_{\mathrm{imp}})$ (Shoemaker 1961). In graze-and-merge collisions, the two bodies spiral inward (potentially in a series of off-axis recollisions) to form the final remnant, producing a more globally distributed provenance of thermodynamically processed material than we present in Figure 3. On the other hand, in the case of stirring, say from 1.4 v_esc to 2 v_esc, collisions are ∼10%–20% more likely to result in a "hit-and-run" (depending on impactor-to-target mass ratio and the core-mass ratio; Kokubo & Genda 2010; Gabriel et al. 2020). ⁶ We predict that these higher-energy, hit-and-run collisions will exhibit velocity-scaled heating, with subsequent overprinting by a trail of reimpacting debris. Together, the diverse outcomes enabled by common geometries challenges the assumption of full reequilibration post-giant-impact, a common feature in most planetary evolution models.

We outline important new functional dependencies of impact heating that are unique to planet-scale collisions. Because of post-impact resettling, simulations suggest that merging collisions can involve additional heating beyond that predicted by pure velocity scaling. Furthermore, due to the dependence of impact velocity on mass of the colliding bodies and stirring by dynamical neighbors, even in the gentlest collisions between Mars-sized bodies and larger, melting and vaporization are unavoidable. This effect has important consequences for evaluating terrestrial planet formation models and interpreting debris disks observations and exoplanet diversity. In self-stirred dynamical systems and in pebble accretion paradigms debris is generally less abundant than in dynamically excited systems for a given impact scenario. This overall smaller amount of debris also experiences a greater degree of thermodynamic processing, together helping to provide a solution to the long-standing missing mantle paradox for the main asteroid belt. We stress the importance of the timing and size distribution of pairwise accretion in different formation models, initial conditions (temperature and depth/pressure), thermodynamic paths (e.g., whether material unloads into high- or low-pressure environments), impact geometry, and numerical effects for accurately understanding the production of melt and vapor in planetary collisions. The latter of which we demonstrate has important implications for the prospect of developing accurate impact heating models from modern simulations that are currently considered to have sufficiently high numerical resolution.

The authors acknowledge Dr. W. Benz (University of Bern), Dr. M. Jutzi (University of Bern), Dr. A. Reufer, Dr. A. Emsenhuber (University of Arizona), and H. Ballantyne (University of Bern) for SPHLATCH code development, maintenance, support, and fruitful collaborations, as well as Dr. H. Piet, Dr. J. A. Tyburczy, Dr. C. Agnor, and Dr. C. Raskin for their helpful remarks. The authors thank Dr. B. Ivanov for his insightful review. Dr. T. S. J. Gabriel acknowledges support from the Arizona State University Space Technology and Science ("NewSpace") Initiative and from NASA grant 19-PLANET19R-0046, as well as thanks Dr. E. Asphaug for his insight over the years, as well as the late Dr. H. J. Melosh for his inspiring work and stimulating conversations. H. Allen-Sutter acknowledges support from NASA grant 80NSSC18K0353.

Software: SPHLATCH (Reufer 2011; Emsenhuber et al. 2018).

Following IOP Publishing's Name Change Policy, one of the authors has changed their name since this article was published. Author Harrison Allen-Sutter has been corrected to Harrison W. Horn, as of 2023 May 17.

SPH is a convenient Lagrangian methodology that allows us to track the thermodynamic history of material nodes through time and space, as opposed to Eulerian (grid-based methods). However, when interpreting results from SPH simulations, it is critical to understand potential sensitivities to numerical effects that can influence thermal and dynamical outcomes. Resolution, for example, is known to have a well-documented effect on the mass of the largest remnant in a collision, with higher-resolution simulations producing less disruptive outcomes under identical conditions to lower-resolution simulations (Genda et al. 2015). To assess resolution effects, we simulated lower-resolution impacts with 1 × 10⁴, 5 × 10⁴, and 3 × 10⁵ SPH nodes in each body, in addition to our nominal resolution of 1 × 10⁵ SPH nodes (see Figure 4). Overheating is consistently associated with lower resolution. The final mantle fraction that reaches incipient vaporization of each M_tar = 1 M_⊕ simulation approaches ∼100%, but shows distinct differences at earlier times. In the M_tar = 10⁻¹ M_⊕ simulations (gray lines), the final heat evolution shows a difference of over ∼10% between the nominal and high-resolution case; the deviation in lower-resolution simulations is even more egregious. A similar giant impact heating study reported that the overall budget is resolution invariant up to 1 × 10⁵ (Nakajima et al. 2021); however, our results demonstrate that the time evolution is not resolution invariant and impact heating in lower-energy simulations is not resolution invariant at any time after impact. We generally conclude that SPH simulations may seem converged in certain scenarios or diagnostics (e.g., disruption threshold), but can be otherwise unconverged elsewhere (e.g., total heating budget), even at widely accepted levels of numerical resolution.

Another important consideration is the choice of artificial viscosity. Shocks are spatially discrete phenomena that are resolved in SPH through the use of "artificial viscosity," which ensures the momentum field is differentiable—a requirement of the SPH conservation equations. Entropy in our simulations is only produced through the application of artificial viscosity, which is triggered in shocks, and our simulations show energy is conserved to better than 1%. However, we recognize the effect of subsonic dissipation in SPH, which can inhibit subsonic turbulence due to the overapplication of artificial viscosity (overheating) in certain types of flows (Deng et al. 2019). In light of this, we also show the result of using a more robust form of artificial viscosity that implements a time-dependent term designed to limit the application of artificial viscosity in the vicinity of a shock (Rosswog 2009; dashed lines in Figure 4). The simulation that is both high resolution and with the time-dependent artificial viscosity term expectedly shows the lowest amount of heating.

We find that the fraction of heated silicates in the escaping debris can also be sensitive to resolution and artificial viscosity, depending on the scale of the collision. For example, the amount of escaping debris that reaches the entropy thresholds only decreases by a few percent in the M_⊕, 1.1 v_esc head-on cases that utilize high resolution and/or time-dependent artificial viscosity. At the roughly Mars scale, however, utilizing time-dependent artificial viscosity and high resolution in the head-on, 1.1 v_esc case decreases the escaping debris that reaches incipient vaporization by roughly half, to 43% and 29% at the low and high release pressure, respectively.

Zoom In Zoom Out Reset image size

There is likely a complicated interplay between heating, resolution, and impact scale. Low-resolution simulations tend to produce more escaping debris leaving behind a smaller intact remnant (Genda et al. 2015). High-resolution simulations on the other hand may produce less escaping debris and initial heating, but the larger remaining remnant would have a greater escape velocity and a better populated debris size distribution, resulting in more energetic (potentially supersonic) secondary impacts of more coherent bodies. We expect this effect to be further complicated by the fact that the phase and nature of debris will transition from being solid at small scales and low velocities, to melt and vapor dominated at larger scales and higher impact velocities. This will result in reimpacting remnants of intact bodies in the former case and more gentle plume-like reimpacts in the latter case.

We also point out areas where numerical instability may be reflected in our results. In the small-scale collisions, particles along the margins of the planet reach various entropy thresholds despite the comparatively low, sometimes subsonic velocities (see the M = 10⁻³ M_⊕, 3 v_esc case in Figure 3). This is likely due to well-documented SPH shortcomings that produce numerical instability at interfaces and free surfaces (reviewed in Agertz et al. 2007), which may preferentially affect escaping debris. The toroidal zones of processed silicates near the impact point (see the M = 10⁻² M_⊕ case in Figure 3) may be related to impact jetting phenomena, where material is accelerated in a shear flow that can exceed the impact velocity, driven in part by the collapsing interface on impact (Walsh et al. 1953; Melosh 1989; Gerasimov et al. 1998). Indeed, upon closer inspection we find velocity enhancement is observed in the syn- and post-impact flows. We do not include models for plastic deformation or pore crushing that could serve as heat source of these materials, and we expect that some level of artificial viscosity triggering in these strong shear flows is possible (Rosswog 2009) and may contribute to errant heating in these zones.

The combination of the effects of subsonic dampening and resolution dependence likely plays a role in SPH simulations; however, there are also important physical factors that may help explain melt and vaporization predictions where unexpected. For example, the lowest entropy threshold is reached in some of our simulations below the 8.2 km s⁻¹ threshold reported in Davies et al. (2020). However, the effects of impact jetting , either driven by artificial viscosity or otherwise, play a role in driving up effective mutual velocities. We also note that our planets at this scale are initially warmer by ∼300–350 K than Davies et al. (2020), which can decrease the impact velocity threshold from colder conditions. In any case, caution should be used in the application of our scaling laws as they can can only be strictly applied to a limited number of scenarios and with consideration for the numerical effects discussed above. Instead, we emphasize their utility in describing a relatively understudied thermodynamic transition in large-scale impacts and outlining the role of impact scale and gravitational stirring on the impact heating budgets as produced in common SPH calculations.

Although not an exhaustive numerical study, our tests suggest that the amount of shock-heated silicates predicted in SPH simulations with modest resolution and a common form of artificial viscosity can be overestimated (by at least 10%–20%, depending on the situation) and careful numerical convergence studies and attention to the associated flows regimes should be employed to make robust predictions. We note that the prospect for explicitly reaching numerical convergence with current computing and software may be poor, as simulations with bleeding-edge resolution are not converged for even low-level diagnostics (e.g., mass of the largest remnant; Genda et al. 2015) and can even exhibit unconvergence (Hosono et al. 2017). Numerical aspects are thus a paramount area of future work in order to assess the true fidelity of SPH giant impact simulations.