Collision Chains among the Terrestrial Planets. II. An Asymmetry between Earth and Venus

During the late stage of terrestrial planet formation, hit-and-run collisions are about as common as accretionary mergers, for expected velocities and angles of giant impacts. Average hit-and-runs leave two major remnants plus debris: the target and impactor, somewhat modified through erosion, escaping at lower relative velocity. Here we continue our study of the dynamical effects of such collisions. We compare the dynamical fates of intact runners that start from hit-and-runs with proto-Venus at 0.7 AU and proto-Earth at 1.0 AU. We follow the orbital evolutions of the runners, including the other terrestrial planets, Jupiter, and Saturn, in an N-body code. We find that the accretion of these runners can take $\gtrsim$10 Myr (depending on the egress velocity of the first collision) and can involve successive collisions with the original target planet or with other planets. We treat successive collisions that the runner experiences using surrogate models from machine learning, as in previous work, and evolve subsequent hit-and-runs in a similar fashion. We identify asymmetries in the capture, loss, and interchange of runners in the growth of Venus and Earth. Hit-and-run is a more probable outcome at proto-Venus, being smaller and faster orbiting than proto-Earth. But Venus acts as a sink, eventually accreting most of its runners, assuming typical events, whereas proto-Earth loses about half, many of those continuing to Venus. This leads to a disparity in the style of late-stage accretion that could have led to significant differences in geology, composition, and satellite formation at Earth and Venus.


INTRODUCTION
Earth and Venus are referred to as "sister planets," as they have a similar mass (Venus being only 15 % less massive) and bulk density. Numerical simulations of the formation of the solar system's terrestrial planets reproduce the formation of Earth and Venus analogs under a wide variety of initial conditions (Chambers 2001;O'Brien et al. 2006;Raymond et al. 2006Raymond et al. , 2009Hansen 2009;Fischer & Ciesla 2014). Yet Venus somehow ended up in a completely different dynamical state, rotating retrograde compared to the other planets and with a rotation period of 243 days (Campbell et al. 2019), with no known satellites (Sheppard & Trujillo 2009). The Earth-Moon system has more than 1000 times the angular momentum per unit mass than Venus.
Various explanations have been proposed for these differences, such as a despinning by tidal torques due to the planet's interior and its atmosphere (Correia & Laskar 2001;Auclair-Desrotour et al. 2017), or that it did not undergo giant impacts (Jacobson et al. 2017). The slow rotation of Venus and its lack of an internal magnetic dynamo may be consistent with the absence of giant impacts (Jacobson et al. 2017), or could indicate only headon giant impacts with low angular momentum. However, giant impacts are predicted to dominate the late stage of planet formation (Wetherill 1985;Kokubo & Ida 2002). The early accretion of planetesimals alone does not provide much net spin (Ida & Nakazawa 1990;Lissauer & Safronov 1991;Dones & Tremaine 1993), while giant impacts can strongly augment angular momentum (e.g. Agnor et al. 1999). There are also other formation models where there may have only been one giant impact, namely, the one responsible for the formation of the Moon (Johansen et al. 2021).
Giant impacts are collisions between similar-sized planetary embryos occurring at velocities v coll that are comparable to the mutual escape velocity v esc = 2G(m tar + m imp )/(r tar + r imp ). (1) Here m tar ≥ m imp are the masses of the target and the impactor, r tar and r imp the corresponding radii, and G is the gravitation constant. Late-stage collisions therefore occur at a few to 10 km s −1 for nominal terrestrial embryos, which is comparable to the sound speed in ices, oxides, and metals (Asphaug et al. 2015), resulting in global shocks that can generate magma oceans (Tonks & Melosh 1993) and trigger differentiation (see Kaula 1979;Rubie et al. 2015). Nonetheless, simulations show that the gross outcomes of giant impacts are governed primarily by gravitational forces and angular momentum (e.g. Canup 2004;Leinhardt et al. 2010). Late-stage giant impacts have already been invoked to explain features of our solar system as far-reaching as the crustal dichotomy of Mars (Wilhelms & Squyres 1984;Marinova et al. 2008), the high obliquity of Uranus (Slattery et al. 1992;Reinhardt et al. 2020;Ribeiro et al. 2020), and the origins of Mercury and the Moon (for reviews, see, e.g., Chau et al. 2018;Asphaug 2014). While the last giant impact around the Sun was probably 4.3-4.5 billion years ago (perhaps the Moon-forming collision itself), giant impacts in several nearby exoplanetary systems are thought to be recent or ongoing, responsible for observations of late-forming debris rings (Thompson et al. 2019) or planets with high eccentricity (>0.5; Frelikh et al. 2019).

