Forming Giant Planets Around Late-M Dwarfs: Pebble Accretion and Planet-Planet Collision

We propose a pebble-driven core accretion scenario to explain the formation of giant planets around the late-M dwarfs of M ⋆ = 0 . 1 − 0 . 2 M ⊙ . In order to explore the optimal disk conditions for giant planet, we perform N-body simulations to investigate the growth and dynamical evolution of both single and multiple protoplanets in the disks with both inner viscously heated and outer stellar irradiated regions. The initial masses of the protoplanets are either assumed to be equal to 0 . 01 M ⊕ or calculated based on the formula derived from streaming instability simulations. Our ﬁndings indicate that massive planets are more likely to form in disks with longer lifetimes, higher solid masses, moderate to high levels of disk turbulence, and larger initial masses of protoplanets. In the single protoplanet growth cases, the highest planet core mass that can be reached is generally lower than the threshold necessary to trigger rapid gas accretion, which impedes the formation of giant planets. Nonetheless, in multi-protoplanet cases, the cores can exceed the pebble isolation mass barrier aided by frequent planet-planet collisions. This consequently speeds up their gas accretion and promotes giant planet formation, making the optimal parameter space to grow giant planets substantially wider. Taken together, our results suggest that even around very low-mass stellar hosts, the giant planets with orbital periods of . 100 days are still likely to form when lunar-mass protoplanets ﬁrst emerge from planetesimal accretion and then grow rapidly by a combination of pebble accretion and planet-planet collisions in disks with a high supply of pebble reservoir > 50 M ⊕ and turbulent level of α t ∼ 10 − 3 − 10 − 2 .


Introduction
Late M dwarfs are the end tail of low-mass stars with a typical stellar mass of ≈0.1−0.2M ⊙ .Gas-dominated giant planets, on the other hand, represent the upper branch of the planetary populations.The giant planet systems around late M dwarfs serve as a useful benchmark for understanding the planet assembling processes, thereby providing important constraints on planet formation models under extreme environments.Studying how giant planets form around late M-dwarfs is thus a matter of great interest.
The occurrence rate of giant planets (η J ) with a mass higher than 30 M ⊕ has been observed to correlate with stellar mass (Johnson et al. 2010).Exoplanet surveys have found that η J around early M dwarfs is approximately 3−5% for cold giant planets (Johnson et al. 2010;Bonfils et al. 2013;Suzuki et al. 2016;Sabotta et al. 2021) and 0.27% for hot Jupiters (Gan et al. 2023b).These values are even lower around mid-to-late dwarfs.For instance, Pass et al. (2023) indicated an upper occurrence rate of 1.5% for planets more massive than Jupiter at the locations out to the water-ice line.Bryant et al. (2023) inferred the hot Jupiter occurrence rate of 0.14% around stars less massive than 0.26 M ⊙ .
Despite the intrinsically low η J , a few such systems have been discovered around late dwarfs.For instance, the Transiting Exoplanet Survey Satellite (TESS) has confirmed one young, warm giant planet TOI-1227 b with a maximum mass of M p ≈0.5 M J orbiting around a very-low-mass star of 0.17 M ⊙ (Mann et al. 2022).In addition, three gas giants have been detected by radial velocity surveys: GJ 3512 b and c (Morales et al. 2019;Lopez-Santiago et al. 2020;Ribas et al. 2023) and GJ 9066 c (Feng et al. 2020;Quirrenbach et al. 2022), all of which have masses comparable to that of Saturn.Meanwhile, microlensing observations contribute a high number of cold, massive planets/brown dwarfs around these stellar objects.However, the planet's masses and orbital properties are not well-constrained in most circumstances (Suzuki et al. 2016;Zang et al. 2023).
The lack of giant planets around low-mass stars can be naturally attributed to the scarcity of solid material in their protoplanetary disks.Observations of nearby star-forming regions at A&A proofs: manuscript no.aanda (sub)millimeter wavelengths have shown a steeper-than-linear correlation between the solid disk masses and stellar masses (Pascucci et al. 2016;Ansdell et al. 2017), although with a huge intrinsic scatter (Manara et al. 2023).This correlation implies a significant shortage of building blocks for planet formation around very low-mass stars.To form giant planets around such systems, an extremely high conversion efficiency from dust to planets would be required.
Several mechanisms have been proposed for giant planet formation.The disk instability theory suggests that giant planets form directly through the gravitational fragmentation of young, massive protoplanetary disks (Boss 1997;Stamatellos & Whitworth 2009;Deng et al. 2021;Boss & Kanodia 2023).Mercer & Stamatellos (2020) explored the conditions for giant planet formation around stars with masses of 0.2−0.4M ⊙ , and found that a super-massive disk (>30% of the host mass) is necessary to yield disk fragments, which typical results in planets with masses several times that of Jupiter at a few tens of au.Morales et al. (2019) reported that gas giant GJ 3512 b could form through disk instability, with a similar disk mass requirement.However, in contrast to observations indicating an increase in η J with higher stellar metallicity (Santos et al. 2004;Fischer & Valenti 2005;Gan et al. 2022Gan et al. , 2023a)), the giant planets formed through gravitational instability appear to be largely unaffected by the disk metallicity (Boss 2002;Cai et al. 2006;Mercer & Stamatellos 2020).
The core accretion theory suggests that the formation of giant planets requires the growth of massive cores (∼10 M ⊕ ) to initiate rapid gas accretion (Pollack et al. 1996) before the dissipation of disk gas.Two channels for core accretion have been proposed: planetesimal accretion and pebble accretion.In the classical planetesimal-driven core accretion scenario (Ida & Lin 2004;Ogihara & Ida 2009;Zhang & Ji 2009;Mordasini et al. 2012;Coleman & Nelson 2016), protoplanet growth occurs by accreting surrounding planetesimals with characteristic sizes of tens to hundreds of kilometers.Miguel et al. (2020) found that only Earth-analog systems could form around stars of 0.1−0.2M ⊙ with sufficiently massive disks.Burn et al. (2021) used subkilometer-sized planetesimals and explored the formation of giant planets by adjusting the initial planetesimal surface density and the type I migration rate.They concluded that giant planets can only form in massive disks (gas mass >0.007 M ⊙ , solid mass >66 M ⊕ ) and when the type I migration rate is reduced by an order of magnitude (see, e.g., Ogihara et al. (2015Ogihara et al. ( , 2018))).A follow-up study from Schlecker et al. (2022) indicated that with the current planetesimal accretion model, it is difficult to reproduce the observed giant planet population around stars less massive than 0.5 M ⊙ .
In the pebble-driven core accretion scenario (Ormel & Klahr 2010;Lambrechts & Johansen 2012, 2014a;Ida et al. 2016;Ormel 2017a;Liu et al. 2019aLiu et al. , 2020;;Venturini et al. 2020;Chachan & Lee 2023), planets accrete millimeter-to-centimetersized pebbles to increase their masses.Compared to planetesimal accretion, the accretion cross-section of pebbles is largely enhanced by aerodynamic gas drag (Ormel 2017a;Johansen & Lambrechts 2017;Liu & Ji 2020).Ormel et al. (2017b) and Schoonenberg et al. (2019) proposed a pebbledriven planet formation model for the TRAPPIST-1 system, and the resulting low water content and characteristic masses of all seven planets are in good agreement with observations.Coleman et al. (2019) examined both planetesimal and pebble accretion modes for this peculiar system and noted that diverse outcomes may arise by considering pebble ablation and planet atmosphere recycling.Notably, in the context of pebble accretion scenarios, the core masses of the planets are limited by the pebble isolation mass (Lambrechts et al. 2014b), which is defined as the point when the planets reach a mass that opens a shallow gap in the protoplanetary disk and truncates the drifting pebbles.Liu et al. (2019aLiu et al. ( , 2020) ) showed that planets around 0.1 M ⊙ can only reach a maximum mass of 2−3 M ⊕ due to the low pebble isolation mass around these stellar hosts.Moreover, the characteristic masses of forming planets clearly exhibit a dependence on stellar metallicity (Liu et al. 2019a).
On the other hand, planets approaching such isolation masses can undergo substantial orbital migration.The strength and direction of migration are determined by the disk properties (Paardekooper et al. 2011).Planets at different disk radii can migrate convergently towards and become trapped in mean motion resonances (Wang & Ji 2017;Pan et al. 2022) at some special disk locations with zero net torque (also known as the transition radius, see Lyra et al. (2010); Horn et al. (2012); Kretke & Lin (2012); Liu et al. (2015)).Such planet migrationinduced mass concentration is likely to trigger dynamical instability with frequent orbital crossings and close encounters (Zhang et al. 2014).Massive cores can be attained through mutual planet-planet collisions in the gas-rich disk phase, promoting subsequent rapid gas accretion.This planet-planet collision driven core accretion scenario has been mainly explored around solar-type stars (Liu et al. 2015;Wimarsson et al. 2020) within a limited stellar mass range (Liu et al. 2016).Noticeably, inferred from the Juno measurement, Jupiter features a dilute core with an extended heavy element layer (Wahl et al. 2017;Helled & Stevenson 2017).This pattern could be explained by the giant impacts among protoplanets during their final assembling phase (Liu et al. 2019c).
In light of the literature studies presented, we speculate that the rapid growth of planetary cores can be achieved by a combination of the above two growth models.In this regard, we propose a hybrid growth model that utilizes both pebble accretion and planet-planet collisions to explain the formation of giant planets around late dwarfs.In this new scenario, the pebble isolation mass is not a barrier that limits core accretion, unlike in previous single protoplanet growth models (Liu et al. 2019a(Liu et al. , 2020)).Pebble accretion plays a key role in the growth of individual protoplanets.As these protoplanets reach certain masses and migrate towards the transition radius, planet-planet collisions take over to assemble massive cores that can transition to runaway gas accretion.We evaluate the above hypothesis in this paper.
The paper is structured as follows.The model setup and Nbody implementation are described in Section 2. The growth of individual protoplanets with different disk and planet parameters is examined in Section 3. Section 4 explores the formation and evolution of multiple protoplanets, and Section 5 presents an assessment of the model and its implications.Finally, we summarize the key results in Section 6.

