© The Authors 2026

1 Introduction

Over the last decade, giant exoplanets have been discovered orbiting M dwarfs, and their number continues to grow (e.g., Kanodia et al. 2024). Giants with orbital periods of less than 100 days are less common around M dwarfs compared to Sun-like stars, with occurrence rates below 1% (Bryant et al. 2023). In contrast, giants with orbital periods greater than 100 days show higher occurrence rates, ranging between 1% and 5% (Clanton & Gaudi 2014; Meyer et al. 2018). In particular, all giant planets with masses above 0.3 M_Jup discovered around stars between 0.08 and 0.2 M_⊙ have orbital periods longer than 100 days and were detected mostly via microlensing¹. The orbital distances of single giant planets range from 0.4 to 5 au. From the observed sample, 19 correspond to single cold-giant-planet systems, representing about 85% of the population, while the remaining 15% consists of four planets belonging to two-giant-planet systems. Notably, no system has yet been found hosting both a giant planet and an inner rocky companion, unlike several systems around Sun-like stars (e.g., Zhu & Dong 2021).

The mechanisms underlying giant planet formation remain poorly understood and continue to challenge existing theoretical models (e.g., Liu et al. 2019; Burn et al. 2021; Schlecker et al. 2022; Shibata & Helled 2025). Limitations in dust and gas disk mass are thought to restrict the material available for core growth and gas accretion (e.g., Ansdell et al. 2017; Pinilla 2022). For instance, the discovery of a close-in Neptune-mass planet orbiting a 0.1 M_⊙ star could only be explained with core accretion models that assume a disk mass an order of magnitude higher than typically observed (Stefánsson et al. 2023). Planet population synthesis models, when confronted with planets detected around M dwarfs by HARPS and CARMENES radial velocity surveys, failed to reproduce a population of giant planets around stars less massive than 0.5 M_⊙, suggesting a gap in our understanding of giant planet migration around very low-mass stars (Schlecker et al. 2022). Recently, Pan et al. (2024) demonstrated that giant planets can form around late M dwarfs with periods shorter than 100 days. In their models, lunar-mass protoplanets can grow through a combination of pebble accretion and planet–planet collisions provided that the disks are sufficiently massive, with M_dust>50 M_⊕ and M_gas ∼ 15 M_Jup, assuming a low gas surface density regulated by a high accretion viscosity (α_g) of 10⁻² and a local turbulent viscosity (α_t) of 10⁻²−10⁻³. However, observational studies suggest that turbulence in disks is better characterized by lower turbulent viscosities, closer to 10⁻⁴ (e.g., Rosotti 2023). Moreover, the model assumed inward migration for giant planets, which may not be accurate for late M dwarfs. A new study based on hydrodynamical simulations using the FARGO3D code, assuming high planet-to-star mass ratios in low-viscosity disks (α_t=10⁻⁴), shows that planets more massive than 66 M_⊕ orbiting a 0.1 M_⊙ star can instead undergo outward migration (Sanchez et al. 2025).

This observational and theoretical context suggests that cold giant planets around the least massive stars likely require an alternative formation pathway. To address this, we performed a set of N-body simulations that implement the new prescriptions for giant planet migration around late M dwarfs proposed by Sanchez et al. (2025), starting from planetary seeds that grow by pebble and gas accretion. Our findings point to a new pathway for the formation of cold giant planets around the least massive M dwarfs and the survival of inner rocky worlds.

2 Planet formation model

We used a modified version of the N-body code Mercury (Chambers 1999) that includes additional physical processes:

• An evolving gas-disk model. We adopted the disk model proposed by Ida et al. (2016), which describes a viscous disk in the inner region and an irradiated disk in the outer region. We assumed an initial disk age of 1 Myr and a disk lifetime of 10 Myr, consistent with disk observations of M dwarfs (e.g., Downes et al. 2015; Pfalzner et al. 2022).

• Star–planet tidal interactions and general relativistic corrections. We followed the equilibrium tidal model adapted by Bolmont et al. (2011) for the case of planets around brown dwarfs and M dwarfs. We included acceleration corrections due to relativistic effects proposed by Anderson et al. (1975).

• Planet-disk interactions for rocky planets. We implemented the dynamical friction approach proposed by Ida et al. (2020) for the acceleration corrections experienced by planets embedded in the disk, and the non-isothermal torque prescription of Paardekooper et al. (2010, 2011).

• Growth by pebble accretion. We let planets grow via pebble accretion until they reach the isolation mass (Bitsch et al. 2018), following the model of Lambrechts et al. (2014), with efficiencies for planets on eccentric orbits given by Ormel & Liu (2018) and Liu & Ormel (2018), assuming that the planet’s inclination satisfies i<h_peb, with h_peb the pebble aspect ratio.

The detailed prescriptions of these effects are presented in Sánchez et al. (2020, 2022); Sanchez et al. (2024). In this work, we present an updated version of the code that also includes the following modifications (see the appendix for details):