The diversity of giant impact outcomes
The nature of the giant impact stage, and hence the geology and habitability of planets that are created by this process, depends on several factors. The starting conditions are represented by the orbital distribution of early-forming embryos (a.k.a. oligarchs), and their sizes and compositions are determined by disk physics (Kokubo & Ida 2000) and nebula chemistry and transport (Ciesla 2009). Equally important is the velocity distribution of the embryos, which evolves owing to self-stirring or forcing by giant planet perturbations, or Giant impacts occur at velocities that are generally faster than v esc , because in order for two planets to collide, their orbits must intersect. This requires a relative velocity v rel > 0, where v 2 coll = v 2 esc + v 2 rel . In classical self-stirring of a similar-sized population, relative velocities are excited to ∼ v esc so that impact velocities of order 1.4v esc are expected (Safronov 1969). Considering Earth and Venus formation, for an embryo in circular orbit starting at 1 AU to collide with an embryo at 0.7 AU requires v rel of at least 2.5 km s −1 . Lower-velocity collisions near v esc , including the so-called "graze-andmerge" collisions (Leinhardt et al. 2010) that lead to successful scenarios for Moon formation (e.g. Canup & Asphaug 2001), are only one outcome among others with similar or greater likelihood .
Faster giant impacts can also happen, but they are less probable because relative motions are damped by planetesimals and collisions. Also, ejection from the planet formation zone (loss of a participant) becomes more likely as relative velocities increase. Thus, most giant impacts occur in a rather sensitive region of velocity parameter space where accretion may or may not happen. Most N -body simulations of terrestrial planet formation indeed find that giant impacts occur in the velocity range v coll /v esc ∼ 1 − 2 with outliers that are faster (e.g., Agnor et al. 1999;Kokubo & Genda 2010;Quintana et al. 2016). Figure 1 shows the distribution of impact velocities for N -body simulations in Emsenhuber et al. (2020, hereafter E20). It should be noted that this set only contains 16 N -body simulations and that these all produce planetary systems less massive than the solar system's terrestrial planets. Also, the distribution only accounts for cases where the target mass m tar > 0.1 M ⊕ in order to provide a better comparison with the simulations that we present later in this work. The median value of v coll /v esc is 1.6, which is above the hit-and-run velocity threshold across expected impact angles (Kokubo & Genda 2010;Gabriel et al. 2020). These simulations have less dynamical friction than some other studies (Raymond et al. 2006;O'Brien et al. 2006;Chambers 2013;Kobayashi et al. 2019) because debris have been ignored, but they have more dynamical friction than studies that assume perfect merging without planetesimals. They have greater gravitational excitation than perfect merging simulations because they allow impactperiapsis 'runners' to escape, with the resulting large deflection. Moreover, in our calculations there is no inflation factor, so objects must come within the sum of the radii to detect a collision. The latter aspects tend to increase the rate of collisions early on, as well as their velocities, compared to other studies, but make it harder to ultimately accrete.
The N -body simulations by E20 determine the outcome of an individual giant impact from the key nondimensional parameters, which are the impact angle θ coll , scaled velocity v coll /v esc , mass ratio γ = m imp /m tar , and composition (Asphaug 2010;Genda et al. 2012;Gabriel et al. 2020) represented by the metallic core radius relative to the planet radius. For a differentiated chondritic sphere with ∼30 wt% iron and 70 wt% silicate mantle the core is half the planet's radius.
Initial rotation, normalized to the spin-disruption limit, is another important parameter, and while it has been included in prior studies (Canup 2008) and special cases close to the breakup limit (Ćuk & Stewart 2012), the free parameters add up quickly, especially in three dimensions, making it difficult to run a systematic study at adequate numerical resolution. The N-body simulations of E20 are therefore based on the assumption that the bodies are not spinning prior to each collision, and rotation is ignored.
To obtain a sufficiently accurate working model of giant impacts, for the regime of nonrotating chondritic bodies and in the mass range of 10 −3 and 1 M ⊕ , the Nbody simulations in E20 used the surrogate model of gi-ant impacts (Cambioni et al. 2019, hereafter C19). This is a set of machine-learning algorithms trained on 800 smoothed particle hydrodynamics (SPH) simulations of giant impacts, varying the impact angle, velocity, and mass ratio within the limits defined by the ranges explored in building the dataset. Sparsely sampled parameter spaces such as these can be used by machine learning to generalize the underlying model. C19 have used machine learning to classify the outcomes (disruption, accretion, graze-and-merge, hit-and-run) and to develop a neural network that accurately predicts the mass of the largest remnant m lr , or specifically the accretion efficiency within the parameter limits of the database. E20 have also developed neural networks in the case of hit-andruns that predict the masses and velocity vectors of the target and runner. As mentioned above, E20 have not included the fate and influence of debris produced by giant impacts, which are defined as anything smaller than the runner, objects that are generally orders of magnitude lower in mass. Over the course of the late stage, the cumulative debris can represent a substantial fraction of the total mass (E20) and reduce the eccentricities and inclinations of the embryos.