Method
We adopt the planet formation model from Liu et al. (2019a).We summarize the key physical processes and main equations here.A detailed and complete model description is referred to Section 2 of Liu et al. (2019a).a quasi-steady manner.The disk is divided into two distinct components based on different heating mechanisms (Garaud & Lin 2007).The inner optically thick disk region is viscously heated (Ruden & Lin 1986), in which the gas surface density, temperature and disk aspect ratio are given by (1) and Ṁg , M ⋆ , and r are the disk accretion rate, stellar mass, and the distance to the central star, respectively.At 1 au, the initial gas surface density is approximately 100 g cm −2 .Given the low surface densities, the magneto-rotational instability (MRI) might be driven by cosmic-ray ionizationGammie (1996).In this paper, we focus on late M dwarfs; therefore, the stellar mass is specified as M ⋆ =0.1 M ⊙ unless otherwise explored (see Section 5.5).We assume that the disk opacity is κ=κ 0 T g /1 K g cm −2 (Garaud & Lin 2007), where κ 0 = 0.01 is the opacity coefficient, and α g represents the global disk angular momentum transport efficiency Shakura & Sunyaev (1973).We choose α g =10 −2 , inferred from disk observations of Ṁg and Σ g (Hartmann et al. 1998;Andrews et al. 2009).
The outer disk is assumed to be optically thin in the vertical direction and primarily heated by stellar irradiation (Chiang & Goldreich 1997), in which the gas surface density, temperature and disk aspect ratio are given by (Ida et al. 2016) where L ⋆ is the stellar luminosity.We use f s = 1/(1 + r tran ) 4 as a smooth function to combine the inner and outer regions, therefore the global disk quantity can be calculated by The transition radius between the inner and outer regions is given by which decreases as disk dissipation.
Initially the young disks can maintain a continuous supply of infall of material from their parent molecular clouds (Padoan et al. 2014).At later times the infall is quenched and the disks gradually deplete gas by combined effects of viscous accretion and stellar photoevaporation (Hartmann et al. 1998;Alexander et al. 2014;Ercolano et al. 2018).In this work we simply assume that the disk accretion rate remains a constant Ṁg = Ṁg,0 at t≤t 0 , and follows an exponential decay Ṁg = Ṁg,0 exp [−(t − t 0 )/τ dep ] at t>t 0 where t 0 separates the early infall and later dissipation stages, and τ dep is the disk depletion timescale.The total gas disk mass is therefore parameterizedly described by Ṁg,0 , t 0 and τ dep .
Observations of disk accretion rate onto very low-mass stars show a wide spread, ranging from <10 −9 M ⊙ yr −1 up to 1−2 × 10 −8 M ⊙ yr −1 (Hartmann et al. 2016;Pinilla et al. 2021).In this paper we focus on giant planet formation around late dwarfs.The paucity of massive planets around these stars indicates such systems are expected to grow only in relatively massive disks.Hence, in the fiducial model for studying the systems around stars of 0.1 M ⊙ , we choose a relatively high initial disk accretion rate of Ṁg,0 =10 −8 M ⊙ yr −1 , and the disk starts to dissipate at t 0 =1 Myr on a timescale τ dep of 0.5 Myr.Here we define the disk lifetime t disk as the time when Σ g at 1 au drops below 1 g cm −2 .In this circumstance, the initial disk mass is 15% of its stellar mass and t disk =3.7 Myr.
On the other hand, the disk lifetime is inferred to be longer around lower-mass stars (Williams & Cieza 2011;Bayo et al. 2012;Manara et al. 2012;Ribas et al. 2015;Picogna et al. 2021).In order to investigate the influence of disk lifetime on planet formation, we construct a disk with the same initial mass as the fiducial one but vary Ṁg,0 =6 × 10 −9 M ⊙ yr −1 , t 0 =1.5 Myr and τ dep = 1 Myr.In such a circumstance t disk =6.3 Myr.

Initial mass of protoplanet
We start the growth of protoplanet with an initial mass M p0 .There are two considerations regarding the choice of M p0 .First, a canonical value of 0.01 M ⊕ is widely adopted in literature, which can date back to the pioneering numerical N-body simulation study of Kokubo & Ida (1998).They found that a few lunar-mass oligarchs naturally emerge out from a swarm of small planetesimals by mutual collisions.Their study is limited to the circumstances of host stars with a solar mass.The formation of embryos with 0.01−0.1 M ⊕ around M dwarfs after oligarchic growth has also been suggested by Ogihara & Ida (2009) (note that their study was conducted under the assumption of the solar nebular conditions).No further extended numerical work has been conducted to explore how the forming masses of protoplanets around lower-mass dwarfs.Ormel et al. (2010) analytically derived the transition mass between runaway and oligarchic growth such that M 0 ∝M −3/7 ⋆ Σ 6/7 plt R 9/7 plt (their Eq.13), where Σ plt and R plt are the surface density and size of planetesimals.The latter two quantities (Σ plt and R plt ) should be also dependent on the host stellar environment.Therefore, without any further assumptions on the planetesimal formation models, the exact stellar mass dependency of M 0 is unknown.Bearing these uncertainties, we follow similar literature studies (Liu et al. 2019a;Burn et al. 2021) and adopt M p0 =0.01 M ⊕ as a prior in the follow-up explorations.This constant mass assumption can serve as one benchmark, which differs from the secondary con-sideration that specifically assumed one planetesimal formation model.We also note that planets with masses that exceed this value are well in the settling pebble accretion regime, where the planet-pebble interaction is substantially aided by gas drag Ormel & Klahr (2010); Liu & Ormel (2018); Liu et al. (2019b).
On the other hand, streaming instability provides a valuable pathway for the formation of planetesimals (Youdin & Goodman 2005;Johansen & Youdin 2007).It occurs when the volume density of pebbles approaches that of the gas, resulting in significant back-reaction from pebbles onto the gas.As a result, these pebbles concentrate radially into dense clumps and eventually collapse into planetesimals through self-gravity.Defining the protoplanet as the largest planetesimal generated from the streaming instability clumps, Liu et al. (2020) obtain the mass of the protoplanet from the extrapolation of literature numerical investigations (e.g., Johansen et al. (2015); Simon et al. (2016); Schäfer et al. (2017); Abod et al. (2019)).This streaming instability-induced protoplanet mass can be expressed as (see section 2.4 of Liu et al. (2020) for derivations) Where γ=4πGρ g /Ω 2 K is the relative strength between self-gravity and tidal shear, ρ g =Σ g /( √ 2πh g r) is the gas volume density and To summarize, we assume two scenarios for the starting mass of a protoplanet.In the equal-mass scenario, protoplanets form from classical planetesimal accretion (Kokubo & Ida 1998).We assume they all have M p0 =0.01 M ⊕ .In the second scenario the protoplanets are specifically generated by streaming instability, and their birth masses follow Eq. ( 8).In contrast to the equalmass scenario and as can be seen in Eq. ( 8), the mass of protoplanets formed by streaming instability correlates with gas disk density, gas disk aspect ratio and stellar mass.In this respect, we expect a higher M p0 at a larger orbital distance since both γ and h g increase with r.

