Atmosphere loss in planet–planet collisions Free

ABSTRACT

1 INTRODUCTION

There have been a wealth of exoplanet discoveries over recent years, bringing the total number of confirmed planets to more than 4000. About 20 per cent have a mass between that of the Earth and the Neptune (henceforth super-Earth mass). Of these, two thirds orbit at a distance less than that of Mercury from their host stars, and more than three quarters are in systems with at least one other planet (NASA 2019).

Planets in this mass and orbital distance range are very diverse, with measured densities ranging from between 0.03 and 12.7 |$\mathrm{g\, cm}^{-3}$|⁠. Even planets within the same Solar system at similar orbital radii can be vastly dissimilar. The two innermost planets of the Kepler-107 system (Bonomo et al. 2019), for example, both have similar radii (1.5–1.6R_⊕), but Kepler-107b has a density of |$5.3\, \mathrm{g\, cm^{-3}}$| compared to Kepler-107c’s |$12.6\, \mathrm{g\, cm^{-3}}$|⁠. This high density implies that Kepler-107c must have a different composition from 107b.

There are two main theories that explain the diversity of densities observed in Super-Earths and Mini-Neptunes: (1) XUV radiation from the central star (Lopez, Fortney & Miller 2012; Lopez & Fortney 2013; Owen & Wu 2013; Jin et al. 2014) strips close orbiting planets of some of their lighter elements, leaving them denser. (2) Giant impacts eject the lighter outer material (such as crust or atmosphere) from both bodies leaving a denser remnant planet (Inamdar & Schlichting 2016; Bonomo et al. 2019).

XUV radiation cannot always explain large differences in density for planets orbiting at similar semimajor axes in the same planetary system, for example, Kepler-107b and c. These two planets orbit at a similar distance from their parent star and have similar physical radii but differ significantly in planetary mass. The outer Kepler-107c is more than twice as dense as the innermost 107b. Kepler-107c’s density can not be explained by XUV radiation because it is orbiting outside the less massive and also less dense Kepler-107b, which would have lost more material due to irradiation. Thus, a giant impact is the best explanation (Bonomo et al. 2019).

Giant impacts have been suggested as the explanation for density diversity in several planetary systems (Liu et al. 2015; Inamdar & Schlichting 2016). In addition, many of the tightly orbiting high-multiplicity systems detected by Kepler appear to be on the borders of stability (Fang & Margot 2013). Therefore, it is not unreasonable to consider the planets g⊕tationally interacting with one another in such a way that either they eject one or both planets or that they collide with one another. Barnes & Raymond (2004) showed that small perturbations in orbit will lead to such ejections or collisions. Volk & Gladman (2015) even suggested close orbiting groups of planets are common in the formation of inner solar systems (a < 0.5 au), including our own, but that they are not stable long term, and will undergo collisional disruption and consolidation. This suggests that the high multiplicity systems that we observe are either these initial unstable systems, or that the planets are in an arrangement that is stable for long time periods.

Several previous works (e.g. Schlichting, Sari & Yalinewich 2015; Inamdar & Schlichting 2016) have demonstrated that giant impacts have the potential to remove large fractions of these planets’ atmospheres, leading to substantial density enhancement. In this paper, we directly calculate atmosphere stripping via 3D modelling of head-on giant impacts between Mini-Neptune mass planets with significant atmospheres and lower mass bare super-Earth impactors.

1.1 Previous work

Until recently, simulations of atmospheric losses due to such giant collisions were too computationally expensive to run within a reasonable time frame. Recent advances, however, have meant high resolutions are now much more attainable.

Because of computational demands much early work focused on analytical predictions (some of which focus on significantly lower mass atmospheres than we consider here, such as Earth-like atmospheres), for example, Genda & Abe (2003) who discuss how much of an Earth-like protoplanet’s atmosphere is likely to survive the giant impact phase. They used one-dimensional models to calculate the amount of material jettisoned after a collision from the ground velocity beneath it, showing that for the canonical Moon-forming impact only ≈20 per cent of mass would be lost. Genda & Abe (2003) also showed that the ground velocity needs to reach the escape velocity for total atmosphere erosion. Schlichting et al. (2015) expanded upon this by formulating a method of predicting atmospheric loss from a wide range of projectiles colliding with terrestrial planets. This method is the one used by Inamdar & Schlichting (2016) when discussing the density diversity of super-Earth mass exoplanets. Inamdar & Schlichting (2016) showed that a single collision between similarly sized exoplanets can cause a decrease in the mass ratio of atmospheric envelope to central core material by a factor of 2, which in turn leads to a density increase of a factor of 2–3. Inamdar & Schlichting (2016) thus suggested giant impacts as a cause for the observed density diversity.

Despite computational limitations some progress has been made in simulating collisions of planets with gaseous envelopes. Liu et al. (2015) considered a model of the Kepler-36 system (with target planets of 6–8 per cent atmosphere by mass). They found that a collision could cause the density difference observed between the two planets in this system, and suggested that giant collisions might therefore be the cause of the dispersion we observe in the mass–radius relationship for super-Earth mass planets. Hwang et al. (2017a,b) used an N-body code to model the evolution and stability of high multiplicity planetary systems (specifically Kepler-11 and Kepler-36). For the collisions, they used a smoothed particle hydrodynamics (sph) code and assumed atmosphere mass fractions of 5–15 per cent to determine the results of highly grazing collisions (collisions where only the atmospheres overlap), they showed that typically the higher mass core will accrete more of the disrupted gas envelope leading to increasing density contrasts between the two planets. Due to problems with their equation of state, however, they where unable to simulate the results of head-on collisions. Kegerreis et al. (2018) also used sph to model collisions involving targets with atmospheres, to see if they could model the formation of Uranus’ off-axis rotation and unusual magnetic field. Unlike this paper, Kegerreis et al. (2018) focus mostly on higher impact parameter collisions in order to study the change in rotation. They considered ice giants as opposed to the metal and silicate planets with atmospheres studied in this work.

