The evolution of collision debris near the $\nu_6$ secular resonance and its role in the origin of terrestrial water

This work presents novel findings that broadens our understanding of the amount of water that can be transported to Earth. The key innovation lies in the combined usage of Smoothed Particle Hydrodynamics (SPH) and $N$-body codes to assess the role of collision fragments in water delivery. We also present a method for generating initial conditions that enables the projectile to impact at the designated location on the target's surface with the specified velocity. The primary objective of this study is to simulate giant collisions between two Ceres-sized bodies by SPH near the $\nu_6$ secular resonance and follow the evolution of the ejected debris by numerical $N$-body code. With our method 6 different initial conditions for the collision were determined and the corresponding impacts were simulated by SPH. Examining the orbital evolution of the debris ejected after collisions, we measured the amount of water delivered to Earth, which is broadly 0.001 ocean equivalents of water, except in one case where one large body transported 7\% oceans of water to the planet. Based on this, and taking into account the frequency of collisions, the amount of delivered water varies between 1.2 and 8.3 ocean's worth of water, depending on the primordial disk mass. According to our results, the prevailing external pollution model effectively accounts for the assumed water content on Earth, whether it's estimated at 1 or 10 ocean's worth of water.


INTRODUCTION
The inner Solar System planets are generally dry (Abe et al. 2000;Raymond & Morbidelli 2022).However, there is some evidence of water on certain inner rocky planets.Earth is the only planet in the inner Solar System known to have liquid water on its surface.About 70% of Earth's surface is covered in oceans.All of Earth's oceans contain roughly 2.5 • 10 −4 M ⊕ water, where M ⊕ denotes the mass of the Earth.This amount is often referred to as one 'ocean' of water, hereafter its mass is denoted by M o .However, the amount of water in Earth's interior is much more uncertain, with different studies suggesting vastly different quantities of water trapped in hydrated silicates and other minerals.It is generally assumed that the amount of water in the interior and on the surface of Earth is between 1 and 10 'oceans' (e.g.Abe et al. 2000;Hirschmann 2006;Marty 2012).
Overall, studies suggest that the origin of Earth's water is a complex process that involves the accretion of water from multiple sources over a long period of time.Below, we list these potential sources of water on rocky planets.Cosmochemical tracers such as isotopic ratios, mainly the D/H ratio, may constrain the origin of Earth's water (for an overview of the D/H measurements, see e.g.Alexander et al. 2012;Morbidelli et al. 2000;Marty 2006).There are six potential models to explain how water was delivered to Earth.They can ★ E-mail: suli.aron@ttk.elte.hu(ÁS) broadly be divided into two categories based on whether the sources of water are local or external.The local source models include (a) the adsorption of water vapour onto silicate grains and (b) the oxidation of H-rich primordial envelope.The external source models include (i) the "pebble" snow model, (ii) the wide feeding zone, (iii) the external pollution, and (iv) the inward migration model.For a detailed description of these models see e.g.Raymond & Morbidelli (2022).
According to recent research by Raymond & Morbidelli (2022); Morbidelli et al. (2000), the prevailing theory suggests that the main source of water on Earth is external pollution from the outer regions of the asteroid belt, rather than local sources.It is believed that this water originates from a few planetary embryos located in the outer part of the asteroid belt, as proposed by Morbidelli et al. (2000).According to the external pollution model, the presence of giant planets strongly influences the orbits of the remaining planetesimals, especially during the phase of rapid gas accretion.In the course of this period, the mass of the giant planets can increase rapidly, from ∼ 10 − 20 M ⊕ to hundreds of Earth masses on a few 10 5 years timescale.The rapid mass increase excites the orbits of the nearby planetesimals that managed to avoid accretion by the growing giants.A multitude of planetesimals experience close encounters with the growing planet, resulting in scattering in all directions.Bodies that are on elliptical orbits that cross the orbit of the giant planet at apostar distance may experience significant energy loss due to gas drag.As a result, the planetesimal's farthest point moves inward and the body will not suffer any more close encounters with the proto-Jupiter.
The most significant secular resonance in our solar system is the  6 resonance (Ito & Malhotra 2006;Minton & Malhotra 2011), which has been in the focus of research of resonant dynamics for several decades and thus has an abundant literature.The  6 secular resonance involves the apsidal precession of asteroids denoted by  and Saturn denoted by  6 .The resonance condition is defined by | −  6 | < , where  6 ≈ 28.2455 arcsec/year (derived from the LONGSTOP 1B numerical integration (Nobili et al. 1989)), and the half-width of the resonant strip, , is usually taken as 1 arcsec/yr.Froeschle & Scholl (1986) conducted numerical integrations over 1 million years to study the orbital evolution in the  6 at 2.05 au.The results revealed that significant variations in eccentricity led to the presence of direct Earth-crossers originating from this secular resonance.The authors suggest that secular resonances should be considered as potential sources of meteorites.Later on numerical simulations by Scholl & Froeschle (1991) demonstrated that the  6 secular resonance often in synergy with the 4/1 mean motion resonance with Jupiter, creates a wide region where Earth-crossing chaotic orbits are produced.The 4/1 resonance rapidly increases eccentricity to 0.5 earth-crosser, whereas remaining in the  6 resonance requires at least several hundred thousand years to become an earth-crosser.The analytical model of Morbidelli & Henrard (1991) is suitable for a global description of dynamics in secular resonances of order 1.The model is applied to investigate the secular resonances  6 ,  5 , and  16 , and illustrations of the secular motion are provided.The  6 resonance sets the inner boundary of the asteroid belt at approximately 2 au.A significant number of near-Earth asteroids originate from the inner main belt rather than the middle or outer belt (Bottke et al. 2002).Many of these asteroids within the inner main belt are affected by the  6 resonance and eventually acquire Earth-crossing orbits (Morbidelli et al. 1994;Bottke et al. 2000).
In a recent paper by Smallwood et al. (2018) the authors conducted -body simulations of a planetary system with an asteroid belt to study how the architecture of the system affects the asteroid impact rate on Earth.The results revealed the significant role of the  6 resonance in the frequency of asteroid collisions with Earth.In further investigations, they replaced Earth with a 10 M ⊕ super-Earth and measured how the asteroid collision rate with Earth changes.If the super-Earth orbited beyond 0.7 au, a significant increase in the number of collisions was observed.However, when it orbited between 1 and 1.5 au, the collision rate decreased.Additionally, altering Saturn's semi-major axis substantially reduced the asteroid collision rate, while increasing its mass raised it.
The formation of the Hungaria family provided further motivation for selecting the location of the collision.The inner edge of the main asteroid belt close to  6 is home to the Hungaria asteroids, including the Hungaria family, which originated from the catastrophic collision of the (434) Hungaria asteroid (Warner et al. 2009;Forgács-Dajka et al. 2022).This is a direct evidence that collisions occurred in the past within this region.
In a recent paper by Martin & Livio (2021) (hereafter ML) the researchers compared the impact efficiencies of different locations in the asteroid belt, namely the  6 secular resonance at  ≈ 2.05 au, at the 2:1 mean-motion resonance with Jupiter at  ≈ 3.32 au and the outer asteroid belt at  ≈ 4.05 au, where  denotes the semi-major axis.They estimated the potential amount of water delivered to Earth and found that the majority of asteroids that collided with Earth originate from the  6 resonance.About 2% of asteroids from  6 collide with Earth, while the collision probabilities from other resonances are orders of magnitude smaller.If the majority of asteroids in the primordial asteroid belt were displaced to the  6 secular resonance, e.g. by chaotic diffusion (Kővári et al. 2023), gravitational scattering, gas drag or the Yarkovsky effect, it is possible that up to ∼ 8 oceans' worth of water could have been transported to Earth.Thus the delivery of at least 1 M o of water from the asteroid belt is possible, but to explain the delivery of 10 M o is more challenging.
In this study, we perform an analysis similar to that of ML but specifically focus on tracking the evolution of debris ejected during collisions of large protoplanets.Following the oligarchic growth model proposed by Kokubo & Ida (2000), we assume that the terrestrial planet formation has been almost finished.At this stage, protoplanets in the inner disk may have grown to Ceres size and potentially larger beyond the snowline.The impact of two Ceres-sized bodies occurs near the  6 secular resonance, where the ejected debris instantly acquires significant eccentricity and inclination.The resonance further increases the eccentricity and inclination of a portion of the debris.As a result, the majority of the debris is rapidly transported to chaotic orbits and can reach the inner planets.The water content of the debris is determined by the original composition of the colliding bodies, and the smoothed particle hydrodynamics (SPH) code consistently accounts for the composition of individual fragments.Thus, in this aspect, we take a step closer to simulating reality.In this study, the simulations of collisions were conducted using the SPH code, which is described in detail in the paper of Schäfer et al. (2020).
In this research, we take a step forward in the simulation of the origin of Earth water.We focus on estimating the amount of water that could have been transported to Earth via collision fragments which were created by a giant impact between two protoplanets.To this end, the dynamical evolution of all fragments will be followed until they impact a large body or the simulation time ends.This scenario is an extension of the external pollution model, since the sudden eccentricity and inclination growth is principally due to collision and not to gravitational interactions.In Section 2, we present the definition of the collision geometry, a short description of the applied SPH and -body simulations and a detailed description of how the initial conditions were constructed from the parameters defining the impact.Section 3 presents the results of the simulations of both the SPH and the -body codes, and we present estimations of how much water could have been transported to the Earth via fragments resulted form a giant collision.We discuss the consequences of water delivery to Earth and draw our conclusions in Section 4.