Pebble accretion
Pebbles undergo fast radial drift towards the central star.A fraction of these drifting pebbles can be accreted by the planet when they cross the planetary orbit.The pebble accretion rate onto the planet's core is given by where Ṁpeb is the pebble mass flux and ε PA is the total pebble accretion efficiency, the formulas of which are adopted from Liu & Ormel (2018) and Ormel & Liu (2018) that include both 2D and 3D accretion efficiencies (ε PA,2D and ε PA,3D ) taking into account the eccentricity and inclination of the planet.In brief, the pebble accretion efficiency firstly gets boosted when the planets have relatively low eccentricity.It then drops with the further increase of eccentricity due to the fact that high pebble-planet impact is not in the settle regime anymore.On the other hand, the pebble accretion efficiency decreases with inclination since the planets are more likely to lift off the pebble plane when they are on inclined orbits.
In the limit of zero eccentricities and inclinations, the pebble accretion efficiency in the settling regime can be approximated as where v K is the Keplerian velocity, ∆v is the relative velocity between the pebbles and planet (dominated by ηv K in the headwind regime, Ω K R H in the shear regime, and R H =(M p /3M ⋆ ) 3 r is the planet Hill radius), h peb is the pebble disk aspect ratio, η = −h 2 g (∂ ln P/∂ ln r)/2 and P is the gas disk pressure.It is important to recognize that both ε PA,2D and ε PA,3D increase as h g (or equivalently η) decreases.Physically, the pebble accretion efficiency becomes higher when the inward drifting pebbles is slower in the 2D regime and/or the pebble disk is less vertically extended in the 3D regime.
The pebble disk scale height H peb = √ α t /(α t + τ s ) H g (Youdin & Lithwick 2007), where α t is the turbulent diffusion coefficient, approximately equivalent to the local turbulent viscous parameter when the disk is driven by magneto-rotational instability (Johansen & Klahr 2005;Zhu et al. 2015).Physically, α t can differ from α g -the average value of the global disk angular momentum transport efficiency -due to instances of layered accretion (Turner & Sano 2008).The midplane of the disk is quiescent and the high altitude region is turbulent active.We note that α t is more relevant to the midplane of the dead zone while α g represents the vertically average, global disk angular momentum transport efficiency.These two parameter are not always equal.The planet gap opening occurs at the disk midplane, and this process can also be closed by local turbulent diffusion.Hence, α t is more relevant to the planet formation processes such as dust stirring, pebble accretion and gap opening (Xu et al. 2017).
Meanwhile, τ s is the pebble's dimensionless stopping time (termed Stokes number hereafter) that characterizes the aerodynamic size of pebbles.The detailed dust evolution is not modeled here.Advanced dust coagulation studies find that the largest pebbles dominate the total mass of the population and their Stokes number is almost a constant (e.g., in the fragmentation-limited regime).For the sake of simplicity, we ideally treat that all pebbles reach a fixed Stokes number of τ s =0.05.The potential influence of α t on τ s is discussed in Section 5.2.
The pebbles are assumed to be constituted of 50% water ice and 50% silicate.The water-ice line r H 2 O is calculated when the disk temperature is 170 K.When the pebbles drift inside of the water-ice line, their icy component sublimates, and the pebble mass flux decreases accordingly.We neglect the pebbles' Stokes number variation when they cross r H2O .
Same as Liu et al. (2019aLiu et al. ( , 2020)), we assume that the pebble and gas flux ratio remains a constant such that ξ p/g = Ṁpeb / Ṁg .Pebbles are well-coupled to the disk gas when their Stokes number is very low.Thus, pebbles and gas drift at the same speed and the initial disk metallicity is preserved, where disk metallicity Z=Σ peb /Σ g .When the pebbles have a higher Stokes number, they drift faster than disk gas.In this case, in order to maintain a constant flux ratio, Σ peb /Σ g becomes lower than the initial disk metallicity.We assume ξ p/g =0.01 in the fiducial model, corresponding to totally 50 M ⊕ solid in pebbles.It is worth noting that the above constant mass flux ratio is a global concept.The disk metallicity can still be enriched at local places due to various mechanisms.For instance, several studies proposed that the local solid density can be enhanced at the waterice line (Ros & Johansen 2013;Schoonenberg & Ormel 2017;Dr ążkowska & Alibert 2017).We do not take these localized effects into account in this study.The enrichment of disk metallicity in the late gas disk dispersal phase is also not considered.
As the planet grows, it becomes massive enough to perturb the surrounding gas and produce a local pressure bump.The inward drifting pebbles stop at the outer edge of the gap generated by the planet.As such, the planet cannot further accrete pebbles.This onset planet mass is defined as the pebble isolation mass (Lambrechts et al. 2014b).On the other hand, the gap opening mass is typically defined when the planet opens a gap whose surface density drops by 50%.We adopt the gap opening mass based on Kanagawa et al. (2015)'s 2D hydrodynamical simulations, Based on the 1D numerical simulations conducted by Johansen et al. (2019), the pebble isolation mass is approximately 2.3 times lower than the gap opening mass.We apply this scaling and convert Kanagawa et al. (2015)'s gap opening mass into pebble isolation mass, which reads We also demonstrate a comparison of M iso adopted in this work and other literature studies (Ataiee et al. 2018;Bitsch et al. 2018) in Appendix A. Same as Liu et al. (2019b) and Jang et al. ( 2022), we consider the filtering of flux when pebbles drift through different planets in a multi-planetary system (see Eq.12 of Liu et al. (2019b)).That means the pebble mass flux entering the inner disk can be reduced due to the accretion of planets in the outer disk region.We simplified that the pebble accretion of the planets that reside in the interior of its orbit is terminated when the planet reaches M iso .The diffusion of small, fragmented particles through the gap is not considered (Liu et al. 2022;Stammler et al. 2023).

Gas accretion
Gas accretion can be divided into the early hydrostatic phase and later runaway phase (Pollack et al. 1996).During hydrostatic accretion, the planet slowly captures disk gas to form a tiny atmosphere.The envelope hydrostatic equilibrium is established when the gravitational energy is balanced by radiative heating.The heat from solid accretion is quenched when the planet reaches the pebble isolation mass.The envelope is expected to undergo subsequent Kelvin-Helmholtz contraction.For simplicity, we ignore gas accretion in the hydrostatic phase and treat M iso as the onset mass for gas accretion (Ogihara & Hori 2020).We note here that the literature critical core mass is estimated to be 5−15 M ⊕ (Ida & Lin 2004;Alibert & Venturini 2019), higher than M iso (typically 1−2 M ⊕ ) around very low-mass stars.Our simplification remains justified as the gravitational force of planets with isolation mass is insufficient to retain a substantial atmospheric envelope (Alibert & Venturini 2019).
The gas accretion rate in the Kelvin-Helmholtz contraction reads (Ikoma et al. 2000) where κ env is the envelope opacity, a crucial parameter that sets the amount of gas accreted by the planet.A variety of κ env have been tested and it has been found that a very low κ env ≪0.1 cm 2 /g might lead overpopulated massive giant planets, contradicting with observations.On the other hand, the disk opacity is estimated to be ∼1 cm 2 /g close to and beyond r H 2 O based on an ISM-like dust size distribution (Bell & Lin 1994).The envelope opacity is expected to be no higher than the disk opacity.This is because when the planet reaches M iso , large pebbles get completely blocked, whereas only small dust well coupled to the gas can drift across the gap and get accreted onto the planet (Liu et al. 2022;Stammler et al. 2023).This dustsize filtration lowers the opacity in the planet envelope compared to the disk gas.In addition, the envelope opacity could be further reduced by grain sedimentation, coagulation and evaporation (Movshovitz et al. 2010;Ormel 2014;Mordasini et al. 2014).Considering the above reasons, we adopt a moderate κ env =0.1 cm 2 /g and assume it does not vary with the disk metallicity (see Fig. 8 of Mordasini et al. (2014)).
The gas accretion in Equation ( 13) decreases dramatically with the lowering of the planet mass.For instance, Ṁp,g ∼1.5 × 10 −7 M ⊕ yr −1 at M p =2 M ⊕ , indicating the gas contraction is very limited over the disk lifetime in this circumstance.However, Ṁp,g ∼4 × 10 −5 M ⊕ yr −1 at M p = 8 M ⊕ and the planet double its mass within a few 10 5 yr.The mass of planetary core plays a critical role in the gas accretion.It significantly affects the accretion rate and the total amount of gas that the planet accumulates before disk dissipation.
Besides, only a fraction of gas within the planet Hill sphere can be accreted (Tanigawa & Watanabe 2002;Machida et al. 2010).We adopt the corresponding accretion rate from Eq. 29 of Liu et al. (2019a): The further gas accretion onto the planet is restricted to the gas flux in the protoplanetary disk.In sum, the total gas accretion rate can be expressed as Ṁp,g = min dM p,g dt KH , dM p,g dt Hill , Ṁg . (15)

Planet migration
Planets embedded in disks exchange angular momentum with the surrounding gas, resulting in their orbital migration, eccentricity and inclination damping.We adopt a combined torque formula including both type I and type II regimes (Kanagawa et al. 2018): where g is the normalized torque strength, f I and f II are the type I and type II migration prefactors.The type II migration coefficient f II = − 1 whereas the type I migration coefficient f I is set by the disk thermal structure and local turbulent α t (see Paardekooper et al. (2011) for details).Differing from the traditional criterion within the viscous accretion disk framework (Lin & Papaloizou 1986;Rafikov 2002;Tanaka et al. 2002), we employ α t instead of α g when addressing the gap opening and migration, motivated by the magnetohydrodynamic effect in the wind-driven disk.(Aoyama & Bai 2023).A smooth function of ] is chosen to avoid discontinuity and ensures that Γ≈Γ I when M p ≪M gap and Γ≈Γ I /(M p /M gap ) 2 when M p ≫M gap (Kanagawa et al. 2018).
We note that the heating torque from gas and pebble accretion (Benítez-Llambay et al. 2015;Masset 2017;Cornejo et al. 2023) as well as other potential planet traps at the opacity transition regions such as ice-lines (Kretke & Lin 2012) are not taken into account in our study.
Planets in the inner viscously heated region can undergo outward migration ( f I > 0) when their masses are comparable to the optimal mass, which reads M opt = 0.9 A notable feature is that even though the protoplanets are distributed widely over the whole disk region, they would migrate convergently towards r tran when reaching such a mass.The inner disk is truncated by the stellar magnetospheric torque (Lin et al. 1996;Liu et al. 2017) and the corresponding cavity radius is around 0.03 au around young T Tauri stars around solar-mass.We set the inner disk boundary as r in =0.01 au around late dwarfs which is assumed to equal to their stellarcorotation radius with spin orbits of ≈3 days.Any planets migrating interior to this radius are immediately stopped.In numerical integrations, we remove these planets inside the cavity radius to save the computational cost.