In this paper, we present the results from a series of head-on collisions between super-Earth mass planets where each of the three different mass targets is a Mini-Neptune with a significant hydrogen envelope (8–33 per cent). We run a large series of numerical simulations with a wide array of atmosphere-less super-Earth projectiles. We provide scaling relations for material loss applicable to a wide range of masses that can be used in N-body simulations and population synthesis models.

2 METHODS

2.1 Numerical code

The simulations presented in this paper were run using the sph code gadget-2 (Springel 2005). Although gadget-2 was initially designed for simulations on cosmological scales, we have used a modified version to model our planets that includes tabulated equations of state (EOS) for the planetary constituents. For further detail on these modifications see Marcus et al. (2009) and Ćuk & Stewart (2012). We further modified gadget-2 to include an ideal gas atmosphere component.

In sph codes, such as gadget-2, the material is split into separate particles each representing an ensemble of material. The continuous fluid properties, such as the density, for each ensemble of particles are calculated using kernel interpolation methods. The g⊕tational forces on the other hand are calculated using hierarchical tree methods (Springel 2005). We ran gadget-2 in ‘Newtonian’ mode with time-step synchronization, and the standard relative cell-opening criterion. We use the standard time-step criterion as described in Springel (2005), where the smaller of the time-step based on the g⊕tational softening and the acceleration, or the courant condition is used (with a courant factor of 0.1). We use the standard artificial viscosity formulation for gadget-2 as described in Springel (2005), with a strength parameter of 0.8.

We modelled the planets as two or three material systems, each planet (both projectile and target) consisting of iron core, forsterite (silicate) mantle, and a hydrogen atmosphere for the target only. In a similar fashion to prior studies (Marcus et al. 2009; Ćuk & Stewart 2012), the mantle was twice the mass of core. We modelled the atmosphere as a monatomic ideal gas for simplicity and ease of comparison (see Section 4.1). Tabulated ANEOS/MANEOS EOS from Melosh (2007) were used for the iron and forsterite (these tables are available from Carter, Lock & Stewart 2019).

2.2 Hardware

The simulations were each run using a full node on the University of Bristol’s Bluecrystal supercomputers, either on phase 3 or phase 4. Phase 3 nodes consist of 16 core 2.6 GHz SandyBridge processors with 59.7 GiB RAM altogether (Gardiner 2015), whilst phase 4 nodes have two 14 core 2.4 GHz Intel E5-2680 v4 (Broadwell) CPUs with the whole node having a combined 128 GiB of RAM (Gardiner 2017). Each collision took between half a day and a day depending on the masses and impact energies involved.

2.3 Initial conditions

To generate the initial planets, we began by creating the central core and mantle. To generate density profiles we used radial temperature profiles (from Valencia, O’Connell & Sasselov 2006) as well as estimates for the average bulk density of each type of material (i.e core or mantle) and the radial range that the core and mantle occupies. From this initial assumption of constant density per material layer, we iteratively generated new density profiles using g⊕tational and hydrostatic pressure calculations, until we obtained a density profile that was consistent with our EOS, the expected hydrostatic pressure, and the input temperature profile.

Once a consistent profile was obtained, it was then used to generate the position of each of the particles by splitting the planet into a number of radial bins based on its mass and randomly positioning a number of particles within each bin proportional to the bin’s mass. The final number of bins being adjusted so that we could reach our desired mass to a tolerance of 1 per cent. The type of material being added to each bin was decided by the material ranges given for the initial density profile. After generating an initial sph planet it was then equilibrated in isolation for 10⁵ s (simulated seconds) using a preliminary run of gadget-2, to ensure that we were running simulations with a planet that was stable. During equilibration we used two ‘artificial cooling’ methods, velocity damping of the particles, i.e applying a restitution factor of 50 per cent each time-step (see Carter et al. 2018), and also entropy forcing, the entropy of each particle is reset to a constant value for each material at every time-step (⁠|$1.3\, \mathrm{kJ\, K^{-1}\, kg^{-1}}$| for the iron core, |$3.2\, \mathrm{kJ\, K^{-1}\, kg^{-1}}$| for the mantle). This entropy forcing ensures that we produce planets with isentropic layers.

We then used this radial mass profile for the atmosphere to determine the position of atmosphere particles, by splitting the profile into radial bins and placing particles proportionally to the mass of the bins at random positions within them, in a similar fashion to the core and mantle. This new body was equilibrated for a longer time of 4–8|$\times 10^5\, \mathrm{s}$| until the radius of the planet had converged to a constant value. The pseudo-entropy of the atmosphere was forced to a value of |$5\times 10^{11}\, \mathrm{Ba\, g^{-\gamma }\, cm^{3\gamma }}$| where the adiabatic index was γ = 5/3, this ensured that the atmosphere would not reach densities where it would sink underneath the mantle and core material, but also that the base of the ideal gas atmosphere was as close as possible to our predicted temperatures.