COLLISION GEOMETRY, SPH AND 𝑁-BODY SIMULATIONS
In this study the Runge-Kutta-Nystrom 7(6) (RKN) Dormand & Prince (1978) algorithm was utilized with adaptive step size control.The method has an acceptable speed and accuracy in most situations.
The integrator was used to model the evolution of a population of asteroids and test particles under the gravitational attraction of Venus, Earth, Mars, Jupiter and Saturn.Uranus and Neptune are not considered in the simulation due to their negligible impact on the dynamics of Jupiter, Saturn and the inner Solar System objects.In the -body code based on the RKN algorithm, all massive bodies interact with all others.Test particles do not influence the motion of the other bodies.In the model we have assumed that all planets have reached their current mass and move on their current orbits.In the simulations, the gas disc had dissipated by this stage and thus gas friction in our simulations is ignored.In the -body simulations, all collisions between bodies are treated as follows: when a close approach with a separation less than the sum of the two radii is detected, the two colliding protoplanets merge to form a new body with a mass equal to the sum of their masses.The orbital elements of the newly formed body are determined based on the position and velocity of the center of mass of the two colliding bodies at the moment of impact.In each collision the composition of the resulting body is determined based on the material composition of the colliding bodies.Most of the previous work on planet formation used this oversimplified treatment of collisions.In the -body runs in order to determine the effect of collisions on the composition of massive bodies, we track the collision history of each individual fragment with mass.

