Issue 
A&A
Volume 659, March 2022



Article Number  A28  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202142180  
Published online  01 March 2022 
The origin of the high metallicity of closein giant exoplanets
II. The nature of the sweet spot for accretion
^{1}
Institute for Computational Science (ICS), University of Zurich,
Zurich,
Switzerland
email: sho.shibata@uzh.ch
^{2}
Department of Earth and Planetary Science, Graduate School of Science, The University of Tokyo,
731 Hongo,
Bunkyoku,
Tokyo
1130033, Japan
^{3}
Division of Science, National Astronomical Observatory of Japan,
2211 Osawa,
Mitaka,
Tokyo
1818588, Japan
^{4}
Department of Astronomical Science, The Graduate University for Advanced Studies (SOKENDAI),
2211 Osawa,
Mitaka,
Tokyo
1818588, Japan
Received:
8
September
2021
Accepted:
23
December
2021
Context. The composition of gas giant planets reflects their formation and evolution history. Revealing the origin of the high heavyelement masses in giant exoplanets is an objective of planet formation theories. Planetesimal accretion during the phase of planetary migration could lead to the delivery of heavy elements into gas giant planets. In our previous paper, we used dynamical simulations and showed that planetesimal accretion during planetary migration occurs in a rather narrow region of the protoplanetary disk, which we refer to as the “sweet spot” for accretion.
Aims. Our understanding of the sweet spot, however, is still limited. The location of the sweet spot within the disk and how it changes as the disk evolves were not investigated in detail. The goal of this paper is to reveal the nature of the sweet spot using analytical calculations and to investigate the role of the sweet spot in determining the composition of gas giant planets.
Methods. We analytically derived the required conditions for the sweet spot. Then, using the numerical integration of the orbits of planetesimals around a migrating planet, we compared the derived equations with the numerical results.
Results. We find that the conditions required for the sweet spot can be expressed by the ratio of the aerodynamic gas damping timescale of the planetesimal orbits to the planetary migration timescale. If the planetary migration timescale depends on the surface density of disk gas inversely, the location of the sweet spot does not change with the disk evolution. We expect that the planets observed inner to the sweet spot include a much greater amount of heavy elements than the planets outer to the sweet spot. The mass of planetesimals accreted by the protoplanet in the sweet spot depends on the amount of planetesimals that are shepherded by mean motion resonances. Our analysis suggests that tens Earthmasses of planetesimals can be shepherded into the sweet spot without planetesimal collisions. However, as more planetesimals are trapped into mean motion resonances, collisional cascade can lead to fragmentation and the production of smaller planetesimals. This could affect the location of the sweet spot and the population of small objects in planetary systems.
Conclusions. We conclude that the composition of gas giant planets depends on whether the planets crossed the sweet spot during their formation. Constraining the metallicity of cold giant planets, which are expected to be beyond the sweet spot, with future observations would reveal key information for understanding the origin of heavy elements in giant planets.
Key words: methods: numerical / planets and satellites: composition / planets and satellites: formation / planets and satellites: gaseous planets
© ESO 2022
1 Introduction
The formation of gas giant planets is a complex process that involves many underlying physical processes, such as core formation, gas accretion, and planetary migration. Recent planetary population synthesis models (e.g. Ida et al. 2018; Mordasini 2018) aim to consider as many processes as possible and to construct a comprehensive planetary formation model that begins at the early stages of planetesimal accretion and continues up to the formation of gas giant planets. However, there are still many uncertainties related to each of these processes and to the initial conditions. Constraining theories using observed physical parameters is crucial to understanding the formation and evolution of gas giant planets. Along with the development of observation and characterisation of exoplanets, many studies have aimed to validate theoretically proposed formation models and to constrain the initial parameters that control giant planet formation. The planetary chemical composition is an important characteristic that can be used to constrain formation models because the formation history is imprinted in the composition of gas giant planets.
Recent observations provide information on the chemical composition of gas giant planets. Typically, the planetary metallicity is used to characterise them and to constrain their formation and evolution histories (e.g. Helled et al. 2021). The composition of gas giant planets was first constrained for the solar system giant planets Jupiter and Saturn (e.g. Guillot 1999; Saumon & Guillot 2004). From the gravitational moments observed by spacecrafts, the total heavyelement mass for Jupiter (M_{Z,Jup}) was estimated to be ~25–45 M_{⊕} (Wahl et al. 2017; Debras & Chabrier 2019) and for Saturn M_{Z,Sat} ~ 16–30 M_{⊕} (Saumon & Guillot 2004; Helled & Guillot 2013; Movshovitz et al. 2020; Mankovich & Fuller 2021). Heavyelement enrichment has also been inferred for giant exoplanets (e.g. Guillot et al. 2006; Miller & Fortney 2011; Thorngren et al. 2016). According to Thorngren et al. (2016), closein gas giant planets contain large amounts of heavy elements of the order of tens of Earthmasses, with bulk planetary metallicities that are significantly higher than stellar metallicities. Surprisingly, some of these planets are estimated to consist of more than ~ 100 M_{⊕} of heavy elements. Revealing the origin of this significant heavyelement enrichment in giant planets is a key objective in giant planet formation theory.
Planetesimal accretion has been proposed to be one of the main sources of heavyelement enrichment of giant planets. After the onset of runaway gas accretion, a large number of planetesimals enter the expanding feeding zone where they can be captured by the protoplanet (Zhou & Lin 2007; Shiraishi & Ida 2008; Shibata & Ikoma 2019; Podolak et al. 2020). During the rapid gas accretion phase, up to 30% of planetesimals inside the feeding zone are captured by the protoplanet (Shibata & Ikoma 2019). For a massive protoplanetary disk, this accretion process is sufficient to explain the estimated heavyelement masses in Jupiter and Saturn. For closein giant exoplanets, however, the heavyelement mass that accreted during rapid gas accretion is rather low because the total mass planetesimals inside the feeding zone of these planets is significantly lower.
During the planetary migration phase, migrating planets encounter planetesimals in wide regions within the protoplanetary disk. Planetesimal accretion during the migration phase was investigated in the context of accelerating the formation of Jupiter’s core (e.g. Tanaka & Ida 1999; Alibert et al. 2005). It was found that the formation timescale of Jupiter’s core can be shortened if planetary migration is included. For the origin of large amounts of heavy elements in closein giant exoplanets, planetesimal accretion bya migrating Jupitersize planet was also investigated by Shibata et al. (2020) and Turrini et al. (2021). Shibata et al. (2020) performed orbital integration calculations of planetesimals around a migrating planet, and found that several tens Earthmasses of planetesimals can be accreted by a gas giant planet if the planet migrates tens of AU in a massive protoplanetary disk. Turrini et al. (2021) also investigated this process and found that planetesimal accretion during planetary migration is a useful tracer of the formation location of giant planets. Thus, planetesimal accretion during the planetary migration phase is a key mechanism in the formation history of gas giant planets.
If the migrating planet is as massive as Jupiter, efficient planetesimal accretion is limited to a narrow ringlike region within the protoplanetary disk (Shibata et al. 2020); in this paper we refer to this region as the sweet spot for planetesimal accretion (hereafter SSP). Even in the case of a migrating Earthmass planet the accretion area is limited to an outer disk region (without outer edge) because of the disk gas drag (Tanaka & Ida 1999). However, such a ringlike region of the SSP is specific to a Jupitermass planet. The shape and location of the accretion region during the planetary migration determine the accretion area and the accretion efficiency. The accretion area of planets with different masses depends on the strength of mean motion resonances accompanied by the migrating planet. Nevertheless, the detailed nature of the SSP is not fully understood.
In this paper we investigate the nature of the SSP in detail and derive analytical equations for determining its location. In Sect. 2, we analytically investigate the orbital evolution of planetesimals around the migrating planet. We derive a condition required for planetesimal accretion during the planetary migration phase. In Sect. 3, we compare the derived condition with the numerical simulations. In Sect. 4, we discuss the location of SSP in a protoplanetary disk and the heavyelement enrichment of giant planets. Our conclusions are summarised in Sect. 5.
Fig. 1 Orbital evolution of a planetesimal initially located at 13.8 AU dictated via the gravitational interaction with a migrating Jupitermass protoplanet. The dots and crosses show the locations of the planetesimal and the protoplanet, respectively, and are colourcoded according to time. In this simulation the surface density of disk gas is set to that of the minimummass solar nebula and aerodynamic gas drag is included. Here, the protoplanet migrates from 15.0 AU with a migration rate of 10^{5} AU yr^{−1}. For the details of the numerical simulation, see Sect. 3. 
2 Analytical expression of the accretion sweet spot
2.1 Overview of planetesimal accretion during planetary migration
Mean motion resonances play an important role in the accretion of planetesimals onto a migrating protoplanet. Here we summarise the orbital evolution of planetesimals dictated via their gravitational interaction with a migrating protoplanet, which was found by the numerical simulations of Shibata et al. (2020). Figure 1 shows the orbital evolution of a planetesimal initially located at 13.8 AU around a migrating Jupitermass protoplanet. The orbital evolution of the planetesimal is divided into three characteristic phases as discussed below.
In the first phase (t = 0 yr–4 × 10^{5} yr; blue symbols) the planetesimal encounters the migrating protoplanet and is trapped in the 2:1 mean motion resonance; this phenomenon is referred to as resonant trapping. The planetesimal trapped in a mean motion resonance is transported inward together with the migrating protoplanet and its eccentricity is highly enhanced. This phenomenon is known as resonant shepherding (Batygin & Laughlin 2015).
In the second phase (t = 4 × 10^{5} yr–7 × 10^{5} yr; green–yellow symbols), as time progresses, the resonantly trapped planetesimal starts to escape from the mean motion resonance, and its eccentricity begins to decrease because of the stronger aerodynamic drag in the inner regions. The break of resonant trapping shown in Fig. 1 is caused by overstable libration (Goldreich & Schlichting 2014). The eccentricity of the trapped planetesimal is excited through gravitational scattering by the migrating protoplanet, but is also damped by the disk gas drag. This overstable equilibrium condition makes the planetesimal orbit unstable and breaks the resonant trapping.
In the third phase (t = 7 × 10^{5} yr; red symbols), in the farther inner region, the disk gas is dense enough that the resonantlytrapped planetesimal is damped faster than the planetary migration, and therefore the planetesimal is eliminated from the feeding zone. We refer to this phenomenon as aerodynamic shepherding, which was found by Tanaka & Ida (1999) in the context of terrestrial planet formation.
According to the numerical simulations preformed in Shibata et al. (2020), planetesimal accretion occurs during the second phase, and the region in which efficient planetesimal accretion occurs is referred to as the SSP. During the first and third phases, planetesimal accretion is negligible. The existence of these three different phases (or regions) is associated with the change in strength of the aerodynamic gas drag which increases as the protoplanet migrates inward.
In the following sections, we derive the conditions for the SSP analytically. In Sect. 2.2, we show that planetesimals are generally trapped in mean motion resonances with an approaching protoplanet in standard protoplanetary disks. In Sect. 2.3, we show that the eccentricity of the resonantly trapped planetesimal reaches an equilibrium value, e_{eq}, that is determined by the ratio of the aerodynamic damping timescale (or the friction timescale) for the planetesimal (τ_{aero, 0}; Eq. (16)) to the tidal damping timescale (or the migration timescale) of the protoplanet (τ_{tide,a}). A necessary condition for a planetesimal to be captured is that it enters the feeding zone of the protoplanet. As for the resonantly trapped planetesimal, e_{eq} must exceed the eccentricity at the feeding zone boundary (e_{cross}; Eq. (27)), which is derived in Sect. 2.4. Even if this condition is satisfied, for the trapped planetesimal to escape the resonance completely the perturbed planetesimal’s orbit must oscillate more widely relative to the resonant width (i.e. the overstable libration), and this oscillation must be amplified more rapidly than the protoplanet’s migration. The conditions are derived in Sect. 2.5. Finally, by putting everything above together, we infer the condition for the SSP in Sect. 2.6.
2.2 Before resonant trapping
When the radial distance between two objects decreases, which is referred to as convergent orbital evolution, the objects are eventually locked in a resonant state. Resonant trapping has been discussed in the contexts of the formation of satellite pairs (e.g. Goldreich 1965; Dermott et al. 1988; Malhotra 1993b), the transport of small bodies (e.g. Yu & Tremaine 2001; Batygin & Laughlin 2015), and the formation of exoplanet pairs (e.g. Fabrycky et al. 2014; Goldreich & Schlichting 2014; Batygin 2015, and references therein). However, not all resonant encounters result in resonant trapping. Resonant trapping requires the following conditions for a planetesimalprotoplanet pair (e.g. Malhotra 1993a):
 (i)