Numerical setup
We used numerical N-body simulations to study the growth and evolution of multi-protoplanets.We have employed the MER-CURY code (Chambers 1999) with the Bulirsch-Stoer integrator.In the code planet-planet collisions are treated as inelastic mergers with conserved angular momentum when the separation of two planets is smaller than the sum of their physical radii.We ignore the influence of the potential energy released during the giant impact on the cooling and gas accretion.The mass of the remnant planet is a sum of both impactor and target.Fragmentation and restitution (Leinhardt & Stewart 2012;Mustill et al. 2018) are not considered in this work (see further discussions in Sect.4.1).A planet with a distance greater than 100 au from the central star is considered to be ejected from the planetary system.
The planet-disk interactions are implemented as accelerations: where v is the velocity vector.Modifying from Cresswell & Nelson (2008), the migration, eccentricity, and inclination damping timescales are given by where In short, the modified version of the code can handle planetplanet interactions and collisions and additionally account for the effects of planet mass growth by pebble accretion, planetgas disk interaction torques, and the corresponding eccentricities/inclinations damping.

-2
10 -1 10 0 10 1 10 2 Semi-major axis (au) This is when the fastest growing planet reaches its isolation mass.The arrow indicates the disk transition radius, and the increasing sizes of the dots denote the disk evolution at one Myr intervals.The planet attains the highest mass at a moderate birth radial distance close to the transition radius.The planet and disk parameters are referred to Table 1.

Growth of a single protoplanet
In this section we explore the growth and migration of a single protoplanet around a star of M ⋆ =0.1 M ⊙ .In our fiducial run we assume the initial mass of the protoplanet to be 0.01 M ⊕ and the local turbulent diffusivity coefficient α t =10 −3 .The initial disk accretion rate is chosen as Ṁg,0 =10 −8 M ⊙ yr −1 , and it starts to dissipate at t 0 =1 Myr with a dispersal timescale τ dep of 0.5 Myr.This corresponds to a disk lifetime t disk =3.7 Myr.The initial pebble flux is Ṁpeb,0 =3.3 × 10 −5 M ⊕ yr −1 , equivalently to ξ p/g = Ṁpeb / Ṁg =1%.
We also investigate the influence of α t in Sect.3.2, pebbleto-gas mass flux ratio in Sect.3.3, disk lifetime in Sect.3.4 and initial mass of protoplanet in Sect.3.5, respectively.The setup of disk and protoplanet parameters are listed in Table 1.

Fiducial case
Figure 1 illustrates the growth of individual protoplanets at various birth locations r 0 .The dashed line refers to M iso at the time when the fastest growing protoplanet reaches (t=1.3Myr), and the vertical arrow indicates the transition radius between two disk heating sources.The increasing size of the dots represents the time evolution, with intervals of a Myr.
The red curve in Figure 1 depicts the growth of a protoplanet at r 0 =1 au.The mass of the protoplanet increases by two orders of magnitudes through pebble accretion during the first Myr.As the planet reaches M opt ∼0.5M ⊕ , it starts to migrate outward to r tran due to a strong, positive corotation torque in the viscously heated disk region (Liu et al. 2019a).However, this corotation torque gradually diminishes as disk dissipating and the planet mass further increasing, causing rapid inward migration.The planet reaches M iso =1.8 M ⊕ at t=1.3 Myr and r=1.2 au.Because of a relatively low core mass, it since then only accretes a limited Notes.t disk is calculated for the timespan when the gas surface density at 1 au drops to 1 g cm −2 .
amount of gas and eventually grows into a close-in, super-Earth planet of M p =2.4 M ⊕ .
The growth differs when the protoplanets are born at different r 0 .Only protoplanets with r 0 ∼r tran can grow sufficiently massive and undergo large-scale radial migration.We term this region that protoplanets can grow beyond 0.5 M ⊕ as the efficient planet growth region.Protoplanets with initial closer-in and furtherout orbits end up as lower-mass planets.This r-dependent mass growth correlates with the gas disk scale height, which governs the efficiency of pebble accretion.In the 2D accretion case, a larger gas scale height indicates a faster headwind speed.The pebbles drift too fast and are less likely to be accreted by the planets.On the other hand, in the 3D accretion case, a larger gas scale height also means more vertically extended pebble layers, leading them less efficient to be attracted by the planet.Taken together, the highest efficient pebble accretion occurs when the disk scale height has the lowest value.This corresponds to the disk location at r tran , since h vis ∝r −1/16 in the inner disk and h irr ∝r 2/7 in the outer disk (see Eqs. 3 and 6).As a result, the planets exhibit a peak growth rate at r∼r tran .However, none of these protoplanets finally grow into massive, gas-dominated planets, due to the fact that their core masses are too low to initiate rapid gas accretion.