• New parameterization of the pebble flux. We defined an evolving pebble flux set by a squared exponential decay to ensure consistency with the dust mass evolution inferred from observations (Appelgren et al. 2023), expressed as $$M_{\mathrm{flux}}=M_{\mathrm{flux},\mathrm{in}}exp{(-t/Myr)}^2,$$(1) with M_flux,in=2.5 × 10⁻⁵ M_⊕ yr⁻¹ the initial pebble flux at 1 Myr, set to half of the flux assumed in Sanchez et al. (2024). We assumed a compact disk with a radius of 20 au and a pebble density profile consistent with millimeter-flux observations (Sanchez et al. 2024; Guerra-Alvarado et al. 2025).

• New inputs for star-planet tidal interactions. We included the necessary inputs to also account for the interaction of giant planets around M dwarfs (Bolmont et al. 2015).

• Gas accretion onto the planets. We let planets accrete gas once they reach the isolation mass, as in Pan et al. (2024). We highlight that the mass needed for a planet to achieve an accretion rate of 10⁻⁵ M_⊕ yr⁻¹ is ∼ 5 M_⊕. At this rate, the planet can double its mass within ∼ 0.5 Myr, thereby triggering the onset of runaway gas accretion.

• Updated planet-disk interactions. We implemented updated prescriptions to calculate the torques exerted by the disk onto the planets. For gap-opening planets with planet-to-star mass ratios (q) satisfying q<0.002 we adopted the prescription of Kanagawa et al. (2018), while for planets with q>0.002, we applied the new torque prescription proposed by Sanchez et al. (2025), which accounts for outward migration.

We characterized the code following Sanchez et al. (2024), adopting the hybrid integrator and modeling collisions between bodies as perfect mergers. We ran a set of ten simulations for 50 Myr, with a central object of 0.1 M_⊙ and 25 lunar-mass planetary seeds, a mass value commonly adopted in previous formation models around M dwarfs (e.g., Coleman et al. 2019; Sanchez et al. 2024; Pan et al. 2024). We note that even in the inner disk, a seed with only 1% of the lunar mass can still undergo efficient pebble accretion (Ormel & Liu 2018; Liu & Ormel 2018). The seeds are initially distributed between 0.1 and 2 au, with separations of 15–30 mutual Hill radii, slightly larger than those typically used (e.g., Kokubo & Ida 2000). Previous simulations with more compact configurations and larger samples (Sanchez et al. 2024) showed no significant differences in the early system evolution compared to this work. We used a twoα model: one, α_g, to characterize the global angular momentum transport that regulates the gas surface density evolution, and another, α_t, associated with local disk turbulence, which controls the vertical stirring of pebbles and planet-disk interactions. For our main set of simulations, we adopted α_t=α_g=10⁻⁴. The initial pebble mass (M_peb) was 6 M_⊕ and the gas mass (M_gas) was 10 M_Jup, corresponding to 10% of the stellar mass. These values were obtained by integrating the initial pebble surface density profile and the gas surface density profile, respectively, over the disk radius of 20 au (see profiles in Sanchez et al. 2024). The resulting pebble-to-gas ratio is 10⁻³, which represents a limiting case since not all dust is assumed to be in pebble form (e.g., Lambrechts & Johansen 2014). The Stokes numbers within 2 au range between 0.02 and 0.05, corresponding to centimeter-sized particles at the mid-plane due to the high gas density. This remains fully consistent with ALMA observations since in dense disks, a vertically settled and coagulated midplane population of centimeter-sized pebbles can still be entirely compatible with the observed millimeter-sized emission. This is because the observed flux depends primarily on the integrated dust mass rather than on the local midplane size distribution. To explore the effect of lower disk masses, we used the same initial conditions of the simulations that formed a giant planet but reduced the gas masses to 0.05 M_* and 0.01 M_*, by lowering the gas disk surface density profiles assuming α_g=10⁻³ and α_g=10⁻², respectively. Additionally, we extended the disk lifetime to 15 Myr to assess its impact on the final mass that the giant planet could reach.

3 From seed to cold giant

We investigated the conditions required for a planetary seed to grow into a cold giant planet more massive than 0.3 M_Jup around a star of 0.1 M_⊙ using our numerical model. Here we describe the dynamical evolution of the cold giant planet within the system in which it is formed, and the impact on the survival of an inner rocky planet companion.

Fig. 1 Evolutionary pathway of cold giant planet formation. Top: evolution of the planetary mass (black line) across different growth regimes: the pebble-accretion–dominated phase (orange area), limited by the isolation mass (dotted orange line); the phase dominated by planet–planet collisions (pink area), with the critical mass required to trigger runaway gas accretion indicated (dotted pink line); and the gas-accretion–dominated phase (green area), during which the planet grows until it reaches its final mass (dotted green line) before disk dispersal (dash-dotted black line). Bottom: evolution of the semi-major axis (black line) across the different migration regimes: type I migration (light purple area), transition regime where migration slows (light blue area), and outward migration (orange area).

3.1 Conditions needed for the formation of a cold giant planet