For the core and mantle of our targets, we used a resolution of 10⁵ particles. All other particles in each simulation were made the same mass as the core and mantle particles, resulting in a total resolution between 1.2 × 10⁵ and 2.5 × 10⁵ particles. This results in an atmosphere ‘thickness’ of between 5 and 10 layers of particles. Our total particle number is smaller than that suggested by Kegerreis et al. (2019), however, we are interested in large changes in mass of the largest remnant so our resolution should be sufficient. We have run resolution tests of head-on collisions of 5 M_⊕ planets with 0.9 M_⊕ atmospheres against one another at both half and double our standard resolution (of 10⁵ particles in the target core and mantle) at impact velocities of 20 and 40 |$\mathrm{km\, s^{-1}}$|⁠. The key quantities we measure are mass of the largest remnant, atmosphere loss fraction, and core and mantle loss fraction, each vary from the mean value at that velocity by less than 5 per cent except for the core and mantle loss fraction at low velocity, where only a very small number of particles are being lost. We did not consider losses or remnants of only a few hundred particles resolvable with sph methods.

2.3.1 The point of impact

In collisional studies, impact parameters such as velocity, impact energy, etc. are normally measured in terms of the point of first physical contact between the two planetary bodies (e.g. Leinhardt & Stewart 2012). For planets with atmospheres this becomes more complex, however, as atmosphere densities decay approximately exponentially so there is no clear boundary at which the atmosphere ends. Tidal forces between planets also have a stronger distorting effect on the atmosphere than the core and mantle. Thus, we use the time when mantle surfaces touched as our point of impact because this can be clearly defined.

To obtain an initial start position from our desired collision parameters (for example, velocity and impact angle), the following process was used: (1) the two planets were represented by point masses at their centres, and these two planets were placed at their position at the point of impact and set at the predicted velocity; 2) Time reversal symmetry was then used to trace the path of the projectile back to a separation of five times the sum of the projectile and target mantle radii (excluding atmosphere). To determine the projectiles path we used a simple Verlet integrator, set the target planet to be centred on the origin, and calculated the acceleration as given by the relative g⊕tational force between two-point particles. The choice of starting separation was a somewhat arbitrary one intended to reduce the tidal forces compared to the starting distances of previous similar studies due to us equilibrating planets in isolation and atmospheres being more easy to tidally distort.

Note using the point where mantles touch as our point of impact is not without its drawbacks. As can be observed in Fig. 1, the presence of the atmosphere causes a significant observable slowing of the relative impact velocity. This slowing is due to atmospheric drag that distorts the projectile planet. For some collisions, this distortion is even present when the leading edge of the projectile is greater than an atmospheric scale height from the target. In further calculations, which use impact energy, we used the measured velocity instead of the predicted input value. To measure this collision point precisely, we re-simulated the point of impact for our collisions with a higher output frequency (snapshots were taken every 3 s as opposed to every 100).

Figure 1.

Comparisons of the predicted velocity for the three targets and the ratio between measured and predicted impact velocities. Impact velocities are given in terms of the mutual escape velocity (⁠|$v_{\mathrm{esc}}=(2GM_{\mathrm{tot}}/R^{^{\prime }})^{\frac{1}{2}}$| where |$R^{^{\prime }}=\left(3M_{\mathrm{tot}}/4\pi \rho _{\mathrm{bulk}}\right)^{\frac{1}{3}}$| and ρ_bulk is the bulk density of the simulated target). Edge colours show target mass, while central colours show projectile mass. The reduced velocity compared to the prediction can be considered a measure of the drag caused by the atmosphere. As might be expected, the denser, higher mass atmospheres around the larger targets tend to cause more drag, the lower mass and lower velocity projectiles tend to experience proportionally more drag as well due to their lower momenta.

Open in new tabDownload slide

2.3.2 Input parameters

2.4 Run parameters

2.5 Analysis methods

2.5.1 Determining bound material

To determine what material was bound in the largest post-collision remnant we used the same iterative method as Marcus et al. (2009) and Carter et al. (2018). We began by locating the particle closest to the potential minimum and using kinetic energy and g⊕tational potential to determine which other particles were g⊕tationally bound to this deepest particle. Using the total mass and centre of mass of these particles as our new seed, we then iteratively ran through the process of determining that extra particles were bound to the seed, and then adding them to the seed, until either the change in mass of the seed was below a set tolerance or a maximum number of iterations had been reached. For our simulations we, obtained single particle differences in mass within a few iterations so our tolerance was set to a single particles mass, and the maximum number of iterations was never reached. Material type was also tracked to determine the mass of each type of material, which was unbound or part of the largest remnant.

2.5.2 Collision categorization

To categorize the collision outcomes for our data, we consider separately atmosphere (A) and core and mantle material (C). For atmospheres if there was greater than 95 per cent of the initial atmosphere mass in the final remnant we considered there to be a merger (AM), on the other hand if there was less than 5 per cent remaining we categorized this as total loss (TAL), the rest were considered to undergo partial loss (AL). Although hydrogen atmospheres as small as |$0.1\small{--}1{{\ \rm per\ cent}}$| of a planet’s total mass can have a significant effect on its radius, we do not have the resolution to accurately probe atmosphere mass-losses that small. For core and mantle, we define mergers (CM) for >95 per cent of the total mass of core and mantle from both projectile and target remaining in the largest remnant. If the mass of core and mantle in the largest remnant is greater than that initially in the target we have an accretion event (CA), if it is less we have an erosion event (CE). If there was less than 10 per cent of the total mass remaining in the largest remnant we define it as a supercatastrophic disruption (SCD). Tables 1, A1, and A2 categorize our final results for each collision using this system.