The perfect merging assumption
From a perspective of planetary formation, the variety of giant impact outcomes described in the previous section would not matter so much if the leftover material is reaccreted later by the same target. This assumption can seem justifiable, because whenever two planets emerge from a hit-and-run collision they may be expected to experience a follow-on collision, so that merger seems to be a foregone conclusion. Most N -body simulations of terrestrial planet formation (e.g. Chambers 1999;O'Brien et al. 2006;Fischer & Ciesla 2014) have traditionally used this approach and treated impacts as perfect mergers, ξ = 1 for collisions that are not catastrophically disruptive.
However, Emsenhuber & Asphaug (hereafter Paper I 2019a) showed that this is not generally the case. They studied the fate of the runner following hit-and-runs into proto-Earths at 1 AU, for thousands of geometries, and found that, contrary to expectation, only about half the time (depending on the runner's egress velocity, which depends on the impact velocity and angle) do they return to collide again with proto-Earth. When they do, the return collision happens on a timescale of thousands to millions of years.
In this sense, the assumption of perfect merging would lead to a poor estimate of the accretion timescale, and the inferred thermodynamic (e.g., lack of cooling between collisions in a chain) and differentiation history ) would be unrealistic. Furthermore, Paper I found that a majority of those runners that do not return to proto-Earth are likely to collide with Venus, assuming the present masses and orbits of the planets. They also showed that for returning runners, the impact velocity of the second collision is usually similar to the egress velocity following the hit-and-run, which is slower than the original impact owing to momentum loss. So, the follow-on collisions tend to be comparatively slower than the previous hit-and-run encounter. The offset angle between impacts, however, is uncorrelated and random, meaning that the returning collision is off-axis by about 90°on average.

Fate of the runner
In most hit-and-runs, a large part of the projectile survives the giant impact. The velocity of the runner can be strongly reduced in magnitude and deflected in direction relative to the inbound velocity (Genda et al. 2012). If the runner is slowed so that the two bodies are bound gravitationally after the collision, then instead of a hit-and-run it is a graze-and-merge collision (e.g., Leinhardt et al. 2010), the archetype being the impact origin of the Moon (Canup & Asphaug 2001). The boundary between hit-and-run and graze-and-merge occurs when the runner gets to about two-thirds of the Hill radius, as studied by Emsenhuber & Asphaug (2019b). This leads to a sensitivity around the boundary between hit-and-run and merger.
Hit-and-runs can be further subdivided into lowvelocity hit-and-runs (that is, not much faster than the escape velocity), where the runner is an identifiable remnant of the projectile (e.g., a mantle-stripped core, sometimes barely so) and high-velocity hit-and-runs where the runner is disrupted or dispersed (Asphaug et al. 2006). Even in dispersive hit-and-runs the target remains mostly intact; disruption of the target by giant impact requires much greater energies than are considered here, 3 to 5 v esc or more. We ignore these faster collisions for now, for several reasons: they are much less probable (Figure 1 of this work and, e.g., Agnor et al. 1999), their debris reaccumulation is not a single event and thus harder to constrain, and the surrogate model can be used with highest fidelity at velocities lower than 3 v esc .
Paper I modeled a series of dynamical evolutions of the two largest remnants (target and runner) after relatively low velocity hit-and-runs with the proto-Earth, an 0.9 M ⊕ planet at 1 AU. The authors found the following: • Between one-third and two-thirds of the runners collide back onto proto-Earth, depending on the egress velocity, after an interlude of thousands to millions of years.
• Of those that do collide with another body, most end up on Venus. Indeed, it was found that the destination of moderate-velocity hit-and-runs with proto-Earth is almost as likely to be Venus as proto-Earth itself.
• Venus-bound runners following hit-and-runs with proto-Earth have impact velocities that are faster than their counterparts having return collisions with proto-Earth, implying that accretion might not end with this event, leading to longer collision chains.
This motivates us to explore more generally the fate of runners and their exchanges between Venus and the Earth, where our aim is to identify systematic differences in the kinds of late-stage collisions they each experience, which might account for the major differences in the geophysical and dynamical states of these otherwise-similar "sister planets."

This work
Here we explore how the capture, loss, and interchange of runners affects the growth of Venus compared with the Earth. Collision chains start with a hit-and-run, potentially involve multiple subsequent hit-and-runs, and end in terminal mergers (Paper I). To study the general process, we need to model each collision, especially in the case where there may be one or more intermediate hit-and-runs. To make this analysis tractable, we cannot simulate every giant impact explicitly, and thus we utilize the work of E20 who extended the machinelearning approach, training on the same simulation data as C19, Reufer (2011), andGabriel et al. (2020) to retrieve not only the accretion efficiency (mass added to or removed from the target) but also the properties of the two main remnants, in the case of hit-and-run, and their relative orbits (i.e. egress velocities). The result is a machine-learning-derived surrogate model for giant impacts, providing a functional map of input parameters into outputs that allows us to realistically model collisions as they happen during the N -body evolution.
Using this machine-learning surrogate model, it is possible to perform N -body evolution that can continue from the initial hit-and-run through the subsequent collisions, while improving the fidelity in the dynamics of the system compared to other models. This allows us to extend the procedure of Paper I to continue the dynamical evolution past the first subsequent collision, to model collision chains lasting one, two, three, or more hit-andruns until the terminal merger, using the procedure described in C19 and E20. We are able to determine the ultimate destination of the runner and the number of intermediate collisions required until final accretion.

METHODS
Our procedure is similar to Paper I, but several improvements have been implemented. The first collision of the chain (C 1 ) is modeled directly with SPH to simulate several candidate low-velocity hit-and-runs into proto-Earth and proto-Venus at high resolution. The results (masses and velocities of the two largest remnants) are then transferred into an N -body code to track the dynamical evolution of proto-Earth or proto-Venus and the runner.
For each first collision, we perform 1000 N -body evolutions for random collision orientations to obtain a statistical description of the fate of the runner. During the N -body evolution, subsequent collisions (C 2 . . . C F , for collisions 2 through F) are treated using surrogate models following the methodology of E20. Each N -body evolution ends when the runner is lost (by accretion onto another body or ejection) or the maximum time elapses (generally 50 Myr, but extended to 400 Myr for some cases in Section 3.5).

Initial collision simulations
To model the initial collisions, we use an SPH code that has been developed for modeling similar-sized collisions in the size and velocity regime of late-stage terrestrial planet formation (Reufer et al. 2012;Emsenhuber et al. 2018; Paper I), with targets of ∼0.1 M ⊕ and larger. For these size bodies, strength effects can be ignored, so the bodies are treated as fluids. SPH is a Lagrangian technique with material subdivided into mass particles. A kernel interpolation is used to compute hydrodynamic quantities at any location. Spatial derivatives are computed using an interpolation with the derivatives of the kernel. Time evolution is provided by the Euler equations, except for the density, which is retrieved using the kernel interpolation with a correction term for particles close to the surface (Reinhardt & Stadel 2017). The pressure p(ρ, S) and the other physical quantities necessary for the hydrodynamical equations are retrieved using a tabulated form of the M-ANEOS equation of state (Thompson & Lauson 1972;Melosh 2007). To compute self-gravity and retrieve neighboring particles, a hierarchical spatial tree is used (Barnes & Hut 1986). For reviews on SPH, the reader is referred to, e.g., Monaghan (1992) and Rosswog (2009).

Thermodynamic equation
Energy conservation in the Euler equations is provided by the energy equation. In its standard form, the equation is simply the adiabatic compression or expansion of the material. In SPH, it is modified (as is the momentum equation) for the addition of the artificial viscosity to handle shocks (see below). Nevertheless, in fluid SPH (i.e., without additional solid forces) the largest contribution to the changes in internal energy comes from adiabatic compression or expansion.
To avoid the numerical integration of the energy equation for an adiabatic process, which can lead to nonphysical states (Reinhardt & Stadel 2017), we use the first law of thermodynamics to determine an entropy equation, where S is the specific entropy, T is the temperature, u is the specific internal energy, p is the pressure, V = 1/ρ is the specific volume, and the ∂u AV /∂t term denotes the contribution of artificial viscosity to the energy equation.
A drawback of this formulation is that the integration of the entropy equation introduces additional imprecision for the artificial viscosity. This can lead to the total energy suffering from a slightly larger drift overall, which is not expected to affect the results substantially.

Artificial viscosity
Artificial viscosity is a nearly ubiquitous numerical tool for resolving shocks in SPH. The standard formulation that was used in previous studies suffers from drawbacks, such as being triggered in cases where no shocks occur, for instance, in shear motions. Here, we use a form that is inspired by Riemann solvers (Monaghan 1997). The acceleration term is computed as where i and j are particle indexes, m j is the particle mass, ρ ij = (ρ i + ρ j )/2 is the averaged particle density, v ij = v i − v j the relative velocity of the particles,r ij = r ij /r ij with r ij = r i − r j the unit vector along the relative position of the two particles, and ∇ i W ij is the gradient of the kernel W (r ij , h ij ) with respect to r i . v sig represents the maximum signal velocity between particles i and j and can be estimated as where c s is the sound speed from the equation of state (Monaghan 1997). The corresponding energy change is given by In SPH simulations, artificial viscosity can trigger in unwarranted circumstances, e.g., during shear motion in accretion disks and potentially in hit-and-runs. A simple method to reduce the artificial viscosity is the "Balsara switch" (Balsara 1991) to apply artificial viscosity only during compression. This can, however, cause problems with shocks in accretion disks (Owen 2004) and artificial particle alignments. So, in our case, we have decided instead to use time-dependent viscosity factors following Morris & Monaghan (1997). Here the α parameter of the artificial viscosity is treated as another time-integrated quantity of each particle. Its derivative is given by with being the timescale over which α i is reduced and ξ = 0.1. This is set so that the α i factor is reduced after shock passed by about 10 smoothing length. In our case, we set α min = 0.1 and α max = 1.5.

Initialization
We construct the projectile and target bodies using an improved methodology over that of Paper I. The first step is to obtain a 1D radial hydrostatic profile using the scheme of Benz (1991). With this profile an initial SPH body is generated using the methodology described in Reinhardt & Stadel (2017). This body is made of successive layers of SPH particles, whose properties are taken from the 1D radial profile. The locations of the particles on each layer are obtained using the HEALPix (Górski et al. 2005) software package. The location of SPH particles in each layer is found iteratively with the goal to minimize the difference of horizontal and vertical spacing between the particles. The resolution is chosen so that the number of particles is proportional to the body's mass, with a 1 M ⊕ body being represented by 500 000 SPH particles. This results in initial SPH bodies whose particles have a lower residual velocity. Nevertheless, we further evolve these bodies for 6 h of physical time, which further reduces the rms of the residual velocities to about 2.5 m s −1 , or 2.5 × 10 −4 of the escape velocity in our Earth-scale planets.

Evolution until Subsequent Collision
The largest remnants of the hit-and-runs are identified following the SPH simulation using the methodology of E20. These bodies, the target and runner, are then mapped into the mercury6 N -body code (Chambers 1999) to follow their dynamical evolution. We in-clude the evolution of the other major planets as well, assuming their present orbits, which is important considering our result that planets can often exchange runners. The presence of Jupiter and Saturn on their presentday orbits is consistent with them either forming eccentric (e.g., Woo et al. 2021) or having migrated early in the process of terrestrial planet formation (e.g., Clement et al. 2018).
For now, the lesser debris produced by the hit-andrun is ignored in the further evolution, and in any case for our simulations at most 2.4 % of the material is in objects smaller than the runner (m lost as reported in Tables 1 and 2). We proceed as in Paper I: for each hydrodynamical simulation, we perform 1000 realizations of the dynamical evolution, each with a different orientation of the pre-impact orbit of the impactor. The target is assumed to be on a circular orbit, while the orbit of the impactor is computed according to Jackson et al. (2018). To obtain the orientation, we use a Monte Carlo approach, and we assume that the orientation of the relative orbit follows a uniform distribution in space (Paper I).
The criterion for the end of the dynamical evolution has been updated following the improvement in our treatment of collision handling in the N -body evolution, described below. Using surrogate models for giant impacts, we are now able to continue the dynamical evolution until the runner has been fully accreted, which can require successive collisions, or until a predetermined time that we set to 50 Myr, as justified in Section 3.5, with several cases extended further to 400 Myr.

Collision model
The underlying methodology of the collisional handling during N -body evolution is described in C19. A data set comprising about 800 SPH simulations that were obtained by Reufer (2011) was used to performed machine learning (data available in Gabriel et al. 2020). The simulations span a range of target masses m tar from 10 −2 to 1 M ⊕ , impactor mass ratio γ = m imp /m tar from 0.1 to 0.7, and impact velocity ratio v coll /v esc from 1 to 4. This dataset is well suited for this study, as it encompasses the parameter range of targets and runners from our SPH simulations. A classifier was trained on the dataset to determine whether a collision is in the hit-and-run regime (two remnants) or not (single remnant). Then, two neural networks were also trained to obtain one regressor that provides the mass of the two largest remnants and another regressor that provides the orbital characteristics of the two remnants in the case of a hit-and-run (E20). The classifier and neural networks were implemented in the collresolve 1 library (Emsenhuber & Cambioni 2019) that we are using in this work. The full procedure for the treatment of collisions and its adaptation to the mercury6 code are discussed in detail in E20.
Additionally, we determine the mass exchange between the two bodies in a similar way that is computed for the Earth-disk equilibration in simulation of potential Moon-forming giant impacts (Reufer et al. 2012): where f l←t mant and f s←t mant are the mass fraction of the mantle of largest and second remnants that are coming from the target, respectively, and f l←i mant is the mantle mass fraction of the largest remnant coming from the impactor. These results for mass exchange are not important for the present study but are applied to the third paper in this series, on the origin of the Moon (Asphaug et al. 2021).
A value of −100 % indicates that the mantles do not mix at all during the collision, i.e., the mantle of the largest remnant is made exclusively of target material, while the mantle of the second remnant is made exclusively of impactor material. A value of 0 % indicates that each remnant has a mantle made of the same proportion of target and impactor material. We assume that this exchanged material becomes well mixed after the collision, within each core and mantle, respectively. The collisions where the two remnants remain gravitationally bound after the collision (graze-and-merge) are analyzed using the procedure outlined in Emsenhuber & Asphaug (2019b).

Initial Hit-and-run Collisions with Proto-Earth
For the SPH simulations of the initial hit-and-runs with the proto-Earth, we start with bodies whose masses are m tar = 0.9 M ⊕ and m imp = 0.15 M ⊕ . The impactor is somewhat smaller than in Paper I; this is to end up with an approximately Mars-sized runner, which would be the suitable size for the scenarios of Moon formation (Piet et al. 2017), which we consider in the next paper in this series (Asphaug et al. 2021). In addition, we perform new SPH simulations of hit-and-runs more relevant to Venus, with m tar = 0.7 M ⊕ and m imp = 0.15 M ⊕ .
The results are provided in Table 1 for the collisions with the proto-Earth and in Table 2 for those with proto-Venus. All the initial bodies are nonrotating, which was 1 https://github.com/aemsenhuber/collresolve chosen for two reasons. First, the most common orientation is for impact orientation and spin to be perpendicular, which provide little angular moment. The second reason is that pre-impact spin affects the outcomes less than the other effects we study here (Timpe et al. 2020). We provide several new outputs of the simulations compared to Paper I that are related to the SPH modeling of possible return collisions. P lr is the rotation period of the largest remnant, computed as P lr = 2π/ω lr , with ω lr = L lr /I lr , L lr the spin angular momentum, and I lr the moment of inertia along the direction of L lr . Also, Z lr and Z sr are the core mass fractions of the largest and second remnants.
For the new analyzed quantities from the SPH collisions (i.e., P lr and δf T ), we see that the rotation period and mixing factor of the largest remnant are mostly dependent on the impact angle, with limited effect of the velocity. It should be noted that all the rotation is induced by the collision as the initial bodies are nonrotating. The core mass fraction of the largest remnant is not much affected by the collision, compared to the original target. It usually slightly decreases because there is often some accretion of impactor mantle material in a hit-and-run (Gabriel et al. 2020). In some cases, though, for common impact angles (42.5°and 43.0°at v coll /v esc = 1.20 in the case of collisions with proto-Earth) it increases by a very small amount. The runner always ends up with a higher core mass fraction than the impactor, which is consistent with previous simulations (Asphaug & Reufer 2014;Chau et al. 2018) and planetary differentiation studies ). Usually, steeper impact angles result in the largest core mass fractions, but this is not always the case. For instance, for collisions with the proto-Earth at v coll /v esc = 1.20, the lowest core mass fraction is found for the collision at 45°; the core mass fraction increases again in the simulations with 50°and 55°.