Initial condition
The process to create collision fragments that originate from impacting bodies with given orbital elements and collision parameters is complex.A description of this procedure is given below and can be divided into the following four main steps: (i) Calculate the radii of the target  t and the projectile  p depending on the composition of the objects.
(ii) Generate the position and velocity vectors for the planets as well as for the projectile and target.Using the specified orbital elements of the target and the collision-defining parameters one has to calculate the initial conditions for projectile.With these initial conditions, the system is integrated backwards in time for approximately  back ≈ 10 4 seconds.
(iii) Compute the collision with the SPH code.From the system obtained at the time  =  back , the initial conditions for the SPH code have to be calculated.With these initial conditions defining the barycentric position and velocity of the Sun, the 5 planets, and the target and projectile, the SPH code is used to simulate the collision for a duration of approximately  SPH ≈ 8 • 10 4 seconds.The end state of the SPH run defines the initial state for the next step.
(iv) Calculate the initial conditions for the -body code from the SPH run's end state.In this step, it is necessary to identify the cohesive debris fragments, which are bodies composed of one or more SPH particles.Once this is accomplished, their position and velocity vectors must be computed which in conjunction with the data of the major planets constitute the initial conditions for the Nbody simulation.The -body simulator is then integrates the system for a duration of  = 10 6 years.
Typically, it takes about 1 million years for these fragments to be transported into Earth-crossing orbits when their eccentricity is directly increased to around ≈ 0.6 by the resonance (Morbidelli et al. 1994;Smallwood et al. 2018).Below, we provide a detailed description of the above four steps.
To complete the first step, we based our assumptions on the observations of Thomas et al. (2005).The authors have analyzed the Ceres' spectral features in reflected sunlight and demonstrated that the presence of a rocky core and a water-ice mantle aligns with the data assuming the densest materials for the core with a mean density of 2.077 kg m −3 .They estimated that the ice mantles have a thickness ranging from 110 to 124 km, accounting for 24-26% of the total mass of the object.In accordance with these results, in the present paper, the composition of the impacting bodies consists of a basaltic core and a mantle composed of water ice.This structural composition of bodies were utilized by Maindl et al. (2014) when studying the fragmentation of icy objects.
Therefore, our scenarios include two Ceres-mass  Ceres = 9.38 • 10 20 kg protoplanets, where both bodies consist of a basalt core (75% of mass) and a mantle/shell of water ice (25% of mass).The  two protoplanets together contain approximately 0.31 M o of water.In Fig. 1, we depict the moment of collision between two bodies in the BC barycentric coordinate system.The target is the larger of the two colliding bodies, and if the masses are equal, then the object with the smaller identifier number becomes the target.The barycentric position vectors of the target and projectile are denoted as r  and r  , respectively.The relative position vector of the projectile with respect to the center of the target is d, given by d = r  − r  , and the relative velocity vector v is defined as v = v  − v  , where v  and v  represent the barycentric velocity vectors of the target and projectile, respectively.The  vector corresponds to the angular velocity vector of the target.
In Fig. 2 the two impacting bodies are shown in the plane of the collision defined by the d and v vectors.The basalt core is represented by a light grey coloured circle, while the ice mantle is denoted by a sky blue coloured one.The spatial arrangement of the SPH particles composing the target and projectile is governed by a hexagonal lattice structure.This lattice structure was determined based on the given core and shell masses.The radii of the target and the projectile are equal and denoted by  t =  p = 4.974 • 10 5 m, resulting in a bulk density of  = 1.82 g cm −3 .The radii of the target's and projectile's core are equal too and denoted by  t,c =  p,c = 3.956 • 10 5 m.
To execute the second step, multiple consecutive rotational transformations need to be performed.It is convenient to perform these rotational transformations in the coordinate system fixed to the target's center denoted as  t .We denote this coordinate system as  t , where the  t  axis corresponds to the rotation axis  of the target, and the plane perpendicular to the  t  axis and passing through the centre is the  −  plane, see Fig. 3.
The geometry of the collision is completely determined by the impact position vector d and the velocity vector v relative to the origin  t .In previous studies, the d vector is usually not provided, only the collision velocity  = |v| and the impact angle  are given, where  is defined as the smaller angle between the vectors d and v at the moment of first contact, see Fig. 2. Therefore, only a "collision cone" is defined, with the apex at the centre of the target  t and the lateral surface composed of possible collision velocity vectors.The impact velocity is given in mutual escape velocity units which is defined as follows (Genda et al. 2012): where  t and  p denote the mass of the target and projectile, respectively and  G is the Gaussian constant of gravity.
When  = 0 • , the collision is referred to as head-on, while  = 90 • corresponds to a grazing collision.In terms of geometric shapes,  represents the half-angle of the cone.For the sake of generality, in the following discussion, we also take into account the d vector.Given r t and v t of the target, as well as the  t  axis, , and , we aim to determine the barycentric position vector r p and velocity vector v p of the projectile, see Fig. 1.Let us express the components of the d and v vectors in the  t  "targetgraphic" coordinate system where  1 is the longitude and  1 is the latitude of vector d, and similarly  2 is the longitude and  2 is the latitude of vector v (since v is a free vector, thus we have imaginarily shifted it to the origin  t ).Let f (1) and f (2) denote the unit vectors in the directions of d and v, respectively Rotate the  t  coordinate system around the  t  axis with  1 angle and then the resulting the  t  ′  ′  ′ coordinate system around the  t  ′ axis with  1 .Let us denote the rotation matrix around axis  t  by R  ( 1 ), and around the  t  ′ axis by R  ′ ( 1 ).They are defined as (5) The rotated system is shown in  possible collision vector, which we denoted as ṽ in Fig. 4.An -angle rotation about the  t  ′′ axis connects the original v and ṽ.Then if we denote by f (2) the unit vector in the direction of ṽ: and then rotate it around the axis  t  ′′ by  it yields Let us transform f (2) back to the original  t  coordinate system.Applying R  ′ (− 1 ) and then R  (− 1 ), we get where the components of f (2) are given in Eq. ( 7).
As a last step, one has to transform the vectors into the inertial barycentric coordinate system BC as shown in Fig. 1.The connection between BC and  t  systems is defined by the  angular velocity vector of the target.The components of  in the BC coordinate system may be given as where  is the longitude and  is the latitude in the "barygraphic" coordinate system.The transformation is defined by two subsequent rotations by applying R x (−) and then R y ′ (−).The computation of these matrices is similarly done as before.
After performing the aforementioned transformations, we obtain the barycentric coordinates and velocities of the projectile.From these, one can calculate the orbital elements associated with the given collision.These computed orbital elements together with the those of the other bodies are listed in Table 1 for all scenarios and serve as the initial conditions for the backward integration.For all the scenarios the initial orbital elements of the bodies are the same except for the projectile, which are given in the table for all the six scenarios.
According to the last objective in step 2, the backward integration is performed for  back ≈ 10 4 seconds.As a result, the separation between the colliding bodies will be at least 5 times the sum of their radii.The SPH code requires this initial separation distance in order for the SPH particles to arrange themselves properly due to mutual tidal and gravitational forces.The end state of this backward simulation is used to generate the initial condition for the SPH code.In this article, for simplicity, we assume that initially neither the target nor the projectile rotates around its axis, and that the collision occurs in a plane parallel to the BC plane, thus  1 = 0 and  1 = 0 therefore d = (, 0, 0).
In step 3, the SPH code developed by Schäfer et al. (2020) was utilized.This code has proven to be highly effective in simulating highresolution surface impacts (Maindl & Haghighipour 2014;Maindl et al. 2015) as well as investigating collisions between planetesimals and protoplanets (Maindl et al. 2014;Maindl & Dvorak 2014).It has been successfully utilized in both scenarios, yielding valuable insights and accurate results.In SPH-based simulations, continuous bodies are represented by discrete SPH particles (e.g.Gingold & Monaghan 1977;Benz & Asphaug 1995;Maindl et al. 2013;Schäfer 2005).These particles carry the specific physical properties of interest and contribute to the overall physical properties of their respective volume elements.In this study, we employ a relatively low resolution of  SPH ≈ 20 • 10 3 SPH particles, which means that each SPH particle has a mass of  SPH =  t / SPH = 4.695 • 10 16 kg.
Due to a collision, a large amount of debris is generated and scattered.If we were to consider every debris as a massive body in the -body code, the simulation would take a very long time to run.In order to strike a balance between the statistically significant number of bodies and the runtime of the simulation, we only considered debris with masses denoted as  frag above a certain threshold as massive bodies and treated the rest as test particles.A fragment is considered as protoplanet if  frag ∈ [10 −10 , 10 −7 ] M ⊕ , if  frag ∈ [10 −13 , 10 −10 ) M ⊕ then it was treated as planetesimal, and any remaining mass below 10 −13 M ⊕ = 1.87 • 10 17 kg as test particle.This lower mass limit corresponds to approximately 4 SPH particles.This compromise ensures to get statistically meaningful results and to keep the number of massive bodies manageable while guaranteeing that the debris is not only present in the form of test particles in the -body simulator.
In Fig. 5, a series of four snapshots are presented for the collision with  = 30 • and  = 3.0 [ esc ].The spherical objects are composed of basaltic rock represented by grey and water ice represented by sky blue.Neither of the bodies is rotating at the start of the simulation.The collision happens at  ≈ 10 4 seconds, and the evolution of the colliding bodies is followed via the SPH code for  = 8 • 10 4 seconds.The lower right frame in Fig. 5 illustrates the final outcome of the SPH simulation, which will serve as the input for the subsequent body simulation.The region occupied by the fragments in the final snapshot is  ∈ [1.981, 1.982] au,  ∈ [8.638, 10.61] × 10 −3 and  ∈ [−4.737, 182.6] × 10 −5 au, which fits in a rectangular cuboid with sizes of Δ = 10 −3 , Δ = 2 • 10 −3 , and Δ = 2 • 10 −3 au.
From the barycentric position and velocity vectors of the fragments at  = 8 • 10 4 seconds, the orbital elements were computed and the eccentricities  and inclinations  versus the semi-major axis  are plotted in Figs. 6 and 7, respectively for all six runs.The large black plus sign indicates the location of the collision at  = 2.08 au,  = 0.05, and  = 5 • (the values of the other orbital elements were set to zero, c.f. the data for the Target in Table 1).The blue-filled circles denote the protoplanets, the smaller red circles are the planetesimals, while the small grey dots are the test particles.
In the ( − ) planes, an interesting, checkmark-shaped structure forms in all cases for  = 0.Both the massive and the massless bodies are confined to this shape.On the contrary, no well-defined shape appears on the ( −) plane, where the fragments are scattered around the impact point.However, as velocity increases, distinct arcs become visible, particularly noticeable in the planetesimal distribution.The mechanism and theory behind the formation of the checkmark-shaped structure will be addressed in another article, it is out of the scope of the present paper.We also note that the semi-major axes of the smaller fragments span a wide range, e.g. for run1 it is approximately  ∈ [1.75, 2.5] au, which is a difference of Δ ≈ 0.75 au, and about 1000 times larger than in the BC space.The higher the collision velocity, the wider the range over which the debris scatters.This fact plays a crucial role in the frequency of collisions with Earth, as discussed in Section 3 and demonstrated in Table 3.For all six runs the initial  falls in the interval of [0, 0.25] while the  is in the range of [0, 20] degree, c.f. Figs. 6 and 7.With an increase in , the upper limits of the initial  and  grows.
The simulation parameters and initial orbital elements for the backward integrations are provided in Table 1.The planets were initially placed in their current orbits.Now we turn our focus to the determination of the fraction of collisions with an inflated Earth.The collision parameters ,  were selected based on the results of Süli (2021), where the author performed 2D -body simulations of planetary accretion.From the compiled statistics, the impact parameter  = sin  has a uniform distribution and the minimum of the impact speed is 1.Based on preliminary results from 3D planetary accretion simulations, the collision angle has a maximum of around 40 degrees, and the collision velocity reaches its maximum between 3 and 10.In light of these results, we selected the collision angles of 30 and 45 degrees, as well as the collision velocities of 3, 5, and 7, see Table 1 for details.
At the beginning, both the target and the projectile are initially in the  6 resonance, raising the question of why they have not been ejected from there, given the relatively short timescale required for a significant increase in eccentricity (Morbidelli et al. 1994;Smallwood et al. 2018).On the other hand, Yoshikawa (1987) and Knezevic et al. (1991), have shown that the  6 resonance is an efficient collector Sim.
capable of capturing objects from the region between 2.4 and 2.7 au with moderate inclination (15 • to 20 • ).Therefore, we may assume that the target and the projectile were formed further away, and when they reached the size of Ceres, they entered the  6 resonance and eventually collided.

