© The Authors 2024

1 Introduction

The formation history of giant planets is still not fully understood. Many theoretical studies investigate the accretion of heavy elements onto the planet at various stages since the observed atmospheric composition of giant planets is thought to be related to the formation history. Heavy-element accretion is now known to occur mainly before the dissipation of the protoplanetary disk. The accretion of dust (Morbidelli et al. 2023), planetesimals (Shibata et al. 2020, 2022a; Hands & Helled 2021; Turrini et al. 2021; Knierim et al. 2022; Pacetti et al. 2022), and disk gas accretion (Madhusudhan et al. 2014; Booth et al. 2017; Booth & Ilee 2019; Notsu et al. 2020; Schneider & Bitsch 2021a; Schneider & Bitsch 2021b) are all regarded to be critical processes that can affect the planetary composition. The above studies suggest that the ratio of molecular species in the planetary atmosphere could be a tracer of the initial formation location of giant planets. However, these studies focus on the accretion of heavy elements before the disk dissipates. In reality, accretion could also occur at later stages and change the composition of the planetary atmosphere. Therefore, it is important to investigate the effects of solid accretion after disk dissipation to link the initial formation location of giant planets to their atmospheric composition.

Heavy-element accretion post-disk dissipation depends on the orbital evolution of giant planets and the remnant solid disk. By the end of the disk dissipation, solid particles like planetesimals are eliminated from a feeding zone due to the aerodynamical drag from the protoplanetary disk gas. This means that no solid accretion occurs without the scattering process triggered by the orbital instability of planets. Planetesimal scattering by giant planets is typically investigated using N-body simulations (Matter et al. 2009; Frewen & Hansen 2014; Mustill et al. 2018; Seligman et al. 2022; Rodet & Lai 2023). Seligman et al. (2022) consider the accretion of comets, which are planetesimals scattered by the other giant planets. Their analytical model suggests that a close-in giant planet can accrete a significant fraction of comets before they are scattered from the planetary system. However, they assume two groups of giant planets: (i) giant planets that scatter the remnant solid disk in the outer disk region, and (ii) giant planets that accrete the scattered solid disk in close-in orbits. Before being accreted by the inner giant planets, some solids can be accreted by the outer giant planets.

Rodet & Lai (2023) investigate the scattering and accretion of planetesimals around an eccentric giant planet. They find that the planet can accrete some scattered planetesimals and that the accretion probability changes with the planet’s semi-major axis, eccentricity, mass, and radius. However, this study focuses on the planetesimals that fall onto the central star. The initial inclination of the planetesimal disk and the mutual inclination between the planetesimal disk and the planet should affect the accretion probability, but this has not been investigated in previous studies. Also, they start the simulations with the planet whose orbit was already crossed to the planetesimal disk. Since the planetary eccentricity is gradually excited by the mutual interactions with the other planets after the disk dissipation (Carrera et al. 2019; Anderson et al. 2019; Garzón et al. 2022), in reality, the scattering of the remnant solid disk starts before the orbits are crossed. It is, therefore, still unknown how much solid material can be accreted by the planet during the scattering process.

In order to investigate the possible accretion of a remnant solid disk, we perform N-body simulations and analyze how the orbits of planetesimals and embryos evolve around an eccentric giant planet. First, we derive the analytical equation for the accretion probability of solid particles onto an eccentric giant planet in Sec. 2. In Sec. 3, we compare the derived formulae to the N-body simulations. Following Rodet & Lai (2023), we begin the simulations with a planet with a finite eccentricity. In order to model the evolution of giant planets in more detail, we investigate the accretion probability when the planetary eccentricity increases with time in Sec. 4. We discuss our results, the limitations of our analytical model, and the effects of solid accretion on the planetary composition in Sec. 5. We summarize our results and present the key conclusions in Sec. 6.

2 Accretion probability to an eccentric planet: Analytical investigations

We start with a simple situation where a giant planet crosses the orbits of surrounding planetesimals, following (Rodet & Lai 2023). We neglect the mutual interaction between planetesimals; therefore, our analysis can be applied to the smaller solid particles, such as debris and dust. The mass M_p, radius R_p, semimajor axis a_p, eccentricity e_p, and inclination i_p of the planet are fixed values.

2.1 Orbital evolution around an eccentric giant planet

First, we consider the gravitational perturbation from the eccentric planet to the planetesimals at the orbital radius r. In a-particle-in-a-box approximation, a velocity perturbation by the passing planet is given by (e.g. Armitage 2010): $$\delta v\sim \frac{2\mathcal{G}{M}_{\mathrm{p}}}{b{v}_{\mathrm{r}\mathrm{e}\mathrm{l}}},$$(1)

where $\mathcal{G}$ is the gravitational constant, b is the impact parameter, and v_rel is the relative velocity to the passing planet. Here, we assume that the impact parameter is the order of r, and the relative velocity is e_pv_K(r), where v_K is the Kepler velocity. The weak perturbation from the passing planet is given by: $$\delta v_{\text{weak }}\sim v_{\mathrm{K}}(r)\frac{2}{e_{\mathrm{p}}}\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}},$$(2)

and the eccentricity perturbation on a planetesimal is approximated to be: $$\delta e_{\text{weak }}\sim C_{\mathrm{v}}\frac{\delta v_{\text{weak }}}{v_{\mathrm{K}}(r)}=\frac{2C_{\mathrm{v}}}{e_{\mathrm{p}}}\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}},$$(3)

where C_v is a fitting parameter. The eccentricity excitation occurs in every orbit of the eccentric planet, and the excitation rate is independent of the radial distance from the central star. Due to the large mass difference between the planet and planetesimals, the perturbations increase the mean eccentricity of the planetesimals ⟨e⟩ linearly. The eccentricity evolution can be written as: $$\frac{\mathrm{d}⟨e⟩}{\mathrm{d}t}\sim \frac{\mathrm{\Delta }e}{\mathrm{\Delta }t}=\frac{\delta e_{\mathrm{w}\mathrm{e}\mathrm{a}\mathrm{k}}}{P_{\mathrm{o}\mathrm{r}\mathrm{b},\mathrm{p}}},$$(4) $$⟨e⟩\sim ⟨e_0⟩+\delta e_{\text{weak }}\frac{t}{P_{\text{orb,p }}},$$(5)

where P_orb,p is the orbital period of the eccentric planet and ⟨e₀⟩ is the initial mean eccentricity of planetesimals.

Different from eccentricity excitation, inclination excitation requires close encounters with the eccentric planet (e.g. Goldreich et al. 2004). When the eccentric planet passes through the planetesimal disk with sin i_p < sin i, a z-component of the velocity perturbation is given by: $$\delta v_{\perp }=\delta v\frac{r\sin i}{b}.$$(6)

We define the distance of the close encounter b_close as the distance at which the velocity perturbation doubles the planetesimal’s inclination. Therefore, $$C_{\mathrm{V}}\frac{\delta v_{\perp }}{v_{\mathrm{K}}}=\sin i,$$(7) $$\frac{b_{\text{close }}}{r}=\sqrt{\frac{2C_{\mathrm{V}}}{e_{\mathrm{p}}}\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}}},$$(8) $$=0.062{\left(\frac{C_{\mathrm{v}}}{2}\right)}^{1/2}{\left(\frac{e_{\mathrm{p}}}{0.99}\right)}^{-1/2}{\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{J}\mathrm{u}\mathrm{p}}}\right)}^{1/2}{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^{-1/2}.$$(9)

At the beginning of the simulation, the thickness of the planetesimal disk a sin i is smaller than b_close. The typical timescale of the close encounter τ_close,2D is given by: $${\tau }_{\text{close },2\mathrm{D}}=\frac{2\pi r}{2b_{\text{close }}}{P}_{\text{orb},\mathrm{p}},$$(10) $$\sim 51P_{\text{orb,p }}{\left(\frac{C_{\mathrm{v}}}{2}\right)}^{-1/2}{\left(\frac{e_{\mathrm{p}}}{0.99}\right)}^{1/2}{\left(\frac{M_{\mathrm{p}}}{M_{\text{Jup }}}\right)}^{-1/2}{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^{1/2}.$$(11)

Close encounters occur in three dimensions once the disk thickness exceeds b_close. In this case, the typical timescale of the close encounter τ_close,3D is given by: $${\tau }_{\text{close,}3\mathrm{D}}=\frac{2\pi r\cdot r\sin i}{\pi b_{\text{close }}{}^2}{P}_{\text{orb,p }}.$$(12)

The time evolution of mean inclination ⟨sin i⟩ can be written as follows: $$\frac{\mathrm{d}⟨\sin i⟩}{\mathrm{d}t}\sim \frac{\mathrm{\Delta }\sin i}{\mathrm{\Delta }t}=\{\begin{array}{l}\frac{⟨\sin i⟩}{{\tau }_{\text{close,2D}}}\text{ for }⟨\sin i⟩<\sin i_{\mathrm{b}}\\ \frac{⟨(\sin i⟩}{{\tau }_{\text{close,3D }}}\text{ for }\sin i_{\mathrm{b}}<⟨\sin i⟩\end{array},$$(13) $$⟨\sin i⟩\sim \{\begin{array}{l}⟨\sin i_0⟩\exp \left(\frac{t}{{\tau }_{\text{close},2\mathrm{D}}}\right)\text{ for }⟨\sin i⟩<\sin i_{\mathrm{b}}\\ \frac{b_{{\text{clase}}^2}}{4r^2}\frac{t-t_{\mathrm{b}}}{P_{\text{orb,},\mathrm{p}}}+\sin i_{\mathrm{b}}\text{ for }\sin i_{\mathrm{b}}<⟨\sin i⟩\end{array}$$(14)

with $$\sin i_{\mathrm{b}}=\frac{b_{\text{close }}}{r},$$(15) $$t_{\mathrm{b}}={\tau }_{\text{close,2D }}\log \frac{\sin i_{\mathrm{b}}}{⟨\sin i_0⟩}.$$(16)

Here, ⟨sin i₀⟩ is the mean initial inclination of planetesimals.

2.2 Accretion to close encounter ratio

The planetesimal accretion probability depends on the geometry of the planetesimal disk. The capture radius of an eccentric planet is given by: $$R_{\text{cap }}=R_{\mathrm{p}}{F}_{\text{grav }}=R_{\text{cap }}{(1+\frac{v_{{\text{esc }}^2}}{2v_{{\text{rel}}^2}})}^{1/2}\text{, }$$(17)

where F_grav is the gravitational focusing factor, v_esc and v_re1 are the escape velocity of the planet and the relative velocity between the planet and planetesimals. When the half thickness of the planetesimal disk r sin i is smaller than the capture radius R_cap, particle accretion occurs two-dimensionally (2D). The accretion probability per orbit of the eccentric planet is given by: $$f_{\mathrm{a}\mathrm{c}\mathrm{c},2\mathrm{D}}=\frac{R_{\mathrm{c}\mathrm{a}\mathrm{p}}}{\pi r}.$$(18)

Once the thickness of the particle disk exceeds R_cap, the particle collision occurs in a 3D manner. The accretion probability is then given by: $$f_{\mathrm{a}\mathrm{c}\mathrm{c},3\mathrm{D}}=\frac{R_{\text{cap}}{}^2}{4\pi r^2\sin i}\text{. }$$(19)

The close encounters with the eccentric planet increase the planetesimals’ inclination. The close encounter probability per planetary orbit is given by: $$f_{\text{close},2\mathrm{D}}=\frac{b_{\text{close}}}{\pi r}.$$(20)

When the thickness of the planetesimal disk exceeds b_close, the close encounter probability is: $$f_{\text{close},3\mathrm{D}}=\frac{b_{\text{close}}{}^2}{4\pi r^2\sin i}.$$(21)

Using Eqs. (18), (19), (20), and (21), we obtain the accretion fraction per single close encounter f_a/_c to be: $$f_{\mathrm{a}/\mathrm{c}}=\{\begin{array}{ll}\frac{f_{\text{acc},2\mathrm{D}}}{f_{\text{close},2\mathrm{D}}}=\frac{R_{\text{cap }}}{b_{\text{close }}}& \text{ for }\sin i<\frac{R_{\text{cap }}}{r},\\ \frac{f_{\text{acc},3\mathrm{D}}}{f_{\text{close},2\mathrm{D}}}=\frac{R_{\text{cap}}}{b_{\text{close}}}\frac{R_{\text{cap}}}{r\sin i}& \text{ for }\frac{R_{\text{cap}}}{r}<\sin i<\frac{b_{\text{close}}}{r},\\ \frac{f_{\text{acc},3\mathrm{D}}}{f_{\text{close},3\mathrm{D}}}={\left(\frac{R_{\text{cap}}}{b_{\text{close}}}\right)}^2& \text{ for }\frac{b_{\text{close }}}{r}<\sin i.\end{array}$$(22)

2.3 Accretion probability

The accretion probability of planetesimals onto an eccentric planet $P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }}$ is given by: $$P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }}=min(1,P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }\mathrm{\prime }}),$$(23)