The ultimate destination of the runner
Here, we analyze the N -body evolutions for the fasterthan-v esc remnants of the collisions described in Section 3.1. As in Paper I, we have included the other planets on their present orbits, out to Saturn. In these calculations we use the surrogate model of Cambioni et al. (2019) to extend the calculation until the runner is lost, either by accretion or by ejection, whereas in Paper I the N -body evolutions were only integrated up to the first encounter. To compare the demographics of the results of the two methods, we provide in Tables 3 and 4 the number of evolutions comparing the bodies encountered on the next collision (C 2 ; shown with variables  refer to the orbits of the bodies after the collision. P lr is the rotation period of the largest remnant, Z lr is its core mass fraction, while Zsr is the core mass fraction of the second remnant. f l←i core and f s←t core are the core mass fraction of the largest and second remnants coming from the impactor and target, respectively. f l←i mant and f s←t mant are similar, but for the mantle. δfT is the mantle equilibration factor, Eq. (9).   Table 1. N 2 ) versus the bodies onto which the runner was finally accreted (C F ; shown with variables N F ).
For the hit-and-runs that we modeled, the runners are more massive than Mercury and Mars, so accretion with those planets is actually a case of Mercury or Mars being accreted by the runner. In this kind of situation the final result is therefore counted as a survival of the runner (unless it has subsequently been accreted by another more massive planet). In a smaller number of cases the runner is accreted by these planets after it has had an intermediate hit-and-run with another planet, which resulted in mass loss such that the successor of the runner is less massive than Mars or Mercury at the moment of the collision. To show the evolution of the runner's mass by intermediate collisions, we present in Figure 2 the cumulative distribution of their masses at the moment of the final collision (C F ) for the simulations. To focus on the processing by intermediate collisions, only Note-The N 2 columns denote the number of evolutions where the runner collides with the specified body; N 2 None indicate evolutions where the runner did not collide with any body during the entire evolution. The N F columns denote the number of evolutions where the runner has been accreted by the specified body; N F None indicate the number of evolutions where the runner survived until the end of the evolution.  Figure 3) and their mass is identical to that of the final bodies from the SPH models. The results show that most runners undergoing intermediate collision lose mass, which is consistent with hit-and-run collisions generally resulting in the transfer of some material from the impactor onto the target. As seen in Paper I, we obtain that the outcomes are principally determined by the relative velocity after the initial collision. To simplify the presentation of the results, we now select four sets of dynamical evolutions that we will discuss in more detail during the remain-der of this work, two from each set of collisions with proto-Earth and proto-Venus. For each target planet we select one hit-and-run with a slow runner (with v coll /v esc = 1.15 and θ coll = 47.5°) and one with a faster runner (with v coll /v esc = 1.20 and θ coll = 55.0°). The former case was selected because it is the lowest-velocity collision (for the given angle range) where hit-and-run occurs (at lower velocities the runner has an egress velocity less than the mutual escape velocity, resulting in a merger). The latter case is the fastest runner that we modeled. The results of the other sets of dynamical evolution usually lie in between the two that were selected. A graphical representation of the results is provided in Figure 3, where the inner ring corresponds to the destination of the first collision (C 2 ) following the hit-andrun (which is C 1 ), and the outer ring is the destination of the final accretion C F (the last collision of the chain). Black indicates a runner with no subsequent accretion during the 50 Myr of our simulations.
The outcome of the first collision is a relatively good predictor of the destination of the final accretion, because most chains end with C 2 . Another common outcome is that the percentage of runners that survive through the end of our simulations (the black section of the outer circle) is larger than the percentage of those that did not undergo any subsequent collision (the black section of the inner circle). This is due to hit-and-runs whose runner survives to the end of evolution, as it happens for those that do not have any collision during the dynamical evolutions.
Comparing similar runners emerging from collisions with proto-Earth and proto-Venus, we find that the ones emerging from proto-Venus are more likely to return and be finally accreted by Venus than vice versa for proto-Earth. Consider our slow runner examples; in these cases about 3/4 of runners are reaccreted by Venus, while 11.3 % of runners emerging from proto-Venus are accreted by proto-Earth. For the reverse scenario, the percentage of runners from proto-Earth that  . For this sample, Venus accretes nearly all of its slowest runners, and Earth accretes a bit more than half. Venus accretes about half of its faster runners, and a quarter go to Earth. Earth accretes only a quarter of its faster runners and Venus another quarter, for this example.
are accreted by the Earth is significantly smaller, while the probability of ending up at Venus is significantly larger (18.8 %). For the faster-runner case, the fraction that are exchanged between the planets increases to 24.5 % and 24.1 % respectively, as discussed further below. We show four systems from the series of evolutions with the "faster" runner emerging from proto-Venus (v coll /v esc = 1.20 and θ coll = 55.0°) in Figure 4. Each panel represents a common outcome of our evolution, with a runner accretion coming either directly (panel (a)) or after a hit-and-run onto the proto-Earth (panel (c)). Panel (b) shows a case where the runner gets accreted by proto-Earth. Finally, panel (d) shows a runner surviving through the end of the N -body evolution, close to the position of Mars after a hit-and-run with the proto-Earth. We see that in most cases of early return on the same planet (akin to panel (a)), the runner remains in the orbital vicinity of the planet between the collisions. However, the runner may wander to a more distant location than what could be implied by only looking at the collisions, with the example of panel (d).