The evolution of the mass and semi-major axis of a planetary seed that becomes a cold giant planet is shown in Fig. 1. In the early stages, growth is dominated by pebble accretion, until the outermost planet reaches the isolation mass within ∼ 300 000 yr after migrating inward on a short timescale. Thereafter, collisions between Earth-mass embryos drive further growth, allowing the planet to reach ∼ 5 M_⊕. Within the following million year, the planet’s mass grows to ∼ 10 M_⊕, the critical mass required to trigger runaway gas accretion. Subsequently, the planet grows to ∼ 66 M_⊕ within ∼ 2 Myr, with its growth regulated by the gas supply set by the stellar accretion rate. Once the planet reaches this mass, it transitions to a different migration regime (see the appendix for details), reverses its migration direction, and begins to migrate outward. In the late stages of disk evolution, the planet accretes in the Hill regime, but once it migrates beyond ∼ 1 au, the stellar accretion rate again limits the gas supply, slowing its growth until disk dispersal. The final mass of the planet depends on how quickly it reaches the ∼ 5 M_⊕ threshold. Including gas accretion at such low core masses is essential rather than waiting until ∼ 10 M_⊕ as in classical models. Notably, the 5 M_⊕ threshold is consistent with the minimum mass required for nucleated instability to begin (Rafikov 2006). Outward migration is crucial in this scenario, as without it giant planets could not reach wide orbits. Pebble accretion efficiencies are high (∼ 50%) because seeds grow mainly within 0.5 au. While we do not include pebble fragmentation inside the snow line (e.g., Venturini et al. 2020), which could reduce accretion, most inner seeds either survive as one close-in planet or collide with the star. In contrast, giant planets form from seeds beyond the snow line. Thus, including fragmentation would likely have little impact on their formation but might require a higher dust mass to reach similar growth.

3.2 Inner rocky planet survival in a system with a cold giant

We analyzed the evolutionary pathways of two simulations in which a giant planet was formed: one leading to a system with a single cold giant, and the other to a system hosting both a cold giant and a close-in rocky companion. Figure 2 shows the evolution of the planetary mass, M_p, semi-major axis, a, eccentricity, e, and inclination, i, for both cases. In the first 0.5 Myr, planet growth is dominated by pebble accretion, and shortly thereafter by collisions among embryos. By 1 Myr, both systems consist of several Earth-like planets and one protoplanet that reaches ∼ 5 M_⊕, triggering rapid gas accretion and the formation of a giant planet. As the giant grows, it accretes the outer planets, either while still migrating inward (two-planet case) or immediately after beginning outward migration (single-planet case). In the single-planet system, the innermost planet is also accreted by the giant. In contrast, in the two-planet system the innermost planet survives because it enters the inner disk cavity early in its evolution and settles into a stable orbit at ∼ 0.015 au. By the end of the disk’s lifetime, the giant in the single-planet system has a mass of ∼ 0.6 M_Jup at 2.6 au, while in the two-planet system the giant’s mass is ∼ 0.3 M_Jup, as it crosses the 5 M_⊕ threshold later in its evolution. In both scenarios, the surviving planets remain on quasi-circular orbits: the innermost planet due to strong star-planet tidal interactions, and the giant planet due to eccentricity damping during the disk lifetime. During the pebble accretion phase, most planets satisfy i<h_peb, and for most of the integration time i<h_g, where h_peb and h_g are the pebble and gas aspect ratios, respectively (see Sanchez et al. 2024). Close encounters can temporarily raise the inclinations up to ∼ 5 h_g. In both systems, the giant planet remains near the midplane, while in the system with an inner rocky companion, the rocky planet retains a high inclination due to its location inside the disk cavity.

In our explored scenarios, the survival of an inner rocky planet in systems with a cold giant primarily depends on whether it enters the inner disk early, preventing collisions with the growing giant and allowing coexistence. Under these conditions, the rocky planet remains atmosphere-free, but further simulations and longer integrations are needed to explore its composition and long-term evolution, which lie beyond the scope of this work.

4 Discussions

In our simulations, giant planet formation occurs in lowviscosity environments characterized by α_t=10⁻⁴. This appears to differ from the findings of Pan et al. (2024), who identify giant planet formation as most favorable for higher turbulent viscosities, α_t=10⁻³. However, this distinction arises primarily from the different gas surface density profiles adopted in the two studies. Although Pan et al. (2024) assumed a more massive disk (15% of the stellar mass, compared to our 10%), their choice of a higher accretion viscosity of 10⁻² yields a significantly lower gas surface density than our simulations. In contrast, our model uses α_g=10⁻⁴, resulting in gas densities that are a factor 30 higher than those in Pan et al. (2024). This difference in disk structure has direct dynamical consequences. In our denser disk, type I migration is substantially faster, up to two orders of magnitude, than in lower-density disks, causing embryos to converge toward the inner disk on timescales of less than 0.1−0.2 Myr (see Fig. 2), compared with the 1−2 Myr in Pan et al. (2024). Such rapid migration strongly increases mutual collision rates, allowing planetary cores of 5 M_⊕ to form very early (within 0.5 Myr). Provided that the disk lifetime is sufficiently long (up to 10 Myr), these cores can subsequently undergo gas accretion and grow to Saturn–Jupiter masses.