RESULTS
We numerically integrate the trajectory and simulate the outcomes of fragments.Each simulation incorporates the presence of the three terrestrial planets Venus, Earth and Mars and two giants Jupiter and Saturn.The initial distribution of the eccentricities and inclinations are shown in the first and third rows of Figs. 6 and 7, respectively.The large black plus marks the impact location, blue circles indicate protoplanets, red ones represent planetesimals, and grey dots correspond to test particles.In the ( − ) plane, a distinct checkmark shape is clearly visible in the initial distribution of the generated debris.Additionally, it can be observed that collisions at higher velocities scatter the debris over a wider range along the  axis, widening the arms of the checkmark.Increasing the collision angle from 30 to 45 degrees roughly halves the number of fragments.It's also worth noting that the distribution of planetesimals concentrates on the right arm of the checkmark, and as the velocity () increases, it clusters towards the lower and upper boundaries of the checkmark shape.In general, after each collision, two protoplanets remain, with the exception of run5, where collision parameters are such that both the target and the projectile completely disintegrate into smaller fragments, leaving no protoplanets.Fig. 7 similarly illustrates the initial and final states of collision remnants but in the ( − ) plane.With increasing velocity, arcs become apparent in the distribution of planetesimals.
At the end of the simulation ( = 10 6 years), we also depicted the distribution of debris in both the ( − ) plane, (Fig. 6) and the ( − ) plane (Fig. 7).For better comparison, the displayed range is the same as the one used in the initial state.Consequently, some of the debris is not visible in the figures, even if it did not collide.The most prominent feature in the final state depicted in the figures is the depletion of the region around the  6 secular resonance, between ≈ 2.05 and ≈ 2.15 au.This phenomenon is most pronounced in the ( − ) plane, as it is more sensitive to secular perturbations (see, for example, Figure 11 in Forgács-Dajka et al. ( 2022)) than in the ( − ) plane.Generally, dynamical maps are usually constructed using the maximum eccentricity method.In the ( − ) plane, the  6 secular frequency is also beginning to manifest itself and is expected to become more pronounced with longer integration times.
In Figs. 8 and 9 you can see the initial positions of celestial bodies that collided with the Sun or one of the rocky planets.These figures clearly illustrate the extent and shape of the  6 secular resonance.The  6 resonance becomes more prominent as the initial fragment distribution covers a larger area in the ( − ) plane.This inspires us to conduct further research in which the entire area will be mapped with test particles to precisely determine the shape and extent of the  6 resonance in both the ( − ) and ( − ) planes.In Fig. 9 the colored region delineated by colliding bodies extends further to larger -values as inclination increases.This observation aligns with Figure 1a in the study by Froeschle & Scholl (1986), although in our figure, the vertical axis represents the osculating inclination while in Froeschle & Scholl (1986)'s figure, it represents the proper inclination.
According to the results presented in ML, there were no collisions observed when considering the actual radius of the Earth ( ⊕ = 6371 km).Therefore, based on their results, we used a 10 times larger radius for the Earth  Earth = 10 ⊕ in the simulations to ensure that collisions occur.The radii of the other planets were taken as their actual values.We considered a particle to hit the Sun if its astrocentric distance became smaller than 0.3 au.
We will use the estimation employed by ML to determine the number of objects colliding with a real-sized Earth based on the number of objects colliding with an inflated Earth of size 10 ⊕ .The estimation resulted in a factor  ≈ 20 for the position of the  6 secular resonance.
According to Table 2 it can be observed that at  = 30 • , approximately twice as many fragments are produced compared to the case of  = 45 • collision for all three impact velocities.A significant portion of the debris collides with larger fragments, primarily with the remnants of the collided bodies, which are the two largest objects Basalt material is represented in grey, while water ice is depicted in sky blue.Both bodies have a layer of water ice (25 percent of its total mass).We note that the scale increases on successive snapshots.In the top left panel, it is a few times ∼ 10 3 km, while in the bottom right panel, it is a few times ∼ 10 5 km kilometers.
in the debris cloud.The sum of these collisions can be found in the table under the  TP header (the subscript TP stands for target and projectile).This process, known as reaccretion, returns a substantial fraction of the mass back to the target and projectile, the specific mass can be found in Table 5.The reaccretion occurs within a few years following the collision, and after that, debris rarely impacts on the remnants of the impacted bodies.However, in the following years collisions still occur between smaller fragments and test particles, but the rate of these collisions decline with time.
Overall, based on the table, it can be concluded that a smaller impact angle, representing a more frontal collision, generates more debris than a larger impact angle, which corresponds to a more tangential collision.Of course, for both impact angles, higher velocities further increase the amount of debris, and they scatter more widely in both the ( − ) and ( −) planes.The wider the range in which the fragments are scattered along the semi-major axis, the proportionally fewer of them end up in the  6 resonance.Therefore, in the case of higher-velocity collisions, fewer fragments enter the  6 resonance region.
The outcomes for  = 3,  = 30 • and  = 45 • are shown in the left panels of Fig. 11 and the data are listed in the first two rows of Table 2.It is clear that the majority of collisions occur with the Sun and with Earth.The number of collisions with the Sun significantly exceeds the number of impacts with Earth.Within the simulation time span of 1 Myr, there were a total of 47 impacts on the Earth, for  = 30 • , while 33 impacts on Earth for  = 45 • .In the case of  = 5 there were 95 and 52 Earth impacts for 30 • and 45 • impact angles, respectively.For the case of  = 7, the measured values are 116 and 58, respectively.For the Sun, these numbers are as follows: in the case of run1 and 2, 178 and 181; in the case of run3 and 4, 427 and 307; while in the last two runs, 600 and 245.The observed trend is that with more oblique impacts, a larger number of fragments are  2).The initial eccentricities were computed using the barycentric position and velocity vectors taken from the last frame of the SPH simulation (see Fig. 5).The big black plus denotes the position of the impact.Blue circles represent protoplanets, red ones planetesimals while grey dots correspond to test particles.
produced and as a result of this the number of impacts also increases.The number of impacts into the giant planets are in close agreement with the results of ML and are in the order of unity.
We have calculated two ratios using two different benchmarks to measure the chance of collision with Earth.In the first case  10 =  E /( pp +  pl +  tp ), while in the second one  ★ 10 =  E / event In order to directly compare the results with those of ML, the latter definition matches what was used in the ML paper.Both  10 and  ★ 10 characterize the likelihood of objects impacting the inflated Earth, thus these ratios can be considered as a loosely defined probability of impact with Earth.The calculated values of  10 and  ★ 10 are listed in Table 3.The measured  10 decreases with increasing , the maximum is ≈ 0.0258 ( = 3,  = 45 • ), the minimum value is ≈ 0.0111 ( = 7,  = 30 • ), thus both the minimum and the maximum fall into the same order of magnitude.The decreasing trend is a consequence of the fact that at higher velocities, the fragments making up the debris scatter over a larger range along the semi-major axis (as mentioned above).Consequently, proportionally fewer fragments enter the  6 resonance, resulting in fewer fragments on Earth-crossing orbits.
The ratio  ★ 10 does not follow a similar trend, it fluctuates around a mean value of 0.0876 with a standard deviation of 0.0112.To adjust for the inflated Earth, these values need to be divided by a factor of  = 20, as described above.The computed ratios  and  ★ are given in the last two columns of Table 3.Since  and  ★ differs by a factor of  they follow the same pattern as  10 and  ★ 10 .In general, the histograms in Fig. 11 are similar.In the time follow- ing the first collision, approximately equal amounts of debris collide with both the Earth and the Sun.However, over a period of a few hundred thousand years, this ratio significantly shifts in favour of collisions with the Sun.The earliest event occurs between 10 5 and 2 • 10 5 years, indicating that it takes roughly this amount of time for the first debris fragments from the  6 resonance to reach the inner Solar System.Interestingly, in the run1 simulation, the first collision occurs with Venus at  ≈ 1.107 • 10 5 years, followed by a longer period without any event until around  ≈ 4 • 10 5 years.From this point on, debris fragments collide with both the Earth and the Sun.
We note that a few impacts on the other rocky planets were also detected, 1 impact into Venus (green bar) for run1 and 2 collisions with Mars both in run5 and run6 (red bar).Mars collision occurs only for the  = 7 cases.
In Fig. 10 a part of the collision history of an individual protoplanet is visualised using a collision tree.The collision events are plotted only for  ≤ 30 years.The gray-black lines show the mass of fragments as a function of time, while the blue line shows the increase in mass of the target body left after the SPH simulation.In the figure, for clarity, only one collision tree has been drawn.The smaller fragments merge together during the first 10 years, forming a larger body with mass ≈ 3 • 10 19 kg, which hits the target at  ≈ 23 years.Notice that initially several fragments have the same mass, thus gray-black lines are overlap in the beginning of the simulation.