Debris production
In reality, all collisions generate some level of escaping material, which we do not explicitly track dynamically, and the total amount can be a substantial fraction of the initial planetesimal population mass. To estimate the amount of mass that is lost in this way, we provide in Figure 5 the cumulative distributions of the total amount of debris produced for the same four dynamical evolution sets presented in Figure 3. These include the debris produced in the initial collision, which are shown at the left end of each curve.
In all four sets, the mass distributions show that some evolutions do not produce debris at all. For the evolutions with a low relative velocity after the collision (labeled as slow on the figure), this includes cases where  the runner is reaccreted without the production of any debris. This happens most for collisions with a velocity of v coll ≈ v esc and impact angle around 45°. In this region of the parameter space, the surrogate collision model indicates that no debris are produced (right panel of Figure 5 in Emsenhuber et al. 2020). For the fasterthan-v esc evolutions, this is due in part to evolutions that do not produce any subsequent collisions.
For the other cases, the amount of debris produced is generally between 10 −2 and 10 −1 M ⊕ , that is, a few lunar masses. In most cases, the masses are much smaller than the planetary masses; hence, debris reaccretion is generally not a problem here. Only in a few cases do we obtain a total debris mass of 10 −1 M ⊕ or more, and this would be greater for more energetic collision chains.
Neglecting the debris in our simulations can affect the simulations in two ways. First, debris affect the orbits of the planets (Kobayashi et al. 2019). Small debris can provide dynamical friction, resulting in lower collision velocities. Second, debris reaccretion affects the boundaries of the different collision regimes. This effect depends on how reaccreted debris are distributed between the remnants. Proportionally more accretion on the largest remnant results in more dissimilar bodies, which would reduce the amount of hit-and-run collisions, and vice versa. It should nevertheless be noted that reimpact by debris need not result in net accretion. There is no simple prescription to determine how debris affects the evolutions; taking all their effects into account is a difficult task that needs to be treated carefully.