Disk turbulence
It is expected that disks are turbulent, which dynamically stirs up solid particles and affects the pebble accretion efficiency of planets (Johansen & Lambrechts 2017), as well as the planetesimal formation through the streaming instability (Johansen et al. 2014;Dr ążkowska et al. 2023).One direct method to assess the strength of the turbulence is by deriving the turbulence-induced broadening observed in molecular line emissions (Najita et al. 1996).Instead of a universal turbulent viscosity, the value of α t varies from one disk to another (Flaherty et al. 2018(Flaherty et al. , 2020;;Teague et al. 2018).The inner disk region and upper layer of the source SVS 13 shows supersonic turbulence (Carr et al. 2004), while HD 163296 demonstrates moderate levels of turbulence with α t < 2.5 × 10 −3 (Flaherty et al. 2015(Flaherty et al. , 2017)).
Another approach to constrain turbulence is through geometric considerations (Rosotti 2023), for example the dust vertical extent or the radial width of disks influenced by settling/radial drift and turbulence diffusion (Whipple 1972;Pinte et al. 2016;Rosotti et al. 2020).Observations of Oph 163131 (Villenave et al. 2022) and the DSHARP survey (Andrews et al. 2018;Dullemond et al. 2018) suggest a preference for low turbulent viscosity of α t 10 −4 to moderate values of α t ∼ 10 −3 .Overall, turbulent viscosity typically ranges from α t = 10 −2 to 10 −4 in different disks, leading us to investigate the influence of disk turbulence on planet growth and migration within this parameter space.
We maintain a constant global disk angular momentum transport efficiency α g , and the results are depicted in Figure 2. In highly turbulent disks with α t =10 −2 , pebbles are vertically extended over the gas scale height, leading to a suppression of pebble accretion compared to the fiducial run.In addition, the pebble isolation mass increases with disk turbulence (Eq.12), making planets more challenging to reach M iso before disk dissipation.The maximum mass that a planet can attain is approximately Venus-mass in Fig. 2a.On the other hand, in weakly turbulent disks with α t =10 −4 , planet growth speeds up due to efficient pebble accretion.Planets can grow up to M iso within 1 Myr at r 0 ∼r tran .However, in such a case M opt is lower, and the effect of outward migration is insignificant.Planets migrate rapidly into the inner disk region.Moreover, M iso is also lower, and planets are prevented from accreting substantial gas to become gas giants.Only small planets with the highest mass of ∼1 M ⊕ form in the end.
In brief, the growth of massive planets from a single protoplanet is largely impeded, in the disks with either very high or very low turbulent levels.

Pebble-to-gas mass flux ratio
The total solid mass in disks, quantified by ξ p/g , the pebble-togas mass flux, crucially determines the planet growth timescale.A higher pebble mass flux facilitates the formation of massive planets.
The right panel of Figure 2 illustrates the growth of protoplanets in metal-rich disks with a pebble-to-gas flux ratio ξ p/g of 2%.The efficient planet growth region is wider, and protoplanets at further-out disk regions can grow more quickly and substantially in this situation compared to the case of the nominal ξ p/g =1% (left panel of Figure 2).For instance, super-Earth can form in disks with α t = 10 −3 at r 0 =10 au and α t =10 −4 disks at r 0 =20 au in Figure 2e and f, respectively, because of their considerable core masses.
Notably, protoplanets in highly turbulent disks experience even more significant mass growth due to the fact that less efficient pebble accretion is largely compensated by a large supply of the pebble reservoir (Figure 2d).In such a case, the optimal mass is higher, allowing them to retain at r tran for a longer time to proceed the mass growth.The pebble isolation mass is also higher.Once they reach such a massive core, they are more likely to initiate rapid gas accretion.Evidently, in Figure 2d protoplanets born at r<5 au reach M iso ∼3 M ⊕ at a relatively early time and eventually grow into Neptune-mass planets.
To conclude, in metal-rich disks, high disk turbulence may no longer pose a threat to the formation of massive planets.The Fig. 2: Growth and migration of single protoplanets born with lunar masses at different turbulent levels and pebble-to-gas mass flux ratios around stars of M ⋆ =0.1 M ⊙ .Three turbulent coefficients of α t = 10 −2 (upper), 10 −3 (middle) and 10 −4 (lower) and two pebble-to-gas flux ratios of ξ p/g =1% (left) and 2% (right) are shown.The planet and disk parameters are listed in Table 1.The dashed line represents the pebble isolation mass at the time when the planet born at 1 au reaches this value.Note that in panel (a) planets never approach the isolation mass.We instead adopt the isolation mass using the stellar irradiation model.The increasing sizes of the dots denote the disk evolution at one Myr intervals.Massive planets prefer to form in mental-rich and highly turbulent disks.
true barrier is M iso , which determines the ability of runaway gas accretion.

Disk lifetime
Giant planets assemble gas and solids within a finite protoplanetary disk lifetime.Therefore, the survival time of the gaseous disk is expected to play a decisive role.Here we keep the total gas disk mass the same as the fiducial run and investigate the influence of a longer disk lifetime on planet growth and migration.This long-lived disk is characterized by a lower initial disk accretion rate Ṁg,0 =6 ×10 −9 M ⊙ yr −1 and an extended disk dissipation such that t 0 =1.5 Myr and τ dep =1.0 Myr.The disk lifetime is 6.3 Myr in this case.Owing to the slow disk dissipation Ṁg becomes higher than the fiducial case after t=1.3 Myr.The results are presented in Figure 3. Since the pebble flux is attached to the gas flux, protoplanets grow slowly in the early stage, and they approach M iso at later times (generally later than 2 Myr) in Figure 3 compared to Figure 2. The disk mass also dissipates slower.After t=1.3 Myr, both M iso and Ṁpeb are higher in long-lived disks.As such, protoplanets speed up their growth at these advanced phases and reach a higher core mass.Meanwhile, M opt is also higher in disks with a relatively high Ṁg .Outward migration is also more profound in Figure 3d when the planet grows beyond a few Earth masses.
Importantly giant planet formation is significantly promoted in long-lived, highly turbulent, and metal-rich disks, as shown in Figure 3d.

Protoplanets formed by streaming instability
We assume that the protoplanets form by the streaming instability mechanism.Unlike previous circumstances where protoplanets had an equal lunar mass, the mass derived from Eq. 8 correlates with M ⋆ , Σ g , and h g , increasing with r within a range of 10 −4 M ⊕ to approximately 0.1 M ⊕ (Liu et al. 2020).We note that the streaming instability triggering condition also Fig. 3: Growth and migration of single protoplanets in disks of a long lifetime at different turbulent levels and pebble-to-gas mass flux ratios around stars of M ⋆ =0.1 M ⊙ .Three turbulent coefficients of α t = 10 −2 (upper), 10 −3 (middle) and 10 −4 (lower) and two pebble-to-gas flux ratios of ξ p/g =1% (left) and 2% (right) are shown.The planet and disk parameters are listed in Table 1.The dashed line represents the pebble isolation mass at the time when the planet born at 1 au reaches this value.The increasing sizes of the dots denote the disk evolution at one Myr intervals.Compared to Figure 2, more massive planets from in disk with a longer lifetime.
correlates with disk properties such as the local disk metallicity (Yang et al. 2017;Li & Youdin 2021).We assume that even though the global disk metallicity is below the threshold value, the local solid density can still be enhanced to fulfill the streaming instability by various of hydrodynamical and magnetic instabilities (see references in Lenz et al. (2019)).Following this, the mass of forming planetesimal is degenerate from the global ξ p/g .We perform simulations with disk parameters identical to the fiducial run, and the results are demonstrated in Fig. 4. Compared to Figure 2, protoplanets formed by the streaming instability face significant challenges in growing masses in highly turbulent disks.This situation holds both true for disks with ξ p/g =1% and 2% (Fig. 4a and d).In a disk with α t =10 −3 , the region of efficient planet growth spans approximately 2−5 au (Fig. 4e), narrower than that in equal-mass cases.The masses of the planets increase by two orders of magnitude at ξ p/g =1% (Fig. 4b), whereas the formation of Earth-sized planets becomes feasible at ξ p/g =2% (Fig. 4e).The outcome is natural to understand from the difference in initial protoplanet masses.The mass from streaming instability is generally lower than lunar mass within 20 au.Therefore, protoplanets with lower masses have lower gravitational potential to capture pebbles, resulting in a longer growth time.
Massive planets are more likely to grow in disks with low α t and high ξ p/g .We find that the general growth pattern is similar, but protoplanets located at the outer disk region in Fig. 4f can attain slightly higher masses than those in Fig. 2f since M p0 there are higher than 0.01 M ⊕ .But the growth is still limited and only cold, super-Earth planets form eventually.
In conclusion, the formation of massive planets is more difficult when protoplanets are assumed to form by the streaming instability rather than being born with an equal lunar mass.

Growth of multi-protoplanets
In the previous section we present the growth of a single protoplanet under various disk and planet parameters.However, multi-planetary systems are commonly observed.It is essential to understand how the formation and evolution of the planetary Fig. 4: Growth and migration of single protoplanets formed by streaming instability at different turbulent levels and pebble-to-gas mass flux ratios around stars of M ⋆ =0.1 M ⊙ .Three turbulent coefficients of α t = 10 −2 (upper), 10 −3 (middle) and 10 −4 (lower) and two pebble-to-gas flux ratios of ξ p/g =1% (left) and 2% (right) are shown.The planet and disk parameters are listed in Table 1.Compared to Figure 2, the growth of the planets is significantly impeded unless in disks with low turbulence and high pebble flux.system from more realistic configurations started from multiprotoplanets.
In order to investigate this, we conduct N-body numerical simulations that account for the gravitational interactions among multiple protoplanets.We test whether the growth pattern differs in single and multi-protoplanet cases.The illustrations of the Nbody simulations are presented in Section 4.1, and we discuss their outcomes in a parameterized manner in Section 4.2.
For the multi-protoplanet cases, we start with N=20 protoplanets initially.Since our goal is to explore the possibility of giant planet formation, we place these protoplanets within the efficient planet growth zone that we explored in our single-planet growth study (Section 3).Note that the radial width of the zone varies among different parameter setups.We randomly select the separation between protoplanets from 10 to 50 mutual Hill radius to fill all bodies within this zone.We test that the final outcome is not sensitive to the choice of their mutual separations, as long as they are well separated at the beginning.For comparison, we also plot the single protoplanet growth by optimizing r 0 to let the planet reach the highest mass.
We also explore a few cases where N=30.However, due to the gravitational interactions and orbital excitations, there is always a limited number of planets that can grow sufficiently massive and dominate the subsequent dynamical evolution.As a result, the final masses and numbers of giant planets are not strongly dependent on the adoption of N from 20 to 30 (Emsenhuber et al. 2021).However, increasing N would significantly increase the computational time.We thus limit our multiprotoplanet explorations to N=20.
The initial eccentricities and inclinations of the protoplanets follow the Rayleigh distributions, with a scaled eccentricity and inclination of e 0 =2i 0 =0.01.We also randomize the initial phase angles of these bodies.The simulations are terminated when the disks are fully dissipated (5 Myr for the fiducial disks and 10 Myr for the disks with a longer lifetime).Initially, protoplanets are widely separated and their masses increase through pebble accretion.Due to the heterogeneity in growth rates at different r 0 , protoplanets near r tran (dashed line in the left panel of Fig. 5) acquire higher masses than others.Once their masses exceed ∼0.1 M ⊕ , they undergo inward migration.The continued mass increase and convergent migration lead to the compression of these bodies' orbits, triggering dynamical instabilities and frequent planet-planet collisions (denoted by dots in Fig. 5) at t∼1 Myr.

Illustration runs
The outcome of a close encounter between two planets can be determined by the ratio of the surface escape velocity v esc and the escape velocity of the planetary system (= √ 2v K where v K is the Keplerian velocity at the planet location).This can be calculated by (Goldreich et al. 2004) The planet bulk density ρ planet ∼2 g cm −3 since they form beyond the water-ice line.We also find that more massive planets collide at closer-in orbits with generally Λ<1.As such, collisions between close-encounters are favored rather than ejections.At this stage the planets of M p ∼M ⊕ at r∼1 au have moderate eccentricities of 0.1.So the impact velocity among these planets approximates ∼ev K ∼1 km/s, lower than their escape velocity.In this regime the accretion efficiency is very high (see Fig. 6 of Cambioni et al. (2019)) and the perfect merger treatment is therefore appropriate (Asphaug 2010).The collisional timescale is given by τ col ∼(nσ col ∆v) −1 , where n = N/(2πr∆r∆z) is the planet number density, r∼3 au and ∆r∼3 au and ∆z∼i × r, σ col ∼πR p (1 + v 2 esc /∆v 2 ) is the collisional cross section, and ∆v∼ev K is the relative velocities among planets.We can estimate τ col ∼0.3 Myr for the planets with R p =1 R ⊕ and e∼2i∼0.1, in agreement with the results shown in Fig. 5.
After a series of mergers, a few protoplanets double their masses, which boosts their subsequent pebble accretion.These bodies, with masses around M p ∼M opt , begin outward migration, leading to chaotic orbits of planets close to r tran .This triggers a second phase of strong perturbations and planet-planet collisions at t∼1.5−2 Myr.Collisions in this phase are not as intense as the first one, since the number of planets in the system has been reduced.
As the gas disk gradually dissipates, the pebble isolation mass drops below ∼3 M ⊕ at t∼1.5 Myr close to the transition radius.Only a few massive bodies remain in the system after multiple planet-planet collisions and scatterings.The largest body reaches the core mass of ∼11 M ⊕ after a collision at ∼ 1.6 Myr.It initaites runaway gas accretion and quick becomes a giant planet with a mass of M p =100 M ⊕ and an orbital period of 30 days.The other three lower mass bodies attain M iso at later times, acquiring the residual disk gas and growing into super-Earth planets.All these planets undergo inward migration, and end up in final orbits of r∼0.1−0.4 au.The outermost three planets are trapped into 4:2:1 mean motion resonances.

Parameter survey
In order to validate our hypothesis that the presence of multiple planet-planet collisions promotes giant planet formation, we conduct an extensive parameter study using N-body simulations by varying two key disk parameters: turbulent level and total solid mass.The solid disk mass correlates with ξ p/g , which is calculated by integrating the pebble mass flux over the disk's lifetime.
We employ a 5 × 5 grids to explore the ranges of α t and ξ p/g .The simulations are conducted at the boundaries where α t = 10 −4 , 5 × 10 −4 , 10 −3 , 5 × 10 −3 , 10 −2 and solid disk mass of 50, 63, 75, 88 and 100 M ⊕ .The intermediate region is populated using linear interpolation.In order to account for the statistical nature of multi-planet interactions, we perform five sets of N-body simulations by randomizing their initial mutual separations and orbital phase angles at each point.The final planet mass is adopted from the largest forming planets over these five realizations.
We discuss the results and implications of the runs with fiducial parameters in Sect.4.2.1, longer disk lifetime in Sect.4.2.2 and initial mass of protoplanet from streaming instability in Sect.4.2.3, respectively.The planet and disk parameters are provided in table 2.

Fiducial case
Fig. 6 displays the highest masses that planets can attain through the growth and evolution of either a single protoplanet (left) or multiple protoplanets (right).The color denotes the final planet mass, with yellow representing the most massive planets and blue representing the lightest ones.The white lines correspond to planet masses of 10, 30, and 100 M ⊕ , respectively.
In the simulations of single protoplanets, gas giant planets exclusively form in disks with α t ∼0.5×10 −3 −10 −2 and total solid mass of 100 M ⊕ (ξ p/g ∼2%).This is because M opt is higher in the moderate and high-turbulent disks, allowing the planet to retain a long time outside before rapid inward migration.Giant planets with orbital period up to 100 days form in our model (an illustration simulation is shown in Fig. A.2). Besides, a higher M iso in these higher α t cases means that planets have the ability to reach more massive solid cores, facilitating subsequent gas accumulation.
Yet, the drawback is that pebble accretion efficiency is lower in such high-turbulent disks.Therefore, giant planet formation only succeeds when disks have a massive supplier of pebble reservoir (high ξ p/g ).This is why massive planets only occur in the yellow region on the right corner of Fig. 6a.Alternatively, in cases of low turbulent disks, planets are incapable of growing massive, as M iso is too low to initiate rapid gas accretion.
In the simulations starting with multi-protoplanets, we find that the giant planet formation zone becomes wider in α t −ξ p/g space and shifts towards lower turbulence and less massive disks (yellow region in Fig. 6b).For disks with α t =10 −3 , M iso at r=1 au is approximately 3 M ⊕ .However, planets with a core mass of M iso fail to grow massive within the disk lifetime (see Fig. 2).Nevertheless, their core mass can be further increased by planetplanet collisions.After a few giant impacts, planets with a core mass of 6−8 M iso can rapidly accrete gas, leading to giant planet formation (see Fig. 5).We find that the growth of giant planets of M p >50 M ⊕ is feasible when disks have more than 70 M ⊕ pebbles at α t 10 −3 (Fig. 6b).
To conclude, compared to the growth of a single protoplanet, giant planet formation is more pronounced in the presence of multiple protoplanets by considering their subsequent convergent migration and planet-planet collisions.Table 2: Disk and planet parameters in Section 4 and Section 5

Disk lifetime
We explore the impact of longer disk lifetime on planet growth and the parameter map is given in Figure 7. Compared to previous runs with a short disk lifetime, we find in Fig. 7a that the giant plant formation zone becomes narrower and only peaks at higher α t and ξ p/g in the single protoplanet case.This can also be understood by comparing Fig. 2 and Fig 3.
In the multi-protoplanet case (Figure 7b), the giant planet formation zone gets extended to lower-mass and less turbulent disks.This is because, first, disk mass decreases slowly in disks with long lifetimes, leading to a protracted supply of gas and pebbles.Second, planets with a higher M opt undergo more pronounced convergent migration (Fig. 3d), enhancing the probability of planet-planet collisions.Third, M iso is higher at later times.All these factors facilitate the giant planet's growth in disks with long disk lifetime.
Therefore, the longer dissipation timescale of turbulent disks promotes the formation of massive solid planetary cores and the accumulation of substantial gas envelopes, resulting in the formation of gas giant planets with a mass approximately 0.7 times that of Jupiter.While in low-turbulence disks, only super-Earths with a few Earth masses can be formed.

Mass of protoplanets
We explore the influence of M p0 by assuming that protoplanets form by streaming instability (equation 8).The parameter map is given in Figure 8.
In single protoplanet cases, planet growth is most efficient at the lowest α t and highest pebble disk mass.This trend is also shown in Figure 4. Since M p0 is much lower than lunar mass in most region of the planetary disk, the planets take a longer time to grow their core masses, resulting in final planets with lower masses compared to those start with equal lunar mass in Figure 6.
In the case of multiple protoplanets, the optimal zone for massive planet formation shifts to α t ∼10 −3 .This is because protoplanets with shorter orbital distances can grow beyond a few Earth masses and migrate quickly into the inner disk region (Figure 4e).They accumulate in a compact configuration, leading to late-phase giant impacts.It is important to note that by this stage, the disk gas has been substantially dissipated, leaving planets The left and right panels illustrate the single protoplanet and multi-protoplanet cases.The protoplanets are assumed to be equal lunar mass and the disk lifetime is 3.7 Myr.Other model parameters are listed in Table 1.The colorbar gives the resulting planet with the highest mass, while the contour lines indicate masses of 10, 30, and 100, respectively.Compared to the single protoplanet case, the parameter ranges for the formation of giant planets are wider in the multi-protoplanet case.  1.The colorbar gives the resulting planet with the highest mass, while the contour lines indicate masses of 10, 30, and 100, respectively.Giant planets are more likely to occur in disks with longer lifetimes.with limited gas envelopes to accrete.Thus, Neptune-mass planets can form in moderately turbulent disks with M solid 85 M ⊕ .

Discussion
We discuss the tension between the observed low dust masses and the model required high solid disk masses for giant planet formation in Sect.5.1.A few aspects of our model are also assessed, including varying Stokes number with disk turbulence, different disk structures, stellar luminosity, and stellar masses in Sect.5.2, Sect.5.3, Sect.5.4 and Sect.5.5.The limitations and caveats are discussed in Sect.5.7.

Disk mass budget for planet formation
Recent observations of protoplanetary disks in various starforming regions have shown that their dust masses typically range from a few hundred to a few Earth masses with a huge scatter (Pascucci et al. 2016;Long et al. 2018;Tobin et al. 2020;Tychoniec et al. 2020;Miotello et al. 2022;Manara et al. 2023).8) the disk lifetime is 3.7 Myr.Other parameters listed in Table 1.The colorbar gives the resulting planet with the highest mass, while the contour lines indicate 10, 30, and 100, respectively.planets around stars M ⋆ =0.1 M ⊙ are difficult to form when the protoplanets from by streaming instability.
These observations indicate that the solid mass of the disk decreases over time, with the highest values in Class 0 disks and a gradual depletion towards the Class II and III phases (e.g., see Figure 2 in Dr ążkowska et al. ( 2023)).Furthermore, dust mass is lower around M-dwarfs compared to their solar-mass counterparts Andrews et al. (2013); Pascucci et al. (2016).For instance, the average solid mass is only ∼1 M ⊕ in ∼2 Myr old Lupus disks around stars of 0.1 M ⊙ , probably due to the radial drift of pebbles (Appelgren et al. 2023).This poses a serious challenge for the formation of giant planets, as it raises the question of whether such disks contain enough solids to form sufficiently massive planetary cores.
It is worth noting that in previous studies, in order to derive the solid disk mass from dust continuum measurements, two assumptions were made: the dust emission is optically thin, and the opacity is mainly due to absorption rather than scattering.However, both of these assumptions have been called into question.Zhu et al. (2019) and Liu (2019) pointed out that optically thick disks with scattering can be misinterpreted as optically thin disks, leading to an underestimation of the disk mass in the literature.
In a recent study by Macías et al. (2021), both scattering and absorption in dust opacity are considered, without making any underlying assumptions on the optical depth.The authors found that in the TW Hydrae disk, the dust mass is ∼300 M ⊕ , a factor of 5 or higher than what would be estimated using typical assumptions.Similar findings are also obtained for the study of the disk of low-mass star ZZ Tau IRS (Hashimoto et al. 2022).These results highlight the importance of considering more realistic dust opacity models in estimating the mass of protoplanetary disks.
On the other hand, the occurrence rate of giant planets around early M dwarfs has been estimated to be less than 5% (Bonfils et al. Sabotta et al. 2021), and an even lower occurrence rate is anticipated around stars of 0.1 M ⊙ .Therefore, the disk conditions that are preferred for the growth of giant planets cannot be considered typical, but rather represent out-liers.Based on the above discussions, it is still likely that early protoplanetary disks around such low-mass stars could contain pebbles with a total mass of >50 M ⊕ .