Water delivery by fragments
As previously stated, the inner region of the original planetesimal disk is dry, therefore several explanations for the presence of water on Earth were devised.According to a group of models explaining the origin of Earth's water, the source is located in the outer regions of the Solar System.Among these models, the external pollution  2).The upper panels depict simulations with  = 30 • while the lower panels for  = 45 • .The colors indicate the fragments colliding with the different celestial bodies, while the initial position of all bodies involved in the  -body simulation is indicated in grey in the background.model is currently widely accepted, which successfully accounts for the quantity and chemical composition of water on Earth and is consistent with various theories of Solar System formation.The present study seamlessly fits into this and the other models too, assuming that Earth's water originate from the outer Solar System, as collisions between protoplanets are inevitable and frequent during the early stages of planetary evolution.
The composition of the bodies in the asteroid belt depends on the radial distance, as indicated by studies e.g.Kerridge (1985); Gradie & Tedesco (1982); DeMeo & Carry (2014).The snowline  SL , which is approximately 2.7 au from the Sun is a boundary beyond which volatile substances, such as water can condense and freeze.The inner regions,  <  SL au are predominantly populated by S-type asteroids, their fraction of water is typically quite low, with estimates ranging from less than 0.1% to 1% by mass.In contrast, beyond  SL the asteroid belt are dominated by C-type asteroids, which generally have a water fraction  water,C of approximately 10% by mass (Kerridge 1985).The current mass of water in Earth's oceans is 2.5 • 10 −4 M ⊕ , while the mass of the asteroid belt is only 4.5 • 10 −4 M ⊕ despite covering a vast region from 2 to 4 au (Krasinsky et al. 2002;Kuchynka & Folkner 2013).This mass is orders of magnitude smaller than what would be predicted by models of planet-forming disks, such as the Table 2.The simulation outcomes with  Earth = 10 ⊕ for  = 10 6 years.The first column contains the names of the simulations.Columns 2 and 3 show the impact velocity  and angle , respectively.Column 4 represents the total number of bodies (including the Sun and 5 planets)  .Columns 5,6 and 7 display the number of protoplanets  pp , planetesimals  pl and test particles  tp respectively.Column 8 -13 contains the number of events (outcomes)  event , fragments that collided with the Sun  Sun , Venus  V , Earth  E , Mars  M and Jupiter or Saturn  JS , respectively.Column 13 lists the number of collisions  TP with the remnants of the target and projectile.