Length of collision chains
To understand the number of hit-and-runs that precede an eventual accretion, we provide histograms for the length N of the collision chains in Figure 6 for the same four scenarios we presented in Figure 3. In all cases, the most common outcome of these slow hit-andruns is a follow-up collision that accretes the runner (C 2 ). We observe that fast runners tend to have longer chains (more follow-up collisions). The reason for that is linked to our findings from Paper I, that is, the median impact velocity of a follow-up collision is similar to its egress velocity from the previous collision. Thus, fast runners will naturally lead to faster C 2 collisions, which are more likely to result in hit-and-runs. As noted, the preponderance of single-link chains (N = 2) observed in this analysis accounts for the similar distribution between the initial collision and final accretion observed in Figure 3.
We observe that the fraction of runners that are accreted at each link of the chain decreases. For example, the proportion of surviving runners in the second chain (C 2 ) in Figure 6 (left) is small (21 % for the slow chain and 46 % for the faster one). Survival becomes the dominant outcome by the fourth collision in the chain (C 4 ) (61 % and 63 %, respectively). This is, in part, due to the natural bias that for any subsequent collision to result in hit-and-run, it must happen at sufficiently high velocity to overcome the merging velocity threshold.
Another, less straightforward effect arises from the gravitational stirring that occurs between each chain link (as the runner evolves dynamically before reimpacting). As described in Paper I, the median followup impact velocity after a hit-and-run is near the egress velocity; however, a tail end of the distribution extends to much higher velocities from N -body dynamics. For example, in the low-velocity hit-and-run chains for proto-Venus ( Figure 6, right), the median value for impact velocities in C 2 is v coll /v esc = 1.01. However, the survivors (hit-and-runs) of C 2 have a median impact velocity of v coll /v esc = 1.21, biasing the subsequent collisions to higher velocities.
This bias is compounded at each link of the chain, with the median survivor velocity increasing each time.
Longer chains thus increasingly sample less probable, high-velocity scenarios. The mutual dynamical evolution between chains, which is inherently chaotic, will introduce variability in impact angle and velocity, producing more diverse outcomes in longer chains.  Tables 3 and 4. For each collision, the left and right bar indicates the 'slow' and 'fast' runner case respectively. Bars that have at least 1 but less than 10 items have the total number shown above them.