We also note a key difference in the pebble reservoir required to form giant planets. In Pan et al. (2024), giant planet formation required pebble masses greater than 50 M_⊕, whereas our model begins with a much smaller reservoir of 6 M_⊕ (total of 10 M_⊕). This could be related to the fact that their slower-migration environment produces fewer embryo-embryo collisions, making growth more dependent on reaching the pebble isolation mass, a process less efficient at large orbital distances. In contrast, the much faster migration in our dense disk makes early collisional growth the dominant pathway, allowing some embryos to reach the critical core mass even with a modest pebble reservoir. This strong dependence on planet–planet collisions also explains why only a subset of our simulations produce cores massive enough to trigger runaway gas accretion. Overall, the high gas surface density in our model plays a central role by accelerating migration, increasing early dynamical interactions, and enabling the rapid formation of massive planetary cores in disks with relatively small pebble reservoirs.

The pebble accretion efficiencies in our main simulations did not include the inclination-dependent reduction factor. However, to assess its impact, we repeated two runs that formed a giant planet, now incorporating the efficiency reduction for cases where i > h_peb, as in Ormel & Liu (2018). As expected, pebble accretion efficiency decreased by about a factor of 2 during the first 10⁵ yr, when some embryos briefly reached inclinations comparable to the pebble aspect ratio. However, by 3 × 10⁵ yr, a similar number of embryos reached the pebble isolation mass in the two sets of simulations. This is because planet–planet collisions consistently occurred after (1.4 ± 0.05) × 10⁵ yr in all runs, demonstrating that collisional growth is the dominant mechanism driving embryos in the inner disk to reach isolation masses. Therefore, neglecting the early-time reduction does not alter when or where embryos reach this critical stage. Regarding giant planet formation, one revised simulation still produced a giant planet, while the other did not. This divergence did not arise from differences in pebble accretion, but rather from stochastic variations in the dynamical histories, which in one case prevented the sequence of collisions needed to assemble a sufficiently massive core. This underscores that, in our model, collisional growth and the system’s intrinsic dynamical sensitivity primarily determine whether giant planet formation is ultimately achieved.

Fig. 2 Evolution of the planetary mass, semi-major axis, eccentricity, and inclination of the planetary seeds (gray lines) in systems where a single giant planet forms (top row), and where a close-in rocky planet forms together with a cold giant planet (bottom row). We show the evolution of the giant planet (black line) and the evolutions the other planets that survived in the system (colored lines) before the more massive planet became a giant planet. The inner edge of the disk, rin, is indicated (horizontal dashed black line). The location of the innermost planet (pink line) is marked with an arrow for comparison with the inner edge of the disk. The gas aspect ratio, hgas (dotted black lines), and the pebble aspect ratio, hpeb (dotted-dashed black lines), are also overplotted, evaluated at a=0.02 au and a=2 au.

5 Conclusions

We performed N-body simulations to explore the formation of cold giant planets around a 0.1 M_⊙ star in a low-viscosity disk (α_t=10⁻⁴). We find that giant planet formation does not require extreme disk conditions in low-viscosity environments where pebble drift is efficient. For instance:

• Saturn- and Jupiter-like planets can form in disks with a pebble supply sustained over ∼ 3 × 10⁵ yr and a total pebble mass of ∼ 10 M_⊕. The total disk mass should be ∼ 10% of the stellar mass for a compact disk with a high gas surface density (20 au, α_g=10⁻⁴);

• Additional simulations suggest that cold-giant formation is unlikely in disks less massive than 5% of the stellar mass, which are associated with lower gas density profiles and higher accretion viscosities (α_g > 10⁻⁴);

• A combination of planet–planet collisions and efficient pebble accretion before 1 Myr, together with a long disk lifetime of at least 10 Myr, plays a key role in enabling the formation of a cold giant planet. A long disk lifetime has also been found to be crucial for giant planet formation in recent studies (e.g., Shibata & Helled 2025). However, extending the disk lifetime to 15 Myr results in a 5% difference in the final planetary mass, indicating that extending the disk lifetime to more than 10 Myr has a negligible impact on the giant planet formation process.

We note that three of our ten simulations produced a giant planet (two with a single giant and one with a giant and an inner rocky companion), while the others formed Earth- and super-Earth-like planets. A larger set of simulations, currently in progress, will better capture the low occurrence rate of cold giants and further clarify how different initial pebble masses and gas density profiles shape planetary diversity.

Acknowledgements

This work was performed using the compute resources from the Academic Leiden Interdisciplinary Cluster Environment (ALICE) provided by Leiden University. This work made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The authors thank Gudmundur Stefánsson for his helpful comments, and the referee for their detailed report, which has improved the clarity of our manuscript.

References

Appendix A Model of gas accretion onto planets