with $$P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }\mathrm{\prime }}={\int }_0^{{\tau }_{\text{eject }}}\frac{f_{a/c}}{{\tau }_{\text{close }}}\mathrm{d}t,$$(24)

where τ_eject is the mean duration time of the planetesimals before they are eliminated. If the eccentricity exceeds e_fall = 1 − R_star/a, where R_star is the radius of the central star and a is the semi-major axis of the planetesimals, the planetesimal will be removed from the planetary system by falling to the central star or ejection from the planetary system. By equating the mean eccentricity given by Eq. (5) to e_fall, we obtain τ_eject as follows: $$\frac{{\tau }_{\text{eject }}}{P_{\text{orb},\mathrm{p}}}=\frac{e_{\mathrm{p}}}{2C_{\mathrm{v}}}\frac{M_{\mathrm{s}}}{M_{\mathrm{p}}}\{(1-\frac{R_{\text{star}}}{a})-⟨e_0⟩\}.$$(25)

Using the initial semi-major axis of a planetesimal a₀, we integrate Eq. (24) and obtain $$P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }\mathrm{\prime }}=P_{2\mathrm{D}}+P_{3\mathrm{D}2\mathrm{D}}+P_{3\mathrm{D}},$$(26)

with $$P_{2\mathrm{D}}=\{\begin{array}{ll}\frac{R_{\text{cap}}}{b_{\text{close}}}\log \frac{R_{\text{cap}}}{a_0⟨\sin i_0⟩}& \text{ for }⟨\sin i_0⟩<\frac{R_{\text{cap}}}{a_0},\\ 0& \text{ other }.\end{array}$$(27) $$\begin{array}{rl}& P_{3\mathrm{D}2\mathrm{D}}=\{\begin{array}{ll}\frac{R_{\text{cap}}}{b_{\text{close}}}(1-\frac{R_{\text{cap }}}{b_{\text{close }}})& \text{ for }\sin i_0<\frac{R_{\text{cap }}}{a_0},\\ \frac{R_{\text{cap }}}{b_{\text{close }}}(\frac{R_{\text{cap }}}{a_0\sin i_0}-\frac{R_{\text{cap }}}{b_{\text{close }}})& \text{ for }\frac{R_{\text{cap }}}{a_0}<\sin i_0<\frac{b_{\text{close }}}{a_0},\\ 0& \text{ other. }\end{array}\end{array}$$(28) $$\begin{array}{rl}& P_{3\mathrm{D}}=\{\begin{array}{ll}{\left(\frac{R_{\text{cap }}}{b_{\text{close }}}\right)}^2& \log \{\frac{1}{4}\frac{b_{\text{close }}}{a_0}(\frac{{\tau }_{\text{eject }}}{P_{\text{orb, p}}}-\pi \log \frac{\sin i_{\text{b}}}{\sin i_0})+1\}\\ {\left(\frac{R_{\text{cap }}}{b_{\text{close }}}\right)}^2& \log \{\frac{1}{4\sin i_0}{\left(\frac{b_{\text{close }}}{a_0}\right)}^2\frac{{\tau }_{\text{eject }}}{P_{\text{orb, p}}}+1\}\\ & \text{ other. }\end{array}\end{array}$$(29)

We assume that the orbital distance of planetesimals r does not change from a₀ during the scattering by the eccentric planet. This is not true for each planetesimal, but numerical simulations in Sec. 3 show that this could be a reasonable assumption by taking an average for a large number of planetesimals.

The inferred accretion probability depends on R_cap/b_close, that is $$R_{\text{cap }}/b_{\text{close }}=\frac{R_{\mathrm{p}}}{r}\sqrt{\frac{e_{\mathrm{p}}}{2C_{\mathrm{v}}}\frac{M_{\mathrm{s}}}{M_{\mathrm{p}}}}{(1+\frac{\mathrm{\Theta }}{e_{\mathrm{p}}{}^2})}^{1/2},$$(30) $$\propto \{\begin{array}{ll}{e}_{\mathrm{p}}^{-1/2}{\left(\frac{R_{\mathrm{p}}}{r}\right)}^{1/2}& \text{ for }{e}_{\mathrm{p}}\ll \sqrt{\mathrm{\Theta }},\\ e_{\mathrm{p}}{}^{1/2}\left(\frac{R_{\mathrm{p}}}{r}\right){\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}}\right)}^{-1/2}& \text{ for }{e}_{\mathrm{p}}\gg \sqrt{\mathrm{\Theta }},\end{array}$$(31)

where Θ is the Safronov number and given as $$\mathrm{\Theta }=v_{\mathrm{e}\mathrm{s}\mathrm{c}}^2/2v_{\mathrm{K}}^2.$$(32)

The Safronov number is usually used as an indicator for describing the fate of planetary objects in an unstable system. An encounter grazing the planetary surface puts a velocity kick of ~v_esc to a planetesimal, which leads to $\delta e\sim \sqrt{2\mathrm{\Theta }}$. If the Safronov number is larger than one, the gravitational scattering by the planet is likely to lead to ejection. Therefore, it is predicted that accretion is dominant in the region of Θ < 1 while ejection or falling into the star is dominant with Θ > 1. Our model suggests that the accretion probability is affected by not only the Safronov number but also the planet’s mass, radius, eccentricity, and planetesimals’ initial inclination.

3 Accretion probability to an eccentric planet: Numerical simulations

In the previous section, we derived an analytical formula for the accretion probability of planetesimals. We found that the accretion probablity depends on not only M_p, R_p, e_p as found in previous studies, but also i₀. To assess the obtained accretion probability formula, we investigate the orbital evolution of planetesimals using N-body simulations.

3.1 Numerical Setup

We integrate orbits of planetesimals around a central star with a mass of M_s and with an eccentric protoplanet with a mass of M_p. The planet’s semi-major axis a_p and the eccentricity e_p are set as input parameters. In this case, the perihelion and aphelion of the planet are a_p(1 − e_p) and a_p(1 + e_p). The inclination sin i_p is also an input parameter. The other orbital angles, longitude of ascending node Ω_p, longitude of pericentre ϖ_p, and true anomaly τ_p, are set as (Ω_p, ϖ_p, τ_p) = (0, 0, π) or (0, 1/2π, π). We distribute planetesimals around the eccentric planet and integrate the orbits up to a time t = 10⁶ yr, or until all the particles other than the planet have been eliminated from the simulations. In our simulations, we use the IAS 15 integrator (Everhart 1985; Rein & Spiegel 2015), whose accuracy is sufficiently high to calculate the orbital evolution of eccentric objects. The benchmark results of our simulations are described in Appendix A.

We count the number of particles that collide with the eccentric protoplanet N_col, fall to the central star N_fall, and are ejected from the planetary system N_eject. We judge that a particle collides with the eccentric protoplanet once the relative distance becomes smaller than the physical radius of the planet R_p. A particle whose orbital distance from the central star r becomes smaller than the central star’s radius R_star = 0.0046 au (comparable to the solar radius) is assumed to fall into the central star. Once r > 10³ au is achieved, the particle is assumed to be eliminated from the planetary system.

We assume that the planetesimals are distributed uniformly in logarithmic space in the radial direction. We distribute planetesimals from a_pl,in to a_pl,out, where a_pl,in is set as smaller than a_p(1 − e_p), and a_pl,out is set as larger than a_p(1 + e_p). The eccentricity and inclination are then given by Rayleigh distribution with the eccentricity dispersion ${⟨e_0^2⟩}^{1/2}$ and inclination dispersion ⟨sin² i₀⟩^1/2. The other orbital angles of planetesimals, the longitude of ascending node Ω, the longitude of pericentre ϖ, and true anomaly τ, are chosen randomly.

To reduce the computational cost, the planetesimals are treated as test particles. Since planetesimals are represented by test particles, our results can also be applied to solid disks of smaller particles like debris and dust. The number of test particles N_tp is set so that more than 400 particles are in each bin with a width of 0.1 in logarithmic space in the radial direction. The surface density of the planetesimals goes as ∝ r⁻². Note that the density profile of the solid disk is not necessarily the same as our planetesimal disk. The goal of our numerical simulations is to infer the accretion probability, which does not depend on the solid disk profile as long as the mutual interactions are negligible. In this studuy, we assume that the planetesimal disk has the same number density in each bin. In any case, the accretion probablity we derive is independent of the assumed disk structure. However, the total accreted heavy-element mass would change for different assumed disk profiles and could be investigated in future research. At the end of simulations, some planetesimals remain in the system. However, we stop the simulations at t = 10⁶ yr because planetesimal accretions with the protoplanet mainly occur around t ≲ 10⁵ yr. Therefore, our numerical results are unexpected to change if longer simulation times are considered.

We perform parameter studies where we change M_p, R_p, a_p, and e_p, which are also parameterised in Rodet & Lai (2023). We also perform parameter studies where we vary sin i_p, and ⟨sin² i₀⟩^1/2. The values of the different parameters used in our simulations, as well as the range of parameters, are listed in Table 1.

Fig. 1 Time evolution of planetesimals’ mean eccentricity (cross points) and inclination (square points). The different colors correspond to different assumed masses for the eccentric planet. We also plot Eq. (5) with solid lines and Eq. (14) with dashed lines using Cv = 2. Here, we show the results obtained in the standard cases with ap = 10 au.

3.2 Orbital evolution around an eccentric giant planet

Figure 1 shows the time evolution of the geometric mean eccentricity and inclination of planetesimals for various assumed planetary masses. We present the case with a planet of a_p = 1 au, e_p = 0.99, and sin i_p = 0. Eccentricity excitation occurs every planetary orbit and the mean eccentricity increases linearly with time. The inclination excitation is much slower since the typical timescale of close encounters τ_close,2D is 20–100P_orb,p in these cases. We also plot Eqs. (5) and (14) with C_v = 2. Equations (5) and (14) reproduce the numerical results until the mean eccentricity gets closer to the planetary eccentricity ⟨e²⟩^1/2 ~ e_p. This is because secular resonances become effective when ⟨e²⟩^1/2 ~ e_p. We will discuss the secular resonances’ effect on planetesimals’ orbital evolution in Sec. 5.1.

Fig. 2 The fate of planetesimals around the eccentric planet of ep = 0.99. The fractions of collision (blue), falling (green), and ejection (orange) are plotted as a function of the initial semi-major axis of planetesimals. We take the width of bins as 0.2 in logarithm space. Here, we show the results obtained in the standard cases with ap = 1 au (dashed line) and ap = 10 au (solid lines). The vertical dashed line shows the orbit of Θ = 1.

3.3 Fate of planetesimals

Figure 2 shows the final fate of planetesimals as a function of their initial semi-major axes. The blue, green, and orange lines show the fraction of planetesimals that accrete on the eccentric planet F_acc, fall to the central star F_fall, and are ejected from the planetary system F_eject, respectively. We present the case with a planet of M_p = 1 M_Jup, R_p = 1 R_Jup, e_p = 0.99, and sin i_p = 0. The solid and the dashed lines show the cases with a_p = 10 au and 1 au, respectively. F_acc is larger in the inner orbit, which is consistent with the derived analytical model Eq. (23). For a Jupiter mass planet with Jupiter’s radius, f_a/c = R_cap/b_close is ~ 0.76, 0.082 and 0.013 at r = 0.01 au, 0.1 au and 1 au, respectively. We find that F_acc is larger than f_a/c, which means that planetesimals experience several close encounters with the eccentric planet before they are eliminated from the planetary system.

We also find that F_fall is increasing with the increasing semi-major axis of planetesimals. Beyond a radial distance of a ~ 0.1 au, planetesimals mostly fall into the central star. On the other hand, F_eject increases in the outer disk region due to the weak gravitational bound from the central star. However, the ejection from the planetary system is rare, and F_eject does not exceed 0.2. The vertical dotted line in Fig. 2 shows the location of Θ = 1. Our numerical results show that while a single close encounter at Θ < 1 is insufficient to eject planetesimals, a series of weak encounters accumulate perturbations on the planetesimal’s orbit as shown in Fig. 1, which eventually leads to falling or ejection.

Figure 3 shows the geometric mean of the ejection time τ_eject inferred in the numerical simulations. We find that τ_eject is shorter in the inner orbit because e_fall is smaller in this region. τ_eject is shorter around the more massive planet since planetesimal excitation occurs more quickly. The dotted lines show our analytical formula of τ_eject given by Eq. (25). Around the aphelion of the eccentric planet, τ_stay gets higher than Eq. (25) because some planetesimals stay in mean motion resonances, which delay the eccentricity excitation of the planetesimals. Section 5.1 discusses this effect in more detail.

Fig. 3 Averaged time before planetesimals fall to the central star or are ejected from the planetary system. Different colors correspond to the different masses of the eccentric planet. Square and cross plots show the cases with ap = 1 au and ap = 10 au, respectively. The dotted lines show the analytical equation Eq. (25) with a fitting parameter Cv = 2.

3.4 Fitting the numerical results

Figure 4 shows the dependence of the accretion fraction of planetesimals F_acc (as a function of the initial semi-major axis of planetesimals) on the assumed planetary radius R_p, mass M_p, eccentricity e_p, and the initial inclination of the planetesimal disk sin i₀. Parameters other than those labeled in the right boxes of each panel are set to the standard values we list in Table 1. Note that some parameter sets cover the unrealistic region, e.g., R_p = 4 R_Jup with M_p = 1 M_Jup. Still, we include these cases to assess our analytical model’s validity. The accretion probability increases with increasing R_p and decreases with increasing M_p and e_p, which are consistent with the results by Rodet & Lai (2023). We also find that the accretion probability decreases when the initial semi-major axis of planetesimal a₀ and the initial inclination sin i₀ increase because the geometrical picture of accretion changes from 2D to 3D.

As shown in Fig. 1, the eccentricity excitation becomes slower than the prediction by Eq. (5) as the eccentricity reaches values closer to the planetary eccentricity. This effect increases the probability of close encounters before being eliminated from the planetary system. We introduce the following fitting parameters to Eq. (26) to account for this effect: $$P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }\mathrm{\prime }}=\alpha P_{2\mathrm{D}}+\beta P_{3\mathrm{D}2\mathrm{D}}+\gamma P_{3\mathrm{D}}.$$(33)