Time until loss of the runner
The previous section discussed the differences in the collision chains starting at proto-Earth and proto-Venus, beginning with two sets of hit-and-runs. Earth forma-tion involves somewhat longer chains, for example. It is therefore natural to also compare the timescales for hit-and-run chains to resolve either by accretion or by ejection. In Figure 7 we show the cumulative resolution time for collision chains. For this figure, the four sets of dynamical evolutions were extended beyond where we computed our collision/chain statistics (at 50 Myr), out to 400 Myr. These evolutions are much longer than those in Paper I, as the concern of that work was only about a single follow-up collision.
Through this comparison of the cumulative results of the four sets of simulations, we find that there are distinct differences in the temporal length of hit-and-run chains between the Venus and Earth dynamical zones. At any given epoch, ∼10-40 % more of the hit-and-run chains at ∼0.7 AU (our proto-Venus scenario) tend to have concluded as compared to those at Earth. This is an expected consequence of the shorter collisional timescales of the inner solar system, as well as the lower mass of proto-Venus making merger less probable for a given giant impact. We find that the time required for a collision chain to resolve can be comparable to estimates for the duration of late-stage solar system formation (∼200 Myr; Chambers 2013; Raymond et al. 2018). The median time for a chain to conclude in the case of a slow runner at 1 AU (Earth) is 1 Myr. In contrast, the median time for a chain to conclude in the slow-runner case at proto-Venus is only 0.1 Myr. These differences are much larger than the ratio of the orbital periods.
Another important factor is that different collisionaldynamical outcomes of a hit-and-run chain dominate at different epochs. Chain resolution at early times is often due to accretion of the runner, whereas at late times chain resolution becomes increasingly dominated by ejection. For instance, there are no ejections or accretions by the Sun before 2 Myr. We also find that there are no runners that swap from Earth to Venus or Venus to Earth before 10 5 yr, so very early chain resolution is dominated by accretion on the initial body. Despite this systematic behavior as a function of simulation time, we find that by 50 Myr the statistics of chain outcomes do not change considerably, justifying our use of 50 Myr for compiling the probabilities in Tables 3 and 4