3 RESULTS

The main aim of this paper is to determine scaling laws for atmosphere loss as well as total material loss during head-on giant impacts. The wide breadth of impact energies that are simulated in this work have uncovered a broad range of outcomes from near perfect merging events to highly energetic catastrophic disruption (see Fig. 2).

Figure 2.

Cross-sectional snapshots sliced through the mid-plane for a series of head-on collisions showing different collision outcomes from merging to catastrophic disruption as the specific relative kinetic energy Q_R increases from row A to row D. Colour denotes material, red and pink – iron core, orange and yellow – forsterite mantle, and blue – hydrogen atmosphere. The additional colours in the first panel denote material that will not be bound by the end of the simulation, black being atmosphere, and grey core or mantle. All collisions shown have the same target mass (⁠|$M_\mathrm{t} = 6.26\, \mathrm{M}_{\oplus }$|⁠) but differ in projectile mass (M_p) and impact speed (⁠|$\mathbf {v_\mathrm{imp}}$|⁠): (a) Atmosphere and core and mantle merger – |$\mathbf {v_\mathrm{i}} = 20$| km s⁻¹, |$M_\mathrm{p}=0.25\, \mathrm{M}_{\oplus }$| (Table 1, 5–0); (b) Atmosphere loss and core and mantle merger – |$\mathbf {v_\mathrm{i}} = 30$| km s⁻¹, |$M_{\mathrm{p}}=2\, \mathrm{M}_{\oplus }$| (Table 1, 5–17); (c) Total atmosphere loss and core and mantle erosion – |$\mathbf {v_\mathrm{i}} = 50$| km s⁻¹, |$M_\mathrm{p}=4\, \mathrm{M}_{\oplus }$| (Table 1, 5–32); (d) Supercatastrophic disruption |$\mathbf {v_\mathrm{i}} = 60$| km s⁻¹, |$M_\mathrm{p}=5\, \mathrm{M}_{\oplus }$| (Table 1, 5–39). Post-collision remnants were inflated in comparison to the initial planets and the expected radius of the resultant planet. This ‘puffiness’ is because we do not cool our final remnants until they reach equilibrium, we only run the simulations until the mass of bound material converges. Videos of these four collisions can be found online at http://www.star.bris.ac.uk/TDenman/Paper1_Videos/.

Open in new tabDownload slide

In the process of determining the loss scaling laws we have found that the atmosphere loss, core and mantle material loss, and the largest remnant mass are well-behaved functions of specific impact energy: |$Q_\mathrm{R} = \frac{1}{2}\mu V_{\mathrm{imp}}^2 / M_\mathrm{tot}$|⁠, where μ is the reduced mass, V_imp is the impact velocity, and M_tot is the sum of the projectile and target masses (see Leinhardt & Stewart 2012). We find that atmosphere dominates the material lost until the impact energy is large enough to remove more than 80 per cent of the atmosphere at which point mantle and core material begin to be removed significantly as well. It is also at this point that there is a break and steepening in slope of the largest post collision remnant mass (M_LR) as a function of Q_R (Fig. 4). In other words, the atmosphere of the target planet cannot be completely removed in one giant impact without significantly eroding the planet.

3.1 Mass-loss

Fig. 3 shows the largest remnant mass from our simulations as a function of specific impact energy for collisions with a target of mass |$M_\mathrm{t}=10.5\, \mathrm{M}_{\oplus }$|⁠. Each different colour represents a different projectile mass (M_p). Instead of a straight line (as in Leinhardt & Stewart 2012), our results appear to show two separate linear regimes: a shallower low-energy one and a steeper high-energy one.

Figure 3.

A comparison of the mass of the largest remnant for each collision compared with its specific impact energy, overlaid with a graphical representation of the process used to determine |$Q_{\mathrm{RD}}^*$|⁠. The circles represent the fraction of the total mass which remains in the largest remnant after a collision with a target mass of |$M_\mathrm{t} = 10.5 \, \mathrm{M}_{\oplus }$| as a function of specific impact energy Q_R. Each colour indicates one of seven M_p values. A filled circle represents a point used for fits, whilst the open circles are points with |$\frac{M_{\mathrm{LR}}}{M_{\mathrm{tot}}}\lt 0.2$| that we considered too close to the supercatastrophic disruption regime, which our fit is not designed to cover, for all collisions where <100 particles were observed in the largest remnant, we also considered the results to be below the resolution limit of the simulation. The solid lines represent our fit to the data for each projectile-target mass ratio, following equation (8). From this fit, the empirically determined value of |$Q_{\mathrm{RD}}^*$| is shown on the horizontal axis by the intersection of a coloured dotted line matching the M_p colour and M_LR/M_tot = 0.5 on the y-axis.

Open in new tabDownload slide

Figure 4.