Sim.
After a collision is simulated with SPH, the orbital evolution of the fragments is tracked using our  -body integrator and if a collision is detected then the two bodies merge.This figure shows only a part of the whole collision tree, where the mass evolution of some fragments that eventually impact the target are depicted by gray-black lines.Every jump in the lines correspond to a collision.(Weidenschilling 1977;Hayashi 1981).We note that actually the MMSN is not a nebula, but represents a protoplanetary disc with solar composition containing the requisite amount of metals for the formation of the eight planets within the solar system, including the asteroid belts.In order to determine the original total mass of solids in the asteroid belt the MMSN were used in which the surface density of solids is where  n is a nebular mass scaling factor,  ice is the ice condensation factor, Σ 1 is the surface density of solids at 1 au and  is the distance from the Sun.According to Hayashi (1981) Σ 1 = 7.1 gcm −2 ,  ice = 1 for  ≤  SL and  ice = 4.22 for  >  SL .By integrating Eq. ( 10) from 1.55 (3 times the Mars' Hill radius) to 4.1 au (3 Hill radii from Jupiter) we obtain the mass of solid material in M ⊕ units.Despite its widespread use, the MMSN model has several limitations.It assumes planets gathered all nearby solids, but only some of the solid material were in planetesimals, while others were in the form of dust which was later lost as the solar nebula's dissipated.Adjusting for this lost mass relies on the fraction of solids in planetesimals.However, the model assumes that all planetesimals were accreted, which may not hold true for terrestrial planets.In the region between 0.4 and 4 au, where terrestrial planets were forged, simulations show significant, ≥ 50% mass loss and redistribution during the final assembly stage (Raymond et al. 2005).The asteroid belt and the region near Mercury have also been depleted of solids at some point (Weidenschilling 1977).This suggests that the MMSN model doesn't fully capture the complexities of planet formation in these areas.Thus, for the MMSN  n = 1, so Eq. ( 11) results in 6.69 M ⊕ which is the minimum mass required to forge the terrestrial planets.In light of the above, it is evident that  n > 1, and simulations examining the formation of the solar system suggest that 2 ≤  n ≤ 3 (see e.g.Fogg & Nelson 2005, 2007;Raymond et al. 2005;Raymond & Izidoro 2017).Since S-type asteroids are very dry, we are interested only in C-type asteroids which are primarily located beyond the snowline.The mass  C contained in  SL <  < 4.1 au is  C = 5.38  n M ⊕ and it is composed predominantly of C-type asteroids.
If we focus only on C-type asteroids and assume that they entered the  6 resonance via some mechanisms (e.g.chaotic diffusion, gravitational scattering, etc.), the estimation for the maximum amount of water that could be delivered to Earth can be given by the following formula (see Eq. (2) of ML) The computed results for  wd are listed in  values were determined without considering any "delivery" and/or "impact" losses.It was assumed that the the total water inventory of the debris was deposited on Earth.Clearly, this provides a generous upper bound estimate on Earth's water inventory.These results are in good agreement with estimates for the total water inventory on Earth may range up to 10  o mantel water (see e.g.Abe et al. 2000;Hirschmann 2006;Marty 2012).The value of  ★ wd listed in Table 4 was calculated similarly to  wd , but using  ★ .The values of  ★ wd with  n = 3 range from 23.76 to 32.67, which are significantly larger than  wd .These results greatly exceed the estimate for the Earth's total water supply.
In Figure 12 the chronological record of collisions between fragments and Earth is shown.On the left hand the vertical axes shows the cumulative count of these collisions (blue line), while on the right hand axis the cumulative mass of the fragments are plotted.From the figure it is evident that initially smaller debris composed solely of pure water ice collide with the Earth.This typically begins around ∼ 2 • 10 5 − 4 • 10 5 years after the impact event.Larger debris carrying both basalt and water ice, on the other hand, collide with the Earth later, around ∼ 4 • 10 5 − 6 • 10 5 years after the impact.
According to Table 4 in most cases, the quantity of water delivered to Earth by the debris fluctuates around 0.001 ocean units, except for the run3 where more than 0.07 'ocean' of water is delivered.This means that as a result of a single distant collision, over 7% of the surface water is delivered to Earth.In this particular case, the remnants of the projectile collided with Earth, bringing along a significant amount of water (initially the projectile contained ∼ 0.16 'ocean' of water).According to Morbidelli et al. (2000), the majority of the water currently present on Earth was delivered by a small number of planetary embryos that were initially formed in the outer asteroid belt and subsequently accreted by the Earth during its final stage of formation.The results of run3 are consistent with this finding.The special case of run3 is very well visible in Fig. 13 (green solid line), where the quantity of water transported to Earth, measured in Table 4. Maximum amount of water that may deliver to Earth from the asteroid belt.Columns 2 and 3 list  wd as computed from Eq. ( 12) using  from Table 3 for  n = 1 and 3, respectively.Columns 4 and 5 are the same as Columns 2 and 3 but for  ★ .Column 6 presents Δ ice as determined from impacts of fragments into the Earth using the results of the simulations (cf.Table 5).

Sim.
terms of ocean units versus time is plotted for the 6 simulations.In general, the amount of water delivered to Earth initially increases rapidly, then the growth rate slows down and stabilizes at around 0.1%.An exception to this trend is observed in run3, as discussed earlier, when the remnants of the projectile hits the Earth at about  ≈ 8.7 • 10 5 years, and depositing large quantity of water ice.
In Table 5 we have summarized the total mass delivered by impacts in all six simulations, along with its distribution specifically between basalt and water ice.Only the bodies that have participated in at least one collision are listed in the table.The second column provides information about the mass increase resulting from impacts.We can see that different bodies experience varying degrees of mass increase, with values spanning several orders of magnitude.Among the terrestrial planets Earth receives the highest amount of material, while Venus and Mars acquire negligible mass.The higher mass increase on the Earth could be attributed to various factors, predominantly its inflated size and gravitational attraction.
In all six runs the mass increase for the Sun is significantly higher  Generally the mass increase resulting from impacts is roughly the same for the Target and the Projectile in all scenarios, except two cases: In run1, the target accretes five times more mass than the projectile, while in run5, the post-projectile (pl7) accumulates three times more mass than the target.The number of collisions varies between the Target and Projectile.In some scenarios, the Target experiences a higher number of collisions, while in others, the Projectile has more collisions.
Overall, the data in the table provides insights into the effects of impacts on mass increase and collision occurrences for various bodies.

CONCLUSIONS
The focus of this study is to examine the fate of debris ejected during collisions among celestial bodies in the asteroid belt and to study their role in the delivery of water to the Earth.In the applied model all major planets in the Solar System have already formed and are revolving around the Sun on their current orbits.To achieve our objective, we utilized SPH and -body simulations to quantify the water transported to Earth through debris ejected from collisions between two Ceres-sized bodies near the  6 secular resonance.We successfully combined the SPH code with the -body simulator and provided a detailed description of generating initial conditions that lead to collisions at the desired location and velocity.
As a result of -body simulations in the ( − ) planes (see Fig. 6), an intriguing checkmark-shaped pattern emerges at  = 0 in all instances, encompassing both massive and massless bodies.In contract, in the distribution of bodies on the ( − ) plane (see Fig. 7) such a distinct structure is not as apparent.However, at higher collision velocities, two arc-like structures become evident in the distribution of planetesimals.The other type of fragments are dispersed around the impact site.It is clear that in the vicinity of the secular resonance, the ( − ) plane clears out, and the bodies end up on highly eccentric orbits (which is not visible in the figures because  < 0.25).The  6 secular resonance pumps up the eccentricities, and the highly elongated orbits intersect with the orbits of Mars and even Earth, allowing for impacts on these rocky planets.In Figs. 8 and 9, the extent and shape of the  6 secular resonance are clearly visible from the distribution of objects that collided with Earth during the simulations.
After the collision, the created the check-mark-shaped pattern's extension along the semi-major axis, increases with velocity .At  = 3, it spans more than 0.5 au, while at  = 7, it covers nearly 1 au in width.Consequently, if a collision deviates from the location of the  6 secular resonance by no more than 0.5 au, due to the scattering effect, debris will enter the resonance.The quantity of debris that enters depends on the distance and velocity of the collision.Since collisions are frequent in the early stages of planetary formation, it's likely that the  6 resonance receives a continuous supply, which it then "transmits" to the inner regions through the mechanism outlined Table 5. Summary of the collisions with  Earth = 10 ⊕ for  = 10 6 years for all six giant impacts.A specific object is only listed if it has suffered from a collision.The first column contains the names of the bodies.In the second column, the mass increase resulting from impacts is displayed in kg.In the third and fourth columns, this mass increase is specified for different materials, namely basalt and water ice.The fifth column shows the number of collisions (cf. it with the data in Table 2).In run5 both the target and the projectile disintegrate into smaller fragments, leaving no protoplanets but planetesimals pl6 and pl7.

Name
Δ above.From the perspective of the inner planets, this represents a continuous water delivery system that persists as long as there are collisions and debris enters the resonance.
Based on the data from the -body simulations, we calculated the probability of the ejected debris impacting Earth.Using this information, we provided an upper estimate for the amount of water that can be transported by the debris.Considering the significant uncertainties and strong model dependencies, our estimate is in reasonable agreement with previous results, such as ML and Morbidelli et al. (2000).According to our findings, a minimum of ≈ 1.2 and a maximum of ≈ 2.8 ocean equivalents of water can be delivered if the MMSN model is considered.However, if we assume an equivalent mass of 3 × MMSN of solids and calculate the water quantity, the range increases to 3.6 to 8.3 ocean equivalents, indicating the potential delivery of larger amounts of water to Earth.Therefore, if the total water inventory on Earth is estimated to be between 1 and 10 ocean equivalents, the external pollution model can provide an explanation for Earth's total water inventory.If this is true, it would not require assuming local sources to explain the origin of water, although they cannot be ruled out either.
During the -body simulations, we were able to directly measure the amount of water brought to Earth by the ejected debris during individual collisions.Typically, this amounted to 0.001 ocean equivalents of water.An interesting exception was observed when a single large body, the remnants of the projectile, collided with Earth and deposited approximately 0.07 ocean equivalents of water on the Earth's surface.The results of this study support the widely accepted external pollution model.The maximum estimated amount of water that was transported to the Earth via fragments varies between 1.2 and 8.3 ocean's worth, depending on the initial mass of the protoplanetary disk.

Figure 1 .
Figure 1.The collision of two bodies in the barycentric coordinate system.The geometry of the collision is defined by the d relative position and v relative velocity vectors.

Figure 2 .
Figure 2. The plane of  ′′ −  ′′ is defined by the d and v vectors.The grey-coloured circle denotes the basalt core, while the sky blue denotes the ice mantel.
Fig 4, denoted by  t  ′′  ′′  ′′ .The intersection of the collision cone and the  t  ′′  ′′ plane specifies a

Figure 3 .
Figure 3.The  t   targetgraphic coordinate system, where  is the angular velocity of the target,  denotes the longitude and  the latitude, d is the position vector of the projectile, v is the impact velocity vector.The d and v vectors determine the plane of collision,  is the impact angle.

Figure 4 .
Figure 4.The projection of the collision cone onto the plane  t  ′′  ′′ .The ṽ is a possible impact vector with the right magnitude  and collision angle , which is the half-angle of the cone.

Figure 5 .
Figure 5. Collision snapshots for  = 3.0 [ esc ] and  = 30• .The simulation time for each frame is indicated at the top.The actual collision takes place at  ≈ 10 4 seconds.Basalt material is represented in grey, while water ice is depicted in sky blue.Both bodies have a layer of water ice (25 percent of its total mass).We note that the scale increases on successive snapshots.In the top left panel, it is a few times ∼ 10 3 km, while in the bottom right panel, it is a few times ∼ 10 5 km kilometers.

Figure 6 .
Figure 6.The ( − ) plane of the fragments in the  -body simulations.The top two rows show the initial (first row) and final (second row) eccentricities of the fragments with  = 30 • for different impact velocities, while the third and fourth rows for  = 45 • (see Table2).The initial eccentricities were computed using the barycentric position and velocity vectors taken from the last frame of the SPH simulation (see Fig.5).The big black plus denotes the position of the impact.Blue circles represent protoplanets, red ones planetesimals while grey dots correspond to test particles.

Figure 8 .
Figure 8.The positions of the fragments in the ( − ) plane that collided with the Sun or one of the rocky planets for different impact velocities and angles (see Table2).The upper panels depict simulations with  = 30 • while the lower panels for  = 45 • .The colors indicate the fragments colliding with the different celestial bodies, while the initial position of all bodies involved in the  -body simulation is indicated in grey in the background.

Figure 10 .
Figure 10.An example of a collision tree, showing the mass of the fragments (gray-black) and of the remnants of the target (solid blue line) vs. time in run5.After a collision is simulated with SPH, the orbital evolution of the fragments is tracked using our  -body integrator and if a collision is detected then the two bodies merge.This figure shows only a part of the whole collision tree, where the mass evolution of some fragments that eventually impact the target are depicted by gray-black lines.Every jump in the lines correspond to a collision.

Figure 11 .
Figure 11.The event history of the  -body simulations are binned into 50 kyr time intervals for the different impact speed and angle, see Table2.The upper panels depict simulations with  = 30 • while the lower panels for  = 45 • .The Earth radius is 10 ⊕ for all cases, the other planets have their actual radius.In the plots, the yellow bars show the number of fragments that hit the Sun, while the green, blue and red bars represent asteroids colliding with Venus, Earth and Mars, respectively.

Figure 12 .
Figure12.The event history of bodies colliding with the Earth.On the left axis, the cumulative number of bodies that collided with the Earth (blue line).On the right axis, the cumulative mass of bodies that collided with the Earth (brown and blue bars) in units of Earth-mass.The brown colour represents the amount of basalt in the total mass, while the blue colour represents water ice.Note that the horizontal axis begins at  = 10 5 yr.

Figure 13 .
Figure13.The amount of water delivered to Earth in ocean units over time in the 6 simulations.The run3 is a special case where the water is primarily delivered by the remnants of a single object, the remnant of the projectile, which collided with Earth at  ≈ 8.7 • 10 5 years.

Table 1 .
Initial conditions for the  -body simulation scenarios for backward integration.The first column contains the name of the simulation.Columns 2 and 3 show the impact velocity  and angle , respectively.Column 4 displays the name of the body, Column 5 represents its mass in Earth-mass units.Columns 6-11 show the orbital elements of the bodies,  is the semi-major axis in au,  is the eccentricity,  is the inclination,  is the argument of pericenter, Ω is the longitude of the ascending node, and  is the mean anomaly.The orbital elements of planets listed above the double line are the same for all six runs.

Table 3 .
The simulation results with  Earth = 10 ⊕ for  = 10 6 years.The first column contains the names of the simulations.Columns 2 and 3 show the collision probability with an inflated Earth  10 =  E /(  pp +  pl +  tp ) and  ★ 10 =  E / event .Columns 4 and 5 list the approximate rate with actual size Earth  =  10 /  and  ★ =  ★ 10 /  , where  = 20.

Table 4 .
The results range between ∼ 1.19 and ∼ 2.78 'ocean' quantities, which is in terms of magnitude consistent with the previous value of ML's estimation of 8 M o .If we take  n = 3, we obtain that  wd falls within the [3.57, 8.34] M o interval.It should be noted that these