We then fit the accretion probablity given by Eq. (23) to the numerical results in logarithmic space. The inferred fitting parameters are listed in Table 2 with the coefficient of determination being 0.951. Equation (23) reproduces the numerical results except in the cases of low planetary eccentricity e_p ≲ 0.3. Our model underestimates the accretion probability of e_p ≲ 0.3. We discuss this point further in Sec. 5.1.

3.5 Dependence on the planetary inclination

We considered the cases where the inclination of the eccentric planet is aligned to the planetesimal disk sin i_p = 0. If the eccentric planet is more inclined than the planetesimal disk sin i_p > ⟨sin i⟩, the distance of the close encounter b_close′ is: $$C_{\mathrm{v}}\frac{\delta v}{v_{\mathrm{K}}}\frac{r\sin i_{\mathrm{p}}}{b_{{\text{close}}^{}}{}^{\mathrm{\prime }}}=\sin i,$$(34) $$b_{\text{close}}{}^{\mathrm{\prime }}=b_{\text{close }}\sqrt{\frac{\sin i_{\mathrm{p}}}{\sin i}}.$$(35)

The inclination evolution accelerates by $\sqrt{\sin i_{\mathrm{p}}/\sin i}$. However, the planetesimal disk reaches sin i_p < ⟨sin i⟩ within τ_eject, and the following orbital evolution would be the same as the cases with sin i_p < ⟨sin i⟩. Figure 5 shows the numerical results with sin i_p = 10⁻³ and 10⁻². For comparison, we also plot the results with sin i_p = 0, but with ⟨sin² i₀⟩^1/2 = 10⁻³ and 10⁻². We find that the collision fraction with different sin i_p is consistent with the cases of corresponding ⟨sin² i₀⟩^1/2 within a factor of ~2. We also change the longitude of the pericenter ϖ_p and find that the collision fraction is mostly independent of ϖ_p. For the cases with sin i_p > ⟨sin i₀⟩ the collision probability is given by $P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }}$ with sin i_p substituted for ⟨sin i₀⟩.

Fig. 4 The fraction of collided planetesimals Facc as a function of the initial semi-major axis of planetesimals. We plot the results of parameter studies for different assumed planetary radius Rp, the planetary mass Mp, the planetary eccentricity ep, and the initial disk thickness ⟨sin2 i0⟩1/2 from top to bottom. The left and right columns show the cases with ap = 1 au and ap = 10 au, respectively. The solid lines show the fitting formula given by Eq. (23).

3.6 Effects of self-gravity

In the simulations discussed above, we treated planetesimals as test particles. We perform additional simulations using five equal-sized embryos instead of the planetesimal disk to investigate the effects of mutual gravitational interactions. We set the mass of each embryo M_emb to 1 M_⊕ and 0.1 M_⊕. The embryos’ mean density is set to 5.5 g/cm³, and the initial semi-major axis of the innermost embryo is an input parameter a_emb,0. The remaining embryos are distributed with the radial orbital separation of 10 mutual Hill radii. The eccentricity and inclination are given by the Rayleigh distribution with ${⟨e_0^2⟩}^{1/2}=2{⟨{\sin }^2{i}_0⟩}^{1/2}={10}^{-3}$, and the orbital angles are set randomly. We perform simulations with 40 different seeds for each a_emb,0 and assume perfect mergers for collisions.

Due to the mutual gravitational interactions, the eccentricity and inclination are excited to ~ h_emb, which is a reduced mutual Hill radius h_emb = (2M_emb/3M_s)^1/3. After the embryos’ inclination reaches ⟨sin i⟩ ~h_emb, it evolves similarly to the planetesimals. Figure 6 shows the time evolution of the geometrical mean of eccentricities and inclinations of embryos. The embryos’ inclination increases with time more rapidly than planetesimals in the early evolution stage t ≲ 10 P_orb,p because the mutual interactions between the embryos dominate the inclination evolution. After that stage, the inclination evolution follows that expected from Eq. (14) with ⟨sin i₀⟩ ~h_emb.