We modified the core of Mercury N-body code Chambers (1999) to include gas accretion onto the planets. We assumed gas accretion starts once the planets reached the pebble isolation mass (Ogihara & Hori 2020). The planets accrete with a gas accretion rate defined as $${\dot{M}}_{\mathrm{p},\mathrm{g}}=min({\dot{M}}_{\mathrm{g},\mathrm{KH}},{\dot{M}}_{\mathrm{g},\text{Hill}},{\dot{M}}_{\mathrm{g}}),$$(A.1) with Ṁ_g, KH the accretion rate following Kelvin–Helmholtz contraction, which is expected to start after the planet reaches the isolation mass (Ikoma et al. 2000) is set as follows: $${\dot{M}}_{\mathrm{g},\mathrm{KH}}=8\times {10}^{-8}{\left(\frac{M_p}{3M_{\oplus }}\right)}^4{\left(\frac{K_{\mathrm{env}}}{1{\mathrm{cm}}^2{\mathrm{g}}^{-1}}\right)}^{-1}{\mathrm{M}}_{\oplus }{\mathrm{yr}}^{-1},$$(A.2) with M_p the mass of the planet and k_env the envelope opacity, that as it needs to be lower than the disk opacity (in our disk model is set at 1 cm² yr⁻¹) we set at 0.1 cm² yr⁻¹ as in Pan et al. (2024). As the planet increases it mass it will change the accretion rate regime. It could start accreting in the Hill regime, where only a fraction of gas within the planet Hill sphere can be accreted (Tanigawa & Watanabe 2002). Following the expression as appears in Liu et al. (2019) the accretion rate in the Hill regime can be express as follows: $$\begin{array}{rl}{\dot{M}}_{\mathrm{g},\text{\hspace{0.17em}Hill\hspace{0.17em}}}=& 0.004{\left(\frac{M_p}{3M_{\oplus }}\right)}^{2/3}{\left(\frac{M_{\star }}{0.1M_{\odot }}\right)}^{-2/3}\left(\frac{{\dot{M}}_g}{{10}^{-8}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}}\right)\times \\ & \times {\left(\frac{{\alpha }_g}{{10}^{-2}}\right)}^{-1}{\left(\frac{h_g}{0.065}\right)}^{-2}{[1+{\left(\frac{M_p}{M_{\mathrm{g}\mathrm{a}\mathrm{p}}}\right)}^2]}^{-1}{\mathrm{M}}_{\oplus }{\mathrm{y}\mathrm{r}}^{-1}\end{array}$$(A.3) with α_g the viscosity related to gas accretion, h_g the height scale of the gas disk (see the disk model in Sánchez et al. 2022; Sanchez et al. 2024), and Ṁ_g the disk accretion rate which in this work is given by a fitting with observations in the Orion Nebular Cluster made by Manara et al. (2012), as used to describe our disk model (see the appendix in Sanchez et al. 2024) which is express as follows: $$\begin{array}{rl}\log \left(\frac{{\dot{M}}_{\mathrm{g}}}{{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}}\right)=& -5.12-0.46\log \left(\frac{t}{\mathrm{y}\mathrm{r}}\right)-5.75\log \left(\frac{M_{\star }}{{\mathrm{M}}_{\odot }}\right)\\ & +1.17\log \left(\frac{t}{\mathrm{y}\mathrm{r}}\right)\log \left(\frac{M_{\star }}{{\mathrm{M}}_{\odot }}\right)\end{array}$$(A.4) We estimate the necessary mass for a planet to open a gap in the disk following Kanagawa et al. (2015): $$M_{\mathrm{gap}}=17{(\frac{1}{{\Sigma }_{\mathrm{p}}/{\Sigma }_0}-1)}^{0.5}{\left(\frac{h_{\mathrm{g}}}{0.1}\right)}^{2.5}{\left(\frac{{\alpha }_{\mathrm{t}}}{{10}^{-3}}\right)}^{0.5}\left(\frac{M_{\star }}{0.1M_{\odot }}\right){\mathrm{M}}_{\oplus },$$(A.5) with Σ_p/Σ₀ the gas depth associated with gap-opening planets, which we set at 0.08 to assure that it will always be higher that the pebble isolation mass assumed in this work that follows Bitsch et al. (2018, see details in Sanchez et al. 2024).

As an example, in Fig. A.1 we show the accretion rates in the different regimes described above for planets with different masses around a 0.1 M_⊙ star, located close to the inner edge of the disk, at 0.02 au. We assumed a gas-disk of 1 Myr (see the gas-disk model in Sanchez et al. 2024). We highlight the critical mass of around 5 M_⊕ associated with the accretion rate of 10⁻⁵ M_⊕ yr⁻¹, which can trigger the runaway gas accretion.

Fig. A.1 Gas accretion rates as a function of planetary mass. The curves show the Kelvin–Helmholtz contraction regime (blue line), the Hill regime (orange line), and disk-limited accretion (green line). The effective accretion rate onto the planet, defined as the minimum of the three regimes, is shown (dashed black line). The critical planetary masses at which the accretion rate reaches 10−5 M⊕ yr−1 are highlighted (dotted gray lines).

Appendix B Planet-disk interactions

