Importance of Giant Impact Ejecta for Orbits of Planets Formed during the Giant Impact Era

Terrestrial planets are believed to be formed via giant impacts of Mars-sized protoplanets. Planets formed via giant impacts have highly eccentric orbits. A swarm of planetesimals around the planets may lead to eccentricity damping for the planets via the equipartition of random energies (dynamical friction). However, dynamical friction increases eccentricities of planetesimals, resulting in high velocity collisions between planetesimals. The collisional cascade grinds planetesimals to dust until dust grains are blown out due to radiation pressure. Therefore, the total mass of planetesimals decreases due to collisional fragmentation, which weakens dynamical friction. We investigate the orbital evolution of protoplanets in a planetesimal disk, taking into account collisional fragmentation of planetesimals. For 100 km sized or smaller planetesimals, dynamical friction is insignificant for eccentricity damping of planets because of collisional fragmentation. On the other hand, giant impacts eject collisional fragments. Although the total mass of giant impact ejecta is 0.1–0.3 Earth masses, the largest impact ejecta are ∼1000 km in size. We also investigate the orbital evolution of single planets with initial eccentricities of 0.1 in a swarm of such giant impact ejecta. Although the total mass of giant impact ejecta decreases by a factor of 3 in 30 Myr, eccentricities of planets are damped down to the Earth level (∼0.01) due to interaction with giant impact ejecta. Therefore, giant impact ejecta play an important role for determination of terrestrial planet orbits.