Turbulence-induced fragmentation-limited Stokes number
In the previous sections, we assume pebbles with a constant Stokes number.Nevertheless, disk turbulence could raise the relative motion between solid particles.When the pebbles' relative velocity is dominated by turbulence, their maximum Stokes number can be expressed as τ s ≈ v 2 F α t c 2 s (Birnstiel et al. 2012), where c s is the gas sound speed, v F is the fragmentation threshold velocity (Blum & Wurm 2008).This means that the Stokes number of the largest pebbles, constrained by the fragmentation limit, decreases as turbulence increases1 .We consider this α t dependence on τ s in this subsection and assume τ s =0.05 ×(10 −3 /α t ), with the adoption of v F =7 m/s.We explore this circumstance with parameters presented in Table 2, and the result is demonstrated in Fig. 9. Compared to Fig. 6, the notable difference occurs for α t 2 × 10 −3 .When pebbles' Stokes number is dependent on α t , the growth of pebbles is strongly suppressed in highly turbulent disks.These small particles cannot settle effectively and therefore pebble accretion is largely impeded.As a result, massive giant planets can only form in massive disks with moderately turbulent level (α t ∼10 −3 ).

Pure stellar irradiation disk
We test a case with a purely stellar-irradiated disk, keeping all other parameters the same as in Figure 2f (see Table 2).The results are shown in Figure 10a.We observe that the disk's scale  2. The planet formation region is significantly narrower compared to the one depicted in Fig. 6, primarily because the growth of pebble size is suppressed in high-turbulence disks.
height is lower in the inner disk region, leading to the formation of planets with lower pebble mass and rapid inward migration.Super-Earth planets are only favored to form in the outer disk region of 10−20 au.We do not find that giant planets form in disks in the absence of viscously heated regions.The effect of multiple protoplanets is expected to be limited because there is only inward migration, and inner planets reach lower M iso at earlier times.The migration is rather divergent, so we do not anticipate a significant difference in the final mass of the planet between the multiple protoplanet case and the single protoplanet case.