The two orbits converge with each other;
 (ii)
The eccentricity of the planetesimal, e, before it is trapped in the resonance is smaller than the critical value, e_{crit} (see Eq. (2));
 (iii)
The timescale for the planetesimal to cross the resonant width, τ_{cross}, is longer than the libration timescale, τ_{lib} (see Eqs. (3) and (8)).
Here, we consider a trap of a planetesimal in the inner firstorder mean motion resonances (j : j − 1 resonance) with a protoplanet migrating inward, and we assume the mass of the planetesimal is negligibly small relative to that of the protoplanet.
The first condition can be obtained by comparing the timescales of change in the semimajor axis of the planetesimal and the protoplanet. In a protoplanetary gaseous disk the orbits of objects shrink by the aerodynamic gas drag and gravitational tidal drag from the disk gas. The former is dominant for planetesimalsize objects (≲ 10^{24} g) and the latter is for planetsize objects (≳10^{24} g) (Zhou & Lin 2007). Using the aerodynamic damping timescale for semimajor axis τ_{aero,a} for the planetesimal and the tidal damping timescale for semimajor axis τ_{tide,a} for the protoplanet, we can write the first condition as (1)
For the second condition, the critical eccentricity for the j : j − 1 resonances is given by (e.g. Murray & Dermott 1999) (2)
where M_{p} and M_{s} are the masses of the planet and central star, respectively, and f_{d} is the interaction coefficient whose values are summarised in Table 1.
When τ_{aero,a} ≫ τ_{tide,a}, the timescale for the planetesimal to cross the resonant width τ_{cross} is given by (3)
where a_{c} is the semimajor axis of the resonance centre, which is related to the protoplanet’s semimajor axis a_{p} as (4)
ȧ_{c} is the temporal change of a_{c}; and Δa_{res} is the width of the resonance given by (e.g. Murray & Dermott 1999) (5)
The width of the resonance depends on the eccentricity and takes the minimum value at (6)
By substituting Eq. (6) into Eq. (5), we obtain the minimum width of the resonance as (7)
The libration timescale at internal resonances is given by (e.g. Murray & Dermott 1999; Batygin & Laughlin 2015) (8)
where n is the planetesimal mean motion. Substituting Eq. (6) into Eq. (8), we obtain the corresponding libration timescale as (9)
where T_{K} (= 2π∕n) is the Kepler period at the resonance centre. Using Eqs. (3), (7), and (9), we obtain the third (sufficient) condition as (10)
where T_{K,p} is the Kepler period of the migrating protoplanet, specifically T_{K, p} = T_{K} j∕(j − 1).
Far from the protoplanet and outside mean motion resonances, the planetesimal’s eccentricity is determined by the balance between the viscous stirring from the planetesimal swarm and the aerodynamic gas drag. The mean value of r.m.s. eccentricities of the planetesimal swarm is of the order of ≲ 10^{−2} (Ohtsuki et al. 2002). For kilometresize planetesimals with e = 10^{−2}, the radial drift timescale τ_{damp,a} is longer than the planetary migration timescale in the type II regime (τ_{tide,a} ≲ 10^{7} yr) over a wide region of a protoplanetary disk (Adachi et al. 1976). In addition, the critical eccentricity is e_{crit} ~ 0.15 for j = 2 and M_{p} ∕M_{s} = 10^{−3}. Thus, the first and second conditions are expected to be achieved in protoplanetary disks. Figure 2 shows the timescales ratio τ_{cross}∕τ_{lib} as a function of semimajor axis, calculated from Eq. (10). Here, we show the cases of 2 : 1 (thick lines) and 3 : 2 (thin lines) resonances, and consider the type II regime for the migration timescale τ_{tide,a} obtained in Kanagawa et al. (2018). To calculate τ_{tide,a}, we adopt the selfsimilar solution of LyndenBell & Pringle (1974) for the disk structure with a disk accretion rate of 10^{−8} M_{⊙} yr^{−1} and viscosity of 10^{−3}. The disk aspect ratio is set to 0.03r^{1∕4}. The timescales ratio τ_{cross}∕τ_{lib} is higher than unity except for the cases of M_{p}∕M_{s} = 10^{−4} in outer disk a ≫ 10 AU. Thus, we conclude that the resonant trapping of planetesimals usually occurs during the migration of a giant protoplanet.
Values of the secular term f_{s} and the directterm f_{d} for the firstorder mean motion resonances.
Fig. 2 Ratio of the crossing timescale τ_{cross} (Eq. (3)) to the libration timescale τ_{lib} (Eq. (9)) as a function of the planetesimal’s semimajor axis. Here, the type II regime obtained in Kanagawa et al. (2018) is considered for the migration timescale τ_{tide,a}. The solid, dashdotted, and dotted lines show the cases of the planettostar mass ratio M_{p} ∕M_{s}= 10^{−2}, 10^{−3}, and 10^{−4}, respectively. Shown are the inner firstorder mean motion resonances (or j:j − 1 resonances): thick lines are for j = 2 and thin lines are for j = 3. To calculate the migration timescale, the same disk model was used as in Ida et al. (2018). The disk accretion rate, disk viscosity parameter, and disk aspect ratio is set to 10^{−8} M_{⊙} yr^{−1}, 10^{−3}, and 0.03r^{1∕4}, respectively. 
2.3 Orbital evolution under resonant trapping
Trapped planetesimals suffer from gravitational perturbation from the migrating protoplanet and aerodynamic gas drag from the disk gas. Lagrange’s equations under the aerodynamic gas drag are given by (e.g. Murray & Dermott 1999) (11) (12)
where φ is the resonant angle; (13)
and τ_{aero,a} and τ_{aero,e} are the aerodynamic damping timescales for semimajor axis and eccentricity, which are given by (Adachi et al. 1976) (14) (15)
Here C_{d} is the nondimensional drag coefficient, m_{pl} is the planetesimal’s mass, R_{pl} is the planetesimal’s radius, ρ_{gas} is the gas density, and v_{K} is the Kepler velocity. The planetesimal’s mass m_{pl} is calculated as , where ρ_{pl} is the planetesimal’s mean density and η_{gas} is a parameter that defines the subKeplerian velocity of disk gas as (17)
We note that we have neglected the highorder terms in eccentricity in Eqs. (11), (12), (14), and (15). Although the eccentricities of the trapped planetesimals are easily excited to e^{2} ~ 0.1, the location of the sweet spot analytically derived with these equations is found to be consistent with the numerical results, as shown in Sect. 3.
During the resonant trapping, the period ratio of the protoplanet to the trapped planetesimal remains constant, by definition, and (18)
When τ_{aero,a} ≫ τ_{tide,a}, using Eq. (18), we can rewrite Eq. (11) as (19)
and using Eq. (19) in Eq. (12) we obtain (20)
Just after being trapped in the resonances, e is small and ė is positive. Equation (20) indicates that e increases on a timescale ~e^{2} τ_{tide,a} (≪ τ_{tide,a}) and reaches an equilibrium value. Solving Eqs. (19) and (20) with ė = 0, we obtain the equilibrium eccentricity and resonant argument given by (21) (22)
Unlike the eccentricity, the planetesimals’ inclinations are not excited by the migrating protoplanet because the planetesimals are far from the protoplanet. When e ≫ i and e ≫ η_{gas}, from Eq. (15), and thus e_{eq} and sinφ_{eq} are written approximately as (23) (24)
We note that if the highorder terms of the eccentricity are considered, the equilibrium eccentricity is smaller than that given by Eq. (23).
2.4 Aerodynamic shepherding
As shown in Shibata et al. (2020), mean motion resonances play a role like narrow flow channels for planetesimals to enter the planet’s feeding zone; the channels are called the accretion bands. Mathematically speaking, planetesimals with negative values of the Jacobi energy, E_{Jacobi}, can enter the feeding zone once their eccentricities are sufficiently enhanced so that their Jacobi energy becomes positive. The Jacobi energy is given by (Hayashi et al. 1977) (25)
where h is the reduced Hill radius defined as (26)
Substituting E_{jacobi} = 0 and into Eq. (25) and assuming cos^{2}i ~ 1, we obtain the eccentricity of the accretion band, e_{cross}: (27)
The condition required for entering the feeding zone is that the equilibrium eccentricity, e_{eq} (see Eq. (23)), is larger than e_{cross}, which is given by (28)
Figure 3 shows e_{cross} as a function of M_{p}∕M_{s} for different values of j from Eq. (27). As indicated in this figure, e_{cross} decreases with increasing planetary mass; as the protoplanet mass increases, the feeding zone expands, and consequently the crossing point itself (i.e. the accretion band) is located inside the feeding zone once the planetary mass exceeds a certain value. For example, there are three crossing points (j = 2, 3, 4) when M_{p} ∕M_{s} = 10^{−3.5}, but there is only one crossing point (j = 2) when M_{p}∕M_{s} = 10^{−2.5} and no points when M_{p}∕M_{s} = 10^{−2}.
Fig. 3 Dependence of the crossing point eccentricity e_{cross} (see Eq. (27)) on the planettostar mass ratio M_{p}∕M_{s} for four different cases: j = 2 (solid), j = 3 (dashed), j = 4 (dashdotted), and j = 5 (dotted). 
2.5 Resonant breaking
Even when entering the feeding zone, planetesimals trapped in mean motion resonances are not always accreted by the migrating protoplanet.The conjunction point of a trapped planetesimal with the protoplanet is given in terms of the mean longitude as (29)
where ϖ is the longitude of pericenter. The resonant argument, φ, which is the displacement of the longitude of conjunction from the pericenter, takes φ_{eq}~ 0 during the converging orbital evolution. This means that every conjunction occurs at the pericenter of the planetesimal orbit; specifically, the planetesimal’s orbit is far from the protoplanet’s at the conjunction. Thus, trapped planetesimals cannot enter the Hill sphere of the protoplanet and planetesimal accretion does not occur (see Figs. C.1 and C.2 in Shibata et al. 2020).
A planetesimal is captured by the migrating protoplanet only if the resonant trapping is broken. Under the forces damping and enhancing the planetesimal’s eccentricity, the orbit of the trapped planetesimal becomes unstable and starts to oscillate. Once the amplitude of oscillation exceeds the resonant width, the resonant trapping is broken. This phenomenon is called the overstable libration and Goldreich & Schlichting (2014) derived the condition required for the instability to grow. Substituting Eqs. (15) and (23) into Eq. (28) of Goldreich & Schlichting (2014), we obtain the breaking condition for resonant trapping due to overstable libration under the aerodynamic gas drag as (30)
Even when this condition is satisfied it takes a time of ~ τ_{aero,e} for the instability to grow (Goldreich & Schlichting 2014). If τ_{aero,e} > τ_{tide,a}, the trapped planetesimal is shepherded into the inner disk before escaping the resonance. The second condition required for the overstable libration is τ_{aero,e} < τ_{tide,a}, which we can write as (31)
2.6 Required condition for the accretion sweet spot
The required condition for the SSP is summarised as (32)
The inner edge of the SSP is regulated by the aerodynamic shepherding and resonant shepherding. The dominating process depends on the planetary mass. On the other hand, the outer edge of SSP is regulated by the resonant shepherding alone.
To derive these equations, we assume that τ_{aero,a} ≫ τ_{tide,a}. Using the equilibrium eccentricity, we transform this condition as (33)
In the nominal disk model, η_{gas} is of the order of 10^{−3}, and this condition is easily achieved in the SSP. Thus, Eq. (32) can be used for small planetesimals as long as the gas drag strength is scaled by Eq. (16). We find that the location of the SSP is determined by the ratio of the gas drag damping timescale, τ_{aero,0}, to the planetary migration timescale, τ_{tide,a}. We therefore conclude that our theoretical prediction is robust (and rather general) since the derivation of this condition does not depend on the assumed disk and migration models. As a result, Eq. (32) can be generally applied to various models of planetesimal accretion.
3 Comparison with numerical results
In this section we validate the derived equation for the SSP by comparison with numerical results. To this end, we perform Nbody simulations similar to those presented in Shibata et al. (2020). We adopt different values of model parameters from those in Shibata et al. (2020) in order to identify whether the SSP is indeed regulated by the ratio of the two timescales, τ_{aero,0} and τ_{tide,a}.
3.1 Numerical model and settings
We considerthe following scenario. Gas accretion has terminated and the young giant planet no longer grows in mass. The planet then migrates radially inward from a given semimajor axis within a protoplanetary disk. Initially there are many equalsized planetesimals interior to the planet’s orbit. The migrating planet then encounters these planetesimals and captures some of them. The planetesimals are represented by test particles, and therefore are affected only by the gravitational forces from the central star and planet, and the drag force by the disk gas. Here, we summarise our numerical model. The details of our numerical model are described in Appendix A.
We assume that the planetary migration timescale is given by (34)
where τ_{tide,0} is a scaling parameter. For the drag force by the disk gas, we adopt a model by Adachi et al. (1976). The gas drag force depends on the size of planetesimal R_{pl} and the profile of ambient disk gas. To focus on the dynamic process, we adopt simply the minimummass solar nebula model (Hayashi 1981) as our gas disk model. The surface density, Σ_{gas}, and the disk temperature, T_{disk}, are given by (35) (36)
where r is the radial distance from the initial mass centre of the starplanet system, Σ_{gas,0} = 1.7 × 10^{3} g cm^{−2}, α_{disk} = 3∕2, T_{disk,0} = 280 K, and β_{disk} = 1∕4. The protoplanetary disk is assumed to be vertically isothermal.
For the planetesimals we adopt a simple surface density profile given by (37)
where Z_{s} is the solidtogas ratio (or the metallicity) and . To speed up the numerical integration, we follow the orbital motion of superparticles, each of which contains several equalsized planetesimals. The superparticles are initially distributed in a_{pl,in} < r < a_{pl,out} with a_{pl,in} = 0.3 AU and a_{pl,out} = 20 AU. The planet is initially located in such a way that the inner boundary of the feeding zone is consistent with the outer edge of the planetesimals disk: (38)
The calculation is artificially stopped once the planet reaches the orbit of a_{p,f}= 0.5 AU. First, we investigate the planetesimal accretion in a reference case, where τ_{tide,0} = 1 × 10^{5} yr, R_{pl}= 1 × 10^{7} cm, and M_{p}= 1 × 10^{−3}M_{s}. In the reference case the total number of planetesimals is M_{tot} ~ 43 M_{⊕}. The choices of the parameter values for the reference model are summarised in Table 2. We perform a parameter study for various values of τ_{tide,0}, R_{pl}, and M_{p}. The dynamical integration for the planetesimals is performed using the numerical simulation code developed by Shibata & Ikoma (2019). In this code we integrate the equation of motion using the forthorder Hermite integration scheme (Makino & Aarseth 1992). For the timesteps we adopt the method of Aarseth (1985).
Parameters used in the reference model.
Fig. 4 Planetesimal accretion rate vs. the semimajor axis of the migrating planet. The vertical solid lines showthe boundaries of the SSP given by Eq. (32). The red lines are for the case of j = 2 and the green lines are for j = 3. The filled areas show the SSP. 
3.2 Comparison with the numerical result
Figure 4 shows the planetesimal accretion rate at each semimajor axis of the migrating protoplanet, obtained in the simulation for the reference model. Planetesimal accretion is found to occur efficiently in a region of 2 AU ≲ a_{p} ≲ 10 AU. Beyond 10 AU, planetesimals are trapped in mean motion resonances of j = 2 and 3 with the migrating protoplanet. Our analysis of the numerical results shows that the condition of overstable libration is achieved, but the timescale of the instability growth is longer than the migration timescale; thus, the trapped planetesimals are shepherded without being captured by the protoplanet. Around ~ 10 AU, planetesimals trapped in the 2:1 mean motion resonance start to escape due to the overstable libration. The large number of shepherded planetesimals are supplied into the planetary feeding zone, so that the accretion rate peaks at around ~ 8 AU. Aerodynamic shepherding starts around 5 AU for the planetesimals in the 2:1 mean motion resonance, halting the supply of planetesimals through the 2:1 mean motion resonance and reducing the accretion rate significantly. Instead, planetesimal supply resumes through the 3:2 mean motion resonance around 7 AU and the accretion rate increases again. More specifically, planetesimals that remain to be brought into the feeding zonethrough 2:1 mean motion resonance move to the 3:2 resonance (i.e. the accretion band shifts from j = 2 to j = 3). The accretion rate peaks again around 3 AU and then rapidly decreases around 2 AU due to aerodynamic shepherding.
In the simulation setting the timescales ratio is given by (39) (40)
The redshaded and greenshaded areas in Fig. 4 show the predicted regions of the SSP, 5 AU ≲ a_{p} ≲ 10 AU for j = 2 and 1 AU ≲ a_{p} ≲ 7 AU for j = 3, which are obtained using Eqs. (40) and (41) in Eq. (32). These regions turn out to cover the areas where the accretion rate is higher than 10^{−6} M_{⊕} yr^{−1}. Thus, we conclude that the derived equation for the SSP reproduces well the numerical results in the reference model.
Fig. 5 Results of the parameter study regarding the size of planetesimals R_{pl}. Left: planetesimal accretion rate shown as a colour contour on a plane of the semimajor axis of the migrating protoplanets and the radius of the planetesimals, The analytically predicted inner and outer boundaries of the sweet spot are indicated with solid and dashed lines, respectively. The red lines are for j = 2 and the green lines are for j = 3. Right: total mass of captured planetesimals as a function of the planetesimal radius. 
3.3 Results of parameter studies
As seen in Eq. (32), the location of the SSP depends on the ratio of the gas drag damping timescale, τ_{aero,0}, to the planetary migration timescale, τ_{tide,a}, and on the planetary mass, M_{p}. To investigate this dependence, we perform parameter studies for the size of planetesimals, R_{pl}, which determines τ_{aero,0}; the scaling factor of migration timescale, τ_{tide,0}; and the planetary mass, M_{p}.
3.3.1 Dependence on the planetesimal size
Since the gas drag damping timescale, τ_{aero,0}, increases with the size of planetesimals, the location of the SSP is predicted to be closer to the central star for larger planetesimals (see e.g. Eqs. (40) and (41)). We perform the numerical simulations, changing the size of planetesimals from 1 × 10^{5} cm to 1 × 10^{8} cm (or ~ 0.01M_{⊕}). We note that the assumption that mutual gravitational interaction among planetesimals is negligible is not valid for large planetesimals with sizes of R_{pl} ~ 10^{8} cm. Nevertheless, a study for a wide parameter range is useful to deepen the physical understanding of the effect of the damping timescale onplanetesimal accretion.
Figure 5 shows the results of the parameter study when changing the planetesimal’s radius, R_{pl}. In the left panel the planetesimal accretion rate is shown as a colour contour on a plane of the semimajor axis of the migrating protoplanet and the planetesimal’s radius. The analytically predicted inner and outer boundaries of the SSP derived in Sect. 2.3 are shown with the solid and dashed lines, respectively. The red lines are for j = 2 and the green lines are for j = 3. In the right panel, the total mass of captured planetesimals M_{cap, tot} is plotted asa function of the planetesimal’s radius.
The location of the SSP (i.e. the region of high accretion rates) is found to be well reproduced by the analytical expressions derived in Sect. 2.3. The maximum accretion rate and the total captured heavyelement mass are higher for larger planetesimals. This occurs because as the SSP moves inward, the number of planetesimals shepherded into the SSP increases. More planetesimals are shepherded into the SSP when the planetesimals are larger, where the SSP locates farther from the initial planet location (or closer to the central star).
Fig. 6 Same as Fig. 5, but for the results of the parameter study regarding the migration timescale τ_{tide,0}. The grey circles in the right panel show the results of the parameter study regarding the size of planetesimals, R_{pl}, for comparison using the fraction of timescales τ_{aero,0}∕τ_{tide,a} on the vertical axis. 
3.3.2 Dependence on the planetary migration timescale
Next, we perform a parameter study when changing the planetary migration timescale. Our analytic consideration above indicates that the location of the SSP is closer to the central star for faster planetary migration. We perform numerical simulations with various migration timescales by changing the scaling factor of the migration timescale, τ_{tide,0}, from 1 × 10^{4} yr to 1 × 10^{6} yr.
Figure 6 shows the dependence of the results on the migration timescale, τ_{tide,0}. In the left panel the SSP is found to be closer to the star for faster planetary migration (or smaller τ_{tide,0}). Since the SSP moves inward as τ_{tide,0} decreases, for the same reason described above, the total mass of captured planetesimals, M_{cap, tot}, increases with decreasing τ_{tide,0}, which is confirmed in the right panel. The right axis shows τ_{aero,0}∕τ_{tide,a} at 1 AU calculated with Eq. (41). For comparison, we also plot with grey circles the M_{cap, tot} versus τ_{aero,0}∕τ_{tide,a} relation by converting the result shown in the right panel of Fig. 5 when using Eq. (41). Except for τ_{tide,0} ≲ 10^{4.6} yr, the two results are similar. This confirms that the location of the SSP scales with the ratio of the timescales, τ_{aero,0}∕τ_{tide,a}.
The effects of planetary migration, however, seem somewhat more complicated relative to the dependence on the planetesimal’s radius. First, as shown in the left panel, for relatively fast migration (τ_{tide,0} ≲ 10^{4.6} yr), the outer boundary of the SSP expands beyond the one inferred analytically. Second, the M_{cap, tot} versus τ_{aero,0}∕τ_{tide,a} relation differs between the cases when we vary R_{pl} and varying τ_{tide,0} for τ_{tide,0} ≲ 10^{4.6} yr. This is related to the stability of resonant trapping. As shown in Sect. 2.2, the condition required for resonant trapping is that the timescale for planetesimals to cross the mean motion resonance is longer than the libration timescale. Using Eq. (10) we find that stable resonant trapping is difficult for fast planetary migration like τ_{tide,0} ≲ 10^{4} yr.
Figure 7 shows the changes in the phase angle of a planetesimal initially located at a_{0} = 13.8 AU in the cases of (a) τ_{tide,0} = 10^{4} yr and (b) τ_{tide,0} = 10^{5} yr. In panel a, the phase angle continues libration even after trapped into the resonance. This occurs because fast planetary migration perturbs the orbit of the trapped planetesimal and accelerates the overstable libration. The resonant trapping is thus broken earlier in the case of τ_{tide,0} = 10^{4} yr than in the case of τ_{tide,0} = 10^{5} yr, which results in the outward expansion of the SSP. However, such an outward shift of the outer edge of the SSP does not lead to an increase in the mass of captured planetesimals because fast planetary migration breaks the resonant trapping as soon as the trapped planetesimal reaches the equilibrium condition. The planetesimal in panel a escapes from the resonant trapping with e ~ 0.7 and that in panel b with e ~ 0.4. The capture probability of planetesimals decreases with increasing eccentricity (Ida & Nakazawa 1989); thus, the total mass of captured planetesimals does not increase with the expansion of the SSP in the cases of τ_{tide,0} ≲ 10^{4.6} yr.
Fig. 7 Change in the phase angle of a planetesimal initially located at a_{0} = 13.8 AU. Panel a shows the case of τ_{tide,0} = 1 × 10^{4} yr and panel b shows the case of τ_{tide,0} = 1 × 10^{5} yr. In panel a, the phase angle oscillates indicating that the resonant trapping is unstable due to fast planetary migration. Then, the resonant trapping is broken earlier in the case of panel a than that shown in panel b. 
3.3.3 Dependence on the planetary mass
Finally, we perform a parameter study of the simulations when varying the mass of the migrating protoplanet, M_{p}, in the range of M_{p}∕M_{s} between 1 × 10^{−4} and 1 × 10^{−2}. The results are summarised in Fig. 8. We find that the semimajor axis of the outer boundary of the SSP is insensitive to the protoplanet mass, which is consistent with Eq. (32). In this figure, we indicate the locations of the inner edge of the SSP for j = 4 and 5 in addition to those for j = 2 and 3 given by Eq. (32), confirming that Eq. (32) accurately reproduces the numerical result. As shown in Fig. 3, the crossing points between the resonance centre and the feeding zone boundary disappear once the planetary mass exceeds a specific value. In Fig. 8, we show the M_{p} –a_{p} relation corresponding to the disappearance of the crossing points with grey solid lines. This is also consistent with the numerical result.
As found in Shibata et al. (2020), the total mass of captured planetesimals M_{cap,tot} changes with the planetary mass in a nonmonotonic manner and takes local maxima at M_{p} ∕M_{s} ~ 10^{−3.6}, 10^{−3.3}, 10^{−2.7}, and 10^{−2.2}. These periodical changes come from the corresponding shifts of the inner boundary of the SSP. The shepherding process that regulates the inner boundary of the SSP changes from the aerodynamic shepherding to the resonant shepherding with increasing planetary mass. The accretion bands also change with the protoplanet mass. Because of these effects the location of the inner boundary of the SSP shifts nonmonotonically which results in the nonmonotonic change of M_{cap,tot}. The global decreasing trend with increasing protoplanet mass comes from the outward shift of the inner boundary.
Fig. 8 Same as Fig. 5, but for the results of the parameter study when varying the protoplanet’s mass, M_{p}. The red, green, sky blue, and blue lines represent the SSP inner (solid lines) and outer (dashed lines) boundaries determined by the resonant trapping, which is predicted by Eq. (32) for j = 2, 3, 4, and 5, respectively. The coloured lines change into the grey lines once the protoplanet mass exceeds a threshold value, beyond which the accretion bands disappear (see Fig. 3). 
4 Discussion
4.1 Collision and ablation of highly excited planetesimals
Our model neglects the effect of planetesimal collisions. During planetary migration, planetesimals are swept by the resonances, and the enhanced local density of planetesimals increases the frequency of planetesimal collisions, which might initiate a collisional cascade (Batygin & Laughlin 2015). A collisional cascade (e.g. Dohnanyi 1969) changes the size distribution of planetesimals and the perturbation adding on the planetesimal orbits by the collisions would also affect the resonant configuration (Malhotra 1993b). In this section we consider the effect of planetesimal collisions on the location of the SSP and the planetesimal accretion.
4.1.1 Collision timescale of shepherded planetesimals
First, we estimate the collision timescale using the number density of planetesimals, n_{pl} ; the collision cross section, σ; and the relative velocity, u. We considera swarm of planetesimals that is trapped in mean motion resonances with total mass of M_{shep}. For simplicity, we consider that these planetesimals have the same mass. The planetesimals are distributed in the ringlike region with a width of and thickness of . The collision timescale between these planetesimals τ_{col} is given as (42) (43) (44)
Here we neglect the effect of gravitational focusing because of the high relative velocities of planetesimals. The value of M_{shep} increases as the protoplanet migrates and sweeps planetesimals, but decreases after reaching the SSP because the trapped planetesimals start to escape from mean motion resonances. As suggested by Eq. (44), τ_{col} rapidly decreases with the inward planetary migration. Planetesimals start to collide with each other where the collision timescale, τ_{col}, is comparable to the migration timescale, τ_{tide,a}. The maximum mass of planetesimals shepherded into the SSP without active planetesimal collisions, M_{shep,max}, is estimated as (45)
Figure 9 shows M_{shep,max} as a function of the radial distance from the central star. During the planetary migration, despite the increase in M_{shep}, M_{shep,max} rapidly decreases. We conclude that the trapped planetesimals collide with each other during the planetary migration phase if the protoplanet migrates a long distance in the radial direction before entering the SSP.
Fig. 9 Maximum mass of planetesimals that can be shepherded in mean motion resonances without effective mutual collisions. We show the cases of R_{pl} = 10^{7} cm and <i> = 10^{−3}. The solid, dashed, and dotted lines correspond to τ_{tide,a} = 10^{6} yr, 10^{5} yr, and 10^{4} yr, respectively. 
4.1.2 Collisional acscade in the resonant shepherding
The outcome of planetesimal collision is determined by the specific energy, which is the kinetic energy of impacting planetesimal normalised by the mass of impacted planetesimal, given by (e.g. Armitage 2010) (46)
where R_{im} is the radius of the impacting planetesimal. The critical value for shattering planetesimals can be scaled as (47)
Benz & Erik (1999) and Leinhardt & Stewart (2009) determined the values of fitting parameters in this equation. For icy planetesimals with an impact velocity of 3 km s^{−1} they obtained q_{s} = 1.6 × 10^{7} erg g^{−1}, q_{g} = 1.2 erg cm^{3} g^{−2}, a = −0.39, and b = 1.26. Using Eqs. (46) and (47), we find that kilometresize planetesimals are shattered by the collisions of planetesimals larger than R_{im,th}, which is given as (48)
We find that collisions of trapped planetesimals easily result in the dispersal and formation of many smaller planetesimals. The generated smaller planetesimals shortens the collisional timescale. This means that during the planetary migration phase M_{shep} increases, and once it exceeds the threshold value M_{shep,max} the trapped planetesimals start to collide with each other and initiate collisional cascade. The planetesimal size distribution is largely changed and the major source of solid material would be the smaller planetesimals. As shown by our simulations, the decrease in planetesimal size shifts the location of the SSP outward. Therefore, the possibility of collisional cascade occurring during planetary migration could affect the efficiency of planetesimal accretion.
The value of M_{shep,max} after the initiation of the collisional cascade would be affected by the amount of dispersed planetesimals left in mean motion resonances. Nondispersed planetesimals with radii of R_{pl} can be supplied as the planet migrates; however, if a large number of dispersed planetesimals larger than R_{im,th} are left in the resonance, M_{shep,max} is decreased by a factor of ~R_{im,th}∕R_{pl} at most. As we discuss below, the impacted planetesimals might be removed from the mean motion resonances, and the small fragments generated by the collisional cascades would drift into the inner disk region faster than the planet. Thus, after the initiation of the collisional cascade, the size distribution of planetesimals in resonant shepherding would be affected by various physical processes. We suggest that the mutual collisions of planetesimals in resonant shepherding could be important for determining the population of small objects in planetary systems.
In the above discussion, we assume that the relative velocity of planetesimals, u, is approximated by ~e v_{K}. However, the alignment of planetesimal orbits via gravitational perturbation from the protoplanet is known to reduce the relative velocity of planetesimals. Guo & Kokubo (2021) investigated the effect of orbital alignment due to the secular perturbation which converges the longitude of pericenter into a certain value, namely that planetesimal orbits are aligned on the inertial plane. Using Nbody simulations, they found that the relative velocity of planetesimals are reduced outside the mean motion resonances. According to their results, orbital alignment is not effective inside the mean motion resonances because of the resonant perturbation. In the planetary migration phase, however, the resonance angle φ converges to zero, and planetesimal orbits can be aligned on the plane rotating with the protoplanet (see Figs. C1 and C2 in Shibata et al. 2020). In the ideal case, the relative velocity of these planetesimals is smaller than ~ e v_{K} and some collisions would result in sticking. In addition, the mean inclination of planetesimals is difficult to estimate, but it is important for the collision timescale. The viscous stirring between trapped planetesimals might affect the mean inclination of planetesimals , which is not considered in our simulations. Under effective viscous stirring, the velocity dispersion of the planetesimal swarm keeps a simple relation of (e.g. Ohtsuki et al. 2002). In this case, the collision timescale increases and M_{shep,max} would be larger than indicated in Fig. 9. However, it is still unclear whether the energy equipartition works in the resonant trapping where the planetesimal’s eccentricity is excited by the protoplanet. It is clear that both M_{shep,max} and the planetesimal size distribution depends on the interactions between planetesimals in mean motion resonances. We suggest that this topic should be investigated in detail in future work.
Fig. 10 Same as Fig. 7, but (a) for the case with artificial perturbation given by Eq. (50) and (b) without artificial perturbation. Here, we use our reference model with R_{pl} = 1 × 10^{7}cm, τ_{tide,0} = 1 × 10^{5} yr and M_{p}∕M_{s} = 1 × 10^{−3}. 
4.1.3 Breakup of resonant trap
Highvelocity collisions of planetesimals perturb the orbit of trapped planetesimals and break the resonant trapping. Using numerical simulations, Malhotra (1993b) found that for a planet of 10 M_{⊕}, resonant trapping is broken with a velocity kick of ≳0.003v_{K}. Therefore, once M_{shep} exceeds the threshold value and the trapped planetesimals start to collide with each other, the resonant trapping can be broken by planetesimal collisions. However, the resonant width increases with increasing planetary mass (i.e. breaking the resonant trapping is more difficult for more massive planets; see Eq. (5)) and not all highvelocity collisions lead to the breaking of resonant trapping (see Fig. 6b of Malhotra 1993b). In addition, as mentioned in the previous section, the alignment of the orbits by resonant trapping might reduce the relative velocity of planetesimals. The possibility of resonant breaking by planetesimal collisions need to be investigated in future work. Here we consider the cases where the resonant trapping is broken by planetesimal collisions.
As shown in Sect. 2, planetesimal accretion occurs only in a limited ringlike region within the protoplanetary disk. This accretion process seems different from that found for Earthmass planets in Tanaka & Ida (1999). In their dynamical simulations Tanaka & Ida (1999) broke the resonant traps instantaneously by exerting artificial perturbations on the trapped planetesimals (see below), supposing that the planetesimals have their motion changed via interactions with other planetesimals. Then, comparing the timescale for the protoplanet’s migration with that for scattering of planetesimals by the migrating protoplanet, they derived a condition required for planetesimal accretion as (49)
By contrast, when deriving Eq. (32) and also performing the dynamical simulations, we considered an extreme case without modification to the planetesimal’s motion.
Using the artificial perturbation developed in Tanaka & Ida (1997), we perform additional numerical simulations.We randomise the phase angle by adding displacements in orbital elements except a, e, and i. The size of the displacement is set as (50)
where f_{p} is a scaling factor of perturbation strength andrandomly given from − 0.01 to 0.01. We add Δθ to orbital angles (longitude of ascendingnode Ω, longitude of pericentre ϖ and mean longitude at epoch ϵ) every conjunction time. Figure 10 shows the temporal change in the phase angle in the case with (a) and without (b) the artificial perturbations. As shown in panel a, even under strong perturbations (corresponding to up to 1% displacement in the phase angle at every conjunction), mean motion resonances with a Jupitermass planet are efficient in trapping planetesimals, and thereby resonance shepherding occurs. This occurs because the resonance width (see Eq. (5)) is larger for Jupitermass planets than for Earthmass planets. The perturbations is, however, found to hasten the breakup of resonant trapping because resonant trapping is weakened and overstable libration is accelerated.
Figure 11 shows the results of the parameter study when varying the planetary mass M_{p} ∕M_{s} including the artificial perturbations given by Eq. (50). For comparison, we also show the results without artificial perturbations with open circles in the right panel. Due to the artificial perturbations, the SSP shifts outward (see the left panel) and the total mass of captured planetesimals is reduced (see the right panel). In the left panel the black solid line shows the condition of shepherding obtained in Tanaka & Ida (1999) (see Eq. (49)); planetesimal accretion occurs in regions exterior to it. This black line corresponds to the inner boundary of the SSP for the case where resonant traps are perfectly broken. Although it might seem counterintuitive, we conclude that without perturbations, resonant trapping shifts the SSP inward and enhances planetesimal accretion due to the effect of accretion bands.
Fig. 11 Same as Fig. 5, but for the results of a parameter study for planetary mass M_{p} with the artificial perturbations (see Eq. (50)). In the right panel the results of a parameter study for planetary mass M_{p} without the artificial perturbations are plotted as grey circles for comparison. 
4.1.4 Ablation of eccentric planetesimals
A highly excited planetesimal has a high relative velocity compared to that of the disk gas. Strong gas drag from the gaseous disk (e.g. Pollack et al. 1986) and generated bowshock (e.g. Tanaka et al. 2013) heats the planetesimal’s surface and triggers the ablation of volatile materials. Planetesimals trapped in mean motion resonances have eccentric orbits and feel strong gas drag, and therefore ablation is expected to be important. The ablation rate depends on the rate of heat transfer from the ambient disk gas, Λ. We set Λ to range between 0.003 and 0.6 for gas drag heating, and between ~0.01 and ~0.1 for bow shock heating (see Tanaka et al. 2013, and references therein). Recently, Eriksson et al. (2021) explored the possibility of perfect ablation of planetesimals around a massive protoplanet using Λ =C_{d}∕4, which is the upper limit of the heat transfer rate (D’Angelo & Podolak 2015). According to their nominal model result, ablation of highly excited planetesimals is effective in the inner disk region ≲ 10AU and a large amount of solid material is ablated. The heating process strongly depends on the eccentricity of planetesimals. Using Eqs. (24) and (32), we find that the eccentricity of planetesimals inside the SSP decreases from ~ 0.5 to 0.1. To avoid effective ablation planetesimals must enter the SSP in the outer disk region ≳ 10AU. If the secular resonance triggered by the disk’s gravity is effective, the ablation would be effective in a disk region that is farther out (Nagasawa et al. 2019). On the other hand, if Λ is much smaller than C_{d}∕4, the ablationof planetesimals would be inefficient.
The investigation presented above concerning planetesimal collision and ablation implies that planetesimal accretion is favoured in the outer disk region ≳10 AU. Detailed future studies that focus on collision and ablation processes are clearly desirable.
4.2 Type II migration with shallow gap
In the classical picture of type II migration (e.g. Lin & Papaloizou 1993), a giant planet migrates with a deep gap opened in the vicinity of its orbit in a gaseous protoplanetary disk. For an extremely deep gap, the radial gas flow toward the central star is stemmed completely by the gap, and consequently the planet migrates with the gaseous disk contracting via viscous diffusion. Recent hydrodynamic simulations, however, reveal that the gap bottom is not so deep and the gas can cross the gap (Kanagawa et al. 2015, 2016, 2017). In this case the planetary migration is not in the classical regime; instead, the torque exerted on the planet similar to that in the type I regime (Kanagawa et al. 2018). Thus, the migration timescale in the new type II regime is written as (51)
where c is a constant of the order of unity and Σ_{gap} is the surface density of disk gas at the gap bottom; Σ_{gap} depends on the unperturbed surface density of disk gas Σ_{gas} and is given from the hydrodynamic simulations of Kanagawa et al. (2017) as (52)
where α_{vis} is the alpha parameter for disk gas turbulent viscosity (Shakura & Sunyaev 1973). From Eqs. (51) and (52), we obtain (54)
Here we discuss the location of SSP and planetesimal accretion onto the protoplanet migrating with Eq. (54).
4.2.1 Location of the accretion sweet spot
The location of the SSP depends on the ratio of the planetary migration to aerodynamic gas drag timescales. The timescale of damping by aerodynamic gas drag τ_{aero,0} in the disk midplane also depends on the disk structure (see Eq. (16)) as (56)
It should be noted that τ_{tide,a,II} is a functionof r_{p}, whereas τ_{aero,0} is a function of r_{pl}; for a planetesimal in the j:j − 1 mean motion resonance with the planet. In this case Eq. (56) can be written as a function of r_{p} as (57)
Here, we assume that the planetesimals trapped in mean motion resonances are outside the gap. If the disk viscosity is low, however, the gap slope reaches the locations of the resonances; in that case the damping timescale increases a few times. For simplicity, we neglect this effect.
From Eqs. (54) and (57) the ratio of the two timescales is given by (59)
Equation (59) indicates that the surface gas density profile only weakly affects the timescale ratio, which does not explicitly depend on Σ_{gas} and weakly depends on α_{disk}: . This means that as the protoplanetary disk evolves, the timescales ratio hardly changes, suggesting that the location of the SSP is fixed.
Assuming α_{disk} = 1, which is valid in a steady viscous accretion disk except for its outermost region (e.g. LyndenBell & Pringle 1974), and β_{disk} = 1∕4, which is valid in an optically thin disk, and substituting Eq. (59) into Eq. (32), we obtain the semimajor axes of the inner and outer edge of SSP,a_{SS,in} and a_{SS,out}, respectively,as (61) (62)
Figure 12 shows exoplanets identified so far and the predicted location of the SSP in the a_{p} –M_{p} plane for C = 10. The SSP locates from ~ 1 AU to ~ 30 AU and shifts outward with increasing planetary mass. The red and black symbols show exoplanets observed using the transit method and the radialvelocity method, respectively (http://exoplanet.eu). Given that planetary growth occurs via runaway gas accretion, followed by inward migration, the evolution path of a gas giant planet can be drawn from the bottom right to the top left in the a_{p} –M_{p} plane (e.g. Ida & Lin 2004; Mordasini et al. 2009). If runaway gas accretion begins in the outer disk region, which might be favoured for planetary core formation because of the increasing isolation mass of protoplanet (e.g. Armitage 2010), the exoplanets currently observed in the region interior to the SSP must have crossed the SSP during their formation stages.
Fig. 12 Distribution of confirmed exoplanets in the semimajor axis vs. planetary mass plane (http://exoplanet.eu). The black points are exoplanets observed by the radialvelocity method. The red points are exoplanets observed by the radialvelocity + transit methods. As for masses of exoplanets indicated with the black points, the observed values of M_{p} sin i are used. The grey area shows the theoretical sweet spot for planetesimal capture given by Eqs. (61) and (62) with C = 10. The red dashed line shows the theoretical sweet spot for planetesimal capture with C = 10^{3}, which corresponds to the case where collisional cascade have broken planetesimals down to R_{pl} ≲ 10^{5} cm. 
4.2.2 Maximum mass of planetesimals shepherded into the SSP
The total mass of accreted planetesimals depends on the number of planetesimals that are shepherded into the SSP. By substituting Eqs. (4), (54), and (62) into Eq. (45), we can estimate M_{shep,max} when the protoplanet reaches the outer edge of the SSP. In this case, M_{shep,max} weakly depends on α_{vis}, ρ_{pl}, R_{pl}, and (− 1∕4 powers of each value), but linearly depends on , Σ_{gas,0} and . Here we adopt the optically thick disk of T_{0} = 140 K (Sasselov & Lecar 2000) and consider two disks with Σ_{gas,0} = 1400 g cm^{−2} and 800 g cm^{−2}. The former disk corresponds to the massive disk of total mass 0.1 M_{s} and typical radius 100 AU. Figure 13 shows M_{shep,max} as a function of planetary mass M_{p}∕M_{s}. We find that a few tens of M_{⊕} of planetesimals can be shepherded into the SSP before initiating collisional cascade even with . For the capture of several tens or more planetesimals, the effect enhancing needs to be considered.
It is important to note that M_{shep,max} only provides the maximum mass of planetesimals shepherded into the SSP and not the mass of planetesimals that can be captured by the protoplanet. The numerical simulations we presented in Sect. 3 suggest that ~ 30% of the shepherded planetesimals are captured by the protoplanet. This accretion efficiency, however, depends on various parameters such as the inclination of planetesimals (Tanaka & Ida 1999) and the planetary capture radius (Inaba & Ikoma 2003). In addition, the initial semimajor axis of planetesimals also affects the capture efficiency because of the mean motion resonances. A detailed investigation of the accretion efficiency under various conditions should be conducted in future research.
Fig. 13 Maximum mass of planetesimals shepherded into the SSP without active planetesimal collisions M_{max,shep} as a functionof planetary mass M_{p}. Here, thecase is considered where the protoplanet migrates with the type II regime obtained by Kanagawa et al. (2018) (see Eq. (51)). The solid and dashed lines show the cases of <i> = 10^{−3} and 10^{−2}, respectively. The black and red lines are cases of massive disk with Σ_{gas,0} = 1400 g cm^{−2} and 800 g cm^{−2}, respectively. The parameters used in this plot are the optically thick disk of T_{0} = 140 K, the disk viscosity α_{vis} = 10^{−3}, the size of planetesimal R_{pl} = 10^{7} cm, and the planetesimal density ρ_{pl} = 1 g cm^{−3}. The scaling parameter of the location of the SSP C is obtained as C = 17. 
4.3 Heavyelement enrichment of giant planets
If giant planets cross the sweet spot during the migration phase, planetesimal accretion could take place. If several tens Earthmasses of planetesimals were shepherded into the SSP as shown in Fig. 13, we would expect a clear difference in the total heavyelement mass of planets depending on whether they migrated interior or exterior to the SSP. If this difference between planetary populations is observed, it could confirm that planetesimal accretion is a significant process that takes place during planetary migration. However, once the collisional cascade occurs, the location of the SSP moves to the outer disk region and the metallicity difference between the populations would be reduced. It is therefore still unclear whether such a difference is detectable.
Currently, there are no estimated metallicities for giant exoplanets within the outer disk region (i.e. at distances ≳ 10 AU). To estimate the metallicity of giant exoplanets, measurements of the radius, mass, and age of the planet are required. At the moment, data of the masses and radii of cold giant planets are still unavailable. Nevertheless, the upcoming decades are expected to provide vast observations thanks to various missions. Wide field observations, such as TESS (Transiting Exoplanet Survey Satellite, Ricker et al. 2014) and Plato (Rauer et al. 2014), will increase the number of cold giant planets. Followup observations using precise radialvelocity measurement, such as HARPS (HighAccuracy Radial velocity PlanetSearcher, Mayor et al. 2003) and ESPRESSO (Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observation, Pepe et al. 2010), and precise transit measurement, such as NGTS (Wheatley et al. 2018) and WASP (Wide Angle Search for Planets, Pollacco et al. 2006), are suitable for obtaining the masses and radii of giant planets at large radial distances. In addition, the highdispersion coronagraphy technique (e.g. Kotani et al. 2020) could constrain the atmospheric composition of cold giant planets. Next generation observations will constrain the main source of heavy elements.
In this study, we focus on planetesimal accretion and consider it as the main source of heavyelement enrichment in gas giant planets. However, other various processes could also lead to heavyelement enrichment. The metallicity of the gaseous disk can be enriched by the sublimation of dust grains and pebbles crossing each ice lines (Booth et al. 2017; Booth & Ilee 2019). The accretion of these enriched disk gas brings massive heavy elements into gas giant planets (Schneider & Bitsch 2021). Giant impacts have also been investigated as a source of enrichment (Ikoma et al. 2006; Liu et al. 2015, 2019; Ginzburg & Chiang 2020; Ogihara et al. 2021). By combining the effects of pebble accretion and giant impacts, Ogihara et al. (2021) found that up to ~ 100 M_{⊕} of heavy elements can be accreted by a giant planet through multigiant impacts of Neptunemass planets.
In order to determine what is the main source of the heavyelement enrichment in giant exoplanets, other observational information such as elemental ratios (e.g. refractorytovolatile, C/O, O/H) is valuable. In recent observations, the existence of TiO and VO in hot Jupiters is suggested by the transmission spectrum analysis (Evans et al. 2016; Chen et al. 2021). The amount of these refractory materials is useful for constraining the dominant process of heavyelement accretion. The accretion of massive refractory materials is only possible through solid accretion; thus, using the information of refractory materials, we can constrain whether the heavy elements are brought by solid accretion or gas accretion.
The elemental ratio of volatile materials also provides information on the accretion process of heavy elements (e.g. Madhusudhan et al. 2014; Notsu et al. 2020; Turrini et al. 2021; Schneider & Bitsch 2021). Turrini et al. (2021) investigated the effect of planetesimal accretion on elemental ratios of volatile materials such as C/O, N/O, or S/N. They considered a Jupitermass planet in a typical disk model; however, the relative position between the SSP and each ice lines changes with the planetary mass and disk conditions. The elemental ratio of planetesimals shepherded into the SSP would depend on the evolution pathways of the giant planets on the a_{p} –M_{p} plane and disk conditions. We hope to further investigate the link between elemental ratios and the formation history in followup studies. This is highly relevant since future missions like JWST and Ariel will provide accurate measurements of the atmospheric composition of gaseous planets that can then be used to better understand the planetary origin.
5 Summary and conclusions
Planetesimal accretion during planetary migration could lead to a significant enrichment of giant planets with heavy elements. As revealed in Shibata et al. (2020), planetesimal accretion occurs in a sweet spot for planetesimal accretion (SSP), where planetesimals can be efficiently captured by the migrating planet. By performing dynamical simulations for planetesimals with a migrating planet, Shibata et al. (2020) demonstrated that the SSP emerges due to shepherding processes caused by aerodynamic gas drag and mean motion resonances. However, they conducted no detailed analysis of the nature (especially, the location) of the SSP, and therefore in this paper we focused on the nature of the accretion sweet spot.
In Sect. 2, we analytically derived an equation describing the location of the SSP and found the following that the SSP is regulated by the ratio of the damping timescale for the planetesimal eccentricity due to aerodynamic gas drag to the migration timescale of the planet, and that the SSP also depends on the planetary mass because the accretion bands change with the relative positions between the mean motion resonances and the feeding zone.
In Sect. 3, we performed numerical simulations and confirmed that the derived conditions reproduce the numerical results well. In Sect. 4, we discussed the effect of planetesimal collisions and ablation during the shepherding process. These effects would be important for inner disk regions ≲ 10 AU and put other constraint for planetesimal accretion process. If the migration timescale is inversely proportional to the disk gas surface density (e.g. Kanagawa et al. 2018), the timescales ratio is rather insensitive to the structure of the protoplanetary disk. Therefore, the location of the SSP is fixed during the evolution of the protoplanetary disk.
Finally, we discussed the effect of the SSP on the heavyelement enrichment of gas giant planets. We suggest that the existence of the SSP makes a difference in the heavyelement mass between planets observed in the regions interior and exterior to the SSP. The detailed chemical composition of gas giant planets, such as the refractorytovolatile ratio or C/O ratio, would be also affected by the location of the SSP. If such compositional features can be observed, it could be a piece of the evidence of planetesimal accretion during planetary migration. Future observations in the upcoming decades will reveal whether planetesimal accretion is a main source of heavy elements in giant planets, and therefore advance our understanding of their origin.
Acknowledgements
S.S. and R.H. acknowledge support from the Swiss National Science Foundation (SNSF) under grant 200020_188460. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Ts 17/2–1 and by JSPS CoretoCore Program “International Network of Planetary Sciences (Planet^{2})” and JSPS KAKENHI Grant Numbers 17H01153 and 18H05439.
Appendix A Model descriptions
A.1 Equations of motion
The equation of motion is given by (A.1)
where r_{i} is the position vector relative to the initial (i.e. t = 0) mass centre of the starplanetplanetesimals system; f_{grav,i,j} is the mutual gravity between particles i and j given by (A.2)
with r_{i,j} being the position vector of particle i relative to particle j (r_{i,j} ≡r_{i,j}); M_{j} is the mass of particle j; is the gravitational constant; f_{aero} is the aerodynamic gas drag; and f_{tide} is the gravitational tidal drag from the protoplanetary disk gas. The central star, planet, and planetesimals are denoted by the subscripts i (or j) = 1, 2, and ≥ 3, respectively. The planetesimals are treated as test particles; therefore, f_{grav,i,j} = 0 in Eq. (A.1) for j ≥ 3. The aerodynamic gas drag force and the tidal drag force are given respectively by (A.3) (A.4)
Here u = v_{pl} −v_{gas} (u = u) is the planetesimal’s velocity (v_{pl}) relative to the ambient gas (v_{gas}), Given the range of the planetesimal mass (~ 10^{16}–10^{22}g) and planet (~ 10^{30}g), we assume f_{tide} = 0 for the former and f_{aero} = 0 for the latter. The central star is not affected by f_{aero} and f_{tide}.
A.2 Gas disk
The protoplanetary disk is assumed to be vertically isothermal, and the gas density ρ_{gas} is expressed as (A.6)
where z is the height from the disk midplane and h_{s} is the disk’s scale height. The aspect ratio of the protoplanetary disk is given by (A.7)
where Ω_{K} is the Kepler angular velocity and c_{s} is the isothermal sound speed of the disk gas. The sound speed is given as (A.8)
where k_{B} is the Boltzmann constant and μ_{disk} is the mean molecular weight in units of proton mass, m_{H}. In this study μ_{disk} is set as 2.3 and the aspect ratio of the disk gas at 1AU is ~ 0.03.
The gasin the protoplanetary disk rotates at subKeplerian velocity because of pressure gradient, namely v_{gas} = (1 − η_{gas})rΩ_{K}, where η_{gas} is ≪ 1 and is given as (A.9) (A.10)
where P_{gas} is the gas pressure.
A.2 Planetesimals
The surface number density of planetesimals is given as Σ_{solid}∕m_{pl}, which gives a maximum value of ~ 10^{6} ∕AU^{2}. To speed up the numerical integration, we follow the orbital motion of superparticles. The surface number density of superparticles n_{s} is given as (A.11)
here N_{sp} is the total number of superparticles used in a given simulation. In our simulation we set N_{sp} = 12, 000 and α_{sp} = 1 to distribute superparticles uniformly in the radial direction. The mass per superparticle M_{sp} is given by (A.13)
where a_{0} is the initial semimajor axis of the superparticle.
Although planetesimals are treated as test particles in our simulations, assuming that planetesimals are in reality scattered by their mutual gravitational interaction, for eccentricity and inclination we adopt the Rayleigh distribution as the initial eccentricities e and inclinations i of planetesimals. We set . The orbital angles Ω, ϖ, and ϵ are distributed uniformly.
During the orbital integration, we judge that a superparticle has been captured by the planet once (i) the superparticle enters the planet’s envelope or (ii) its Jacobi energy becomes negative in the Hill sphere. The planet’s radius R_{p} is calculated as (A.14)
where ρ_{p} is the planet’s mean density. The planet’s radius is assumed to be extended shortly after formation and the density is set to ρ_{p} = 0.125 g cm^{−3}, corresponding to a planet radius twice as large as Jupiter’s current radius.
References
 Aarseth, S. J. 1985, in IAU Symposium, 113, Dynamics of Star Clusters, eds. J. Goodman, & P. Hut, 251 [CrossRef] [Google Scholar]
 Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progr. Theor. Phys., 56, 1756 [Google Scholar]
 Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Batygin, K. 2015, MNRAS, 451, 2589 [NASA ADS] [CrossRef] [Google Scholar]
 Batygin, K., & Laughlin, G. 2015, Proc. Natl. Acad. Sci. USA, 112, 4214 [CrossRef] [Google Scholar]
 Benz, W., & Erik, A. 1999, Icarus, 142, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Booth, R. A., & Ilee, J. D. 2019, MNRAS, 487, 3998 [CrossRef] [Google Scholar]
 Booth, R. A., Clarke, C. J., Madhusudhan, N., & Ilee, J. D. 2017, MNRAS, 469, 3994 [Google Scholar]
 Chen, G., Pallé, E., Parviainen, H., Murgas, F., & Yan, F. 2021, ApJ, 913, L16 [NASA ADS] [CrossRef] [Google Scholar]
 D’Angelo, G., & Podolak, M. 2015, ApJ, 806, 203 [Google Scholar]
 Debras, F., & Chabrier, G. 2019, ApJ, 872, 100 [Google Scholar]
 Dermott, S. F., Malhotra, R., & Murray, C. D. 1988, icarus, 76, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [Google Scholar]
 Eriksson, L. E. J., Ronnet, T., & Johansen, A. 2021, A&A, 648, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Evans, T. M., Sing, D. K., Wakeford, H. R., et al. 2016, ApJ, 822, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [Google Scholar]
 Ginzburg, S., & Chiang, E. 2020, MNRAS, 498, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P. 1965, MNRAS, 130, 159 [NASA ADS] [Google Scholar]
 Goldreich, P., & Schlichting, H. E. 2014, ApJ, 147, 32 [CrossRef] [Google Scholar]
 Guillot, T. 1999, Planet. Space Sci., 47, 1183 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., Santos, N. C., Pont, F., et al. 2006, A&A, 453, L21 [CrossRef] [EDP Sciences] [Google Scholar]
 Guo, K., & Kokubo, E. 2021, AJ, 162, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 [Google Scholar]
 Hayashi, C., Nakazawa, K., & Adachi, I. 1977, Publ. Astron. Soc. Jpn., 29, 163 [Google Scholar]
 Helled, R., & Guillot, T. 2013, ApJ, 767, 113 [CrossRef] [Google Scholar]
 Helled, R., Werner, S., Dorn, C., et al. 2021, Exp. Astron., in press [arXiv:2103.08481] [Google Scholar]
 Ida, S., & Lin,D. N. C. 2004, ApJ, 604, 388 [Google Scholar]
 Ida, S., & Nakazawa, K. 1989, A&A, 220, 293 [NASA ADS] [Google Scholar]
 Ida, S., Tanaka, H., Johansen, A., Kanagawa, K. D., & Tanigawa, T. 2018, ApJ, 864, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Ikoma, M., Guillot, T., Genda, H., Tanigawa, T., & Ida, S. 2006, ApJ, 650, 1150 [NASA ADS] [CrossRef] [Google Scholar]
 Inaba, S., & Ikoma, M. 2003, A&A, 410, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015, MNRAS, 448, 994 [NASA ADS] [CrossRef] [Google Scholar]
 Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, Publ. Astron. Soc. Jpn., 68, 43 [CrossRef] [Google Scholar]
 Kanagawa, K. D., Tanaka, H., Muto, T., & Tanigawa, T. 2017, Publ. Astron. Soc. Jpn., 69 [Google Scholar]
 Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [Google Scholar]
 Kotani, T., Kawahara, H., Ishizuka, M., et al. 2020, in SPIE Conf. Ser., 11448, 1144878 [NASA ADS] [Google Scholar]
 Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199, 542 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. C. B. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 749 [Google Scholar]
 Liu, S. F., Agnor, C. B., Lin, D. N., & Li, S. L. 2015, MNRAS, 446, 1685 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, S.F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
 Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJ, 794, 12 [Google Scholar]
 Makino, J., & Aarseth, S. J. 1992, Publ. Astron. Soc. Jpn., 44, 141 [Google Scholar]
 Malhotra, R. 1993a, Icarus, 106, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Malhotra, R. 1993b, Nature, 365, 819 [CrossRef] [Google Scholar]
 Mankovich, C., & Fuller, J. 2021, Nat. Astron., 5, 1103 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Miller, N., & Fortney, J. J. 2011, ApJ, 736, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Mordasini, C. 2018, in Handbook of Exoplanets, 2425 [CrossRef] [Google Scholar]
 Mordasini, C., Alibert, Y., & Benz, W. 2009, A&A, 501, 1139 [CrossRef] [EDP Sciences] [Google Scholar]
 Movshovitz, N., Fortney, J. J., Mankovich, C., Thorngren, D., & Helled, R. 2020, ApJ, 891, 109 [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Nagasawa, M., Tanaka, K. K., Tanaka, H., et al. 2019, ApJ, 871, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Notsu, S., Eistrup, C., Walsh, C., & Nomura, H. 2020, MNRAS, 499, 2229 [Google Scholar]
 Ogihara, M., Hori, Y., Kunitomo, M., & Kurosaki, K. 2021, A&A, 648, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ohtsuki, K., Stewart, G. R., & Ida, S. 2002, Icarus, 155, 436 [Google Scholar]
 Pepe, F. A., Cristiani, S., Rebolo Lopez, R., et al. 2010, SPIE Conf. Ser., 7735, 77350F [Google Scholar]
 Podolak, M., Haghighipour, N., Bodenheimer, P., Helled, R., & Podolak, E. 2020, ApJ, 899, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Pollacco, D., Skillen, I., Cameron, A., et al. 2006, PASP, 118, 1407 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., Podolak, M., Bodenheimer, P., & Christofferson, B. 1986, Icarus, 67, 409 [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, SPIE Conf. Ser., 9143, 914320 [Google Scholar]
 Sasselov, D. D., & Lecar, M. 2000, ApJ, 528, 995 [Google Scholar]
 Saumon, D., & Guillot, T. 2004, ApJ, 609, 1170 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, A. D., & Bitsch, B. 2021, How drifting and evaporating pebbles shape giant planets I: Heavy element content and atmospheric C/O, Tech. rep. [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Shibata, S., & Ikoma, M. 2019, MNRAS, 487, 4510 [NASA ADS] [CrossRef] [Google Scholar]
 Shibata, S., Helled, R., & Ikoma, M. 2020, A&A, 633, 13 [Google Scholar]
 Shiraishi, M., & Ida, S. 2008, ApJ, 684, 1416 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., & Ida, S. 1997, Icarus, 125, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., & Ida, S. 1999, Icarus, 139, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, K. K., Yamamoto, T., Tanaka, H., et al. 2013, ApJ, 764, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Thorngren, D. P., Fortney, J. J., MurrayClay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Turrini, D., Schisano, E., Fonte, S., et al. 2021, ApJ, 909, 40 [Google Scholar]
 Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [CrossRef] [Google Scholar]
 Wheatley, P. J., West, R. G., Goad, M. R., et al. 2018, MNRAS, 475, 4476 [Google Scholar]
 Yu, Q., & Tremaine, S. 2001, AJ, 121, 1736 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, J., & Lin, D. N. C. 2007, ApJ, 666, 447 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Values of the secular term f_{s} and the directterm f_{d} for the firstorder mean motion resonances.
All Figures
Fig. 1 Orbital evolution of a planetesimal initially located at 13.8 AU dictated via the gravitational interaction with a migrating Jupitermass protoplanet. The dots and crosses show the locations of the planetesimal and the protoplanet, respectively, and are colourcoded according to time. In this simulation the surface density of disk gas is set to that of the minimummass solar nebula and aerodynamic gas drag is included. Here, the protoplanet migrates from 15.0 AU with a migration rate of 10^{5} AU yr^{−1}. For the details of the numerical simulation, see Sect. 3. 

In the text 
Fig. 2 Ratio of the crossing timescale τ_{cross} (Eq. (3)) to the libration timescale τ_{lib} (Eq. (9)) as a function of the planetesimal’s semimajor axis. Here, the type II regime obtained in Kanagawa et al. (2018) is considered for the migration timescale τ_{tide,a}. The solid, dashdotted, and dotted lines show the cases of the planettostar mass ratio M_{p} ∕M_{s}= 10^{−2}, 10^{−3}, and 10^{−4}, respectively. Shown are the inner firstorder mean motion resonances (or j:j − 1 resonances): thick lines are for j = 2 and thin lines are for j = 3. To calculate the migration timescale, the same disk model was used as in Ida et al. (2018). The disk accretion rate, disk viscosity parameter, and disk aspect ratio is set to 10^{−8} M_{⊙} yr^{−1}, 10^{−3}, and 0.03r^{1∕4}, respectively. 

In the text 
Fig. 3 Dependence of the crossing point eccentricity e_{cross} (see Eq. (27)) on the planettostar mass ratio M_{p}∕M_{s} for four different cases: j = 2 (solid), j = 3 (dashed), j = 4 (dashdotted), and j = 5 (dotted). 

In the text 
Fig. 4 Planetesimal accretion rate vs. the semimajor axis of the migrating planet. The vertical solid lines showthe boundaries of the SSP given by Eq. (32). The red lines are for the case of j = 2 and the green lines are for j = 3. The filled areas show the SSP. 

In the text 
Fig. 5 Results of the parameter study regarding the size of planetesimals R_{pl}. Left: planetesimal accretion rate shown as a colour contour on a plane of the semimajor axis of the migrating protoplanets and the radius of the planetesimals, The analytically predicted inner and outer boundaries of the sweet spot are indicated with solid and dashed lines, respectively. The red lines are for j = 2 and the green lines are for j = 3. Right: total mass of captured planetesimals as a function of the planetesimal radius. 

In the text 
Fig. 6 Same as Fig. 5, but for the results of the parameter study regarding the migration timescale τ_{tide,0}. The grey circles in the right panel show the results of the parameter study regarding the size of planetesimals, R_{pl}, for comparison using the fraction of timescales τ_{aero,0}∕τ_{tide,a} on the vertical axis. 

In the text 
Fig. 7 Change in the phase angle of a planetesimal initially located at a_{0} = 13.8 AU. Panel a shows the case of τ_{tide,0} = 1 × 10^{4} yr and panel b shows the case of τ_{tide,0} = 1 × 10^{5} yr. In panel a, the phase angle oscillates indicating that the resonant trapping is unstable due to fast planetary migration. Then, the resonant trapping is broken earlier in the case of panel a than that shown in panel b. 

In the text 
Fig. 8 Same as Fig. 5, but for the results of the parameter study when varying the protoplanet’s mass, M_{p}. The red, green, sky blue, and blue lines represent the SSP inner (solid lines) and outer (dashed lines) boundaries determined by the resonant trapping, which is predicted by Eq. (32) for j = 2, 3, 4, and 5, respectively. The coloured lines change into the grey lines once the protoplanet mass exceeds a threshold value, beyond which the accretion bands disappear (see Fig. 3). 

In the text 
Fig. 9 Maximum mass of planetesimals that can be shepherded in mean motion resonances without effective mutual collisions. We show the cases of R_{pl} = 10^{7} cm and <i> = 10^{−3}. The solid, dashed, and dotted lines correspond to τ_{tide,a} = 10^{6} yr, 10^{5} yr, and 10^{4} yr, respectively. 

In the text 
Fig. 10 Same as Fig. 7, but (a) for the case with artificial perturbation given by Eq. (50) and (b) without artificial perturbation. Here, we use our reference model with R_{pl} = 1 × 10^{7}cm, τ_{tide,0} = 1 × 10^{5} yr and M_{p}∕M_{s} = 1 × 10^{−3}. 

In the text 
Fig. 11 Same as Fig. 5, but for the results of a parameter study for planetary mass M_{p} with the artificial perturbations (see Eq. (50)). In the right panel the results of a parameter study for planetary mass M_{p} without the artificial perturbations are plotted as grey circles for comparison. 

In the text 
Fig. 12 Distribution of confirmed exoplanets in the semimajor axis vs. planetary mass plane (http://exoplanet.eu). The black points are exoplanets observed by the radialvelocity method. The red points are exoplanets observed by the radialvelocity + transit methods. As for masses of exoplanets indicated with the black points, the observed values of M_{p} sin i are used. The grey area shows the theoretical sweet spot for planetesimal capture given by Eqs. (61) and (62) with C = 10. The red dashed line shows the theoretical sweet spot for planetesimal capture with C = 10^{3}, which corresponds to the case where collisional cascade have broken planetesimals down to R_{pl} ≲ 10^{5} cm. 

In the text 
Fig. 13 Maximum mass of planetesimals shepherded into the SSP without active planetesimal collisions M_{max,shep} as a functionof planetary mass M_{p}. Here, thecase is considered where the protoplanet migrates with the type II regime obtained by Kanagawa et al. (2018) (see Eq. (51)). The solid and dashed lines show the cases of <i> = 10^{−3} and 10^{−2}, respectively. The black and red lines are cases of massive disk with Σ_{gas,0} = 1400 g cm^{−2} and 800 g cm^{−2}, respectively. The parameters used in this plot are the optically thick disk of T_{0} = 140 K, the disk viscosity α_{vis} = 10^{−3}, the size of planetesimal R_{pl} = 10^{7} cm, and the planetesimal density ρ_{pl} = 1 g cm^{−3}. The scaling parameter of the location of the SSP C is obtained as C = 17. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.