Figure 7 shows the fraction of collisions to the planet F_col,emb as a function of the semi-major axis. F_col,emb = M_col,tot/M_emb,tot where M_col,tot is the total mass of the embryos colliding with the eccentric planet. M_emb,tot = 200 M_emb is the total mass of embryos used for each a_emb,0. We also plot the accretion fraction F_acc obtained in the simulations using the planetesimal disk with ⟨sin² i₀⟩^1/2 = h_emb. We find that the collision fraction F_col,emb is comparable to the accretion fraction F_acc. Therefore, we can conclude that the collision probability of embryos is also given by $P_{\mathrm{c}\mathrm{o}\mathrm{l}}^{\mathrm{\prime }}$ with ⟨sin i₀⟩ = h_emb.

Fig. 5 Same as Fig. 4 but for different assumed planetary inclinations. The blue and green solid lines show the cases with different longitudes of pericenter ϖp. The orange dashed line shows the result for a non-inclined planet but with an inclination-excited planetesimal disk.

Fig. 6 Same as Fig. 1, but using embryos instead of a planetesimal disk. The blue and green plots show the geometrical mean of eccentricities and inclinations of embryos. The blue solid line shows Eq. (5) with ⟨e0⟩ = 1 × 10−3. The green solid and dashed lines show Eq. (14) with ⟨sin i0⟩ = 5 × 10−4 and ⟨sin i0⟩ = hemb, respectively. Here, we show the results of 1 M⊕ embryos around 0.4 au.

Fig. 7 Same as Fig. 4 but with the blue plots showing the results for embryos instead of the planetesimal disk. We use five embryos for each simulation and sum up the collision fraction. The upper panel shows the case with 1 M⊕ embryos, and the lower panel shows the case with 0.1 M⊕ embryos. For comparison, we also plot the results for the planetesimal disk with green lines. The initial inclination of the planetesimal disk is equal to the reduced hill radius of embryos.

4 The effect of the planetary eccentricity excitation

In the previous section, we inferred the accretion probability of a giant planet with a fixed orbit without considering the planet’s orbital evolution. Here, we investigate how the planet’s eccentricity evolution affects the inferred accretion probability.

Fig. 8 Snapshots of the planetesimal disk when the planet’s eccentricity increases with time. The left and right columns show the cases with τecc = 103 Porb,p and 103 Porb,p, respectively. The upper and lower panels show the snapshots when the planetary eccentricity is ep = 0.3 and 0.81, respectively. The vertical dotted lines are the pericenter and apocenter of the planet. The dashed lines are the analytical models of the planetesimal eccentricity given by Eq. (40). When τecc = 105 Porb,p, the planetesimal eccentricity is smaller than that given by Eq. (40) because the resonances suppress the linear increase of the eccentricity.

4.1 Eccentricity excitation

We consider a case where multiple planets orbit the exterior of the planetesimal disk. The planets initiate orbital instability and increase eccentricities, and then the pericenter of the innermost planet starts to cross the planetesimals’ orbit. We assume that the gravitational interaction with the innermost planet is the dominant perturber onto the planetesimal disk and neglects the gravity from the outer planets since the gravitational force from the outer planets is much smaller than that from the innermost planet.

We start the N-body simulations with e_p = 0 and increase the planetary eccentricity by 1/τ_ecc every orbit where τ_ecc is the excitation timescale. There are various mechanisms for eccentricity excitation, and the excitation timescale differs in each mechanism. Typical timescales of the Kozai-Lidov resonance and secular interactions could be 10⁵ P_orb,p or longer (e.g. Dawson & Johnson 2018). On the other hand, planet-planet scattering can increase the eccentricity more quickly. The eccentricity excitation kicked by a single encounter is given by: $$\mathrm{\Delta }e\lesssim \frac{v_{\mathrm{e}\mathrm{s}\mathrm{c}}}{v_{\mathrm{K}}},$$(36) $$=0.32{\left(\frac{M_{\mathrm{p}}}{0.5M_{\text{Jup }}}\right)}^{1/2}{\left(\frac{R_{\mathrm{p}}}{2R_{\text{Jup }}}\right)}^{-1/2}{\left(\frac{a_{\mathrm{p}}}{0.1\mathrm{a}\mathrm{u}}\right)}^{1/2}{\left(\frac{M_{\mathrm{s}}}{M_{\oplus }}\right)}^{-1/2}.$$(37)

N-body simulations have shown that planet-planet scattering can excite the planetary eccentricity as quickly as 10³ P_orb,p at most (Carrera et al. 2019; Pu & Lai 2021; Garzón et al. 2022). We consider the following values of τ_ecc = 10³, 10⁴, 10⁵ P_orb,p to account for the various excitation mechanisms.

While the secular interactions keep the semi-major axis constant (Teyssandier et al. 2019), the planet-planet scattering changes other orbital parameters, such as the planetary semimajor axis and inclination (Carrera et al. 2019; Pu & Lai 2021; Garzón et al. 2022). However, here, we increase the eccentricity by keeping the semi-major axis and the inclination constant to investigate how the time evolution of the planetary pericenter affects the accretion probability.

4.2 Case with τ_ecc = 10³ P_orb,p

First, we show the result with a planet whose eccentricity increases as rapidly as τ_ecc = 10³ P_orb,p. The left panels in Fig. 8 show snapshots of the planetesimal disk. The upper and lower panels show when the planet’s eccentricity becomes 0.30 and 0.81, respectively. The vertical dotted lines are the pericenter of the eccentric planet. We find that planetesimals are excited to the eccentric orbit before the planet’s pericenter reaches the planetesimals’ orbit. This is because of the indirect forces of the eccentric planet. Kobayashi & Ida (2001) showed that the eccentric object can perturb the inner planetesimal disk without orbital crossing. The eccentricity perturbation put by a passing object is given as (Kobayashi & Ida 2001): $$\mathrm{\Delta }{e}_{\mathrm{p}\mathrm{a}\mathrm{s}\mathrm{s}}\simeq A{\left(\frac{a}{q_{\mathrm{p}}}\right)}^{5/2},$$(38) $$\begin{array}{rl}A=& \frac{15\pi }{32\sqrt{2}}\frac{M_{\mathrm{p}}/M_{\mathrm{s}}}{\sqrt{1+M_{\mathrm{p}}/M_{\mathrm{s}}}}[{\cos }^2{\omega }_{\mathrm{p}}{(1-\frac{5}{4}{\sin }^2{i}_{\mathrm{p}})}^2\\ & {+{\sin }^2{\omega }_{\mathrm{p}}{\cos }^2{i}_{\mathrm{p}}{(1-\frac{15}{4}{\sin }^2{i}_{\mathrm{p}})}^2]}^{1/2}\end{array},$$(39)

where q_p is the pericenter of the planet. We adopt Eq. (38) to the eccentric planet and integrate it along the planetary eccentricity evolution. We obtain the eccentricity of the planetesimal disk, which is given by: $$\begin{array}{rl}⟨e⟩& \simeq {\int }_0^t\mathrm{\Delta }{e}_{\text{pass }}\mathrm{d}{t}^{\mathrm{\prime }}={\int }_0^{e_{\mathrm{p}}}\mathrm{\Delta }{e}_{\text{pass }}\frac{{\tau }_{\text{ecc }}}{P_{\text{orb,p }}}\mathrm{d}{e}_{\mathrm{p}}^{\mathrm{\prime }},\\ & =⟨e_0⟩+\frac{2}{3}A{\tau }_{\text{ecc }}{\left(\frac{a}{a_{\mathrm{p}}}\right)}^{5/2}[{(1-e_{\mathrm{p}})}^{-3/2}-1].\end{array}$$(40)

We plot Eq. (40) with the dashed lines in Fig. 8. We find that Eq. (40) reproduces well the eccentricity evolution of the planetesimals.

The condition where the planet crosses the orbit of planetesimals is given by: $$a(1+e)=a_{\mathrm{p}}(1-e_{\mathrm{p}}).$$(41)

By solving Eqs. (40) and (41) numerically, we can obtain the eccentricity of planetesimals and the planet when the orbits are crossed e_cross, and e_cross,p, separately. Figure 9 shows e_cross as a function of r_cross = a(1 + e_cross). The plots with solid lines are numerical results inferred from the N-body simulations. We obtain e_cross by calculating the mean eccentricity of planetesimals whose orbits are just before crossing to the planet’s orbit: 0.8 < a(1 + e)/a_p(1 − e_p) < 0.9. The dashed lines show the analytical prediction of e_cross obtained from Eqs. (40) and (41). When τ_ecc = 10³ P_orb,p (blue lines), the analytical model reproduces the numerical result well, and e_cross is smaller than 0.1. Therefore, planetesimals have many chances to have close encounters after the orbits are crossed.

Figure 10 shows the accretion fraction of planetesimals F_acc vs. the initial semi-major axis of planetesimals. F_acc exceeds 10% in the inner disk a₀ < 1 au and around the initial orbit of the planet a₀ ~ 10 au. When the planet’s eccentricity is small, the gravitational focusing works effectively, and F_acc can be high. As the planet is excited to the eccentric orbit, its pericenter reaches the inner orbit. The planet accretes planetesimals there, but the accretion probability decreases because of the weaker gravitational focusing effect. This effect leads to a peak on F_acc around a₀ ~ a_p. On the other hand, F_acc gets higher in the further inner disk, too, because R_p/b_close is larger there. The accretion probability can be estimated with Eq. (23) by substituting a₀ = r_cross and e_p = e_p,cross. The green dashed line in Fig. 10 shows Eq. (23). Equation (23) reproduces the numerical results well when τ_ecc = 10³ P_orb,p.

Fig. 9 rcross = ap(1 − ecross,p) vs. ecross obtained by the N-body simulations. ecross and ecross,p are calculated by taking the average of eccentricities of planetesimals whose orbits are just before crossing to the planet’s orbit. The dashed lines show ecross obtained from Eqs. (40) and (41). The thick dashed lines mean Eq. (44) is achieved. When Eq. (44) is achieved, ecross becomes smaller than the analytical models predict. As indicated in the legend, the different colors correspond to the different eccentricity excitation timescales τecc.

Fig. 10 Same as Fig. 4 but with a planet whose eccentricity increases with a given timescale τecc. The blue solid lines show the numerical results and the green dashed lines show Eq. (45).

4.3 Cases with larger τ_ecc