Exchange of runners
In addition to comparing the collisional-dynamical timescales of hit-and-run chains in different dynamical zones, we can also examine the exchange of runners between the two zones. In the pathway diagrams of Figures 8 and 9 we show collision chain networks for the same four scenarios examined previously: low impact velocity hit-and-run chains beginning at Venus and Earth (Figure 8, left and right, respectively), and higher impact velocity hit-and-run chains beginning at Venus and Earth (Figure 9, left and right, respectively). Numbers are per 1000 cases. The patterns that were obtained in the previous sections are reflected here, especially that low impact velocity runners emanating from proto-Venus have a higher rate of return to Venus, with a small fraction transported to Earth, i.e. Venus holds on to its runners. In contrast to this, almost twice as many low impact velocity runners leaving proto-Earth end up colliding with Venus. As expected, the higher-velocity runners show an even greater proportion of transport between the dynamical zones and a greater probability of a second return collision, although, as before, Venus is much more efficient than Earth at retaining its runners.
In all cases, only a small overall fraction of runners ( 10 %) survive past the first follow-up collision, independently of the target of that collision. Importantly, we also notice that the total fraction of second return collisions in the hit-and-run chains originating at Venus changes considerably with initial impact velocity (3 %-10 %), whereas hit-and-run chains originating at Earth are less sensitive to the initial impact velocity (increasing from 6 % to 9 %). This indicates that further study is needed for hit-and-runs in the velocity range ∼ 1.3 − 1.6v esc and faster.
We note an evolution in the dynamical character of runners that survive subsequent collisions. A larger percentage of the runners that survive C 2 do not collide later with either Earth or Venus, compared to the original runners egressing from C 1 . Out of the 42 slow runners that survive C 2 with proto-Venus (left panel of Figure 8), only 22 of them (52 %) were found to have a subsequent collision with proto-Earth or proto-Venus, compared to 91 % of the runners from C 1 . As we discussed in Section 3.4, the median impact velocity of each subsequent collision increases, which results in a greater likelihood of interacting with other planets, decreasing the likelihood of re-impacting the initial target.
For Earth, there is a greater overall transfer of runners compared to Venus and a likelihood of longer chains. We therefore find it more probable that runners from proto-Earth would end up in the inner reaches of the solar system than vice versa. Indeed, as can be seen in Figure 9, higher-velocity hit-and-runs into proto-Earth are less likely to resolve in mergers at Earth than they are to resolve at Venus, the next planet in. The opposite is true for proto-Venus, which tends to resolve chains at Venus in both the slow and faster scenarios.
In summary, we find that in the case of hit-and-runs with Earth at 1.0 AU with a departing runner velocity v dep 1.1v esc , the runner is about as likely to collide with Venus as it is to return to Earth. The converse is not true: runners from hit-and-runs into proto-Venus, at least for the "slow" hit-and-runs we have studied, have a significantly greater likelihood of returning to Venus than making it to Earth. Venus holds onto its runners and is a sink for runners whose planetary parent bodies began farther our in the solar system.

CONCLUSION
In this work, we aim to understand whether the geophysical differences between Venus and Earth can be potentially explained by a systematic difference in the giant impacts that formed them. To study this, we perform dynamical evolutions of remnants of hit-and-run collisions until the runner is finally accreted or ejected. The runners may experience subsequent hit-and-runs, a collision chain C 1 , C 2 , ... The first C 1 is modeled using our updated SPH code, and subsequent collisions in the chain are represented in an N -body routine (Chambers 2012) using the surrogate giant impact model described in Cambioni et al. (2019) and the methodology of Emsenhuber et al. (2020). The starting collisions take place at 1.0 AU and 0.7 AU, respectively.
We arrive at four main conclusions: • The terrestrial planets are not isolated during the late stage of planetary formation. Runners emerging from hit-and-runs with one planet are likely to collide with another, with varying probability Figure 8. Fate of the "slow" runner case, up to the second return collision. The left side is for runners emanating from hit-and-run with proto-Venus while the right side shows runners emanating from hit-and-run with proto-Earth. As an example, among those runners that escape accretion in the initial collisions on Venus, 93 runners are lost, 787 runners re-impact with Venus, and 120 runners impact with Earth. Among those runners that impact with Earth and manage to escape again, 7 runners are lost, 12 runners resolve in mergers at Earth, and 16 runners resolve in mergers at Venus. depending several factors, such as relative velocity and orbital configuration.
• Long chains are less probable because they require larger initial velocity, and higher-velocity runners are less likely to return.
• Earth serves as sort of a vanguard for Venus, capable of slowing late-stage projectiles down in this manner, yet not accreting more than about half of them itself.
• Earth runners end up at Venus with the same likelihood that they return to Earth. Venus, however, retains the majority of its runners in all situations we studied, effectively serving as a sink in these scenarios.
Moderate-velocity hit-and-run collisions are akin to dissipative interactions, acting to slow down impacting planets. If the egress velocity is just above the mutual escape velocity, hit-and-run collisions can lead to subsequent accretionary collisions, as the collision velocity of C 2 is usually close to the egress velocity from C 1 (Paper I). But they can also act like a headwind for the runners, leading to a reduction in energy and an overall inward shift of their orbits.
So, if the terrestrial planets formed in multiple giant impacts, then Venus is significantly more likely than the Earth to have accreted a massive outer solar system body during the late stage of planet formation. The Earth, by contrast, has no terrestrial planet beyond its orbit to act as a vanguard. Mars is about the same mass as the late-stage projectiles (Piet et al. 2017), 0.1 M ⊕ , and thus relatively inconsequential in terms of slowing them down through hit-and-run, so Earth has to do it on its own. These demographic differences have broad implications for Earth and Venus formation specifically, and terrestrial planet formation in general, that require further exploration. The accretion of most runners by Venus, and fewer by Earth, and the handing-off of faster runners from Earth to Venus are a newly identified aspect of their late-stage formation that would influence Venus's bulk composition and could potentially lead to substantial differences in their final spin states, core-mantle dynamics, and satellite formation which is the subject of the next paper in this series (Asphaug et al. 2021).
A.E., E.A., S.C., and S.R.S. acknowledge support from NASA under grant 80NSSC19K0817 and the University of Arizona. T.S.J.G. acknowledges support from the Arizona State University Space Technology and Science ("NewSpace") Initiative. We thank the anonymous reviewers, whose comments helped to improve the manuscript. An allocation of computer time from the UA Research Computing High Performance Computing (HPC) is gratefully acknowledged.