Top: M_LR versus |$Q_\mathrm{R}/Q_{\mathrm{RD}}^*$| for all simulations. Colour indicates target mass, the open circles represent points that were excluded from the fit. Coloured solid lines are best fits to data for a given M_t. The black solid line is the best fit to the entire data set. The dotted black line is the universal law from Leinhardt & Stewart (2012). Middle: Fraction of atmosphere lost versus |$Q_\mathrm{R}/Q_{\mathrm{RD}}^*$|⁠. The dashed coloured lines indicate best fits. The dotted vertical and horizontal lines show the specific energy of the break in the M_LR fit and the respective fraction of atmosphere loss. Bottom: Fraction of core and mantle lost versus |$Q_\mathrm{R}/Q_{\mathrm{RD}}^*$|⁠. The dashed lines are a broken linear fit. The shaded vertical sections show the difference between the specific energy of the break in the M_LR fit and the break in the core loss fit. For all three graphs, the horizontal error bars represent the error in the |$Q_{\mathrm{RD}}^*$| determination method.

Open in new tabDownload slide

To ensure a robust fitting method, despite the low number of data points, we began by generating multiple random sets of input parameters for least-squares fitting. We determined a best fitting set of output parameters for each input set, and chose the fit with the smallest squared residuals as our final best fit. The fits we obtain for each set of projectile and target masses can then be used to determine the catastrophic disruption threshold by measuring the energy at which half the mass of the system is in the largest remnant (as shown in Fig. 3). The errors on |$Q_{\mathrm{RD}}^*$| were determined by propagating errors generated by the least-squares fitting algorithm.

We also investigate the dependence of core and mantle loss on the impact energy (bottom panel Fig. 4). We observe a similar broken linear relation to the broken linear relation seen for the mass of the largest remnant (Figs 3 and 4, top). In this case, we have a shallow gradient below the pivot energy (negligible for the two larger mass targets) indicating very little core and mantle mass-loss, and a much steeper gradient (therefore much greater loss) above. For lower mass targets (lower atmosphere fractions), there is greater loss of core and mantle at low energies.

Parameters for these fits are given in Table 4. As can be seen from the coloured vertical bars in the bottom panel of Fig. 4, there appears to be a correlation between this pivot energy and the pivot energy found for the mass of the largest remnant. The similar energies of these pivot points implies that the break in slope for largest remnant mass is due to the break in slope for loss of core and mantle material. We therefore categorize impacts into two energy loss regimes: below the pivot energy we have atmosphere-dominated loss, and above it we have substantial core and mantle loss.

An important point to note here is that at this pivot energy, the amount of atmosphere remaining is consistently |$20\small{--}30{{\ \rm per\ cent}}$|⁠. This means that a single giant impact cannot remove all of a planet’s atmosphere without also removing core and mantle material.

3.2 Atmospheric loss scaling law

Figure 5.

Top: Fractional atmosphere loss compared to specific impact energy scaled with respect to the pivot energy for each target mass. The amount of energy required to remove more atmosphere increases as the amount of atmosphere removed increases. Fitted to the data is a simple quadratic that has been fixed to cross the origin and have a peak at total atmosphere loss (equation 13). Bottom: Specific energy of the pivot normalized by the total mass of the system compared with projectile-target mass ratio. Mass normalized pivot energy appears to be approximately constant, decreasing slightly as projectile masses increase. The fit is given by equation (14).

Open in new tabDownload slide

Equation (14) has quite a shallow gradient, meaning that the energy of the pivot is strongly dependent on the the total mass, i.e. the impact energy is strongly dependent on the square of the total mass, implying it is related to the total g⊕tational binding energy of the system. The error envelope in Fig. 5 is generated from the variance of the fit parameters. Together equations (13) and (14) can be used to predict the mass fraction of atmosphere lost from the target for any head-on giant impact in the regime tested.

4 DISCUSSION

4.1 Caveats

4.1.1 Atmosphere fraction effects

We chose the initial atmosphere fraction for each target planet based on the most likely atmosphere fractions from the Bern global planet formation model (equation 1). One issue with this method is that it makes distinguishing between the effects of target mass and atmosphere fraction somewhat difficult, as both are directly related. An important result that requires further investigation is the relation between atmosphere fraction and pivot energy. One might expect that less massive atmospheres would require less energy to remove and thus core and mantle material would begin to be removed earlier, decreasing the pivot energy. We observe a decrease in pivot energy but we cannot distinguish whether this is due to a smaller atmosphere fraction or a shallower g⊕tational potential well from a less massive target. In future, we could test the cause of the decrease in pivot energy by simulating collisions involving targets of the same masses with different atmosphere fractions.

4.1.2 Equation of state effects

One simplification we have made in this work is to use ideal gas atmospheres as opposed to a more realistic EOS. For Earth atmospheric pressures this would have a negligible effect, but we are dealing with significantly more massive atmospheres (roughly 7 orders of magnitude). For such massive atmospheres, the pressures and densities at the base of the atmosphere are such that the assumption that molecules do not interact (inherent in using an ideal gas) may not be realistic. The lack of interparticle interactions means that we have unrealistically high densities at high pressures, (the base of our atmospheres post equilibration being 1.2–19.2 × 10⁵ atm or 1.2–19.2 × 10⁹ GPa). One effect of this is a significant compression of the atmosphere when equilibrating a gadget-2 planet. During the equilibration, the radii of atmospheres shrank on average by a factor of 2. Despite these issues, we have used an ideal gas as a starting point, more realistic EOS should be the subject of a future study.