As the eccentricity excitation of the planet gets slower, e_cross becomes larger since the perturbations from the eccentric planet accumulate more before the orbits are crossed. The analytical model of e_cross given by Eqs. (40) and (41) predicts that e_cross becomes as large as 1 when τ_ecc = 10⁵ P_orb,p (the orange dashed line in Fig. 9). This means that planetesimals could be ejected from the planetary system before the orbits of the planetesimals and the planet are crossed, which should decrease the accretion rate. However, numerical results show that e_cross is much smaller than the analytical model predicts (the orange plots in Fig. 9), and the planetesimals can accrete to the eccentric planet as shown in the lower panel of Fig. 10. This is because of the effect of secular resonances. Kobayashi & Ida (2001) found that resonance effects appear around r/q_p ≳ 0.2. This resonance effects suppress the eccentricity excitation and e_cross becomes smaller than the analytical model predicts. Due to the effective secular resonances, the planetesimal disk can accrete onto the eccentric planet even if the eccentricity excitation of the planet is as slow as τ_ecc = 10⁵ P_orb,p.

On the other hand, the accretion fraction is reduced around 1 au and opens a gap in F_acc when τ_ecc = 10⁴ P_orb,p. In the region inner than 1 au, e_cross is small, and the planetesimals can interact with the eccentric planet after the orbits are crossed. However, e_cross is increasing with a₀, and once e_cross exceeds ~ 0.3, F_acc drops to 0. We define a critical value e_cross,crit as the upper limit for planetesimals to accrete onto the eccentric planet. Planetesimal accretion occurs in the region where $$e_{\text{cross }}<e_{\text{cross,crit }}.$$(42)

If the secular resonances become effective, the planetesimals can accrete onto the eccentric planet even if the analytical model given by Eqs. (40) and (41) does not achieve Eq. (42). For the resonances to be effective, the pericenter of the planet needs to move slower than ~ 1/τ_res, where τ_res is the typical timescale of the secular resonances. By analyzing the eccentricity evolution of each planetesimal, we find that the typical timescale of the resonance τ_res is ~ 10³ P_orb,p (see details in Appendix B). The region where the resonances become effective is $$\left|\frac{a_0}{{\dot{q}}_{\mathrm{p}}}\right|>{\tau }_{\text{res }}.$$(43) $$\frac{{\tau }_{\text{res }}}{{\tau }_{\text{ecc }}}{a}_{\mathrm{p}}<a_0.$$(44)

Neither of Eqs. (42) and (44) are not achieved around 1 au when τ_ecc = 10⁴ P_orb,p. Therefore, F_acc has a gap around 1 au. Because of the same reason, F_acc around 0.1 au decreases when τ_ecc = 10⁵ P_orb,p.