We implemented an updated routine to account for the torques exerted by the gas disk on planets across different regimes. For all planetary masses, we retained the same expressions for the acceleration corrections already included in the code (Sánchez et al. 2022), adopting the formulation proposed by Ida et al. (2020): $${\mathbf{f}}_{\text{gas}}=-\frac{v_{\mathrm{r}}}{t_{\mathrm{e}}}{\widehat{e}}_{\mathrm{r}}-\frac{(v_{\theta }-v_{\mathrm{k}})}{t_{\mathrm{e}}}\widehat{e_{\theta }}-\frac{v_{\mathrm{k}}}{2t_{\mathrm{a}}}\widehat{e_{\theta }}-\frac{v_{\mathrm{z}}}{t_{\mathrm{i}}}\widehat{e_{\mathrm{z}}},$$(B.1) with the velocity of the embryos given by v=(v_r, v_θ, v_z), and ê_r, ê_θ and ê_z the versors in each direction associated with a system in cylindrical coordinates (r, θ, z). The gas velocity is given by v_g=(0,(1−η) v_k, 0), with $\eta \sim h_{\mathrm{g}}^2$ and v_k the Keplerian velocity. The variables t_a, t_e and t_i denote the characteristic timescales for the orbital evolution of the semi-major axis a, eccentricity e, and inclination i of the embryos, respectively. In all cases, we calculate t_a as follows: $$t_{\mathrm{a}}=-\frac{t_{\mathrm{wave}}}{2h_{\mathrm{g}}^2}\frac{{\Gamma }_0}{{\Gamma }_{\mathrm{total}}},$$(B.2) with $$t_{\mathrm{wave}}=\left(\frac{M_{\star }}{M_{\mathrm{p}}}\right)\left(\frac{M_{\star }}{{\Sigma }_{\mathrm{g}}{r}^2}\right)h_{\mathrm{g}}^4{\Omega }_{\mathrm{k}}^{-1},$$(B.3) the timescale proposed by Papaloizou & Larwood (2000) and Tanaka & Ward (2004), and the scaling torque factor: $${\Gamma }_0={(M_{\mathrm{p}}/M_{\star })}^2{\Sigma }_{\mathrm{g}}{r}^4{h}_{\mathrm{g}}^{-2}{\Omega }_{\mathrm{k}}^2,$$(B.4) with M_p the mass of the planet, Ω_k the Keplerian angular velocity, and Σ_g the gas surface density.

The calculation of the total torque that the disk exerted onto the planets Γ_total is computed according to the planetary mass regime For planets with masses above 0.01 M_⊕ and below 66 M_⊕ corresponding to a planet-star mass ratio of q=0.002 for a 0.1 M_⊙ and with $K={\left(\frac{M_{\mathrm{p}}}{M_{\star }}\right)}^2{h}_{\mathrm{g}}^{-5}{\alpha }_{\mathrm{t}}^{-1}<{10}^4$, we apply the torque prescription of Kanagawa et al. (2018), which extends the type I migration regime. In this case, the total torque is $${\mathrm{\Gamma }}_{\text{total\hspace{0.17em}}}=\frac{{\mathrm{\Delta }}_{\mathrm{L}}{\mathrm{\Gamma }}_{\mathrm{L}}+{\mathrm{\Delta }}_{\mathrm{C}}{\mathrm{\Gamma }}_{\mathrm{C}}exp(-K/K_{\mathrm{t}})}{1+0.04K},$$(B.5) with Γ_C and Γ_L the corotation and Lindblad torques, respectively, for a circular and co-planar planetary orbit (as in Paardekooper et al. 2010, 2011) and Δ_C and Δ_L the associated reduction factors for an eccentric and inclined orbit (as in Ida et al. 2020). The details of the estimation of the Lindblad and corotation torques, together with the corresponding reduction factors, are described in Sánchez et al. (2022), where they were implemented for the first time. In its original form, the fit from hydrodynamical simulations made by Kanagawa et al. (2018) does not include the reduction factors, as showed in Eq. (B.5). However, by treating their prescription as an extension of type I migration and following the assumptions of Matsumura et al. (2021), we incorporated these factors into the general fitting formula. The parameter K_t is associated with the gap depth at which the corotation torque becomes negligible. Kanagawa et al. (2018) showed that for K ≳ 20 the gap opened by the planet modifies the torque, and therefore adopted K_t=20. For K ≪ K_t and 0.04 K ≪ 1, the total torque reduces to the linear expression, Γ_total=Δ_L Γ_L+Δ_C Γ_C. We note that adopting K_t=20 corresponds to the transition mass M_tran for partial gap-opening planets proposed by Kanagawa et al. (2018), assuming a gas surface density reduced by a factor of 0.5. This transition mass is $$M_{\text{tran}}=2.7{\left(\frac{{\alpha }_t}{0.001}\right)}^{0.5}{\left(\frac{h_{\mathrm{g}}}{0.05}\right)}^{2.5}\left(\frac{M_{\star }}{0.1{\mathrm{M}}_{\odot }}\right){\mathrm{M}}_{\oplus }.$$(B.6)