4.1.3 Impact angle effects

For impact angle we have also elected to investigate the simplest case – head-on collisions. In general, we would expect an average collision angle of 45° (Shoemaker & Hackman 1962). Leinhardt & Stewart (2012) present in their prescription a method to relate the mass-loss of an off angle collision with a head-on one, which uses an interacting mass. In the impact scenarios considered in this work, we would expect the impact angle correction to become more complicated due to the density contrast between atmosphere and mantle. Hwang et al. (2017a,b) have previously done work on grazing collisions where the cores and mantles did not touch. Their results show that mass-loss follows a power-law dependence on impact parameter. We plan to investigate this further in future work.

4.2 Catastrophic disruption threshold scaling

Fig. 6 details how our results compare to Leinhardt & Stewart’s predictions. Our results run parallel to the prediction of equation (15) indicating a larger c* value of 2.52 (compared to the value given in Leinhardt & Stewart 2012). The larger c* means that a greater excess over the binding energy is required to remove material from planets with atmospheres than those without. This is presumably due to the increased compressibility and decreased viscosity of atmosphere compared to mantle and core material.

Figure 6.

Top: Catastrophic disruption threshold |$Q_{\mathrm{RD}}^{*}$| compared with R_C1 the radius a spherical body would have if it had the total system mass and a density of 1000 |$\mathrm{kg\, m^{-3}}$|⁠, with the black and grey lines showing their predicted relationship for different values of the strength parameter c* following equation (15). We observe a value of c* = 2.52, which is a |$32{{\ \rm per\ cent}}$| increase compared to Leinhardt & Stewart (2012)’s value of 1.9 for solely rocky bodies. Coloured lines show the predicted catastrophic disruption threshold for equal mass collisions for particular velocities. The colour of each data point indicates the relative velocity the pair of planets would have for each collision if they were equal mass. Bottom: Comparison between the catastrophic disruption threshold, and the total mass of each set of collisions. The black lines detail predictions of the catastrophic disruption threshold for each total mass for pure energy scaling for different values of c*, the dotted is our measured value of 2.52, the dot–dashed is Leinhardt & Stewart (2012)’s value of 1.9, coloured lines show predictions for pure momentum scaling. Our results (the dots with colours representing target mass as per Fig. 4) appear to follow the prediction for pure energy scaling.

Open in new tabDownload slide

It should be noted that we appear to be observing pure energy scaling for these collisions as no mass ratio correction is required, this is illustrated best in the bottom part of Fig. 6. Here, the black lines represent the predicted catastrophic disruption threshold using energy scaling, whilst coloured lines represent the momentum scaling that Leinhardt & Stewart (2012) predict we should be close to. As can be observed, our results follow the energy scaled prediction closely. This runs contrary to the results in Leinhardt & Stewart (2012) that predict near-perfect momentum scaling. Presumably, this difference is due to the presence of the atmosphere, as the Leinhardt & Stewart (2012) method is based on Holsapple & Schmidt (1987)’s crater scaling, where they show that porous materials tend to follow momentum scaling, while perfect gases tend to follow energy scaling.

4.3 Mass-loss efficiency comparisons

Figure 7.

Comparisons between the mass-loss efficiency for core and mantle, and for atmosphere to the total mass of the system. Atmosphere loss shows a decreased efficiency in comparison to core and mantle loss, as might be expected considering the increased compressibilty of atmosphere material. The filled shapes are gradients where we obtained four or greater data points for that line, open are where we have three or fewer. The grey line shows the loss efficiency predicted by Leinhardt & Stewart (2012) for rocky material, the coloured line beneath them is this value multiplied by 1 plus the atmosphere fraction for that particular target as per equation (18). This appears to show reasonable correlation with our results especially for higher mass targets. The coloured dotted line above this shows a prediction of what the atmosphere gradient must be that uses all previously derived scaling laws, is linear and passes through zero mass-loss for zero input energy (equation 19). Our results show some degree of correlation with this value, but we do not have a high enough density of data in this region to effectively probe the accuracy of this prediction.

Open in new tabDownload slide

We can consider these gradients as a measurements of the efficiency of mass-loss in each regime. The efficiency of mass-loss for the core and mantle loss region is typically at least double that for the atmosphere-dominated loss regime, this is presumably due to the increased compressibility of the atmosphere.

As can be observed, the Leinhardt & Stewart (2012) model (equation 17 above) appears to underpredict the efficiency of mantle and core material loss for planets with substantial envelopes. This is likely due to the pressure of the atmosphere above providing a resisting force to reduce mantle loss. Once the atmosphere is removed, the specific energy of the impact is higher than would have been necessary to remove significant mantle if no atmosphere was present, so more material is being removed per unit of specific energy.

4.4 Implications