Stellar luminosity
We also investigate the impact of stellar luminosity on planet growth.Low-mass young stars typically follow an empirical relation such that L ⋆ ∝ M β ⋆ , where the power-law index β∼1−2.For a star with a mass of 0.1 M ⊙ , its luminosity typically ranges from 0.01 to 0.1 L ⊙ .In our simulations, we adopt a conservative value of 0.01 L ⊙ and assume that the stellar luminosity remains constant over the relatively short disk lifetime of several million years.
We test a higher value of 0.1 L ⊙ while the other model parameters are the same as in Figure 2f.Our simulations show that planet growth becomes more difficult in disks around more luminous stars.In this circumstance, the stellar irradiation region becomes more predominant, and the gas disk scale height is much larger, resulting in extremely low pebble accretion rates.Even in the most metal-rich disks we explored (ξ p/g =2%), growth remains slow and the protoplanet hardly reaches M iso in the outer disk region.Only super-Earth planets of 2 M ⊕ can form at a moderate r 0 ∼5 au (Fig. 10b).In this regard, we expect that massive planets are preferred to form in the late stage of protostellar evolution.
It is important to note that in reality, the stellar luminosity should gradually decline over time (Chabrier & Baraffe 1997;Baraffe et al. 2015), and the more realistic planet growth pattern lies somewhere between the cases of constant low and high lu-minosity (Fig. 2f and Fig. 10b).In future work, we intend to investigate planet formation coupled with a more self-consistent time-dependent evolution of stellar luminosity.

Stellar mass
In addition to the influence of stellar luminosity on the disk profile, the stellar mass is an important factor affecting the availability of supplementary materials within the protoplanetary disk.This, in turn, governs the amount of solid material that protoplanets can accrete.Here we explore the growth and evolution of protoplanets around stars of 0.2 M ⊙ by assuming a linear scaling relationship between the disk mass and stellar mass.Figure 11 displays the maximum planetary mass attained by a single protoplanet as well as multiple protoplanets formed with a lunar mass or through streaming instability.
In the case of a single protoplanet with equal mass, the formation of gaseous planets with mass > 30 M ⊕ is possible in massive disks with M d ≥ 150 M ⊕ or lower mass disks whose α t 10 −3 .The giant planet formation zone extends significantly beyond the fiducial case around 0.1 M ⊕ (see Fig. 6a) for three reasons.Firstly, the pebble isolation mass increases with the stellar mass (Eq.12).Protoplanets orbiting more massive stars can therefore achieve a higher core mass and accumulate a denser atmospheric envelope.Secondly, both the inward migration timescale (Eq.20) and the optimal mass for outward migration (Eq.17) increase with the stellar mass, favoring the growth of protoplanets in the outer regions of the disk.Lastly, the pebble and gas fluxes are also enhanced with stellar mass, providing a greater supply of materials for protoplanet growth.
When considering the mutual interactions of multiple equalmass protoplanets (Fig. 6b), the giant planets are widely available for α t ≥ 5 × 10 −4 due to planet-planet collisions.Whereas the formation for gas giants larger than 0.5 M J is still limited in the turbulent disk with α t > 10 −3 and the solid disk mass > 150 M ⊙ .
Another significant distinction in planet formation around stellar masses of 0.1 M ⊙ and 0.2 M ⊙ is that giant planets can  2).
form from seeds generated by streaming instability (Fig. 6c  and d).Gas giants with masses greater than 100 M ⊕ are exclusively formed in single-protoplanet scenarios when the disk has a turbulent viscosity parameter α t = 5 × 10 −3 , accompanied by a solid disk mass of approximately 200 M ⊙ .In the case of multiple protoplanets, the gas giant formation zone expands to 10 −3 α t 6 × 10 −3 and a solid disk mass exceeding 150 M ⊙ .

Comparison with Observations
We have demonstrated the possibility of forming gas giants with masses ranging from 0.1 to 0.4 M J around stars with a mass of 0.1 M ⊙ , which is consistent with the discovery of four giant planets orbiting host stars with M ⋆ < 0.2 M ⊙ , namely TOI-1227 b (< 0.5 M J ), GJ 3512 b (0.46 M J ), GJ 3512 c (0.45 M J ), and GJ 9066 c (0.21 M J ).We also provide an illustrative example of a giant planet formation (0.3 M J ) promoted by pebble accretion and planet-planet collisions at a location of approximately 0.09 au (see Fig. 5), exhibiting an agreement with the transit observation of planet TOI-1227 b.However, we admit that the smooth disk assumption in our study faces challenges in explaining the orbital characteristics of the other three distant giant planets discovered by radial velocity surveys.The rapid inward migration of planets, predominantly caused by the Lindblad torque, occurs prior to the formation of a surface density gap (Paardekooper et al. 2011).We anticipate that a structured disk with rings and gaps would effectively sup-press the inward migration of planets (Baillié et al. 2016) and facilitate the formation of giant planets in wide orbits.We plan to explore this possibility in our future work.Besides, the interactions between the planets and the gaseous disk result in orbital circularization (Kley & Nelson 2012).In order to reproduce the observed high orbital eccentricities of GJ 9066 c, it is essential to consider the close encounters and mutual scatterings during the later dynamical evolution of planets in a gas-free environment (Ji et al. 2011;Ida et al. 2013).

Caveat
Bell & Lin (1994)'s opacity law is widely adopted in the community, which is based on the assumption of ISM-like grains with compact spherical structures.However, various physical processes have taken place in protoplanetary disk environments (e.g., grain growth and crystallization), and the corresponding disk dust opacity could significantly differ from the ISM's form.A realistic opacity calculation relies on the detailed size distribution, composition and porosity of the dust population, which yet remains poorly understood.In this study we assume a simple opacity law.The opacity variation across the ice lines is also neglected.It is worth pointing out that the migration direction might be reversed and planets get trapped at distinctive regions due to opacity transition (Kretke & Lin 2012).In our study there is only one convergent migration radius.When considering multiple transition radii, although the individual growth pattern differs, the optimal disk condition obtained in this work (e.g., disk mass) for giant planet formation still generally holds.
Besides, we adopt a simplified approach that assumes a constant Stokes number of pebbles and a fixed pebble-to-gas flux ratio under all varied disk and stellar environments (except in Section 5.2).The assumption of a fixed pebble-to-gas flux ratio can be justified if pebbles remain small in the outermost regions of the protoplanetary disc (Johansen et al. 2019), so that the pebble flux follows closely the gas mass flux.In this approach, pebbles may still grow large in the inner region of the disc, even though their total flux remains small and hence the pebble population is not depleted.In more realistic conditions, the coagulation and fragmentation of dust result in them following a powerlaw size distribution (Birnstiel et al. 2011), and the pebble flux is not necessarily always attached to gas flow (Dr ążkowska et al. 2023).In future work, we aim to implement a more sophisticated dust-size population and flux profile (e.g., from Dustpy code, Stammler & Birnstiel (2022)) to gain a better understanding of how these factors influence final planet growth.
We also utilize a simplified gas accretion model that disregards the influence of disk temperature on the gas accretion rate.It has been demonstrated that the elevated temperatures closer to the central star result in higher thermal energy, posing a challenge for the hot gas to be gravitationally captured by and accreted onto the planet (Coleman et al. 2017).Planets located within the innermost disk region require a more extended period to cool down, leading to a reduction in the gas accretion rate (Piso & Youdin 2014).According to Lee et al. (2014)), the runaway accretion timescale follows a power-law relationship with disk temperature (τ run ∝ T 0.34 ).In such circumstances, the giant planet form at a relatively later stage when the inner gaseous disk has cooled slightly but the disk surface density remains high.We also note that in our simulations, some planets initiate gas accretion beyond the water snowline and subsequent migrate inward.These planets accrete a substantial amount of gas before entering the innermost disk region.We note that in Section 5.4 we have demonstrated that the stellar luminosity influences the disk thermal structure and therefore pebble accretion mass growth.The stellar luminosity cannot remain constant during its whole evolution.In particular, in the early stage, the stars are more luminous and pebble accretion may hardly be effective.The core growth can proceed when stars gradually cool.Since we mainly consider the stars with relatively low luminosities, the condition obtained in this paper can be treated as an optimistic perspective.We leave the implementation of a proper stellar evolution in our following studies.