For planets with masses greater than 66 M_⊕ and values of K > 10⁴, we adopted a new prescription for the estimation of the total torque estimation given by Sanchez et al (in review), described as $${\mathrm{\Gamma }}_{\text{total\hspace{0.17em}}}/{\mathrm{\Gamma }}_0=\{\begin{array}{clc}a,& \text{\hspace{0.17em}if\hspace{0.17em}}& {\mathrm{\Sigma }}_0>{10}^{-5}\\ \frac{b}{q},& \text{\hspace{0.17em}if\hspace{0.17em}}& {\mathrm{\Sigma }}_0\le {10}^{-5}\end{array},$$(B.7) with a=5 × 10⁻⁴ and b=10⁻⁵, and Σ₀ the local gas-disk density at the location of the planet r₀ in units of $M_{\star }{r}_0^{-2}$. We emphasize that, in this new regime, planets experience outward migration throughout the disk, in contrast to the previous regime, where migration is predominantly inward except in a narrow region near the inner disk edge (see details on type I migration in Sánchez et al. 2022). We also highlight that this is the first time that such a prescription has been incorporated into an N-body model to study the formation and evolution of giant planets.

We note that all the above equations are evaluated at the semi-major axis of the orbit of the planets. The remaining timescales t_e and t_i are calculated as $$t_{\mathrm{e}}=(1+0.04K)\frac{t_{\mathrm{wave}}}{0.78}[1+\frac{1}{15}{(e_{\mathrm{rat}}^2+i_{\mathrm{rat}}^2)}^{3/2}],$$(B.8) $$t_{\mathrm{i}}=(1+0.04K)\frac{t_{\mathrm{wave}}}{0.544}[1+\frac{1}{21.5}{(e_{\mathrm{rat}}^2+i_{\mathrm{rat}}^2)}^{3/2}],$$(B.9) with e_rat=e/h_g and i_rat=i/h_g. Due to the lack of consensus on eccentricity and inclination damping in the type II regime, we follow Matsumura et al. (2021) and adopt the damping timescales by scaling those of type I migration from Ida et al. (2020) with a factor (1+0.04 K). This ensures that the eccentricity and inclination damping timescales remain shorter than the migration timescale, an essential condition for planets to be resonantly trapped near the inner edge of the disk (Ogihara et al. 2010).

As an example, Fig. B.1 shows the normalized total torques experienced by planets of different masses throughout a 1 Myr disk (see details of the disk model in Sanchez et al. 2024), at various distances from a 0.1 M_⊙ star, assuming coplanar and circular orbits. The figure illustrates type I migration as described by linear theory (Paardekooper et al. 2010, 2011; Ida et al. 2020), with a transition at the mass M_tran (Eq. (B.6)) where type I migration can be extended following Kanagawa et al. (2018). In this regime, the migration slows down due to a reduction in the total normalized torque. In both cases, we adopt the torque prescription given by Eq. (B.5). For more massive planets satisfying a planet-to-star mass ratio q > 0.002 and K > 10⁴, the normalized torque is calculated according to Eq. (B.7). The figure highlights the different regimes of inward and outward migration in each case.

Fig. B.1 Normalized torques as a function of planetary mass and semi-major axis, assuming circular and coplanar orbits. The different migration regimes are indicated, separating inward from outward migration. The mass required to open a gap in the disk (dashed white line), the planet-to-star mass ratio (q) of 0.002, and the gas depth criterion K=104 are over-plotted (dashed black lines).

Appendix C Star-planet tidal interactions

We adopted the equilibrium tidal model proposed by Hut (1981) and reformulated by Bolmont et al. (2011) to describe the star-planet tidal interactions in systems of planets orbiting brown dwarfs and M dwarfs. This model was originally implemented in the N-body code by Sánchez et al. (2020). The tidal forces acting on each planet account for the tides raised both by the host star on the planet and by the planet on the star. The acceleration induced on the planets by tidal forces includes terms that affect the evolution of the argument of periastron, ω(f_tide−ω), as well as terms that influence the evolution of the semi-major axis, a, and eccentricity, e (f_tide−ae), as follows: $${\mathbf{f}}_{\mathrm{tide}-\omega }=-3\frac{\mu }{r^8}[k_{2,\star }\left(\frac{M_{\mathrm{p}}}{M_{\star }}\right)R_{\star }^5+k_{2,\mathrm{p}}\left(\frac{M_{\star }}{M_{\mathrm{p}}}\right)R_{\mathrm{p}}^5]\mathbf{r},$$(C.1) $$\begin{array}{rl}{\mathbf{\text{f}}}_{\text{tide-ae}}=& -3\frac{\mu }{r^{10}}\left[\frac{M_{\text{p}}}{M_{\star }}{k}_{2,\star }\mathrm{\Delta }{\mathrm{t}}_{\star }{R}_{\star }^5(2\mathbf{\text{r}}(\mathbf{\text{r}}\cdot \mathbf{\text{v}})+r^2(\mathbf{\text{r}}\times {\mathrm{\Omega }}_{\star }+\mathbf{\text{v}}))\right]\\ -& 3\frac{\mu }{r^{10}}\left[\frac{M_{\star }}{M_{\text{p}}}{k}_{2,\text{p}}\mathrm{\Delta }{\mathrm{t}}_{\text{p}}{R}_{\text{p}}^5(2\mathbf{\text{r}}(\mathbf{\text{r}}\cdot \mathbf{\text{v}})+r^2(\mathbf{\text{r}}\times {\mathrm{\Omega }}_{\text{p}}+\mathbf{\text{v}}))\right],\end{array}$$(C.2) with r and v the position and velocity vector of each planet, respectively, μ=G(M_*+M_p), M_p and R_p the mass and radius of the planet, k_{2, *} and k_{2, p} the potential Love numbers of degree 2 of the star and the planets, and Δ t_* and Δ t_p the time-lag model constants for the star and each planet, respectively. The factors k_{2, *} Δ t_* and k_{2, p} Δ t_p are related with the dissipation factors by $$k_{2,\mathrm{p}}\Delta {\mathrm{t}}_{\mathrm{p}}=\frac{3R_{\mathrm{p}}^5{\sigma }_{\mathrm{p}}}{2G}{k}_{2,\star }\Delta {\mathrm{t}}_{\star }=\frac{3R_{\star }^5{\sigma }_{\star }}{2G},$$(C.3) with the dissipation factor for each planet σ_p, and the dissipation factor of the host star σ_*=2.006 × 10⁻⁵³ k⁻¹ m⁻² s⁻¹, the same factor for a M dwarf (Hansen 2010).

We include the evolution of the stellar rotational velocity, Ω_*, following the prescription proposed by Bolmont et al. (2011), but neglecting the impact on the planets in its evolution, for simplicity. We also account for the evolution of the stellar radius, R_*, by fitting the models of Baraffe et al. (2015) for a star with a mass of 0.1 M_⊙. We assume that each planet’s rotational velocity, Ω_p, corresponds to the pseudo-synchronous state defined by Hut (1981, see further details in Sánchez et al. 2020).

In this work, we compute the planetary radius assuming a spherical body and adopt three different bulk densities depending on the planetary mass. For planets with M_p<3 M_⊕ we assume a density ρ=5 g cm⁻³, for those with 3<M_p/M_⊕<10, we use ρ=2.5 g cm⁻³; and for planets with M_p>10 M_⊕, we adopt ρ=1 g cm⁻³. These values are consistent with typical densities of planets in the Solar System. Moreover, we assume the same tidal dissipation factor for planets with M_p<10 M_⊕, as that of the Earth, σ_p=8.577 × 10⁻⁴³ k⁻¹ m⁻² s⁻¹, estimated by (Neron de Surgy & Laskar 1997). For planets M_p > 10 M_⊕, we adopt the dissipation factor of Jupiter, σ_p=7.024 × 10⁻⁵² k⁻¹ m⁻² s⁻¹, following Hansen (2010).

All Figures

Fig. 1 Evolutionary pathway of cold giant planet formation. Top: evolution of the planetary mass (black line) across different growth regimes: the pebble-accretion–dominated phase (orange area), limited by the isolation mass (dotted orange line); the phase dominated by planet–planet collisions (pink area), with the critical mass required to trigger runaway gas accretion indicated (dotted pink line); and the gas-accretion–dominated phase (green area), during which the planet grows until it reaches its final mass (dotted green line) before disk dispersal (dash-dotted black line). Bottom: evolution of the semi-major axis (black line) across the different migration regimes: type I migration (light purple area), transition regime where migration slows (light blue area), and outward migration (orange area).

In the text

Fig. 2 Evolution of the planetary mass, semi-major axis, eccentricity, and inclination of the planetary seeds (gray lines) in systems where a single giant planet forms (top row), and where a close-in rocky planet forms together with a cold giant planet (bottom row). We show the evolution of the giant planet (black line) and the evolutions the other planets that survived in the system (colored lines) before the more massive planet became a giant planet. The inner edge of the disk, rin, is indicated (horizontal dashed black line). The location of the innermost planet (pink line) is marked with an arrow for comparison with the inner edge of the disk. The gas aspect ratio, hgas (dotted black lines), and the pebble aspect ratio, hpeb (dotted-dashed black lines), are also overplotted, evaluated at a=0.02 au and a=2 au.

In the text

Fig. A.1 Gas accretion rates as a function of planetary mass. The curves show the Kelvin–Helmholtz contraction regime (blue line), the Hill regime (orange line), and disk-limited accretion (green line). The effective accretion rate onto the planet, defined as the minimum of the three regimes, is shown (dashed black line). The critical planetary masses at which the accretion rate reaches 10−5 M⊕ yr−1 are highlighted (dotted gray lines).

In the text

	Fig. B.1 Normalized torques as a function of planetary mass and semi-major axis, assuming circular and coplanar orbits. The different migration regimes are indicated, separating inward from outward migration. The mass required to open a gap in the disk (dashed white line), the planet-to-star mass ratio (q) of 0.002, and the gas depth criterion K=104 are over-plotted (dashed black lines).
In the text