This study was motivated by the density disparity observed in exoplanet systems. The post-impact bodies in our simulations are hot and inflated, often with a large mass of vaporised silicate, and thus do not represent the structure of planets millions of years after their final giant impacts (Lock & Stewart 2017; Carter, Lock & Stewart 2020; Lock, Stewart & Ćuk 2020). To determine the bulk density of the remnants from their post-collision material composition, we use the approximation suggested in Lopez & Fortney (2014). In their model the total radius of a planet, R_planet, from which we determine the density, can be calculated from summing the radial contributions of the following three components: the core and mantle, which has a power-law relation to mass; the convective envelope, which is dependent on the temperature of the atmosphere, which itself is a function of stellar flux and planet age; and the radiative atmosphere (also dependent on temperature). Here, we ignore the contribution of the radiative atmosphere due to its small effect (∼0.1R_⊕). We used an age of |$5\, \mathrm{Gyr}$| for our comparison planets (see Fig. 8) as this is the most common age for stars in the local Galactic neighbourhood. We used a flux, F_planet, of |$100\, \mathrm{F_{\oplus }}$| as the type of planets we simulated are most commonly observed at ∼0.1 au around the Sun-like stars. Fig. 8 shows contour plots of the radius (top) and density (bottom) as a function of mass and atmosphere fraction according to the Lopez & Fortney (2014) approximation described above.

Figure 8.

Top: Radius as a function of mass following the prescription of Lopez & Fortney (2014), compared against our initial targets (the squares) and also our post-collision largest remnants (the circles), the same colour scheme for target masses is used as for the rest of the paper. Bottom: Density as a function of mass calculated using the radius above, this is similarly compared against our initial targets (the squares) and also our post-collision largest remnants (the circles), collisions with full atmosphere removal have been removed so atmosphere fraction could be plotted logarithmically.

Open in new tabDownload slide

The prediction we obtained for the envelope radius of our initial targets using the Lopez & Fortney (2014) model was within 10 per cent of the radius of our initial thermodynamic profiles. This radius was, however, significantly larger than our gadget-2 targets due to the compression caused using an ideal gas. We also note that we could not reach the 20 mBar pressures which Lopez & Fortney (2014) consider to be the edge of the atmosphere with computationally practical resolutions.

For our collisions, we always observe a decrease in atmosphere fraction. This decrease in atmosphere fraction means that for all the collisions we simulate, except for those resulting in largest remnants below the resolution limit, we observe an increase in density as shown in the bottom panel of Fig. 8.

4.5 New prescriptions for atmosphere loss and largest remnant mass

Here, we summarize the process one would need to use to predict the atmosphere loss from any giant impact in the regime probed by this paper, as well as the modifications to the Leinhardt & Stewart (2012) prescription for largest post-collision remnant mass for an arbitrary head on collision between a mini-Neptune with a significant gaseous envelope and a lower mass super-Earth without an atmosphere. This algorithm can be incorporated into N-body codes and population synthesis models.

• For a given collision scenario (M_p, M_t, and V_imp), calculate the specific relative kinetic energy of the impact,

  $$\begin{eqnarray*} Q_\mathrm{R}=\frac{1}{2}\mu \frac{V_\mathrm{imp}^2}{M_{\mathrm{tot}}}, \end{eqnarray*}$$

  (20)

  where M_tot = M_t + M_p, and μ = M_tM_p/(M_t + M_p) is the reduced mass, M_t is the target mass, M_p the projectile mass, and V_imp is the impact velocity.
• Then, calculate the specific kinetic energy of the transition between the atmosphere loss and core and mantle loss regimes:

  $$\begin{eqnarray*} Q_{\mathrm{piv}}=\left(\frac{M_{\mathrm{tot}}}{M_\oplus }\right)\left(-2.45\frac{M_\mathrm{p}}{M_\mathrm{t}}+14.56\right) \,\,\,\, [10^6\, \mathrm{J\, kg^{-1}}]. \end{eqnarray*}$$

  (21)
• Calculate the catastrophic disruption threshold,

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

  (22)

  using a value of c* = 2.52 for collisions involving planets with atmospheres (⁠|$Q_{\mathrm{RD}}^{*(\mathrm{New})}$|⁠) and c* = 1.9 (⁠|$Q_{\mathrm{RD}}^{*(\mathrm{LS12})}$|⁠, as in Leinhardt & Stewart 2012) for targets with no atmosphere. |$\rho _{\mathrm{1}}=1000\, \mathrm{kg\, m^{-3}}$|⁠, and |$R_{\mathrm{C1}} = \left(\frac{3M_{\mathrm{tot}}}{4\pi \rho _{\mathrm{1}}}\right)^{\frac{1}{3}}$| is the radius a spherical body would have if it had the total system mass and a density of ρ₁.
• Calculate the gradients of each of the linear sections of the largest remnant mass fraction relation. For the core and mantle loss regime this is

  $$\begin{eqnarray*} m_{\mathrm{c\&m}}=\frac{-0.5(1+f_{\mathrm{atmos}})}{Q_{\mathrm{RD}}^{*(\mathrm{LS12})}} , \end{eqnarray*}$$

  (23)

  where f_atmos is the mass fraction of the target which is atmosphere.
• Zero impact energy means zero mass-loss therefore the gradient for the atmosphere loss dominated part of the relation is, from equation (9),

  $$\begin{eqnarray*} m_{\mathrm{atmos}}=\frac{m_{\mathrm{c\&m}}\left(Q_{\mathrm{piv}}-Q_{\mathrm{RD}}^{*(\mathrm{New})}\right)-0.5}{Q_{\mathrm{piv}}}. \end{eqnarray*}$$

  (24)
• Next, calculate the supercatastrophic disruption threshold, taking this to be where <10 per cent of the initial mass ends up in the largest remnant (following Leinhardt & Stewart 2012) we obtain

  $$\begin{eqnarray*} Q_{\mathrm{supercat}}=Q_{\mathrm{RD}}^{*(\mathrm{New})}-\frac{0.4}{m_{\mathrm{c\&m}}}. \end{eqnarray*}$$

  (25)