Conclusions
In this paper we investigate the formation of giant planets around late M dwarf stars with a stellar mass range between 0.1 and 0.2 M ⊙ .Although these planets are rare in exoplanet surveys, their formation mechanism remains unclear.Previous studies have suggested that the core accretion scenario faces difficulties in explaining the existence of such high planet-to-star mass ratio systems around these small stars (Liu et al. 2019a;Coleman et al. 2019;Liu et al. 2020;Miguel et al. 2020;Burn et al. 2021).
To address the issue of whether these giant planets can form through core accretion scenario, we use the pebble-driven planet formation model proposed by Liu et al. (2019a) and perform Nbody simulations to study the growth and migration of single and multiple protoplanets in the protoplanetary disk with inner viscously heated and outer stellar irradiated regions.Our simulations incorporate various physical processes, including pebble accretion onto planet cores, gas accretion onto planet envelopes, planet-planet interactions/collisions, type I and type II planet migration and gas damping.We study the influence of several key disk and planet properties, such as the disk turbulent level, solid disk mass (or flux ratio of pebbles and gas), disk lifetime, birth mass of the protoplanets and stellar mass.
The most favorable region for planet growth is near the transition radius r tran ∼3 au that separates the inner viscously heated and outer stellar irradiated regions.However, in the case of single protoplanet growth, it is difficult for the planet to accrete a massive gaseous atmosphere due to its low pebble isolation mass, which is typically around ∼2−3 M ⊕ in systems around stars of M ⋆ =0.1 M ⊙ (Figure 1).When considering multi-protoplanets with the same disk condition, their core growth is no longer limited by pebble isolation.Planets massive enough can undergo convergent migration and evolve into tightly compact orbits, which likely induces subsequent orbital crossings and planetplanet collisions.In general, this dynamical process can overcome the pebble isolation mass barrier for the single protoplanet and promote the growth of a massive core even in very low-mass stellar host systems (Figure 5).
Two different birth masses of protoplanets are considered.On the one hand we consider the protoplanets form from runaway/oligarchic planetesimal accretion and end up with masses of 0.01 M ⊕ .The parameter space for giant planet formation significantly expands when we take into account the growth of multiple protoplanets.Gaseous planets with masses exceeding 100 M ⊕ can form around stars with mass of 0.1 M ⊙ in disks characterized by α t >10 −3 and solid mass 60 M ⊕ (Figure 6).More massive planets are preferred to grow in disks with a longer lifetime and higher supply of pebble reservoirs (Figures 3 and  7).Meanwhile, the giant planet formation benefits from the increasing stellar mass, due to high pebble isolation mass, massive solid disks and long planet migration timescale in systems around massive stars.
On the other hand, in the streaming instability scenario the birth mass of protoplanets increases with orbital distance and stellar mass.Generally, M p0 is much lower than 0.01 M ⊕ at r<10 au.In the single protoplanet case, super-Earth planets only form in the outer region of low-turbulence disks around 0.1 M ⊙ stars (Figure 4).Neptune-mass planets can form in the multiple protoplanet case in disks with α t ∼10 −3 and solid mass exceeding 80 M ⊕ (Figure 8).For systems around more massive stars of 0.2 M ⊙ , the formation of giant planets takes place in disks with 10 −3 α t 6 × 10 −3 and solid disk masses >150 M ⊙ (Figure 11).
Overall, our study highlights the crucial finding that the formation of giant planets with orbital periods of 100 days is favored in turbulent and massive protoplanetary disks.The extended lifetime of the disk and a higher stellar mass contribute to the formation of more massive planets, despite a narrower formation zone within long-lived disks.If protoplanets arise from streaming instability, they only give rise to the birth of giant planets when the stellar mass exceeds 0.2 M ⊙ .
We propose the formation of giant planets with masses ranging from 0.1 to 0.6 M J around stars with masses of 0.1 M ⊙ .This finding aligns with the observed planetary mass in GJ 3512, GJ 9066, and TOI-1227 systems, which were studied through the CARMENES and TESS programs (Morales et al. 2019;Mann et al. 2022;Quirrenbach et al. 2022).Furthermore, our results suggest an increasing feasibility of giant planet formation as stellar mass increases, indicating a correlation between the occurrence rate of giant planets and stellar mass (Bryant et al. 2023;Gan et al. 2023b;Ribas et al. 2023).We anticipate that ongoing and upcoming exoplanet search projects, such as TESS (Ricker et al. 2015), MEarth (Irwin et al. 2009), TRAP-PIST (Jehin et al. 2011), SPECULOOS (Sebastian et al. 2021), CARMENES (Quirrenbach et al. 2014), EDEN (Gibbs et al. 2020), PLATO (Rauer et al. 2014), ET (Ge et al. 2022), and CHES (Ji et al. 2022), will provide a larger sample of planets with well-constrained mass and orbital properties.As the gas accretion at different place indicates a different gas composition of giant planets, we also expect that JWST observations offer valuable insights into the composition of planetary atmospheres to constrain the birthplaces of giant planets.These datasets will significantly contribute to our understanding of planetary formation around very low-mass stars.3D pebble accretion efficiency v K Keplerian velocity at planet's location v esc Escape velocity of the planetary system ∆v Relative velocity between the pebbles and planet ξ p/g Pebble-to-gas mass flux ratio Fig. 1: Growth and migration of individual protoplanets initiated at different disk locations around stars of M ⋆ =0.1 M ⊙ .The dashed line represents the pebble isolation mass at t=1.3 Myr.This is when the fastest growing planet reaches its isolation mass.The arrow indicates the disk transition radius, and the increasing sizes of the dots denote the disk evolution at one Myr intervals.The planet attains the highest mass at a moderate birth radial distance close to the transition radius.The planet and disk parameters are referred to Table1.

Figure 5
Figure5is an example that demonstrates the growth and migration of multiple protoplanets.The planets that survived after

Fig. 5 :
Fig.5: Semi-major axis (left) and mass (right) evolution from multi-protoplanet growth.The planet and disk parameters are:M p0 = 0.01 M ⊕ , α t =5 × 10 −3 , M d =0.15 M ⋆ , ξ p/g =1.75%and t disk =3.7 Myr.The filled dots indicate the planet-planet collisions.The transition radius is illustrated in the dashed line in the left panel, while the dot-dashed lines in the right panel represent the pebble isolation mass at the transition radius.By a combination of pebble accretion and planet-planet collisions, a system with one gas giant and three super-Earths form in very low-mass stars of M ⋆ =0.1 M ⊙ .

Fig. 6 :
Fig.6: Formation of massive planets in the equal protoplanet mass scenario as a function of disk turbulent level and solid disk mass.The left and right panels illustrate the single protoplanet and multi-protoplanet cases.The protoplanets are assumed to be equal lunar mass and the disk lifetime is 3.7 Myr.Other model parameters are listed in Table1.The colorbar gives the resulting planet with the highest mass, while the contour lines indicate masses of 10, 30, and 100, respectively.Compared to the single protoplanet case, the parameter ranges for the formation of giant planets are wider in the multi-protoplanet case.

Fig. 7 :
Fig.7: Formation of massive planets in the long disk lifetime scenario as a function of disk turbulent level and solid disk mass.The left and right panels illustrate the single protoplanet and multi-protoplanet cases.The protoplanets are assumed to be equal lunar mass and the disk lifetime is 6.3 Myr.Other model parameters are listed in Table1.The colorbar gives the resulting planet with the highest mass, while the contour lines indicate masses of 10, 30, and 100, respectively.Giant planets are more likely to occur in disks with longer lifetimes.

Fig. 8 :
Fig. 8: Formation of massive planets in the streaming instability protoplanet mass scenario as a function of disk turbulent level and solid disk mass.The left and right panels illustrate the single protoplanet and mult-protoplanet cases.The protoplanets are assumed form by streaming instability where their masses follow Equation (8) the disk lifetime is 3.7 Myr.Other parameters listed in Table1.The colorbar gives the resulting planet with the highest mass, while the contour lines indicate 10, 30, and 100, respectively.planets around stars M ⋆ =0.1 M ⊙ are difficult to form when the protoplanets from by streaming instability.

Fig. 9 :
Fig.9: Similar to Fig.6but in the fragmentation-limited disk.The stokes number is inversely related to the turbulence strength, and we adopt τ s =0.05 when α t =10 −3 .parameters are listed in Table2.The planet formation region is significantly narrower compared to the one depicted in Fig.6, primarily because the growth of pebble size is suppressed in high-turbulence disks.

Fig. 10 :
Fig. 10: Growth and migration of individual protoplanets at different disk locations around stars of M ⋆ =0.1 M ⊙ in a pure stellar irradiation disk with L ⋆ =0.01 L ⊙ (top) and in a viscously heated and stellar irradiated disk with L ⋆ =0.1 L ⊙ (bottom).The model parameters are similar to that in Figure 2f (see Table2).

Fig. 11 :
Fig. 11: Similar to Fig. 6 (panel a and b) and Fig. 8 (panel c and d) but for the stellar mass of 0.2 M ⊙ .The disk accretion rate and stellar luminosity increase with stellar following Ṁ⋆ ∝ M ⋆ and L ⋆ ∝ M 2⋆ , respectively.Larger gaseous planets are achieved around more massive stars even in the case that protoplanets form through streaming instability.

Table 1 :
Disk and planet parameter setup in Sect.3.

Table A .
1: List of notations.