Including the above effects, we finally model the accretion probability of planetesimals by $$P_{\mathrm{a}\mathrm{c}\mathrm{c}}=\{\begin{array}{ll}{P}_{\mathrm{a}\mathrm{c}\mathrm{c}}^{\mathrm{\prime }}& \text{ for Eq. (42) or Eq. (44), }\\ 0& \text{ other, }\end{array}$$(45)

We set τ_res = 2 × 10³ P_orb and e_cross,crit = 0.3 from the numerical results. r_cross, e_cross and e_p,cross are obtained with Eqs. (44) and (45). Other parameters such as M_p, R_p, a_p, sin i_p and ⟨sin² i₀⟩^1/2 are input parameters. The green dashed lines in Fig. 10 show the accretion probability given by eq. (45).

Fig. 11 Same as Fig. 1, but for the results of parameter study of the planetary eccentricity ep.

5 Discussion

5.1 The limitations of our orbital evolution model

In Sec. 2, we model the orbital evolution of planetesimals around an eccentric planet. A few effects have been neglected for simplicity, although some of them could be important under certain conditions. Figure 11 shows the time evolution of the mean eccentricities and the mean inclinations of planetesimals for different planetary eccentricity cases. The figure clearly shows that the excitation rates of ⟨e⟩ (cross dots) and ⟨sin i⟩ (square dots) become slower than predicted by our model (solid and dashed lines) once ⟨e⟩ exceeds the planetary eccentricity. Also, our model is invalid when e_p = 0.1. Here, we discuss our model assumptions and the limitations of our orbital evolution model.

5.1.1 A-particle-in-a-box approximation

In a-particle-in-a-box model, we assume that the perturbation is negligibly small so that the trajectory of the planetesimal is unchanged. In the case of a close encounter where the trajectory of the perturbed particle is deflected during the passage of the eccentric planet, a-particle-in-a-box approximation is not valid. By substituting ~ v_rel to δv, we estimate the minimum relative distance available in a-particle-in-a-box approximation b_min as: $$b_{min}\sim \frac{2C_{\mathrm{v}}}{e_{\mathrm{p}}{}^2}\frac{M_{\mathrm{p}}}{M_{\mathrm{s}}}a.$$(46)

Our model is appropriate when: $$b_{\text{close }}>b_{min},$$(47) $$e_{\mathrm{p}}>0.156{\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{J}\mathrm{u}\mathrm{p}}}\right)}^{1/3}{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^{1/3}$$(48)

with C_v = 2. If the planetary eccentricity is lower than Eq. (48), the velocity kick of close encounter becomes $\sim \sqrt{\mathcal{G}{M}_{\mathrm{p}}/b}$ and b_close is redefined as: $$b_{{\text{close }}^{\mathrm{\prime }}}=min(b_{\text{close }},C_{\mathrm{v}}^{2/3}{\left(\frac{2M_{\mathrm{p}}}{M_{\mathrm{s}}}\right)}^{1/3}).$$(49)

Since b_close′ is smaller than b_close, the close encounter timescale gets longer compared to the one we assumed. Equation (14) overestimates the time evolution of the inclination of planetesimals. However, this effect is minor because the difference between b_close and b_close′ is only a few % when e_p = 0.1.

5.1.2 High eccentricity effects on the eccentricity and inclination kick

To transform a given velocity kick to perturbations on eccentricity and inclination, we use the approximated relations of δe ~ δv/δv_κ and δ sin i ~ δv_⊥/δv_κ, where we neglect the eccentricity terms. The exact formulae are given by: $$\delta e=\sqrt{1-e^2}\{\sin f\frac{\delta v_{\mathrm{R}}}{v_{\mathrm{K}}}+(\cos f+\cos E)\frac{\delta v_{\mathrm{T}}}{v_{\mathrm{K}}}\},$$(50) $$\delta \sin i=\sqrt{1-e^2}\frac{\cos (\omega +f)}{1+e\cos f}\frac{\delta v_{\mathrm{Z}}}{v_{\mathrm{K}}},$$(51)

where E is the eccentric anomaly and ω is the argument of perihelion. The term of $\sqrt{1-e^2}$ makes the eccentricity and inclination evolution slower as e increases.

We also assume that orbits between the planetesimals and the planet are crossed, and velocity kicks are put on the planetesimals’ orbit independent of f in every planetary orbit. However, this geometrical picture does not apply when ⟨e⟩ > e_p. Figure 12 shows the histogram of the differences between the pericenters of the planetesimals ϖ and the planet ϖ_p. Each column shows the cases of e_p = 0.1, 0.3, and 0.5. We plot the mean eccentricity of planetesimals divided by the planetary eccentricity in each panel. We find that ϖ − ϖ_p are distributed uniformly or around ~ π/2 if ⟨e⟩ < e_p. In this case, most of the planetesimals cross the planet’s orbit. However, secular resonances align the pericenters between the planetesimals and the planet when ⟨e⟩ > e_p. The alignment of the pericenters detaches the planetesimal’s orbits from the planetary orbit. Since the velocity kicks are weaker after the detachment, the orbital excitation gets slower once the planetesimals’ eccentricities exceed the planetary eccentricity.

5.1.3 Co-rotation resonance

When e_p = 0.1, the mean eccentricities and inclinations of planetesimals are smaller than the other cases even when ⟨e⟩ < e_p. This is because many planetesimals are trapped in the co-rotation resonance (see snapshots in Appendix C). The eccentricities and inclinations of planetesimals in the co-rotation resonances (a ~ a_p) remain small. A small number of planetesimals are trapped in the co-rotation resonances when e_p = 0.3, but no planetesimals are trapped around the more eccentric planets e_p ≥ 0.5. Due to the planetesimals trapped in the co-rotation resonances, the mean eccentricities and inclinations are smaller when e_p = 0.1 than in the other cases. The trapped planetesimals gradually escape from the resonance. By the end of the simulations, most planetesimals escape from the resonance and follow the same orbital evolution as the others.

Fig. 12 Distribution of ϖ − ϖp around an eccentric planet. We plot the cases of ep = 0.1, 0.3, and 0.5 in the left, center, and right columns, respectively. The upper, middle, and lower panels show different times. The mean eccentricities of planetesimals divided by the planetary eccentricity are indicated in the upper left corner of each panel. Once ⟨e⟩ exceeds ep, the planetesimals’ orbits start to align to the planetary orbit.

5.1.4 Effects on our model

For the above reasons, our orbital evolution model is valid only when ⟨e⟩ ≲ e_p. After reaching ⟨e⟩ ~ e_p, orbital excitation gets slower, which makes the mean ejection time of planetesimals τ_eject longer than Eq. (25). Therefore, the accretion probability is larger than that our model expects when e_p ≲ 0.3. Our accretion model underestimates the accretion probability by a factor of 3 (2) when e_p = 0.1 (0.3). This inconsistency gets smaller as the planetary eccentricity gets larger.

5.1.5 Perturbations from other planets

In Sec. 4 we neglect the gravitational perturbation from planets other than the innermost planet. The strength of the eccentricity perturbation Δe_pass depends on ∝ M_p/q_p^5/2 as shown in Eqs. (38) and (39). If two Jupiter-mass planets are separated by 10 Hill radius, the inner planet has ~ 5 times stronger perturbations. If the outer planet is heavier than the inner planet, the perturbation from the outer planet could be comparable to the one from the inner planet. However, in that case, the outer planet is expected to remain in a low eccentricity orbit because of the mass difference. Therefore, only the inner planet enters the highly eccentric orbit and dominates the scattering of the planetesimal disk since Δe_pass strongly depends on q_p. The gravitational perturbation from the outer planets might accelerate the excitation of the planetesimal disk; however, this is expected to be a minor effect.

On the other hand, we found that resonances with the innermost planet play a role in planetesimal accretion. The gravitational perturbation from the other planet might alter the resonance configuration and the accretion probability. If the eccentricity excitation of the planet is as slow as τ_ecc = 10⁵ P_orb,p, almost all planetesimals are affected by the secular resonances as shown in Sec. 4. If the outer planets completely break secular resonances, the planetesimals would be ejected before the innermost planet crosses their orbit, preventing planetesimal accretion. On the other hand, when τ_ecc = 10³ P_orb,p, the accretion probability would not change since the secular resonances do not affect the accretion probablity in this case.

Also, we do not consider the orbital swap between planets. Mustill et al. (2018) ran N-body simulations with three planets and a planetesimal disk, where the planet initially located in the innermost orbit could be scattered into the outermost orbit. Since their main focus was to investigate the accretion onto the white dwarf, it is unclear how many planetesimals accreted onto each planet. Future studies should assess these effects by running the N-body simulations with several planets.

Fig. 13 Semi-major axis vs. eccentricity of observed exoplanets. The data is obtained from Data from Data & Analysis Center for Exoplanets (DACE) in April 2024. Here, we include exoplanets heavier than Neptune, with a measurement uncertainty of less than 20 % in eccentricity. The grey dashed lines show the boundary of the regions where planets could be tidally disrupted (qp < 0.01 au) and where tidal forces with the central star are sufficient to cause inward planetary migration (qp < 0. 05 au).

5.2 Enrichment of eccentric giant planets

Giant planets in exoplanetary systems have a wide range of eccentricities (see Fig. 13). Data from Data & Analysis Center for Exoplanets (DACE) show that Jupiter-mass planets (with M_p > 0.04 M_Jup) have an average eccentricity ≃ 0.30 with a standard deviation of ≃ 0.22, and a median eccentricity of ≃ 0.26. The eccentricity distribution of eccentric Jupiters could be generated by the orbital instability between giant planets after disk dissipation (Carrera et al. 2019; Anderson et al. 2019; Garzón et al. 2022). If solid disks exist, these eccentric Jupiters can accrete solids during their dynamical history. Below, we estimate the solid accretion during the eccentricity excitation process.

Fig. 14 Estimated mass of accreted solid materials Macc normalised with πap2∑sol for various planets. The left panel shows Macc/πap2∑sol as a function of the planetary eccentricty. The different color corresponds to the different semi-major axis of the planet. The planetary mass is set to Jupiter’s mass. The right panel shows Macc/πap2∑sol as a function of the planetary mass. Here, the different color corresponds to the different eccentricity of the planet. The semi-major axis is set to 1 au. In both panels, the solid and dashed lines show the cases of sin ip = 10−2 and 10−1, respectively. We set τecc to 103 Porb,p assuming the planet-planet scattering.

5.2.1 Accretion of remnant solid disk onto cold Jupiters

We consider a planetary system with multiple giant planets that initiates orbital instability and estimate the mass of accreted solids onto the innermost planet M_acc. By the time the planet has reached its current eccentricity e_p, the mass of accreted solids is given by: $$M_{\mathrm{a}\mathrm{c}\mathrm{c}}={\int }_{r_{\text{in }}}^{r_{\text{out }}}2\pi a_0{\mathrm{\Sigma }}_{\text{sol }}{P}_{\mathrm{a}\mathrm{c}\mathrm{c}}\left(a_0\right)\mathrm{d}{a}_0,$$(52)

where ∑_sol is the surface density of solid material. For simplicity, we set ∑_sol as a constant. We set r_in to be the planet’s pericenter a_p(1 − e_p) and r_out to the inner edge of the planetary feeding zone $a_{\mathrm{p}}(1-2\sqrt{3}{\left(M_{\mathrm{p}}/3M_{\mathrm{s}}\right)}^{1/3})$ because the planet scatters planetesimals inside the feeding zone before the disk dissipates (e.g. Shibata & Ikoma 2019; Shibata et al. 2022b). We fix the planet’s semi-major axis during the integration of Eq. (52). However, mutual scattering between the giant planets changes the planetary semi-major axis. If the planet is scattered inward, which is a natural outcome of the planet-planet scattering, the planet sweeps planetesimals orbiting interior to r_in and accretes more solids. Also, the planetesimals orbiting beyond the planet would be scattered by the other giant planets. Some scattered planetesimals could accrete to the innermost planet as shown in Matter et al. (2009), but we neglect this contribution. A comparison between our inferred M_acc and their results is presented below.

Figure 14 shows the inferred mass of accreted solids. We normalise M_acc with πa_p²∑_sol. Here, we adopt the mass-radius (M-R) relation from Otegi et al. (2020) and use Jupiter’s radius as an upper limit of the planetary radius. The planetary radius is set as: $$R_{\mathrm{p}}=min(0.70R_{\oplus }{\left(\frac{M_{\mathrm{p}}}{M_{\oplus }}\right)}^{0.63},R_{\mathrm{J}\mathrm{u}\mathrm{p}}).$$(53)

The solid and dashed lines show the cases with a planetary inclination of sin i_p = 10⁻² and 10⁻¹, respectively. We set the eccentricity evolution timescale to τ_ecc = 10³ P_orb,p assuming planet-planet scattering. The left panel shows M_acc as a function of the planetary eccentricity. When the planetary eccentricity is smaller than ~ 0.2, the eccentric planet cannot accrete planetesimals because the planet’s pericenter is within the planetary feeding zone. When e_p ≳ 0.2, M_acc increases with planetary eccentricity. Solid accretion is more effective for small radial distances and lower inclination. The right panel shows M_acc as a function of planetary mass. The mass of accreted solid rapidly decreases when M_p > M_Jup. This is because e_cross is larger for more massive planets. When M_p > 2 M_Jup, e_cross exceeds 0.3 and P_acc decreases. This means the solids are ejected from the planetary system before the planetary orbit crosses the solid disk.

While it is still unclear how many solids remain around giant planets after disk dissipation, N-body simulations show that many planetesimals (more than 10 M_⊕) could remain around a giant planet (e.g. Zhou & Lin 2007; Shibata & Ikoma 2019; Shibata et al. 2022b). Also, the giant planet itself can generate a new planetesimal disk via a gap opening in the protoplanetary disk (Shibaike & Alibert 2020, 2023; Eriksson et al. 2022). After disk dissipation, the remnant planetesimal disk excites itself through the mutual gravitational interactions and initiates collisional cascades (Wyatt et al. 2007; Löhne et al. 2008). Heng & Tremaine (2010) found that the maximum surface density of solids around an old star of 3 Gyrs exceeds ∑_sol = 1 g/cm² beyond 10 au. Also, the observed dust emission suggests the existence of a massive solid disk at distances of 10–100 au (Krivov & Wyatt 2020). If a solid disk of ∑_sol = 1 g/cm² remains at 10 au, πa_p²∑_sol is 12 M_⊕ and M_acc is 0.01–0.1 M_⊕. This heavy-element mass can’t increase the planetary bulk metallicity but could increase the atmospheric metallicity (Howard et al. 2023; Müller & Helled 2024). Therefore, we predict that the atmosphere of eccentric cold Jupiters will be enriched with heavy elements, but the actual enrichment depends on the planetary eccentricity and mass. Solid accretion does not occur for massive cold Jupiter of M_p > 7 M_Jup.

In this study, we use a single eccentric planet, although we assume that the planet is excited by other giant planets. The existence of the other giant planets would increase the accretion of solids because they scatter planetesimals from the outer orbit inwards. Matter et al. (2009) showed that the accretion of comets scattered by Uranus and Neptune brings ~ 0.1 M_⊕ of heavy elements to Jupiter. Since Jupiter has a low eccentricity <0.1 in their model, comet accretion onto an eccentric cold Jupiter would be smaller than 0.1 M_⊕. Nevertheless, the scattering of the remnant solid disk by multiple giant planets would increase the mass of the accreted solid after disk dissipation. This process should be addressed in future studies.

The available solid surface density would be smaller than ∑_sol = 10⁻³ g/cm² at 1 au (Heng & Tremaine 2010) because the collisional evolution timescale strongly depends on the radial distance from the central star. Despite the higher accretion efficiency, the mass of accreted solids would be smaller than 10⁻⁵ M_⊕ in the inner orbit. Therefore, other solid accretion processes like comet accretion (Seligman et al. 2022) would be required to enrich the planet further.

Fig. 15 Estimated mean collision number of planets onto an eccentric planet. We show the mean collision number as a function of the planetary mass. The left and right panel show the cases with τecc = 103 Porb,p and τecc = 105 Porb,p, respectively. The blue, green, and orange lines correspond to ep = 0.3, 0.6, and 0.9, respectively. The solid and dashed lines show the cases of sin ip = 10−2 and 10−1, respectively.

5.2.2 Collisions of planets

The statistical analysis shows that many planetary systems with a cold Jupiter have inner terrestrial planets (Zhu & Wu 2018; Bryan et al. 2019; Rosenthal et al. 2022). Collisions of such inner terrestrial planets could increase the planetary metallicity of eccentric giant planets. Mustill et al. (2015) used N-body simulations to show that some inner planets collide with the outer giant planets during orbital instability. In this section, we estimate the collision probability between eccentric Jupiters and inner terrestrial planets. Kunimoto & Matthews (2020) derive a fitting function for the occurrence rate of planets with R_p < 4 R_⊕ in the inner orbit, which is given by: $$\begin{array}{rl}{F}_{\text{occ }}& =\frac{\mathrm{d}f}{\mathrm{d}\log P_{\text{orb,p}}}\\ & =C_{\text{KM}20}{P}_{\text{orb,p }}{}^{{\beta }_{\mathrm{K}\mathrm{M}20}}(1-{\gamma }_{\mathrm{K}\mathrm{M}20}\exp (-\frac{P_{\text{orb},\mathrm{p}}}{P_{\text{orb},0}})).\end{array}$$(54)

where C_km20, β_km20, γ_KM20, and P_orb,0 are fitting parameters (see Kunimoto & Matthews 2020 for each value).

Note that Eq. (54) gives the occurrence rate for all planetary systems, but a planetary system with a cold Jupiter would have a larger occurrence rate for the inner planets (Zhu & Wu 2018; Bryan et al. 2019; Rosenthal et al. 2022). The mean collision number of inner planets N_col can be obtained by integrating ${\int }_{r_{\text{in}}}^{r_{\text{out}}}{F}_{\mathrm{o}\mathrm{c}\mathrm{c}}{P}_{\mathrm{a}\mathrm{c}\mathrm{c}}\mathrm{d}{a}_0.$. Figure 15 shows the estimated mean collision number for the inner planets. Here, we consider the case of a_p = 1 au and use the same M-R relation as Sect. 5.2.1. The results are presented for different assumed eccentricity excitation timescales of τ_ecc = 10³ P_orb,p (left panel) and τ_ecc = 10⁵ P_orb,p (right panel), and different planetary inclination sin i_p = 0.01 (solid lines) and sin i_p = 0.1 (dashed lines). We find that N_colincreases with increasing e_p and τ_ecc, and for decreasing sin i_p.

With the averaged eccentricities of detected Jupiter-mass planets ≃ 0.3, N_col is estimated to be smaller than 0.04. Therefore, collisions of inner terrestrial planets onto the eccentric Jupiter would be minor events. On the other hand, hot Jupiters might experience collisions with a higher probability. One of the possible formation pathways of hot Jupiters is high-eccentric migration. A giant planet perturbed to a high-eccentric orbit would dissipate its orbital energy and migrate inward through interaction with the central star. A hot Jupiter formed by high-eccentric migration should have been excited to the highly eccentric orbit of e_p ≃ 0.9 if the initial semi-major axis was ~ 1 au. Using N-body simulations Garzón et al. (2022) found that most highly eccentric planets are formed via slow eccentricity excitation mechanisms like secular interactions. Assuming that hot Jupiters are formed via high-eccentric migration with e_p ≃ 0.9 and τ_ecc > 10⁵ P_orb,p, we predict that roughly 10 % of the hot Jupiters collided with inner terrestrial planets.

The dataset used in Kunimoto & Matthews (2020) shows a peak at R_p ~ 2.5 R_⊕ in the occurrence rate. This corresponds to a mass of M_p~10 M_⊕ from the M-R relation of Otegi et al. (2020). A collision of a 10 M_⊕ planet could increase the planetary core mass and bulk metallicity. However, such collision probability is only of the order of ~ 10% and therefore cannot explain the inferred high metallicity of hot Jupiters (e.g. Guillot et al. 2006). Other heavy element accretion processes such as planetesimal accretion (Hands & Helled 2021; Knierim et al. 2022) and accretion of enriched gas (Schneider & Bitsch 2021a; Schneider & Bitsch 2021b) before disk dissipation are required to explain the enrichment of hot Jupiters.

Fig. 16 Fraction of planetesimals that fall onto the central star as a function of the planetesimals’ initial semi-major axis. The blue, green, and orange lines show the cases of τecc = 103 Porb,p, 104 Porb,p, and 105 Porb,p, respectively.

5.3 Pollution of the central star

Our simulations show that most of the remaining solids fell onto the central star after disk dissipation. Figure 16 presents the fraction of planetesimals that fall onto the star F_fall during the planetary eccentricity evolution (using the results obtained in Sec. 4). Rodet & Lai (2023) find that F_fall is larger for a planet with higher eccentricity. In our simulations, the planet is found to scatter planetesimals between the pericenter and apocenter, and this region expands with the planetary eccentricity. Therefore, F_fall is small around the planet’s semi-major axis and increases as the orbital separation from the planet increases. Rodet & Lai (2023) fixed the planet’s eccentricity, the same as the simulations in Sec. 3. Comparing Fig. 2 and Fig. 16, we find that F_fall decreases when we include the effect of eccentricity evolution and with increasing eccentricity excitation timescale τ_ecc. Our results imply that the pollution of the central star depends on the orbital evolution history of their giant planets. We hope to investigate this in further detail in future research.

The observed enhancement of Li in Sun-like stars implies that a significant amount of solids polluted the stellar atmosphere during their lifetime (Spina 2024, see references there). While the engulfment of planets is suggested as a source of stellar atmospheric pollution (Spina et al. 2021), accretion of scattered planetesimals onto the central star is another possible mechanism for enriching stellar atmospheres with refractory materials. Planetesimal accretion onto a central star has been investigated by various authors using N-body simulations (Beust & Morbidelli 2000; Frewen & Hansen 2014; Mustill et al. 2018; Smallwood et al. 2021; Veras et al. 2022; O’Connor et al. 2022; Rodet & Lai 2023, see also Veras 2021 for more references). Mustill et al. (2018) investigate the accretion of the remnant planetesimal disk onto white dwarfs and find that 13.3% of planetesimals distributed between 5 au and 8.5 au entered the Roche limit of the white dwarf and got accreted. Our simulations suggest a higher accretion probability of 64% for τ_ecc = 10³ P_orb,p and 38% for τ_ecc = 10⁵ P_orb,p in the same orbital region. As shown in Rodet & Lai (2023), the accretion probability is higher for a larger stellar radius and a larger stellar mass. In this study, we use a larger physical radius of the central star than the radius of the Roche limit of the white dwarf and a larger stellar mass. As a result, the accretion probability we infer is larger in comparison to Mustill et al. (2018). This implies that the pollution of the stellar atmosphere is more probable for solar-type stars than white dwarfs. Spina (2024) shows that the accretion of several M_⊕ of solids onto a solar-type star increases the iron fraction in the outer convective cell beyond the detection limit and that this enrichment can be observed because the enriched layer persists for ~ 1 Gyrs (Behmard et al. 2022). Our study also suggests that several M_⊕ of planetesimals could be accreted by the star after the onset of orbital instability of giant planets if a solid disk of ∑_sol = 1 g/cm² remains around 10 au. However, since our orbital evolution model is not valid when ⟨e⟩ > e_p (see Sec. 5.1), we can’t derive an analytical prescription for the accretion process onto the central star. As shown in Fig. 12, the orbital evolution of planetesimals, including the effects of secular perturbations, needs to be modeled, and we hope to address this topic in our future research.

6 Summary and conclusions

We investigate the possibility of solid accretion onto an eccentric gas-giant planet after disk dissipation. We focus on the scenario where a giant planet is excited by other giant planets and investigate the accretion probability of solids onto the excited planet. We find that the orbital evolution of solid particles, such as planetesimals and embryos, is regulated by weak encounters with the eccentric planet rather than strong close encounters. We show that even in the region where the Safronov number is smaller than unity, most of the solids fall onto the central star or are ejected from the planetary system.

We also derive an analytical formula for the accretion probability of solid particles. Our analytical prescription can be applied for both planetesimals and embryos by setting the initial inclination appropriately. We also investigate the solid accretion rate of a giant planet whose eccentricity increases with time. In that case, we find that resonance effects play an important role in the orbital evolution of the solid particles before the planetary orbit crosses the solid disk. We derive the condition determining when the solid particles are ejected from the planetary system before the orbits are crossed. We find that solid accretion mainly occurs in the two orbital regions: the inner disk where the planetary capture radius is larger than the disk’s thickness, and at the planet’s vicinity where the planet’s orbit crosses that of planetesimals when gravitational focusing is efficient.

The main conclusions from our study can be summarized as follows:

1. Eccentric Jupiters around 10 au could accrete 0.01–0.1 M_⊕ of heavy elements from the solid disk remaining around the planet. The accretion probability depends on the planetary eccentricity, inclination, and mass. Therefore the metallicity of the planetary atmosphere should relate to these parameters. While the accretion probability is higher in the inner orbit, the accretion of the solid disk is negligible around 1 au because the available solid disk mass is too small in the inner orbit;

2. Using the occurrence rate of terrestrial planets, we estimate how many terrestrial planets collide with an eccentric planet. The number of planets colliding with the eccentric planet increases with the planetary eccentricity; if hot Jupiters formed via high-eccentric migration, ~10% of hot Jupiters were impacted by a planet of ~ 10 M_⊕

3. We find that several Earth masses of solid material scattered by the eccentric planet could accrete onto the central star. Accretion of the scattered solid disk would be a source of heavy elements of the stellar atmosphere as well as the engulfment of planets.

We conclude that solid accretion onto an eccentric planet would be a minor process for planets in the inner orbit ≲ 1 au. Therefore, a different solid accretion process, such as the accretion of small bodies such as comets and asteroids would be required to enrich the planetary atmosphere with heavy elements. This work represents the first step in investigating how the orbits of planetesimals and embryos evolve around a single eccentric planet. However, since some planetary systems can have more than one giant planet, future studies should explore solid accretion accounting for multiple giant planets. Our study clearly shows that in order to understand the observed composition of giant planets, and even stars, a robust understanding of the early evolution of planetary systems is required.

Acknowledgements

We thank the anonymous reviewer for constructive comments. We thank Hiroshi Kobayashi for the fruitful discussion about the dynamics part of this study. S.S. and R.H. acknowledge support from the Swiss National Science Foundation (SNSF) under grant 200020_188460. Numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A Benchmarks of IAS15

We investigate the orbital evolution of planetesimals and embryos around an eccentric planet. The orbital integration of eccentric objects should be carefully treated because the angular velocity changes orders of magnitude in one orbit. We adopt the IAS15 integrator (Everhart 1985; Rein & Spiegel 2015), which uses the adaptive timestep. Figure A.1 shows the time evolution of numerical errors of an eccentric particle. Here, we solve the two-body problems of a solar-mass star and an eccentric particle with a semi-major axis of 10 au and eccentricities of 0.99 (left) and 0.999 (right), respectively. The numerical error is less than 10⁻⁹ even if we integrate them for 10⁵ Kepler orbits.

Fig. A.1 Time evolution in the errors of orbital elements. Blue, green, and orange lines show the error in the semi-major axis, eccentricity, and longitude of the pericentre, respectively.

Appendix B Resonances of the eccentric planet

As shown in Sec. 4.3, resonances of the eccentric planet play important roles in the orbital evolution of planetesimals. Figure B.1 shows the difference in the longitude of pericenters between planetesimals and the planet ϖ − ϖ_p. cos(ϖ − ϖ_p) = 1 means that the pericenter of the planetesimal and the planet is aligned and the gravitational scattering is weaker than in the case where the pericenters are not aligned. Therefore, the eccentricity excitation rate gets slower than the Eq. (40) expects. The alignment of the pericenter around a giant planet within a gaseous disk is shown by Guo & Kokubo (2021, 2022). The alignment also occurs around a giant planet whose eccentricity is increasing.

Fig. B.1 Same as Fig. 8, but shows the difference in the pericenters between planetesimals and the planet instead of planetesimals’ eccentricity. The left and right columns show the cases with τecc = 103 Porb,p and 105 Porb,p, respectively. The upper and lower panels show the snapshots when the planetary eccentricity is ep = 0.1 and 0.81, respectively. The vertical dotted lines are the pericenter and apocenter of the planet.

In order to investigate the typical timescale of the resonances, we analyze the orbital evolution of planetesimals around an eccentric Jupiter-mass planet with a_p = 10 au and e_p = 0.5. Figure B.2 shows the time evolution of cos(ϖ − ϖ_p) and the amplitude spectrum obtained with the Fast Fourier Transform. In the amplitude spectrum, some planetesimals have a strong peak around ~ 10³ P_orb,p. Therefore, the secular perturbations from the eccentric planet align the planetesimal’s pericenter within a timescale of ~ 10³ P_orb,p.

Fig. B.2 Left columns: Time evolution of sin(ϖ − ϖp) for each planetesimal. Right columns: Amplitude spectrum obtained with Fast Fourier Transform. The amplitude spectrum is scaled with N/2, where N is the number of data points. Each planetesimal has a different initial semi-major axis, and it is shown in the upper right corner of the panels in the right column divided by the pericenter of the planet qp = 5 au.

Appendix C Snapshots of planetesimals

We show the snapshots of our simulations as supporting information for this paper. Figure C.1 shows the snapshots from our simulation. We present the case with a planet of M_p = M_Jup, R_p = R_Jup, and an orbit with a_p = 1 au, e_p = 0.99, sin i_p = 0. The eccentricity of the entire planetesimal disk is excited by the first passage of the eccentric planet. As expected in a-particle-in-a-box approximation, the eccentricity excitation rate is independent of planetesimals’ semi-major axis. Eccentricity excitation occurs in every planet’s orbit and uniformly within the planetesimal disk.

Fig. C.1 Snapshots of orbital integration around an eccentric giant planet. The left and right columns show the eccentricity and inclination, respectively, with the semi-major axis of planetesimals. Time progresses from upper to lower panels. The solid lines in the left column show efall. The small boxes on the right side of each panel show the histogram. The total number of planetesimals N is written in the boxes. Here, we show the case with ap = 1 au, ep = 0.99, sin ip = 0, Mp = MJup and Rp = RJup.

Figure C.2 also shows the snapshots from our simulation but with a different eccentricity of the planet e_p = 0.1. In this case, many planetesimals are trapped in the co-orbital resonance around the planet’s semi-major axis a_p = 10au. The eccentricity and inclination excitation is slower in the co-orbital resonance. Therefore, the mean eccentricity and inclination obtained in the N-body simulations are smaller than the analytical models given by eq. (5) and eq. (14) predict.

Fig. C.2 Same as Fig. C.1, but for ep = 0.1. Unlike in Fig. C.1, here, many planetesimals are trapped in the co-rotation resonance of the planet.

References

All Tables

All Figures

	Fig. 1 Time evolution of planetesimals’ mean eccentricity (cross points) and inclination (square points). The different colors correspond to different assumed masses for the eccentric planet. We also plot Eq. (5) with solid lines and Eq. (14) with dashed lines using Cv = 2. Here, we show the results obtained in the standard cases with ap = 10 au.
In the text

Fig. 2 The fate of planetesimals around the eccentric planet of ep = 0.99. The fractions of collision (blue), falling (green), and ejection (orange) are plotted as a function of the initial semi-major axis of planetesimals. We take the width of bins as 0.2 in logarithm space. Here, we show the results obtained in the standard cases with ap = 1 au (dashed line) and ap = 10 au (solid lines). The vertical dashed line shows the orbit of Θ = 1.

In the text

	Fig. 3 Averaged time before planetesimals fall to the central star or are ejected from the planetary system. Different colors correspond to the different masses of the eccentric planet. Square and cross plots show the cases with ap = 1 au and ap = 10 au, respectively. The dotted lines show the analytical equation Eq. (25) with a fitting parameter Cv = 2.
In the text

Fig. 4 The fraction of collided planetesimals Facc as a function of the initial semi-major axis of planetesimals. We plot the results of parameter studies for different assumed planetary radius Rp, the planetary mass Mp, the planetary eccentricity ep, and the initial disk thickness ⟨sin2 i0⟩1/2 from top to bottom. The left and right columns show the cases with ap = 1 au and ap = 10 au, respectively. The solid lines show the fitting formula given by Eq. (23).

In the text

	Fig. 5 Same as Fig. 4 but for different assumed planetary inclinations. The blue and green solid lines show the cases with different longitudes of pericenter ϖp. The orange dashed line shows the result for a non-inclined planet but with an inclination-excited planetesimal disk.
In the text

	Fig. 6 Same as Fig. 1, but using embryos instead of a planetesimal disk. The blue and green plots show the geometrical mean of eccentricities and inclinations of embryos. The blue solid line shows Eq. (5) with ⟨e0⟩ = 1 × 10−3. The green solid and dashed lines show Eq. (14) with ⟨sin i0⟩ = 5 × 10−4 and ⟨sin i0⟩ = hemb, respectively. Here, we show the results of 1 M⊕ embryos around 0.4 au.
In the text

Fig. 7 Same as Fig. 4 but with the blue plots showing the results for embryos instead of the planetesimal disk. We use five embryos for each simulation and sum up the collision fraction. The upper panel shows the case with 1 M⊕ embryos, and the lower panel shows the case with 0.1 M⊕ embryos. For comparison, we also plot the results for the planetesimal disk with green lines. The initial inclination of the planetesimal disk is equal to the reduced hill radius of embryos.

In the text

Fig. 8 Snapshots of the planetesimal disk when the planet’s eccentricity increases with time. The left and right columns show the cases with τecc = 103 Porb,p and 103 Porb,p, respectively. The upper and lower panels show the snapshots when the planetary eccentricity is ep = 0.3 and 0.81, respectively. The vertical dotted lines are the pericenter and apocenter of the planet. The dashed lines are the analytical models of the planetesimal eccentricity given by Eq. (40). When τecc = 105 Porb,p, the planetesimal eccentricity is smaller than that given by Eq. (40) because the resonances suppress the linear increase of the eccentricity.

In the text

Fig. 9 rcross = ap(1 − ecross,p) vs. ecross obtained by the N-body simulations. ecross and ecross,p are calculated by taking the average of eccentricities of planetesimals whose orbits are just before crossing to the planet’s orbit. The dashed lines show ecross obtained from Eqs. (40) and (41). The thick dashed lines mean Eq. (44) is achieved. When Eq. (44) is achieved, ecross becomes smaller than the analytical models predict. As indicated in the legend, the different colors correspond to the different eccentricity excitation timescales τecc.

In the text

	Fig. 10 Same as Fig. 4 but with a planet whose eccentricity increases with a given timescale τecc. The blue solid lines show the numerical results and the green dashed lines show Eq. (45).
In the text

	Fig. 11 Same as Fig. 1, but for the results of parameter study of the planetary eccentricity ep.
In the text

Fig. 12 Distribution of ϖ − ϖp around an eccentric planet. We plot the cases of ep = 0.1, 0.3, and 0.5 in the left, center, and right columns, respectively. The upper, middle, and lower panels show different times. The mean eccentricities of planetesimals divided by the planetary eccentricity are indicated in the upper left corner of each panel. Once ⟨e⟩ exceeds ep, the planetesimals’ orbits start to align to the planetary orbit.

In the text

Fig. 13 Semi-major axis vs. eccentricity of observed exoplanets. The data is obtained from Data from Data & Analysis Center for Exoplanets (DACE) in April 2024. Here, we include exoplanets heavier than Neptune, with a measurement uncertainty of less than 20 % in eccentricity. The grey dashed lines show the boundary of the regions where planets could be tidally disrupted (qp < 0.01 au) and where tidal forces with the central star are sufficient to cause inward planetary migration (qp < 0. 05 au).

In the text

Fig. 14 Estimated mass of accreted solid materials Macc normalised with πap2∑sol for various planets. The left panel shows Macc/πap2∑sol as a function of the planetary eccentricty. The different color corresponds to the different semi-major axis of the planet. The planetary mass is set to Jupiter’s mass. The right panel shows Macc/πap2∑sol as a function of the planetary mass. Here, the different color corresponds to the different eccentricity of the planet. The semi-major axis is set to 1 au. In both panels, the solid and dashed lines show the cases of sin ip = 10−2 and 10−1, respectively. We set τecc to 103 Porb,p assuming the planet-planet scattering.

In the text

Fig. 15 Estimated mean collision number of planets onto an eccentric planet. We show the mean collision number as a function of the planetary mass. The left and right panel show the cases with τecc = 103 Porb,p and τecc = 105 Porb,p, respectively. The blue, green, and orange lines correspond to ep = 0.3, 0.6, and 0.9, respectively. The solid and dashed lines show the cases of sin ip = 10−2 and 10−1, respectively.

In the text

	Fig. 16 Fraction of planetesimals that fall onto the central star as a function of the planetesimals’ initial semi-major axis. The blue, green, and orange lines show the cases of τecc = 103 Porb,p, 104 Porb,p, and 105 Porb,p, respectively.
In the text

	Fig. A.1 Time evolution in the errors of orbital elements. Blue, green, and orange lines show the error in the semi-major axis, eccentricity, and longitude of the pericentre, respectively.
In the text

Fig. B.1 Same as Fig. 8, but shows the difference in the pericenters between planetesimals and the planet instead of planetesimals’ eccentricity. The left and right columns show the cases with τecc = 103 Porb,p and 105 Porb,p, respectively. The upper and lower panels show the snapshots when the planetary eccentricity is ep = 0.1 and 0.81, respectively. The vertical dotted lines are the pericenter and apocenter of the planet.

In the text

Fig. B.2 Left columns: Time evolution of sin(ϖ − ϖp) for each planetesimal. Right columns: Amplitude spectrum obtained with Fast Fourier Transform. The amplitude spectrum is scaled with N/2, where N is the number of data points. Each planetesimal has a different initial semi-major axis, and it is shown in the upper right corner of the panels in the right column divided by the pericenter of the planet qp = 5 au.

In the text

Fig. C.1 Snapshots of orbital integration around an eccentric giant planet. The left and right columns show the eccentricity and inclination, respectively, with the semi-major axis of planetesimals. Time progresses from upper to lower panels. The solid lines in the left column show efall. The small boxes on the right side of each panel show the histogram. The total number of planetesimals N is written in the boxes. Here, we show the case with ap = 1 au, ep = 0.99, sin ip = 0, Mp = MJup and Rp = RJup.

In the text

	Fig. C.2 Same as Fig. C.1, but for ep = 0.1. Unlike in Fig. C.1, here, many planetesimals are trapped in the co-rotation resonance of the planet.
In the text