• Then, the total mass in the largest remnant is

  $$\begin{eqnarray*} \frac{M_{\mathrm{LR}}}{M_{\mathrm{tot}}}= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}m_{\mathrm{atmos}}Q_{\mathrm{R}}+1 & 0\lt Q_{\mathrm{R}}\lt Q_{\mathrm{piv}} \\ \\ m_{\mathrm{c\&m}}\left(Q_{\mathrm{R}}-Q_{\mathrm{RD}}^{*(\mathrm{New})}\right)+0.5 & Q_{\mathrm{piv}}\lt Q_{\mathrm{R}}\lt Q_{\mathrm{supercat}}.\\ \end{array}\right. \nonumber\\ \end{eqnarray*}$$

  (26)

  We did not probe the supercatastrophic disruption regime in this study due as this would require much higher resolution; for collisions in this energy regime we recommend using the prescription of Leinhardt & Stewart (2012).
• Finally, the atmosphere fraction lost is

  $$\begin{eqnarray*} \chi _{\mathrm{loss}}^{\mathrm{atmos}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{-0.94^2}{4}\left(\frac{Q_{\mathrm{R}}}{Q_{\mathrm{piv}}}\right)^{2}+0.94\frac{Q_{\mathrm{R}}}{Q_{\mathrm{piv}}} & \frac{Q_{\mathrm{R}}}{Q_{\mathrm{piv}}}\lt 2.12 \\ 1 & \frac{Q_{\mathrm{R}}}{Q_{\mathrm{piv}}}\gt 2.12. \end{array}\right. \end{eqnarray*}$$

  (27)

5 SUMMARY

In this paper, we present the results from a series of sph simulations of head-on collisions of planets in which the target has a significant atmosphere. Our findings are summarized next:

• Giant impacts can have sufficient energy to remove large fractions of mass from the target planet; the mass lost is dependent upon the specific kinetic energy of the impact.

• Giant impacts can result in substantial increases in the densities of mini-Neptune planets by ejecting a fraction of their atmospheres.

• The fraction of mass lost splits into two regimes – at low specific impact energies only the outer layers are ejected corresponding to atmosphere dominated loss, at higher energies material deeper in the potential is excavated resulting in significant core and mantle loss.

• Approximately 20 per cent of the initial atmosphere remains at the transition between the two regimes.

• A single collision cannot remove all the atmosphere without also removing a significant percentage of mantle material.

• Mass removal is less efficient in the atmosphere loss dominated regime compared to the core and mantle loss regime.

• The specific energy of this transition (pivot energy, Q_piv) scales linearly with the ratio of projectile to target mass for all projectile-target mass ratios measured:

  |$Q_{\mathrm{piv}}=M_{\mathrm{tot}}\left(-2.45\frac{M_\mathrm{p}}{M_\mathrm{t}}+14.56\right) \,\,\,\,\big[10^6\, \mathrm{J\, kg^{-1}\, M_\oplus ^{-1}}\big]$|⁠.

• The fraction of atmosphere lost is well approximated by a quadratic in terms of the ratio of specific energy to transition energy: |$X_{\mathrm{loss}}^{\mathrm{atmos}} = \frac{-0.94^2}{4}\left(\frac{Q_{\mathrm{R}}}{Q_{\mathrm{piv}}}\right)^{2}+0.94\frac{Q_{\mathrm{R}}}{Q_{\mathrm{piv}}}$|⁠, for |$Q_{\mathrm{R}}\lt 2.12\, Q_{\mathrm{piv}}$|⁠, and total atmosphere loss for energies greater than this.

SUPPORTING INFORMATION

http://www.star.bris.ac.uk/TDenman/Paper1_Videos/

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol – http://www.bristol.ac.uk/acrc/. Thomas R. Denman acknowledges support from an STFC (Science and Technologies Facilities Council) studentship. Phil J. Carter acknowledges support from University of California Office of the President grant LFR-17-449059. Christoph Mordasini acknowledges support from the Swiss National Science Foundation under grant BSSGI0_155816 ‘PlanetsInTime’. Parts of this work have been carried out within the framework of the NCCR PlanetS (National Centre of Competence in Research PlanetS) supported by the Swiss National Science Foundation. This research has used the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

REFERENCES

APPENDIX:

FURTHER IMPACTS

We present results in this paper for collisions against three separate simulated super-Earth targets, a |$3.27\, \mathrm{M}_{\oplus }$| target with a |$0.27\, \mathrm{M_{\oplus }}$| atmosphere, |$1\, \mathrm{M_{\oplus }}$| core, and |$2\, \mathrm{M_{\oplus }}$| mantle (initial conditions and results given in Table A1); a |$6.26\, \mathrm{M}_{\oplus }$| one with a |$1.25\, \mathrm{M_{\oplus }}$| atmosphere, |$1.67\, \mathrm{M_{\oplus }}$| core, and |$3.34\, \mathrm{M_{\oplus }}$| mantle (Table 1); and a |$10.5\, \mathrm{M}_{\oplus }$| one with a |$3.43\, \mathrm{M_{\oplus }}$| atmosphere, |$2.36\, \mathrm{M_{\oplus }}$| core, and |$4.71\, \mathrm{M_{\oplus }}$| mantle (Table A2).