Introduction
In the standard scenario for terrestrial planet formation, Marssized protoplanets are formed prior to the gas depletion of the protoplanetary disk (in several Myr) with large orbital separations ∼10 mutual Hill radii (e.g., Kobayashi & Dauphas 2013). The gas depletion triggers the long-term orbital instability of protoplanets and the collisions between protoplanets induced by the orbital instability result in the formation of Earth-or Venussized planets, which is called the giant impact stage (Chambers & Wetherill 1998;Iwasaki et al. 2001).
Most of the masses of the terrestrial planets in the solar system are in Earth and Venus, which have low eccentricities of 0.017 and 0.007, respectively. The largest planets formed in orbital simulations for giant impact stages have much greater eccentricities and inclinations than those of Earth or Venus (Chambers 2001;Kokubo et al. 2006). Those eccentricities and inclinations are possible to be damped via the equipartition of random energies (dynamical friction) with surrounding planetesimals (O'Brien et al. 2006;Raymond et al. 2009;Morishima et al. 2010). However, the surface density of surrounding planetesimals decreases via the collisional cascade of the planetesimals , which reduces the efficiency of dynamical friction. Therefore, collisional fragmentation plays an important role in this issue.
On the other hand, a series of giant impacts eject fragments with a total mass comparable to Earth, resulting in the increase of infrared emission. Therefore, giant impact ejecta explain infrared excesses caused by warm debris disks around 1 au , while cold debris disks beyond 10 au may be related to collisional fragmentation in planetesimal disks induced by planet formation (Kobayashi & Löhne 2014). Such giantimpact-ejecta disks may affect the orbital evolution of protoplanets. The evolution of total masses of giant impact ejecta is controlled by the collisional cascade. Therefore, we need to consider the orbital evolution and the collisional cascade simultaneously.
The orbital evolution of protoplanets in the giant impact stage is mainly treated by N-body simulation. However, all fragments produced via collisional fragmentation cannot be treated individually by N-body simulation because of computational limitation. Therefore, one applies the superparticle approximation where a superparticle represents a large number of planetesimals and fragments. This method is applied for planet formation (Levison et al. 2012;Chambers 2013;Morishima 2015;Clement et al. 2019;Deienno et al. 2019;Walsh & Levison 2019) and for debris disks (Kral et al. 2013;Nesvold et al. 2013). In this paper, we newly develop the N-body code including the mass evolution of planetesimals via collisional cascade, which allows us to evaluate dynamical friction and collisional fragmentation in the giant impact stage. In Section 2, we explain the method to develop the code. In Section 3, we conduct test calculations for the collisional cascade and validate the method. In Section 4, we perform simulations for the orbital evolution of protoplanets in planetesimal or giant-impact-ejecta disks using the newly developed code. In Section 5, we discuss the effect of remnant planetesimals and giant impact ejecta on the orbital evolution of protoplanets in the giant impact stage. We summarize our findings in Section 6.

Method
In the giant impact stage, the orbital evolution and collisions of protoplanets occur. The gravitational interaction with a planetesimal disk is important for the final orbits of protoplanets. However, the disk mass of planetesimals decreases due to the collisional cascade starting from planetesimal fragmentation. Therefore, we need to treat the evolution of orbits and masses consistently.
We apply the superparticle approximation for planetesimals and smaller bodies ejected by collisional fragmentation, where a superparticle represents planetesimals and collisional fragments. Meanwhile, we apply a single particle for a single protoplanet. We numerically integrate the equations of motion of particles via the fourth-order Hermite scheme (Makino & Aarseth 1992;Kokubo & Makino 2004). The orbital integration allows us to accurately treat dynamical evolution and direct collisions between protoplanets and superparticles. 1 However, the number of superparticles that we apply is much smaller than that of planetesimals and fragments with which we are concerned so that statistical treatment is required for accurate calculation of interactions between superparticles. Morishima (2015) developed a method for the collisions and dynamical interactions between superparticles in N-body simulations. We treat the collisions between superparticles following Morishima (2015), while we ignore the dynamical interaction because the collisional timescale is much shorter than the dynamical interaction timescale for dynamically hot planetesimals.
We consider the jth superparticle at the cylindrical coordinate (r j , θ j , z j ) originated at the host star with mass M * . The surface density around the superparticle, Σ j , is determined by the total mass of superparticles in the area with where N j n, is the number of superparticles in the area and m k is the mass of the kth superparticle in the area. We choose δr and dq according to the accuracy of mass evolution of planetesimals due to the collisional cascade, which we discuss in Section 3.
The relative velocity, which characterizes collisional fragmentation of planetesimals in superparticles, is determined by the orbital elements of superparticles in the area. We calculate the relative velocity v j between planetesimals in the jth superparticle with semimajor axis a j , eccentricity e j , inclination i j , the longitude of pericenter ϖ j , and the longitude of ascending node Ω j , given by The collisional cascade grinds planetesimals down to micron-sized grains, which are blown out due to the radiation pressure of the host star (e.g., Kobayashi et al. 2008Kobayashi et al. , 2009Krivov 2010). Collisional fragmentation of planetesimals thus reduces the surface density of planetesimals due to the collisional cascade. We consider the quasi-steady-state collisional cascade, for which  derived the reduction rate of the surface density of planetesimals, given by where m c is the mass of largest bodies in the collisional cascade, α c is the index of the mass distribution of bodies determined by the collisional cascade, f (α c ) is the dimensionless value dependent on α c , and Q D * is the specific impact energy needed for the ejection of the half-mass of colliders. The mass of largest bodies m c does not change in the collisional cascade, which is valid until the bodies with m c exist. For The mass distribution with the index α c is achieved in the timescale for the collisional cascade. An earlier mass distribution index may be different from α c in Equation (6) . However, the mass loss due to the collisional cascade is negligible in such an early stage, and the mass evolution mainly occurs when the mass distribution given by Equation (6) is achieved. Therefore, Equations (5) and (6) are valid to treat the surface density evolution accurately. According to , f (α c ) is given by where ρ is the density of bodies, b is the power-law index of the mass distribution of fragments, ò is the constant determining the largest mass of fragments, and 2 ò The impact laboratory experiments show b=1.5-1.7 and  0.1 (e.g., Takagi et al. 1984;Nakamura & Fujiwara 1991) so that the choice of b and ò insignificantly change the value of f (α c ). 3 Therefore, we set = b 3 5 and ò=0.2. If v r is fixed, the integration of Equation (5) over time t gives where τ 0 =−Σ(0)/Ṡ (0). In the giant impact stage, Equation (11) is not always valid because of the evolution of v r . Therefore, we use Equation (11) only for the validation of our simulation in Section 3. Based on Equation (5), the mass-loss rate of the jth superparticle due to the collisional cascade is given by a where m j c, is the mass of the largest bodies in the jth superparticle. We calculate the mass evolution of superparticles via the integration of Equation (12) using Equations (1) and (2).
The masses of the largest bodies in superparticles, m c , are set to 10 16 g corresponding to 1 km in radius. Therefore, Q D * of such a body is mainly determined by shuttering and gravitational reaccumulation (Benz & Asphaug 1999;Leinhardt & Stewart 2012;Genda et al. , 2017Jutzi 2015;Suetsugu et al. 2018). Therefore, Q D * has a monotonous increasing function of m c , which is assumed to be (e.g., Benz & Asphaug 1999 where Q 0 and p q are determined from numerical simulations. The collisional cascade is controlled not by individual collisions but by successive collisions. The value averaged over impact angles is applied for Q D *. Taking into account selfgravity and a model of rock fractures, Benz & Asphaug (1999) obtained the averaged values for Q 0 and p q via the smoothed particle hydrodynamics (SPH) simulation. Recently, the dependence of Q D * on the number of SPH particles is argued and the value of Q D * decreases 20% in the limit of highresolution simulation . However, the friction of damaged rock, which was not considered in Benz & Asphaug (1999), increases Q D * (Jutzi 2015). The value of Q D * obtained in Benz & Asphaug (1999) is similar to that given by the high-resolution simulations with the friction (Suetsugu et al. 2018). Therefore, we set Q 0 =9.5×10 8 erg g −1 and p q = 0.453 (Benz & Asphaug 1999).

Validation for Collisional Cascade
We perform a simulation for collisional evolution of a planetesimal disk composed of 2000 superparticles with a= 0.975-1.025 au, δr=0.01au, and m c =10 16 g. The radial distribution of superparticles is put according to Σ(a) ∝ a −1 and e and i have the Rayleigh distributions with mode values e=0.01 (a/1 au) 1/2 and = i a 0.005 1 au 1 2 ( ) . The collisional cascade timescale in this setting is estimated to be t » 0.6 0 yr. We evaluate the accuracy of the method from the comparison with the analytic solution in Equation (11) at t=100 yr.
The accuracy of the simulation depends on the number of particles in the neighbor area, N n . We calculate the relative error of the surface density from comparison of the simulation with the analytical solution. Figure 1 shows the dependence of relative errors on the number of superparticles. The errors are roughly proportional to the number of superparticles, which is proportional to N n for fixed δθ. However, the errors are almost independent of δθ. Although N n decreases with decreasing δθ, the number of particles having the experience in a neighboring area is almost the same in the time much longer than orbital periods if the number of superparticles and δr are fixed. The timescale of the collisional cascade is mainly much longer than orbital periods. Therefore, an accuracy better than 10% is obtained if the number of particles in the annulus with width δr is 10.

Orbital Evolution of Protoplanets in a Swarm of Fragmenting Planetesimals
We carry out simulations for the orbital evolution of three protoplanets with an Earth mass M ⊕ in a planetesimal disk with 30M ⊕ composed of 3000 superparticles around the central star with solar mass M e . The semimajor axis of the intermediate protoplanet is initially set at 1 au, the orbital separation of protoplanets is 10 mutual Hill radii, and their eccentricities and inclinations are 0.03 and 0.015, respectively. The radial distribution of superparticles is initially put according to Σ(a) ∝ a −1 and their e and i have the Rayleigh distributions with mode values e=0.03 and i=0.015, respectively. The intermediate radius and width of the planetesimal disk are 1 au and 30 mutual Hill radii of the protoplanets, respectively. For the treatment of collisional fragmentation, we put δr=0.01 au and δθ=π/8.
The collisional cascade is characterized by planetesimal mass m c . We set m c =10 16 g (≈1 km in radius) for collisional fragmentation, while we have an additional simulation without collisional fragmentation ( show the orbital distribution of protoplanets and superparticles at t=10 3 yr with collisional fragmentation, while Figures 2(c) and (d) are that without fragmentation. Dynamical friction decreases e and i of protoplanets, while e and i of planetesimals increase. Increases in e and i of planetesimals activate their collisional cascade, which reduces the surface density of planetesimals. Collisional fragmentation weakens dynamical friction so that e and i of protoplanets with fragmentation remain higher than those without fragmentation.
Figures 3-5 show the evolution of the rms of eccentricities á ñ e 2 1 2 and inclinations á ñ i 2 1 2 for planetesimals and protoplanets in the same initial setting of planetesimal disks and protoplanets as Figure 2. If we ignore collisional fragmentation, á ñ e 2 1 2 and á ñ i 2 1 2 of planetesimals become much larger than those of protoplanets ( Figure 3). This is caused by dynamical friction. However, collisional fragmentation of planetesimals decreases the surface density of planetesimals (Figures 4 and  5), which suppresses the e and i reduction of protoplanets via dynamical friction. Collisional fragmentation is effective for small m c so that á ñ e 2 1 2 and á ñ i 2 1 2 of protoplanets insignificantly change. On the other hand, the e and i evolution for planetesimals are almost independent of fragmentation. This is because the viscous stirring of protoplanets controls e and i of planetesimals.
If we consider two populations of bodies such as protoplanets with á ñ e 2 1 1 2 , á ñ i 2 1 1 2 , and mass m p,1 and planetesimals with á ñ e 2 2 1 2 , á ñ i 2 2 1 2 , and mass m p,2 , the time differential of á ñ e 2 α and á ñ i 2 α for α=1 or 2 (protoplanets or planetesimals) due to dynamical friction and viscous stirring is analytically given by (Ohtsuki et al. 2002) where a 0 is the mean semimajor axis of bodies, Ω K,0 is the Keplarian frequency at a 0 , N s,β is the surface number density of bodies for β=1 or 2, ] is the reduced Hill radius of bodies with masses a m p, and b m p, , P VS , P DF , Q VS , and Q DF are the efficiencies for viscous stirring and dynamical friction for á ñ e 2 1,2 and á ñ i 2 1,2 , respectively. The analytical formulae of P VS , P DF , Q VS , and Q DF are given as a function of á ñ e 2 1,2 and á ñ i 2 1,2 in Ohtsuki et al. (2002). The e and i variation rates due to dynamical friction are given by the second terms on the right-hand sides of Equations (14) and (15) where Σ is the surface density of planetesimals and we assume á ñ e 2 1 ≈á ñ e 2 2 and á ñ i 2 1 ≈á ñ e 2 2 . The damping rates are mainly    determined by Σ but are almost independent of the choice of planetesimal masses, m p,2 . This means the superparticle approximation is valid. However, dynamical friction leads to the equipartition of the random energies such as á ñm e p,1 2 1 1 2 á ñ m e p,2 2 2 1 2 and á ñ~á ñ m i m i p,1 2 1 1 2 p,2 2 2 1 2 . The values for the equipartition depend on the choice of m p,2 . Therefore, we need to take care with the choice of superparticle masses if we are interested in the equipartition.
To compare the results of simulations with Equations (14) and (15), we set m p,2 to be superparticle mass m j and integrate Equations (14) and (15) over time (see Figure 3). The analytic solution is consistent with the simulation. For t  100 yr, dynamical friction is ineffective because of the achievement of energy equipartition such as á ñ~á ñ m e m e p,1 2 1 1 2 p,2 2 2 1 2 and á ñ~á ñ m i m i p,1 2 1 1 2 p,2 2 2 1 2 . The result depends on the superparticle mass as discussed above. However, collisional fragmentation mainly occurs prior to the achievement of the equipartition (see Figures 4 and 5).
The mass of planetesimals m p,2 in Equations (14) and (15) is initially set to be superparticle mass m j , although m c in Equation (5) is independently given. We then integrate Equations (5), (14), and (15) over time. The analytic solution is roughly in agreement with the results of simulations (see Figures 4 and 5). In the simulations, each superparticle has an independent mass. Superparticles in the inner disk effectively lose their masses via collisional fragmentation because of high collisional speeds, while superparticles tend to have large masses in the outer disk. On the other hand, the mass evolution of planetesimal disks is calculated with the averaged collisional speed in the analytical solution. Therefore, the analytic solution cannot perfectly reproduce the simulations. However, the tendency of orbital interaction and collisional fragmentation is understood from the analytical solution.
It should be noted that the relative increases in protoplanet masses are smaller than 0.2 and 0.02 in the cases without and with collisional fragmentation, respectively. Those are caused by collisions with planetesimals. We estimate the effect of collisional damping due to the collisional accretion according to the model by Kobayashi et al. (2016), resulting in the e damping of 0.006 and 0.0001 for protoplanets without and with collisional fragmentation, respectively. The collisional damping insignificantly affects the orbital evolution of protoplanets. Therefore, the analytical solution without collisional damping is in good agreement with the results of the simulations.
The timescale of dynamical friction is estimated from Equations (14) and (15) where M disk is the total mass of a planetesimal disk, Da is the width of the disk, and = r m M a 3 H p , 1 1 3 0 * ( ) is the Hill radius for protoplanets. Note that P DF is estimated to be ∼10 for m p,1 ∼M ⊕ and á ñ e 2 1 1 2 ≈á ñ e 2 2 1 2 ∼0.1. This estimate implies a planetesimal disk with M disk = M ⊕ leads to eccentricity damping for Earth-sized protoplanets in a long timescale ?10 4 yr via dynamical friction. In addition, dynamical friction increases á ñ e 2 1 2 of planetesimals. The increase in á ñ e 2 2 1 2 makes τ df longer.
On the other hand, the decreasing timescale for M disk due to the collisional cascade, τ cc , is estimated from Equation (5) For 100 km sized planetesimals (m c ≈ 10 22 g) with á ñ e 2 2 1 2 ∼ 0.1, τ cc  τ df . It should be noted that dynamical friction increases á ñ e 2 2 1 2 during the á ñ e 2 1 1 2 damping, which shortens τ cc and elongates τ df as discussed above. Therefore, even if τ cc ∼ τ df initially, τ cc eventually becomes much shorter than τ df .
Planets formed via giant impacts have high eccentricities and inclinations. For Earth-sized planets formed in the giant impact stages, the mean eccentricity is ∼0.1 (Chambers 2001;Kokubo et al. 2006), which is much larger than the current eccentricities of Earth and Venus. The depletion timescale of a planetesimal disk with  m 10 g c 21 is shorter than the e damping timescale via dynamical friction (see Equations (17) and (18)). Therefore, larger planetesimals are required for e damping of Earth-sized planets. We consider such large planetesimals as being produced from a giant impact. The mass ratios of the largest giant impact ejecta to parent protoplanets are ∼0.01 . Therefore, we considerm 10 c 26 g, resulting in τ cc ? τ df .
We perform the simulation for the orbital evolution of a high-eccentricity planet with mass M ⊕ in a swarm of giant impact ejecta. Giant impact ejecta initially have similar orbits to parent planets. However, the orbital phases of giant impact ejecta, ϖ and Ω, are eventually distributed uniformly due to perturbation from other planets. The timescale to achieve a uniform distribution for ϖ and Ω, which is roughly given by the precession rates for ϖ and Ω obtained from the secular perturbation theory (Murray & Dermott 1999), is estimated as ∼10 5 -10 6 yr for giant impact ejecta around 1 au perturbed by a Venus-like planet, which is much shorter than the timescale for dynamical friction caused by giant impact ejecta. Therefore, instead of ignoring perturbation from other planets, we uniformly set ϖ and Ω of giant impact ejecta from the beginning. We initially set a=1 au and = = e i 2 0.1 for the planet and a=0.95-1.05 au and á ñ e 2 1 2 =2 á ñ i 2 1 2 =0.1 for a giant-impact-ejecta disk composed of 120 superparticles with total mass 0.2M ⊕ and m c =10 26 g ( Figure 6). The evolutions of á ñ e 2 1 2 and á ñ i 2 1 2 for giant impact ejecta and the planet differ from those predicted by the analytical solutions given by the time integration of Equations (5), (14), and (15). Once á ñ e 2 1 2 and/or á ñ i 2 1 2 are much larger than 0.1, their orbits are controlled by the higher-order terms for á ñ e 2 1 2 and á ñ i 2 1 2 , which are ignored in Equations (14) and (15). The higher-order terms make dynamical friction less effective. Therefore, the variations of á ñ e 2 1 2 and á ñ i 2 1 2 for the planet and planetesimals are less than those expected by the analytic solution. However, e of the planet becomes comparable to or larger than the analytic estimate until 10 Myr and decreases much greater than the analytic estimate in several 10 Myr. This is caused by the orbital energy damping for the planet rather than its angular momentum variation. If the angular momentum is fixed, the energy damping given by −Δa/a results in D = e -D e a ae 1 2 2 ( ) , where Δe and Δa are the changes in e and a, respectively. Therefore, even a small energy damping of D »a a 0.02 leads to Δe≈−0.1 for the planet. Giant impact ejecta are scattered by the planet and tend to stay in the outer disk so that the planet loses the orbital energy via the scatterings. This is clearly shown in the semimajor-axis evolution for the planet and giant impact ejecta ( Figure 6). The collisional damping between the planet and giant impact ejecta also reduces the orbital energy of the planet. The accretion mass of giant impact ejecta onto the planet is for t=30-100 Myr. Although the e damping is significant after 30 Myr, the collisional damping expected from the accretion mass is slight. Therefore, the energy loss for the planet due to the scattering of giant impact ejecta mainly decreases its eccentricity.
The mass evolution of giant impact ejecta due to collisional fragmentation is less important than that estimated by Equation (5). Increase in e of giant impact ejecta due to planetary perturbation makes the giant-impact-ejecta disk wider. This decreases the surface density of giant impact ejecta, which reduces the efficiency of collisional fragmentation. In addition, collisions between high-eccentricity giant impact ejecta are more frequent around the planetary orbit. Collisional fragmentation between giant impact ejecta with similar ϖ and Ω thus mainly occurs so that those relative velocities are comparatively small (see Equations (3) and (4)). The mass loss by collisional fragmentation is thus less effective due to the increase of ejecta eccentricity. Therefore, the giantimpact-ejecta disk is maintained in a long timescale 10 Myr, resulting in the significant eccentricity damping for the planet.

Discussion
In the giant impact stage, orbital instability of protoplanets occurs in a timescale longer than 10 Myr (Chambers et al. 1996;Chambers & Wetherill 1998;Iwasaki et al. 2001Iwasaki et al. , 2002Kominami & Ida 2002;Kokubo et al. 2006). The protoplanets formed from the long-term evolution mainly have high eccentricities (∼0.1; Chambers 2001; Kokubo et al. 2006), which are much larger than the current eccentricities of Earth and Venus. A swarm of planetesimals may damp the eccentricities of protoplanets due to dynamical friction (O'Brien et al. 2006;Raymond et al. 2009;Morishima et al. 2010). Instead, dynamical friction increases eccentricities of planetesimals, which induces collisional fragmentation between planetesimals. The collisional cascade grinds planetesimals until small fragments are blown out by radiation pressure, which results in the mass loss of the planetesimal disk. Collisional fragmentation may thus suppress the eccentricity damping for protoplanets. Therefore, we have investigated the orbital interaction between protoplanets and planetesimals, taking into account collisional fragmentation.
Eccentricities and inclinations of protoplanets are damped via dynamical friction if we ignore collisional fragmentation (see Figures 2(c) and (d)). However, collisional fragmentation weakens dynamical friction (see Figures 2(a) and (b)). The mass loss due to collisional fragmentation depends on the typical planetesimal size of which planetesimals mainly determine the total mass of the planetesimal disk. Collisional depletion occurs in a short timescale for planetesimal disks with small typical planetesimal sizes. For 100 km sized or smaller planetesimals, the eccentricity damping is ineffective due to collisional fragmentation (Figures 4 and 5).
The primordial planetesimals may remain even in the giant impact stage. Although the typical size of primordial planetesimals is not unknown, the size may be on the order of 100 km, similar to that of Main Belt asteroids (e.g., Kobayashi et al. 2016). Due to collisional fragmentation, dynamical friction for eccentricity damping is ineffective in such a planetesimal disk. On the other hand, the increase in the eccentricities of protoplanets is required for the onset of the orbital instability of protoplanets for giant impacts (Iwasaki et al. 2001). Although primordial planetesimal disks may not explain the small eccentricities of Earth and Venus, such a disk maintains the orbital instability of protoplanets until Earth-and Venus-sized planets are formed (Walsh & Levison 2019). In addition, the direct formation of Earth and Venus via the accretion of planetesimals without the giant impact stage is inhibited due to the depletion of remnant planetesimal disks via collisional fragmentation Kobayashi & Dauphas 2013). Therefore, collisional fragmentation supports the giant impact scenario to form Venus and Earth.
Giant impacts lead to the ejection of fragments as well as collisional growth of protoplanets. The collisional ejecta from single impacts have 0.1-0.3M ⊕ , while the typical size of giant impact ejecta can be as large as 1000 km . The simulation for the interaction between a protoplanet and giant impact ejecta shows significant eccentricity damping of the protoplanet in 30 Myr (see Figure 6). The collisional fragmentation is less effective in an originally narrow planetesimal disk for high-eccentricity planetesimals, while the dynamical friction is also less effective for high eccentricities and inclinations of planetesimals. As a result, eccentricity damping for planets occurs in a long timescale ∼100 Myr.
Giant impact ejecta are produced even in the early giant impact stage . However, the orbital separations of protoplanets are so narrow that giant impact ejecta are distributed widely due to interactions with other protoplanets. In addition, the eccentricity damping timescale due to dynamical friction is much longer than the collisional timescale between protoplanets. Therefore, sequent giant impacts occurs due to insignificance of dynamical friction by giant impact ejecta. On the other hand, stable orbital configurations of protoplanets are achieved after tens of giant impacts. Such orbital separations are wide enough to keep the orbital concentration of giant impact ejecta around the orbits of single protoplanets. Therefore, dynamical friction by giant impact fragments is effective in the late giant impact stage, which may form low-eccentricity planets.
Finally we need to discuss the accuracy of the mass loss due to collisional fragmentation. As we discuss above, the collisional mass loss mainly occurs after the collisional cascade is achieved. In the collisional cascade, the mass loss is mainly determined by the total ejecta mass from single collisions, and is insensitive to the mass distribution of ejecta . Therefore, the uncertainty of collisional fragmentation mainly comes from Q D * at the typical sizes of planetesimals (see Equation (18)). As discussed in the previous studies (Kobayashi et al. 2016;, Q D * of 10 km sized or smaller primordial planetesimals may be much larger than Q D * we set in the simulations. However, taking into account the enhancement of Q D *, 100 km sized or smaller primordial planetesimals insignificantly work for the e damping of protoplanets. On the other hand, Q D * of 1000 km sized bodies are mainly determined by the self-gravity of colliding bodies (Kobayashi et al. , 2011Genda et al. , 2017Suetsugu et al. 2018). The uncertainty is small for Q D * of largest ejecta of giant impacts. Therefore, the e damping of planets formed in the giant impact stage is likely to be caused by giant impact ejecta.

Summary
Terrestrial planets are formed via giant impacts between Mars-sized protoplanets. The resultant planets have larger eccentricities than the current values for Earth and Venus. Dynamical friction with a planetesimal disk is expected to damp the eccentricities of protoplanets. On the other hand, the collisional cascade of planetesimals and blowout of small collisional fragments by radiation pressure decrease the planetesimal disk mass, which weakens dynamical friction with planetesimals. Therefore, we have investigated the orbital evolution of planets with collisional fragmentation. Our findings are as follows.
1. We have developed an N-body simulation code involving the mass loss due to the collisional cascade. We have calculated the mass evolution of a planetesimal disk composed of superparticles using the code, which reproduces the analytical solution for the mass loss due to the collisional cascade with high accuracy. 2. We have investigated the evolution of orbits and masses of protoplanets and planetesimals via gravitational interaction and collisional fragmentation using the N-body code. If collisional fragmentation is ignored, dynamical friction damps eccentricities and inclinations for protoplanets. However, collisional fragmentation suppresses dynamical friction. The timescale ratio of dynamical friction to collisional fragmentation depends on the typical planetesimal size but not on the disk mass. Even a massive planetesimal disk composed of 100 km sized or smaller planetesimals cannot damp eccentricities of planets. Therefore, the disks composed of primordial planetesimals are ineffective for the eccentricity damping for planets. 3. Giant impacts eject collisional fragments. The total masses of giant impact ejecta are several tenths of colliding planets. The typical size of the largest ejecta is ∼1000 km. We have carried out simulations for a planet with an Earth mass in the disk with 0.2 Earth masses composed of giant impact ejecta with the largest ejecta 10 26 g. Collisional mass loss is insignificant for such a large typical size. Eccentricity of the planet is damped from 0.1 to ∼0.01 due to interactions with giant impact ejecta in 30 Myr. Therefore, giant impact ejecta are possible to decrease the eccentricities of planets formed via giant impacts, comparable to those of Earth and Venus.