© The Authors 2025

1 Introduction

The observed exoplanet population provides extraordinary opportunities to constrain planet formation processes. With the detection of ~3000 exoplanets during its lifetime, the Kepler telescope (Borucki et al. 2010) yielded many discoveries that have provided insights into the planet formation process. One such discovery is the apparent enhancement of the occurrence rate of close-in super-Earths for M-dwarfs compared to FGK-stars (Mulders et al. 2015). Most planets are thought to form in protoplanetary discs around young stars through core accretion, where the solid core is assembled first, followed by the accretion of a gaseous envelope. There are multiple mechanisms by which a protoplanet can accrete its solid core, such as through the accretion of planetesimals (kilometre-sized), pebbles (millimetre-sized), or a combination of the two (Alibert et al. 2018; Drążkowska et al. 2023)¹. As the dust inventories of protoplanetary discs are observed to increase with increasing stellar mass (Pascucci et al. 2016; Manara et al. 2023), a first-order expectation would be that the occurrence rates of planets should also increase with increasing stellar masses. Indeed, this has been observed for giant planets (e.g. Johnson et al. 2010; Reffert et al. 2015). The mechanisms that led to a decrease of the super-Earth occurrence rate around high-mass stars are therefore currently uncertain.

As high-mass stars are more likely to form giant planets, it has been proposed that the early formation of an outer giant planet could block the influx of pebbles onto an inner planet, preventing it from growing to a super-Earth (e.g. Mulders et al. 2021). Such a model requires two important assumptions. The first is that the inner planet needs to form inside the water iceline and experience a reduction of the pebble flux due to the sublimation of water ice in the protoplanetary disc. This is required as otherwise the formation of the inner super-Earth is significantly faster than the formation of the outer giant due to the much higher efficiency of accreting drifting pebbles in the inner regions of protoplanetary discs (Danti et al. 2025)². The second assumption is that once the outer planet reaches pebble isolation mass, it blocks the entire pebble flux by creating a gap in the protoplanetary disc, quenching planet formation inside of its orbit. The degree to which giant planets are able to trap inward drifting pebbles for an extended period of time is a topic of intense current research (Haugbølle et al. 2019). Stammler et al. (2023) found that mm-sized pebbles are trapped by the planet and fragment into smaller sizes, which are sufficiently coupled to the gas to flow with the gas through the gap. Similar results were found by Van Clepper et al. (2025) and Huang et al. (2025), who performed 3D simulations around an embedded giant planet and found that small dust particles can indeed pass through the gap. Further, outer giant planets have been found to have either a positive or zero correlation with the presence of inner super-Earths (Bryan & Lee 2024; Van Zandt et al. 2025; Bryan & Lee 2025; Bonomo et al. 2025). Should giant planets be able to block a significant amount of the pebble flux, it would be expected that outer giants are negatively correlated with the presence of inner super-Earths, in contrast with observations. Nevertheless, the presence of super-Jupiters has been found to be negatively correlated with the presence of an inner super-Earth, highlighting that the formation of very massive giant planets might indeed inhibit the formation of close-in super-Earths (Lefèvre-Forján & Mulders 2025). Therefore, it is currently uncertain whether the lack of inner super-Earths around high-mass stars can be explained by the presence of an outer giant planet blocking the inward flux of pebbles.

The presence of binary companions around massive host stars has also been proposed as an explanation for the observed super-Earth population. Binary companions are more prevalent for massive stars and could suppress the formation of planets. However, even when accounting for the presence of binaries, Moe & Kratter (2021) found that the enhanced occurrence rate of super-Earths around M dwarfs is only half of the observed enhancement.

Another proposed explanation is that, depending on the choice of model, the pebble accretion efficiency can be enhanced around low-mass stars (Chachan & Lee 2023). By utilising a planet population synthesis model, Pan et al. (2025) were recently able to explain the observed super-Earth occurrence rate as a function of stellar mass by a reduced migration rate around high-mass stars as well as through an increased pebble accretion efficiency around low-mass stars.

In this work, we aim to understand how the formation and migration of super-Earths are affected by the structure and sizes of protoplanetary discs. In Sect. 2, we present our disc and planet formation model. In Sect. 3 we show the results of our planet population synthesis model, while in Sects. 4 and 5, we discuss and conclude our results, respectively. In Appendix A we show the resulting fraction of close-in super-Earths for different disc photoevaporation models. In Appendix B we show the planet migration coefficient as a function of disc size and stellar mass as well as time, while in Appendix C we show some additional simulation results using different model parameter values.

2 Planet formation model

2.1 Protoplanetary disc model

The viscosity of the disc, v, can be written as (Shakura & Sunyaev 1973) $$v={\alpha }_{\text{v}}{c}_{\text{s}}H$$(1)

where α_v is a dimensionless parameter setting the strength of the viscosity, c_s is the sound speed, and H is the disc scale height. Throughout this work, we assumed a value of α_v = 001. For a detailed discussion on our choice of α_v, we refer to Sect. 7.1 of Nielsen et al. (2023) and Sect. 2 of Nielsen et al. (2025). The sound speed can be written as $$c_{\text{s}}=\sqrt{\frac{k_{\text{b}}T}{{\mu }_{\text{d}}{m}_{\text{H}}}}=\sqrt{ℛT},$$(2)

where k_b is the Boltzmann constant, T is the temperature of the disc, µ_d is the mean molecular weight of the disc, m_H is the mass of a hydrogen atom, and ℛ is the mass-specific gas constant. We use a value of 2.35 for µ_d, consistent with a gas that is composed of mostly H₂ and He. The gas scale height, in turn, can be written as $$H=\frac{c_{\text{s}}}{\Omega },$$(3)

where $\Omega =\sqrt{G{M}_{★}/r^3}$ is the Keplerian frequency, M_★ is the stellar mass, G is the gravitational constant, and r is the distance between the planet and the star. The radial speed of the gas in the inner regions of the disc is $$u_r=-\frac{3}{2}\frac{v}{r}.$$(4)

We used two different models to calculate the temperature in the midplane of the disc. The first is a purely irradiated disc, where the temperature is given as (Ida et al. 2016) $$T_{\text{irr}}=150\text{\hspace{0.17em}K}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{-1/7}{\left(\frac{L_{★}}{L_{\odot }}\right)}^{2/7}{\left(\frac{r}{\text{AU}}\right)}^{-3/7}.$$(5)

We scaled the stellar luminosity, L_★, as $$\frac{L_{★}}{L_{\odot }}={\left(\frac{M_{★}}{M_{\odot }}\right)}^2,$$(6)

which is typical for young stars (Baraffe et al. 2002, 2015).

In our second model, we also included viscous heating³ of the disc, which we took from Ida et al. (2016) as $$T_{\text{visc}}=317\text{\hspace{0.17em}K}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{3/10}{\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-1/5}{\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-7}{\text{M}}_{\odot }/\text{yr\hspace{0.17em}}}\right)}^{2/5}{\left(\frac{r}{1\text{\hspace{0.17em}AU}}\right)}^{-9/10},$$(7)

where Ṁ_g is the gas accretion rate onto the star. When we include viscous heating, the temperature of the disc is given as max(T_irr, T_vise).

From the continuity equation, the surface density of the gas in the accreting regions of the disc, Σ_g, relates to Ṁ_g as $${\Sigma }_{\text{g}}=-\frac{{\dot{M}}_{\text{g}}}{2\pi r{u}_r}=\frac{{\dot{M}}_{\text{g}}\sqrt{G{M}_{★}}}{3\pi {\alpha }_{\text{v}}ℛT{r}^{3/2}}.$$(8)

For the irradiated case, this gives a surface density of $$\begin{array}{c}{\Sigma }_{\text{g,irr}}\approx 2500\text{\hspace{0.17em}g}/{\text{cm}}^2{\left(\frac{L_{★}}{L_{\odot }}\right)}^{-2/7}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{9/14}{\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-1}\\ \times \left(\frac{{\pi }_{\text{d}}}{2.35}\right)\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-7}{\text{M}}_{\odot }/\text{yr}}\right){\left(\frac{r}{1\text{\hspace{0.17em}AU}}\right)}^{-15/14},\end{array}$$(9)

while in the viscously heated regions, the surface density is given by $$\begin{array}{c}{\Sigma }_{\text{g,visc}}\approx 1325\text{\hspace{0.17em}g}/{\text{cm}}^2{\left(\frac{M_{★}}{M_{\odot }}\right)}^{1/5}{\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-4/5}\\ \times \left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-7}{\text{M}}_{\odot }/\text{yr}}\right){\left(\frac{r}{1\text{\hspace{0.17em}AU}}\right)}^{-3/5}.\end{array}$$(10)

When viscous heating is taken into account, the surface density is given by min(Σ_g,irr, Σ_g,visc).

The time evolution of the disc accretion rate can be written as (Hartmann et al. 1998) $${\dot{M}}_{\text{g}}={\dot{M}}_{\text{g,0}}{\left(\frac{t}{t_{\text{s}}}+1\right)}^{\left(\left(2/2\right)-\gamma \right)/\left(2-\gamma \right)},$$(11)

where γ is the power-law index of the viscosity, as a function of distance to the star, r. The characteristic evolution timescale, t_s, is given by $$t_{\text{s}}=\frac{1}{3{\left(2-\gamma \right)}^2}\frac{R_{\text{d,0}}^2}{v\left(R_{\text{d,0}}\right)},$$(12)

where R_d,0 is the initial disc size. The initial gas accretion rate, Ṁ_g,0, at t = 0, is a free parameter. The estimated initial sizes of protoplanetary discs are still quite uncertain and very model dependent. Using observations of CO, gas disc sizes have been estimated to be of the order of ~100 AU (Ansdell et al. 2018), although recent work by Tabone et al. (2025) suggests that the initial disc sizes are much smaller at ~10 AU. Nevertheless, it is an important parameter, as in typical protoplanetary disc models the disc size determines the characteristic timescale for the evolution of the gas in the disc. Further, for discs with a given mass, a smaller disc would indicate higher accretion rates and therefore higher surface densities, which in turn affects the formation and migration of planets (e.g. Drążkowska et al. 2021). We therefore treat the disc size as a free parameter in this work. Throughout this work, we set the initial accretion rates, Ṁ_g,0, according to the observed power law from Almendros-Abad et al. (2024), who found that the accretion rates scale as $${\dot{M}}_{\text{g}}\propto M_{★}^{1.67}$$(13)

for the Ophiuchus star-forming region, which was the youngest region they considered. We normalised the accretion rate such that a solar-mass star has an initial accretion rate of 10⁻⁷ Ṁ_⊙/yr, which is typical for young solar-mass stars (Hartmann et al. 2016).

The power-law index of the viscosity, γ, was determined by the power-law index of the temperature. In the model when we included viscous heating, we used the power-law index for the irradiated case since the length scale for the disc evolution is the initial disc size, which is typically much larger than the transition radius, r_tran, where the disc becomes viscously heated (Liu et al. 2019). The initial mass of the disc in the irradiated case is given by $$\begin{array}{c}{M}_{\text{d,irr,0}}={\displaystyle {\int }_0^{R_{\text{d,0}}}2\pi r{\Sigma }_{\text{g}}dr=0.1}{M}_{\odot }{\left(\frac{R_{d,0}}{71\text{\hspace{0.17em}AU}}\right)}^{13/14}\left(\frac{{\mu }_{\text{v}}}{2.35}\right)\\ \times {\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-1}\left(\frac{{\dot{M}}_{\text{g,0}}}{{10}^{-7}{M}_{\odot }/\text{yr}}\right){\left(\frac{L_{★}}{L_{\odot }}\right)}^{-2/7}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{9/14},\end{array}$$(14)

In the case of pure viscous heating, the initial disc mass is given by $$\begin{array}{c}{M}_{\text{d,visc}}=0.26M_{\odot }{\left(\frac{M_{★}}{M_{\odot }}\right)}^{1/5}{\left(\frac{{\dot{M}}_{\text{g,0}}}{{10}^{-7}{\text{M}}_{\odot }/\text{yr}}\right)}^{3/5}\\ \times {\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-4/5}{\left(\frac{R_{\text{d,0}}}{71\text{\hspace{0.17em}AU}}\right)}^{7/5}.\end{array}$$(15)

The initial disc mass in the viscously heated case is therefore a combination of Eqs. (14) and (15) $$\begin{array}{c}{M}_{\text{d,full}}={\displaystyle {\int }_0^{R_{\text{d,0}}}2}\pi r{\Sigma }_{\text{g}}dr={\displaystyle {\int }_0^{r_{\text{tran}}}2}\pi r{\Sigma }_{\text{visc}}dr+{\displaystyle {\int }_{r_{\text{tran}}}^{{\text{R}}_{\text{d,0}}}2}\pi r{\Sigma }_{\text{irr}}dr\\ =M_{\text{d,visc}}\left(r_{\text{tran}}\right)+M_{\text{d,irr}}-M_{\text{d,irr}}\left(r_{\text{tran}}\right),\end{array}$$(16)

where M_d,visc(r_tran) is the mass enclosed at r_tran in the fully viscous case and M_d,irr(r_tran) is the same but for the irradiated case. The transition radius between the viscously heated region and the irradiated region is given by setting T_irr = T_visc, and yields $$\begin{array}{c}{r}_{\text{tran}}\approx 4.9\text{\hspace{0.17em}AU}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{31/33}{\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-14/33}\\ \times {\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-7}{M}_{\odot }/\text{yr}}\right)}^{28/33}{\left(\frac{L_{★}}{L_{\odot }}\right)}^{20/33}.\end{array}$$(17)

As mentioned, r_tran,0 ≪R_d,0, which means that M_d,full,0 ≈ M_d,irr,0 so for simplicity, we set the disc mass to be M_d,irr,0.

We also took X-ray photoevaporation of the gas disc into account by following the mass loss relation from Owen et al. (2012) $${\dot{M}}_{\text{w}}=6.25\cdot {10}^{-9}{\text{M}}_{\odot }/\text{yr\hspace{0.17em}}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{-0.068}{\left(\frac{L_{\text{X}}}{{10}^{30}\text{\hspace{0.17em}erg}/\text{s}}\right)}^{1.14},$$(18)

where the X-ray luminosity, L_X, is given by Bae et al. (2013) as $$\log \left(L_{\text{X}}\left[\text{erg}/\text{s}\right]\right)=30.37+1.44\left(\frac{M_{★}}{M_{\odot }}\right).$$(19)

When the photoevaporation rate is larger than the accretion rate onto the star, photoevaporation dominates. When this occurred, we followed Liu et al. (2019) and modified the accretion rate to $${\dot{M}}_{\text{g}}={\dot{M}}_{\text{g,0}}{\left(\frac{t}{t_{\text{s}}}+1\right)}^{\left(\left(5/2\right)-\gamma \right)/\left(2-\gamma \right)}\exp \left(-\frac{t-t_{\text{pho}}}{{\tau }_{\text{pho}}}\right),$$(20)

where t_pho is the time t where Ṁ_g = Ṁ_W. The mass-loss timescale, τ_pho, is given by $${\tau }_{\text{pho}}=\frac{M_{\text{d}}\left(t=t_{\text{pho}}\right)}{{\dot{M}}_{\text{w}}}.$$(21)

We then simply set the disc lifetime to be t_pho + τ_pho. For solar-mass stars with an initial accretion rate of 10⁻⁷ M_⊙/yr, and an initial disc mass of 0.1 M_⊙, this results in a lifetime of 4.8 Myr. We also tested the updated photoevaporation model presented in Ercolano et al. (2023), which calculates how the mass-loss rate depends on stellar mass and X-ray luminosity based on the works by Picogna et al. (2021) and Ercolano et al. (2021), respectively. That model predicts higher values for Ṁ_w and thus lower disc lifetimes for different stellar masses, with the disc lifetime around a solar-mass star being ~2 Myr. Nevertheless, we found little difference in the planet population between the two photo-evaporation models, as the accretion rate and consequently gas surface density decreases significantly at these later stages of the protoplanetary disc. This means that extending the disc lifetime with a few megayears does not have a major effect on the overall planet population. We show these results in Appendix A.

We show disc lifetimes as well as the disc-to-star mass ratios for a range of stellar masses and disc sizes in Fig. 1. Large discs around massive host stars are longer-lived. However, discs larger than 100 AU around stars more massive than ~l.5 M_⊙ are typically gravitationally unstable in their outer regions (marked by a hatched pattern). Further, even if the disc lifetime is ~ 10 Myr, the disc is pretty inactive in forming planets after ~3 Myr due to the low gas surface densities at these late times. The ratio between disc masses and stellar masses ranges from ~l% to ~30%.

From this protoplanetary disc model, it is clear that the disc size plays an important role in setting important physical properties of the disc. For a given accretion rate, the disc size determines the characteristic timescale of the disc, as seen in Eq. (12). Further, the disc size also determines the mass of the disc as seen in Eq. (14). The mass of the gas disc also sets the mass budget of solids in the disc, which in turn affects which planets can form. Therefore, it might be more intuitive to have the disc mass as a free parameter instead of the disc size. However, as the disc size also directly affects the disc accretion rate through t_s, we argue that the disc size is the more natural parameter to vary. Moreover, disc sizes can to some extent be observed directly as opposed to disc masses, which have to be inferred indirectly (see e.g. Manara et al. 2023, and references therein).

Finally, we modelled the inner edge of the disc following Mulders et al. (2015), where the inner edge is approximated to the co-rotation radius of the star given by the expression $$r_{\text{in}}=0.04{\left(\frac{M_{★}}{M_{\odot }}\right)}^{1/3}.$$(22)

Fig. 1 Disc lifetimes (top) and disc-to-star mass ratios (bottom) for discs with varying stellar masses and initial disc sizes. For a fixed disc mass, the disc lifetime increases with increasing stellar mass, while a larger disc also results in longer disc lifetime. Large discs (≳100 AU) around stars more massive than ~1.5 M⊙ are gravitationally unstable in some regions of the discs. This region is marked by a hatched pattern.

2.2 Core accretion

Our pebble accretion model is based on the work of Drążkowska et al. (2021) and is detailed in Nielsen et al. (2023). Throughout this work, we assumed that the star has a solar composition, with the elemental composition based on the work by Magg et al. (2022). The resulting dust-to-gas ratio in the disc is therefore 0.014. We calculated the growth of dust in the protoplanetary disc due to sticking by micrometre-sized monomers. The dust size is limited by fragmentation (both turbulent fragmentation and fragmentation from differential drift) as well as radial drift. We set the fragmentation velocity of pebbles to 2 m/s throughout the entire disc so that the Stokes numbers, St, of pebbles in the outer regions are ~0.01. Typical reported fragmentation velocities range between 1 m/s for silicate particles and 10 m/s for icy particles (e.g. Gundlach & Blum 2015). However, as shown by Nielsen et al. (2025) in their appendix A, varying the fragmentation velocity based on the composition of the dust does not have a significant effect on the final super-Earth population. The core then grows following Eqs. (2)–(5) and (20)–(28) in Nielsen et al. (2023).

Pebble accretion continues until the pebble isolation mass, M_iso, is reached. We calculated the pebble isolation mass following Bitsch et al. (2018), $$\begin{array}{c}{M}_{\text{iso}}=25M_{\oplus }{\left(\frac{H/r}{0.05}\right)}^3\left[0.34{\left(\frac{-3}{\log \left({\alpha }_t\right)}\right)}^4+0.66\right]\\ \times \left(1+\frac{\chi -2.5}{6}\right)\left(\frac{M_{★}}{M_{\odot }}\right),\end{array}$$(23)

where α_t is a parameter that defines the turbulence strength and is set in this work to 10⁻⁴. This α-parameter is distinct from the α-parameter that drives the loss of angular momentum in the disc and subsequent gas flux onto the star, α_v, as it has been shown that these two parameters need not be equal (Lesur et al. 2023). The parameter χ is the negative logarithmic pressure gradient of the gas. In the irradiated regime, this is equal to 39/14 ≈ 2.786 while in the viscous regime, it is equal to 51/20 = 2.55.

2.3 Planet migration

The migration of the protoplanet is driven by the excitation of density wakes as a result of gravitational planet-disc interactions (Kley & Nelson 2012). While it has been shown that the dust in the disc can also induce torques onto the migrating planet, which can cause halted inward migration as well as outward migration if the dust-to-gas ratios in the disc are high enough (Benítez-Llambay & Pessah 2018; Guilera et al. 2025). We chose to neglect these effects for simplicity. Neglecting these effects, type-I migration for a fully irradiated disc can be written as $${\dot{r}}_{\text{I,irr}}=f_{\text{I,irr}}\frac{M_{\text{pl}}}{M_{★}}\frac{{\Sigma }_{\text{g}}{r}^2}{M_{★}}{\left(\frac{H}{r}\right)}^{-2}{\upsilon }_{\text{K}},$$(24)

where υ_K is the Keplerian speed at the location of the planet. The factor, f_I,irr, is a constant that has been fit by D’Angelo & Lubow (2010) to be $$f_{\text{I,irr}}=-2\left(1.36+0.62\frac{\partial \ln {\Sigma }_{\text{g}}}{\partial \ln r}+0.43\frac{\partial \ln T}{\partial \ln r}\right)\approx -\mathrm{4.42.}$$(25)

However, when the disc is viscously heated, a more complex migration law needs to be taken into account. Therefore, when we include viscous heating in our model, we follow Paardekooper et al. (2011), specifically their Eqs. (50)–(53), and calculate the total torque, Γ, onto the planet as the sum of the Lindblad torque, Γ_L and the co-rotation torque, Γ_∞. The migration rate can then be written as $${\dot{r}}_{\text{I}}=\frac{2\Gamma }{M_{\text{pl}}{\upsilon }_{\text{K}}}.$$(26)

This can be rewritten to be in the form of Eq. (24)⁴ but with ƒ_I,irr replaced with ƒ_I,visc, which, similarly to ƒ_I,irr, depends on the surface density- and temperature gradients, but also on α_t and the opacity of the disc. For the opacity of the disc, which enters the migration rate in the viscously heated case, we use the opacity law from Bell & Lin (1994). In Fig. 2, we show the value of the migration coefficient, ƒ_I,visc, as a function of planet mass and distance to the star for a solar-mass star with a disc size of 71 AU at t = 0. Clearly, there is a large region of outward migration within r_tran, which prevents super-Earths from migrating inwards significantly. This can be contrasted to discs heated purely through irradiation where type-I migration is stronger than in the viscous case and always in the direction towards the star, as ƒ_I,irr is in that case negative with a value of −4.42 compared to ƒ_I,visc, which varies between −4 and 2 depending on the planet mass and distance to the star. Further, inside of r_tran, the surface density of the star is by definition lower in the viscous case compared to the purely irradiated case, which means that migration is even slower in this region in the viscous case compared to the irradiated case.

Should the planet grow large enough, the planet can start to carve out a gap in the disc (Kanagawa et al. 2018; Paardekooper et al. 2023). This results in a transition towards type-II migration, where the depth of the gap in the disc affects the rate at which the planet migrates inwards. We included type-II migration by modifying our migration equation following Johansen et al. (2019), $$\dot{r}={\dot{r}}_{\text{I}}\frac{1}{1+{\left[M_{\text{pl}}/M_{\text{gap}}\right]}^2},$$(27)

where M_gap is the mass at which the planet starts to carve a gap into the surrounding gas, which we set to be equal to 2.3 times M_iso following the approach of Johansen et al. (2019). This definition of the gap transition mass ensures that it never exceeds the pebble isolation mass, which otherwise might be possible at the low levels of turbulence that we assume in this work. Additionally, how the pebble isolation mass scales with α_t is currently uncertain, see Section 3.1 in Johansen et al. (2019) for more detailed discussions. As such, we chose a more simplified definition of the gap transition mass. Finally, we allowed planets to grow and migrate until they reached the inner edge, r_in, of the disc, after which migration and growth was halted.

Fig. 2 Migration coefficient, ƒI,visc, as a function of distance to the star and planet mass for a solar-mass star with a disc size of 71 AU at t = 0. Inside of the transition radius (~5 AU), the disc is viscously heated, and for certain masses, the planet experiences outward migration. The sharp switch to inward migration at ~0.2 AU is caused by the significant decrease in opacity due to the evaporation of silicate and metal grains, which occurs at ~1000 K. We also show ƒI,visc at different times and for different stellar masses and disc sizes in Appendix B.

2.4 Gas accretion

After the pebble isolation mass was reached, the planet began to cool and accrete gas. We followed Bitsch et al. (2015b), who derived an analytical model based on the work of Piso & Youdin (2014), to yield $$\begin{array}{c}{\dot{M}}_{\text{atm,\hspace{0.17em}KH}}=1.75\times {10}^{-4}{M}_{\oplus }{\text{yr}}^{-1}{f}^{-2}{\left(\frac{{\kappa }_{\text{env}}}{1{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}\\ \times {\left(\frac{{\rho }_c}{5.5\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}}\right)}^{-1/6}{\left(\frac{M_{\text{c}}}{M_{\oplus }}\right)}^{11/3}{\left(\frac{M_{\text{atm}}}{M_{\oplus }}\right)}^{-1}{\left(\frac{T_{\text{d}}}{81\text{\hspace{0.17em}K}}\right)}^{-1/2}\end{array},$$(28)

where κ_env is the envelope opacity, which we took to be 0.05 cm² g⁻¹. The numerical factor, f, is included to match the analytical results to the numerical and is set to 0.2. This analytical gas accretion rate is only valid when the planet is not undergoing runaway gas accretion, i.e. when M_atm < M_c. In the case of M_atm ≥ M_c, we followed Machida et al. (2010) who estimated the rapid gas accretion rate from 3D hydrodynamical simulations $${\dot{M}}_{\text{atm,\hspace{0.17em}rapid}}=\{\begin{array}{ll}0.83\Omega {\Sigma }_{\text{g}}{H}^2{\left(\frac{R_{\text{H}}}{H}\right)}^{9/2}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{R_{\text{H}}}{H}<0.3\hfill \\ 0.14\Omega {\Sigma }_{\text{g}}{H}^2\hfill & \text{if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{R_{\text{H}}}{H}\ge \mathrm{0.3.}\hfill \end{array}$$(29)

These accretion rates are limited by the amount of gas that can be supplied to the planet from the protoplanetary disc. Ida et al. (2018) found that the gas can enter the Hill sphere at a rate of $$\begin{array}{c}{\dot{M}}_{\text{atm,\hspace{0.17em}hill}}=1.5\times {10}^{-3}{M}_{\oplus }{\text{yr}}^{-1}{\left(\frac{H/r}{0.05}\right)}^4\\ \times {\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{-1}\left(\frac{{\dot{M}}_{\text{g}}}{{10}^{-8}{M}_{\odot }{\text{yr}}^{-1}}\right)\frac{1}{1+{\left(M/M_{\text{gap}}\right)}^2}.\end{array}$$(30)

Finally, the gas accretion rate cannot exceed the gas flux that flows through the disc, Ṁ_g. We therefore set the gas accretion rate to be $${\dot{M}}_{\text{atm}}=\{\begin{array}{lll}\min \left({\dot{M}}_{\text{atm,KH}},{\dot{M}}_{\text{atm,Hill}},{\dot{M}}_{\text{g}}\right)\hfill & \text{if}\hfill & M_{\text{atm}}<M_{\text{c}}\hfill \\ \min \left({\dot{M}}_{\text{atm,rapid}},{\dot{M}}_{\text{atm,Hill}},{\dot{M}}_{\text{g}}\right)\hfill & \text{if}\hfill & M_{\text{atm}}\ge M_{\text{c}}.\hfill \end{array}$$(31)

3 Results

We injected our planets by sampling the beginning location from a log-uniform distribution between 0.5 AU and the initial disc size, R_d,0. The injection time was sampled from a uniform distribution between 10³ yr and the disc lifetime, t_d. The initial protoplanet mass was set to be 0.01 M_⊕, which allowed us to ignore Bondi accretion and only work within the 2D and 3D Hill regimes of pebble accretion (Drążkowska et al. 2023). We thereby assumed that embryos form rapidly through pebble accretion (Lyra et al. 2023) or through a combination of planetesimal accretion and pebble accretion (Lorek & Johansen 2022). In all simulations throughout this work, we injected 5000 proto-planets. In other words, for each protoplanetary disc, we sampled 5000 combinations of injection locations and injection times. We defined super-Earths as those with masses between 1 and 10 M_⊕ to facilitate comparison with previous theoretical work (e.g. Pan et al. 2024). We also explored different definitions of close-in planets. The upper semi-major axis limit for planets around M-dwarfs in Mulders et al. (2015) is 0.2 AU. In comparison, Pan et al. (2025) defined close-in planets as planets with orbital periods of less than 100 days. We compared our results using both of these definitions and found no significant difference between them. Therefore, we chose to define close-in planets as those planets with semi-major axes less than 0.2 AU. We show the resulting semi-major axes and masses of 5000 planets injected in two example discs with sizes of 50 AU and 150 AU, respectively, around solar-mass stars in Fig. 3. Here, the disc is only heated through stellar irradiation. We also show our definitions of close-in super-Earths and giants in blue and red, respectively. Even with rapid inwards type-I and type-II migration in the disc, we can form Jupiter-mass planets as far out as 4 AU if the protoplanetary disc is large enough. The formation of giants on wider orbits can also be facilitated if the pebble flux in the outer regions decays significantly, which allows the planet to cool sufficiently and thus starts accreting gas (Gurrutxaga et al. 2024); however, we ignore that possibility here.

3.1 Exploring the parameter space

In this section, we explore the parameter space of disc sizes and stellar masses. For a given disc size and initial accretion rate, we calculated the gas mass of the disc from Eq. (14). In Fig. 4, we show the resulting close-in super-Earth fraction⁵ as a function of stellar mass and disc size.

In the irradiated case, the fraction of close-in super-Earths increases with increasing stellar mass for all disc sizes. For a given stellar mass, disc size has a minor effect on the super-Earth fraction as even small discs form super-Earths. Clearly, the irradiated case is unable to reproduce the observed super-Earth trend, which peaks at ~0.5 M_⊙. When considering viscous heating, however, we find a marked decrease with increasing stellar mass in the resulting super-Earth fraction for nearly all disc sizes. For the largest discs (90–150 AU), the super-Earth fraction peaks at 0.5 M_⊙ while for smaller discs, the peak is shifted towards solar-mass stars.

In order to understand the cause of the reduction of the close-in super-Earth fraction around high-mass stars, we investigated some example growth tracks of individual planets. We injected the protoplanet at r = 3 AU and t = 10⁴ yr around stars with a range of disc sizes and masses and show the growth tracks in Fig. 5. When the disc is heated solely through irradiation, inward migration is very efficient and the planet grows into a super-Earth with a final orbit at the inner edge for all stellar masses and disc sizes. In the viscous case, similar results are found for disc sizes of 30 and 100 AU around a star with a mass of 0.5 M_⊙. In the case of a disc with a size of 150 AU and M_★ = 0.5 M_⊙, the planet has just enough time to grow into a Saturn-mass planet. Around a solar-mass star, only a small disc size yields a super-Earth as now the planet experiences outward migration, delaying the migration to the inner edge, which allows the planet core to become massive enough for the planet to grow into a giant. for a star with a mass of 1.5 M_⊙, the injected embryo manages to grow into a giant planet for all disc sizes, as all embryos experience outward migration and rapid core growth. Further, the final planet mass increases with increasing disc size. For small discs or for low-mass stars, the region of outward migration moves inwards so fast that planets experience very little outward migration and therefore migrate inwards rapidly, halting their growth into a giant. The reduction in the close-in super-Earth fraction reported, for example, by Mulders et al. (2015) and Zink et al. (2023) can thus potentially be explained if some stars host large, viscously heated protoplanetary discs.

We further show the giant planet fraction as a function of stellar mass for different disc sizes in Fig. 6. We include all giant planets here and not just close-in giants, which means that we divide the number of giants with the total number of injected embryos. This is due to type-II migration being much slower and decreasing with increasing planet mass. As a result, close-in super-Earths become giant planets with final orbits farther away from the star when the size of the disc or the mass of the host star is increased. This can be seen, for example, in Figure 5. The giant planet fraction increases with increasing stellar mass and disc size for both disc models, in agreement with observations, which show that the giant planet occurrence rate peaks for stars with mass ~2 M_⊙ (Reffert et al. 2015; Ghezzi et al. 2018).

Estimated giant planet occurrence rates vary significantly in the literature depending on the mass range and semi-major axis range considered. Fulton et al. (2021) found that the giant planet occurrence rate for planets within 1–5 AU and masses above 100 M_⊕ (~0.3 M_J) is ~ 15–20%, while Fernandes et al. (2019) reported an occurrence rate of 26.6% for planets between 0.1 and 100 AU and masses of 0.1–20 M_J. It is important to note that the work of Fernandes et al. (2019) considers the entire stellar mass range of giant hosts observed by Kepler, which includes stars with different masses than the Sun, which can inflate the occurrence rate. Reffert et al. (2015) investigated the occurrence rates of giant planets as a function of stellar mass and metallicity and found that it increases with increasing stellar mass up to 2 M_⊙, with reported occurrence rates ranging from 3–8% between stellar masses of 1 M_⊙ and 2 M_⊙. Ghezzi et al. (2018) found that the giant planet occurrence rate ranges from ~ 10–20% in the same stellar mass range. The discrepancy between these two studies is likely caused by the fact that the smallest confirmed minimum mass used by Reffert et al. (2015) was 2.3 M_J, while Ghezzi et al. (2018) considered planets down to masses of 0.5 M_J. Nevertheless, the observed giant planet occurrence rates for similar mass ranges that we consider seem to vary between 10–20% for solar mass stars, increasing to at least 20% for stars with masses of 2 M_⊙.

In this work, in the case of irradiation heating only, migration is efficient enough to keep the giant planet fraction low, in the range ~0.5–3% for solar-mass stars and increases to 2–10% for stars with masses of 2 M_⊙. Discs smaller than 60 AU for all stellar masses fail to form any giants due to their rapid evolution and low masses. When viscous heating is included and inward migration is delayed, the giant planet fraction increases by a factor of ~5 for most stellar masses and protoplanetary disc sizes. These results are in better agreement with observed occurrence rates for giant planets around both solar mass stars and high-mass stars. However, we caution against comparing our planet fractions directly to occurrence rates due to model uncertainties. We discuss the interpretations of our planet fraction more in detail in Sect. 4.2.

Further, in the viscous case, giant planets can form even around small discs (~30 AU) around solar-mass stars and above, which, in our model, is not possible in the irradiated case. Neither model produces any giants for stellar masses of M_★ ≲ 0.5 M_⊙, which is in rough agreement with previous work modelling giant planet formation around varying stellar masses (Liu et al. 2019; Shibata & Helled 2025).

Fig. 3 Example of the resulting semi-major axes and masses of planets injected around a solar-mass star for two different initial disc sizes, 150 AU (top) and 50 AU (bottom). Here, only irradiation heating is considered. We indicate the planets that we define as close-in super-Earths (blue) and giants (red). The large disc forms significantly more giant planets due to having more available mass and a longer lifetime.

Fig. 4 The fraction of close-in super Earths (r < 0.2 AU) for a range of stellar masses (x-axis) and disc sizes (colours). In all cases, no super-Earths are formed around stars with a mass of 0.1 M⊙. When only considering discs heated by irradiation (left), the fraction of super-Earths increases with increasing stellar mass regardless of disc size. In the viscously heated case (right), the fraction of super-Earths reaches a maximum at 0.5 M⊙ and then declines slightly for large discs (90 AU or larger). For smaller discs (30 and 60 AU), the super-Earth fraction peaks at 1 M⊙, while for larger discs it peaks at 0.5 M⊙.

Fig. 5 Growth tracks of a single planet injected at r = 3 AU and t = 104 yr. In the irradiated case (left), this results in the formation of a super-Earth at the inner disc edge for all stellar masses and disc sizes. In the viscous heating case (right), the protoplanet grows into a giant planet for high stellar masses as well as for solar-mass stars with sufficiently large discs.

Fig. 6 Giant planet fraction at all semi-major axes as a function of stellar mass and disc size. With or without viscous heating, the giant planet fraction increases with increasing stellar mass for all disc sizes where a giant planet can form to begin with. The lower stellar mass limit for giant planet formation is 0.5 M⊙ in both models. In the irradiated case (left), the lower disc size limit for giant planets to form is 60 AU (90 AU for M★ = 0.5 M⊙). In contrast, when viscous heating is included, the lower disc size limit is as low as 30 AU for solar-mass stars and above, while it is 60 AU for stars with M★ = 0.5 M⊙.

3.2 Using a fixed injection time

In our nominal model, we sampled the injection time uniformly between 10³ yr and the disc lifetime. As a consequence, discs with long lifetimes are negatively affected as planets injected late into the disc (after ~2–3 Myr) struggle to grow to pebble isolation mass. Therefore, we also investigated the effects of different sampling schemes for the injection times.

In Fig. 7 we show the close-in super-Earth fraction using a fixed injection time at t = 10⁵ yr. We also tested using a fixed injection time at t = 10³ yr and t = 10⁴ yr and found no significant differences between the three cases. Further, we tested the effects of sampling the injection time from a log-uniform distribution between each of these times and the disc lifetime and found no significant difference from the case with a fixed injection time. In the irradiated case, we find the opposite trend as in our nominal model. As large discs typically have longer life-times and are more massive than small discs, we find that large discs produce fewer super-Earths for a fixed stellar mass as more planets have time to grow into a giant planet instead of a super-Earth. In the case of irradiated heating, the super-Earth fraction decreases beyond 0.5 M_⊙ for all disc sizes with the exception of very small discs (30 AU) where the super-Earth fraction continuously increases with increasing stellar mass as we find that no giant planets form for such small discs. For discs larger than 30 AU, we find that the fraction of giant planets increases enough with increasing stellar mass and disc size that the super-Earth fraction is reduced for increasing stellar masses above 0.5 M_⊙.

For example, for a solar mass star, the giant planet fraction ranges between ~10% and ~25% when the injection is fixed for disc sizes between 60 and 150 AU. When sampling injection times throughout the entire disc lifetime, the giant planet fractions instead range between ~0.5% and ~5%, which is not high enough to impact the super-Earth fractions.

In the viscous case, we find that the super-Earth fraction is heavily suppressed around large discs while for small discs it peaks at 0.5 M_⊙ and decreases significantly for host stars more massive than solar-mass. For discs with sizes of 30 AU, the super-Earth fraction peaks at M_★ = 1.5 M_⊙. As all planets are injected early on in the disc lifetime, the outward migration region is located far from the star when the planets start to grow and more planets are affected by some level of outward migration, allowing them to reach higher pebble isolation masses and grow to become giants. The fraction of giant planets approaches unity at solar masses and above for discs larger than 30 AU, which is the cause of the suppressed super-Earth fraction. Further, even the smallest discs are able to form giant planets for high enough stellar masses.

For a given system with a fixed injection time, the only parameter determining the final mass of the planet is the injection location. In order for a giant planet to form, the planet needs to be injected sufficiently far out in order for it to be able to reach a large enough pebble isolation mass to accrete a substantial amount of gas, but not too far out as the core will then grow too slowly to reach the pebble isolation mass at all. Planets injected outside of this ‘giant planet region’ will thus not have time to migrate inwards fast enough during the disc lifetime. Planets injected interior of the giant planet region, in contrast, grow to become Earths or super-Earths instead with only minimal gas accretion. To illustrate this distinction, we show the final mass as a function of injection location for the model with fixed injection times in Fig. 8. In both models, no super-Earths are formed around host stars with a mass of 0.1 M_⊙. In the irradiated case for high stellar masses, the final planet mass increases with increasing injection location as planets take longer to migrate to the inner edge and reach higher pebble isolation masses. Planets grow into giants if they are injected far out in the disc at r ≳20 AU. This is true for all disc sizes above 30 AU and stellar masses above 0.5 M_⊙. For stars with a mass of 0.5 M_⊙ and small discs, the final planet mass decreases if the planet is injected beyond 10 AU.

With viscous heating included, giants start to form for as low stellar masses as 0.5 M_⊙ for most disc sizes and injection locations due to outward migration, which allows planet cores to become more massive. The giant planet-forming region covers almost the entire parameter space for disc sizes of 120 AU and above for all stellar masses above 0.1 M_⊙, resulting in a severe reduction in super-Earth formation at these disc sizes. For small discs (30 AU), the giant planet-forming region is small for stars with masses of 0.5 M_⊙ and 1 M_⊙, resulting in high super-Earth fractions while for larger stellar masses, the giant planet-forming region starts to cover the full parameter space, reducing the super-Earth fraction significantly.

Fig. 7 Similar to Fig. 4 but using a fixed injection time at t = 105 yr. We also show the giant planet fraction (triangles). In both models, super-Earth formation is suppressed for high stellar masses, as planets have time to grow into giant planets instead. As a reminder, the giant planet fraction is defined as the number of giant planets formed compared to the total number of planets injected, while the super-Earth fraction is the number of close-in super-Earths compared to the number of close-in planets, as described in Section 3.1.

Fig. 8 Final planet mass as a function of injection location for different stellar masses and disc sizes, for a fixed injection time of t = 105 yr. For low-mass stars, neither disc model produces any super-Earths or giant planets. In the viscous model, the region where giant planets form moves close to the star for larger discs and stellar masses, resulting in fewer super-Earths forming. In the irradiated case, the giant-forming region remains far out in the disc and only moves slightly closer to the star for higher stellar masses.

Fig. 9 Fraction of close-in super-Earths as a function of stellar mass with sampled initial accretion rates and disc sizes. The black lines show the median values out of all runs, while the error bars show the 16th and 84th percentiles. Similarly to Fig. 4, the super-Earth fraction increases with increasing stellar mass in the irradiated case (left). When the disc is viscously heated (right), the super-Earth fraction flattens out for host stars more massive than 0.5 M⊙.

3.3 Sampling disc sizes and accretion rates

In the previous section, we scaled the initial accretion rates with stellar mass and explored the (M_★,R_d,0)-space. This meant that the initial accretion rate was assumed to be fixed for a given stellar mass. However, as seen by, for example, Manara et al. (2012) and Almendros-Abad et al. (2024), the observed accretion rate exhibits a large scatter around the linear fit. As discussed by Manara et al. (2012), this scatter might be caused by the observed stars having different ages, as the accretion rate declines over time. Although the scatter remains even when binning the stars by their age, it is therefore possible that the scatter in the accretion rate is either intrinsic or dependent on some other, unknown variable. Further, it is also possible for some of the discs in the parameter space to be gravitationally unstable in some regions of the parameter space. To account for these effects, we turned to a sampling scheme in order to synthesise a more realistic population of discs around stars with different stellar masses and explore the resulting super-Earth fraction.

For a given stellar mass, we sampled the initial accretion rate from a Gaussian centred around the power law in Eq. (13) with a width of σ = 0.2. The disc size was then sampled from a uniform distribution. The minimum disc size was set to 20 AU, while the maximum disc size was set to the smallest value between 150 AU and the largest size at which the entire disc is gravitationally stable. We performed this sampling 50 times. Gravitational stability is determined through the Toomre Q parameter, $$Q=\frac{c_{\text{s}}\Omega }{\pi G{\Sigma }_{\text{g}}}.$$(32)

Here we calculate Σ_g and c_s using an irradiated temperature profile as r_tran ≪ R_d,0. As the disc is gravitationally stable for Q > 1, the maximum possible initial disc size for a gravitationally stable disc is given by $$\begin{array}{c}{R}_{\text{d,max,0}}=173\text{\hspace{0.17em}AU}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{-1/3}{\left(\frac{L_{★}}{L_{\odot }}\right)}^{2/3}{\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{14/9}\\ \times {\left(\frac{{\dot{M}}_{\text{g,0}}}{{10}^{-7}{M}_{\odot }/\text{yr}}\right)}^{-14/9}{\left(\frac{{\mu }_{\text{d}}}{0.01}\right)}^{-7/3}\\ \approx 173\text{\hspace{0.17em}AU}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{-1.6}{\left(\frac{{\alpha }_{\text{v}}}{0.01}\right)}^{14/9}{\left(\frac{{\mu }_{\text{d}}}{2.35}\right)}^{-7/3},\end{array}$$(33)

where the last approximation comes from our specific luminosity and gas accretion rate scalings in Eqs. (6) and (13), respectively. From this, we find that the largest gravitationally stable disc size decreases with increasing stellar mass. Therefore, the range of sampled disc sizes decreases with increasing stellar mass, which means that discs around larger stars are typically smaller than those around low-mass stars when sampled from a uniform distribution. We show the resulting super-Earth fractions in Fig. 9 for both the irradiated and viscous case. As the disc sizes only marginally affect the resulting super-Earth fraction in the irradiated case, we find similar results as when we explored the full parameter space. Due to the sampling of initial accretion rates, however, we find that there is a greater variance in the resulting super-Earth fraction around low-mass stars, with small discs and low initial accretion rates resulting in slightly lower super-Earth fractions.

In the viscous case, we find that the super-Earth fraction peaks at ~0.5 M_⊙ and plateaus for higher stellar masses, although the variance is large for all stellar masses. As the discs around high-mass stars are typically smaller, the super-Earth fraction remains relatively high for high stellar masses (>l M_⊙) since the viscously heated region in small discs moves inwards faster than around large discs, resulting in more efficient inward migration and fewer giant planets forming. Notably, for stellar masses of 0.5 M_⊙ and below, the fraction of super-Earths increases with increasing disc size for a given stellar mass while the opposite relationship is found for higher stellar masses. As shown in Figure 6, giant planets start to form for stellar masses of 0.5 M_⊙ and above. As a result, for higher stellar masses when giant planet formation is possible, increasing the disc size results in even more giant planets forming, reducing the fraction of super-Earths forming. This can also be seen when comparing the right panels of figures 4 and 6. For stellar masses above 0.5 M_⊙, the super-Earth fraction decreases with increasing disc size for a given stellar mass, while the giant planet fraction increases for increasing disc size.

In Fig. 10 we show the super-Earth fraction as a function of stellar mass for the sampled discs but using a fixed injection time at t = 10⁵ yr. When the disc is heated purely through irradiation, migration is efficient enough that only protoplanets injected far out in the disc will grow into a giant planet, even when the protoplanets are injected early on in the disc lifetime. Giant planet formation is nevertheless more efficient than in our nominal case due to the early injection time. Therefore, a minor reduction in the median super-Earth fraction with increasing stellar mass is found.

When viscous heating is included, we find similar results as in Fig. 7 where super-Earth formation is heavily suppressed for most sampled discs, since giant planet formation starts to dominate at solar mass and above. For large discs around low-mass stars, we also find that giant planet formation becomes efficient enough to dominate over super-Earth formation. For most sampled discs, no matter where the planet is injected, it has enough time to grow into a giant planet as outward migration is efficient in the early lifetimes of these discs, resulting in fewer super-Earths formed.

As mentioned, the sample range of disc sizes around high-mass stars is smaller than that of low-mass stars due to gravitational instability, which leads to smaller disc sizes. As a result, the median super-Earth fraction increases slightly around stars with masses of ~2 M_⊙. However, this increase is marginal as even small discs predominately form giant planets when viscous heating is included and a fixed, early, injection time is used.

Fig. 10 Same as Fig. 9 but with a fixed injection time of 105 yr. Large discs form more giants when the embryos are injected early, which results in fewer super-Earths being formed. In contrast, if no giants are able to form, such as when the host star mass is below 0.5 M⊙, the super-Earth fraction increases with increasing disc size. We find that the median super-Earth fraction decreases with increasing stellar mass, both with or without viscous heating, although the effect is significantly stronger when viscous heating is active. Further, we find that the median super-Earth fraction increases slightly above 1 M⊙ in the viscous case as the sampled disc sizes around these high-mass stars are slightly smaller.

4 Discussion

4.1 Comparing with previous work

Considering planetesimal accretion as the main planet formation pathway, Burn et al. (2021) performed planet population synthesis for stars with masses below solar mass. They found that giant planet formation occurs for stellar masses of at least 0.5 M_⊙ and that the occurrence rates of giant planets increase with increasing stellar masses, in line with the results we have found in this work. Further, they found that the super-Earth occurrence rate increased with increasing stellar mass up to masses of ~0.5 M_⊙, after which the trend plateaus, despite including outward migration in their work. These discrepancies compared to our work could be due, for example, to their relatively low disc masses, which result in less efficient giant planet formation and consequently less inhibition of super-Earth formation. Moreover, giant planet formation has been found to be less efficient for planetesimal accretion models compared to pebble accretion, especially in the outer regions of the protoplanetary disc, (e.g. Lorek & Johansen 2022). If giant planet formation is inefficient, super-Earth formation never gets suppressed and continues to increase with increasing stellar mass.

Recently, Pan et al. (2025) performed planet population synthesis using pebble accretion and similar disc models to this work, including viscous heating and outward migration. They found that the occurrence rates of super-Earths increased for increasing stellar masses up to ~0.5 M_⊙, after which it decreases, similar to the models in this work where viscous heating and outward migration is included. Further, they found that giant planets could still form for stellar masses below 0.5 M_⊙, in contrast to this work and the work of Burn et al. (2021), due to the effects of efficient pebble accretion and planet-planet collisions, which aids in the formation of massive planetary cores.

4.2 Interpreting planet fractions

Throughout this work, we determined the fraction of planets by calculating the number of planets in our defined mass ranges compared to the total number of planets injected. While it is tempting to try to compare these values with observed occurrence rates of super-Earths and giant planets, we caution against this due to large uncertainties in the choice of model and simulation setup.

One major reason for this caution is that the final planet fractions are strongly affected by how the injection time and location are sampled. For instance, in our nominal model without any viscous heating, the giant planet fraction varies between ~0.5%–5% and ~2%–10% for stars with masses of 1 M_⊙ and 2 M_⊙, respectively. This can be contrasted to when we fix the injection time to 10³ yr, where the giant planet fraction varies from ~ 10–25% up to ~20–35%.

Further, our super-Earth fractions in particular are sensitive to our model choices, especially the lack of N-body interactions, which we neglected in favour of computational speed. Including multiple planets in our model could cause the planetary system to become unstable and super-Earths to be ejected when, for example, a giant planet migrates inwards to the inner regions of the disc (e.g. Mustill et al. 2015; Bitsch & Izidoro 2023). For further discussion on this, see Sect. 4.5.

As a result, there is no direct translation possible between our resulting planet fraction and observed occurrence rates of both giant planets and super-Earths. Instead, our resulting planet fractions should be interpreted primarily as the probability that a certain type of planet forms in isolation around a given star, given a set of initial conditions and model assumptions, rather than the fraction of stars that host that type of planet.

Despite this, we argue that the trends for the dependence of our planet fractions with respect to stellar mass give useful insight into the underlying causes for the observed occurrence rate. For example, a lower super-Earth fraction does indicate that super-Earths are less likely to form, which in turn would result in a reduced occurrence rate, but the absolute values cannot necessarily be directly converted into each other.

4.3 Viscous heating of discs

Whether or not discs are viscously heated is uncertain. When non-ideal magnetohydrodynamics (MHD) are taken into account, simulations have found that viscous heating mostly takes place at several scale heights above the midplane as the magnetorotational instability (MRI) becomes suppressed in the midplane (Mori et al. 2019). In these cases, the midplane temperature is dominated by stellar irradiation. As shown in this work, the temperature structure of the disc can have a significant impact on the migration of planets in the disc, and understanding to what extent discs are viscously heated is of great significance in order to understand the properties and origins of the resulting planet population.

When we assume that viscous heating takes place in the midplane, the two main parameters that determine the temperature are α_v and Ṁ_g, which are both related to each other. The value of α_v for discs is currently estimated to be between 10⁻⁴ and 10⁻² (Pinte et al. 2016; Najita & Bergin 2018; Trapman et al. 2020; Villenave et al. 2022; Rosotti 2023). For lower values of α_v, the initial accretion rate also needs to be lower in order to achieve realistic disc sizes for a fixed disc mass. For example, for α_v = 10⁻³ and Ṁ_g,0 = 10⁻⁷ M_⊙/yr, the resulting initial disc size is ~6 AU for a fixed disc mass of 0.1 M_⊙, which we consider unrealistic as dust disc sizes of young discs have been observed to be ≳ 10 AU (Ohashi et al. 2023) and are expected to be smaller than the gas disc due to pebble drift (e.g. Ansdell et al. 2018). If instead Ṁ_g,0 = 10⁻⁸ M_⊙/yr, the disc size remains at 71 AU. In order to test the effects of α_v, we ran the same simulations as in Sect. 3.1 but with α_v = 10⁻³ and Ṁ_g,0 = 10⁻⁸ M_⊙/yr. We show the resulting super-Earth fractions in Fig. C.1 in Appendix C and find that the resulting super-Earth fractions are qualitatively similar to our nominal model. In the low viscosity case, the transition radius, r_tran, between the viscously heated region and the irradiated region at t = 0 is closer to the star at 1.8 AU compared to 4.9 AU in our nominal model. This is a consequence of the lower accretion rates, which results in weaker viscous heating. When the transition radius is closer to the star, inward migration is more efficient and the reduction of close-in super-Earths is less significant around high-mass stars. However, the initial surface density at t = 0 is also higher in the viscous regime, meaning outward migration is more efficient. This slightly reduces the effects of the efficient inward migration and we find that the fraction of close-in super-Earths still decreases with increasing stellar masses for large, viscously heated discs. Further, migration (both inwards and outwards) in both the irradiated and viscous case remains fast for a longer time for low α_v, since the surface density of the low viscosity disc remains high over a longer timescale. As a result, inward migration is more efficient in irradiated discs and the fraction of close-in super-Earths peaks at a stellar mass of ~0.5 M_⊙ and decreases only slightly for higher stellar masses as migration is efficient enough that only a few planets grow beyond the super-Earth range. This is found for all disc sizes as disc size does not have a significant effect in the irradiated case.

4.4 Varying fragmentation velocities

Throughout this work, we used a constant fragmentation velocity of 2 m/s throughout the entire disc. This resulted in St ~ 0.01 in the outer disc, which is consistent with observations (Jiang et al. 2024). However, fragmentation velocities of dust in protoplanetary discs are highly uncertain, with experiments finding fragmentation velocities to vary between 1 m/s and 10 m/s depending on the composition of dust (e.g. Gundlach & Blum 2015). Previous work has therefore often chosen to vary the fragmentation velocity such that dry pebbles inside the water iceline have fragmentation velocities of 1 m/s while icy pebbles outside the water iceline have fragmentation velocities of 10 m/s (e.g. Drążkowska & Alibert 2017; Venturini et al. 2020). In previous work, Nielsen et al. (2025) showed that using a varying fragmentation velocity had little to no effect on the super-Earth population around a solar-mass star. However, considering that the location of the iceline is dependent on the temperature of the disc and therefore on the stellar mass, it is possible that varying the fragmentation velocity might play a bigger role around stars with different stellar masses. We therefore tested the effects of different fragmentation velocities on the close-in super-Earth population. We ran the same set of simulations as in Sect. 3.1 but with fragmentation velocities of 1 m/s and 10 m/s for dry and icy particles, respectively. We show the results in Fig. C.2 in Appendix C. In the irradiated case, the results are relatively unchanged while in the viscous case, the super-Earth fraction is even more reduced for higher stellar masses for all disc sizes. This is a result of the water iceline moving outwards for higher stellar masses, which in turn results in the inefficient, dry region, with small pebbles, extending to larger parts of the disc, resulting in fewer super-Earths being formed. Further, the giant planet formation beyond the water iceline is enhanced due to the larger pebbles, which in turn reduces the super-Earth fraction farther.

4.5 Including N-body effects

In this work we only included one protoplanet per disc, and the injected planets do not interact with each other. Including N-body effects could have a significant impact on the final planet population, as planet-planet collisions have shown to be able, for example, to produce giant planets around low-mass stars (Pan et al. 2024). As multiple super-Earths migrate inwards through the disc, they can get trapped in resonant chains, which become unstable over time and break, causing planet-planet collisions (e.g. Izidoro et al. 2021; Esteves et al. 2022; Kajtazi et al. 2023). Given that the Hill radii of planets of a given mass and semi-major axis are higher around low-mass stars, it is therefore natural to expect that planet-planet collisions can be more prevalent around low-mass stars. Further, as inward migration is efficient around low-mass stars, especially in irradiated discs, the enhanced close-in super-Earth fraction could be explained by efficient inward migration and subsequent collisions (e.g. Pan et al. 2025).

Due to their shorter lifetimes on the main-sequence, most high-mass stars that are observed today are young. From galactic chemical evolution, young stars are typically more metal-rich (Feuillet et al. 2018; Delgado Mena et al. 2019). As a result, observed high-mass stars are on average more metal-rich (see e.g. Fig. 2 in Nielsen et al. 2023), and are more likely to host giant planets. The presence of inner super-Earths has been found to be positively correlated with the presence of an outer giant for stars with super-solar mass and super-solar metallicity (Bryan & Lee 2024, 2025; Bonomo et al. 2025). This indicates that if a single giant planet forms and remains on a wide orbit, any potential inner planets typically remain undisturbed by the giant. However, if instead multiple giants form in the system, these can cause instabilities, which in turn can eject inner super-Earths from the system, provided that the giants form or migrate close enough to the super-Earths (Mustill et al. 2015; Bitsch & Izidoro 2023; Kong et al. 2024). Therefore, in order to understand the effects of stellar mass on the inner super-Earth population, it will be necessary to consider the formation of multiple planets and their gravitational interactions in a given system.

5 Conclusions

In this work, we simulated the formation of planets around different stellar masses using two different disc models with either purely irradiated heating or both viscous and irradiated heating. We find the following:

1. When the disc is heated purely through irradiation, type-I migration is exclusively inwards and very efficient, resulting in rapid inward planet migration and an increase in the close-in super-Earth fraction for higher stellar masses These results are roughly independent of disc size;

2. When viscous heating is included, the direction of type-I migration is flipped in the inner disc for a range of relevant planet masses, which delays inward migration slightly. This delay results in a significant fraction of planets growing into giants instead. The transition radius between inwards and outward migration is dependent on stellar mass and moves closer to the star over time. Around small discs or low-mass stars, this region of outward migration moves close to the star quickly enough that planets experience rapid inward migration, which halts their growth into giants. In contrast, larger discs evolve over a longer timescale, which results in inward migration being delayed for longer, allowing more giants planets to form;

3. As discs with larger sizes or those around high-mass stars have longer lifetimes, giant planet formation is negatively affected when sampling injection times over the entire disc lifetime. This is a result of planet growth being inefficient in the later stages of the disc. When fixing the injection time to 10⁵ yr, we find that the super-Earth fraction is reduced around high-mass stars in both the irradiated and viscous case. In the irradiated case, more planets around high-mass stars grow into giants compared to when sampling injection times, resulting in a slightly reduced super-Earth fraction around high-mass stars. In the viscous case, a majority of planets around solar-mass stars and above grow into giants when injected this early due to outward migration, resulting in a highly suppressed super-Earth fraction around high-mass stars and in large discs. We also tested sampling the injection time from a log-uniform distribution and found little difference compared to using a fixed injection time;

4. When sampling disc parameters from a distribution, discs are typically smaller around high-mass stars. As a result, slightly fewer giants form around high-mass stars when sampling disc parameters compared to the case when disc sizes increase with increasing stellar mass. However, we still find that the super-Earth fraction is reduced around high-mass stars when viscous heating is included.

5. We tested our results for lower values of α_v and different fragmentation velocities and find that the fraction of super-Earths around high-mass stars is still reduced when viscous heating is active, similar to our nominal model.

This work shows how stellar masses and the structure of protoplanetary discs can shape the resulting planetary systems. In particular, we find that the reduction of close-in super-Earths around high-mass stars can be explained by an increase in giant planet formation efficiency around more massive stars. Future statistics of planetary systems around a wide range of stellar masses as well as better constraints on the structure of protoplanetary discs could help confirm these predictions.

Acknowledgements

We thank the anonymous referee for their comments that helped improve this manuscript. A.J. acknowledges funding from, the Carls-berg Foundation (Semper Ardens: Advance grant FIRSTATMO), and the Göran Gustafsson Foundation.

Appendix A Different photoevaporation models

In fig. A.1, we show the close-in super-Earth fraction using for both our nominal photoevaporation model from Owen et al. (2012) and the photoevaporation from Ercolano et al. (2023). The photoevaporation model from Ercolano et al. (2023) results in shorter disc lifetimes due to the overall higher photoevaporation rates. However, as the surface density of the disc is low at later stages of its lifetime, extending the lifetime beyond a few megayears, as our nominal model does, does not affect the planet population.

Fig. A.1 Close in super-Earth fraction as a function of stellar mass for different disc sizes. The circles show the results using our nominal photoevaporation model (same results as in fig. 4), while the triangles show the results using the photoevaporation model from Ercolano et al. (2023). We find no significant differences between the two models.

Appendix B Migration coefficient as a function time and disc size

Figs. B.1–B.3 show the migration coefficient as a function of distance to the star and planet mass but for different stellar masses. There are no significant differences for different stellar masses with the exception of r_tran extending slightly farther out for larger stellar masses. As larger discs evolve slower than smaller disc, the region where outward migration is prevalent moves inwards faster for smaller discs, which results in more efficient inward migration around smaller discs.

Fig. B.1 Same as Fig. 2 but at different times and for different disc sizes.

Fig. B.2 Same as Fig. B.1 but for a star with a mass of 1.5 M⊙.

Fig. B.3 Same as Fig. B.1 but for a star with a mass of 0.5 M⊙.

Appendix C Additional simulations

In Figs. C.1 and C.2, we show the close-in super-Earth fraction for a different value of α_v and different fragmentation velocities, respectively. We find minimal differences between these results and our nominal model.

Fig. C.1 Same as Fig. 4 but with αv = 10−3 and Ṁg,0 = 10−8 M⊙/yr. We find similar results in both the irradiated case (left) and the viscously heated case (right) compared to our nominal model. In the irradiated case, the fraction of super-Earths saturates at a lower stellar mass due to efficient inward migration. In the viscously heated case, we find that the fraction of close-in super-Earth still decreases with increasing stellar mass for large discs but the reduction is slightly lower compared to our nominal model.

Fig. C.2 Same as Fig. 4 but with pebble fragmentation velocities set to 1 m/s for dry pebble inside the water iceline and 10 m/s for icy pebbles outside the water iceline. The reduction of super-Earths for higher stellar masses is slightly stronger in the viscous case, as the dry region of the disc – where pebbles are small and planet formation is inefficient - extends farther away from the star for higher stellar masses. In the irradiated case, the results are relatively unchanged.

References

All Figures

Fig. 1 Disc lifetimes (top) and disc-to-star mass ratios (bottom) for discs with varying stellar masses and initial disc sizes. For a fixed disc mass, the disc lifetime increases with increasing stellar mass, while a larger disc also results in longer disc lifetime. Large discs (≳100 AU) around stars more massive than ~1.5 M⊙ are gravitationally unstable in some regions of the discs. This region is marked by a hatched pattern.

In the text

Fig. 2 Migration coefficient, ƒI,visc, as a function of distance to the star and planet mass for a solar-mass star with a disc size of 71 AU at t = 0. Inside of the transition radius (~5 AU), the disc is viscously heated, and for certain masses, the planet experiences outward migration. The sharp switch to inward migration at ~0.2 AU is caused by the significant decrease in opacity due to the evaporation of silicate and metal grains, which occurs at ~1000 K. We also show ƒI,visc at different times and for different stellar masses and disc sizes in Appendix B.

In the text

Fig. 3 Example of the resulting semi-major axes and masses of planets injected around a solar-mass star for two different initial disc sizes, 150 AU (top) and 50 AU (bottom). Here, only irradiation heating is considered. We indicate the planets that we define as close-in super-Earths (blue) and giants (red). The large disc forms significantly more giant planets due to having more available mass and a longer lifetime.

In the text

Fig. 4 The fraction of close-in super Earths (r < 0.2 AU) for a range of stellar masses (x-axis) and disc sizes (colours). In all cases, no super-Earths are formed around stars with a mass of 0.1 M⊙. When only considering discs heated by irradiation (left), the fraction of super-Earths increases with increasing stellar mass regardless of disc size. In the viscously heated case (right), the fraction of super-Earths reaches a maximum at 0.5 M⊙ and then declines slightly for large discs (90 AU or larger). For smaller discs (30 and 60 AU), the super-Earth fraction peaks at 1 M⊙, while for larger discs it peaks at 0.5 M⊙.

In the text

	Fig. 5 Growth tracks of a single planet injected at r = 3 AU and t = 104 yr. In the irradiated case (left), this results in the formation of a super-Earth at the inner disc edge for all stellar masses and disc sizes. In the viscous heating case (right), the protoplanet grows into a giant planet for high stellar masses as well as for solar-mass stars with sufficiently large discs.
In the text

Fig. 6 Giant planet fraction at all semi-major axes as a function of stellar mass and disc size. With or without viscous heating, the giant planet fraction increases with increasing stellar mass for all disc sizes where a giant planet can form to begin with. The lower stellar mass limit for giant planet formation is 0.5 M⊙ in both models. In the irradiated case (left), the lower disc size limit for giant planets to form is 60 AU (90 AU for M★ = 0.5 M⊙). In contrast, when viscous heating is included, the lower disc size limit is as low as 30 AU for solar-mass stars and above, while it is 60 AU for stars with M★ = 0.5 M⊙.

In the text

Fig. 7 Similar to Fig. 4 but using a fixed injection time at t = 105 yr. We also show the giant planet fraction (triangles). In both models, super-Earth formation is suppressed for high stellar masses, as planets have time to grow into giant planets instead. As a reminder, the giant planet fraction is defined as the number of giant planets formed compared to the total number of planets injected, while the super-Earth fraction is the number of close-in super-Earths compared to the number of close-in planets, as described in Section 3.1.

In the text

Fig. 8 Final planet mass as a function of injection location for different stellar masses and disc sizes, for a fixed injection time of t = 105 yr. For low-mass stars, neither disc model produces any super-Earths or giant planets. In the viscous model, the region where giant planets form moves close to the star for larger discs and stellar masses, resulting in fewer super-Earths forming. In the irradiated case, the giant-forming region remains far out in the disc and only moves slightly closer to the star for higher stellar masses.

In the text

Fig. 9 Fraction of close-in super-Earths as a function of stellar mass with sampled initial accretion rates and disc sizes. The black lines show the median values out of all runs, while the error bars show the 16th and 84th percentiles. Similarly to Fig. 4, the super-Earth fraction increases with increasing stellar mass in the irradiated case (left). When the disc is viscously heated (right), the super-Earth fraction flattens out for host stars more massive than 0.5 M⊙.

In the text

Fig. 10 Same as Fig. 9 but with a fixed injection time of 105 yr. Large discs form more giants when the embryos are injected early, which results in fewer super-Earths being formed. In contrast, if no giants are able to form, such as when the host star mass is below 0.5 M⊙, the super-Earth fraction increases with increasing disc size. We find that the median super-Earth fraction decreases with increasing stellar mass, both with or without viscous heating, although the effect is significantly stronger when viscous heating is active. Further, we find that the median super-Earth fraction increases slightly above 1 M⊙ in the viscous case as the sampled disc sizes around these high-mass stars are slightly smaller.

In the text

	Fig. A.1 Close in super-Earth fraction as a function of stellar mass for different disc sizes. The circles show the results using our nominal photoevaporation model (same results as in fig. 4), while the triangles show the results using the photoevaporation model from Ercolano et al. (2023). We find no significant differences between the two models.
In the text

	Fig. B.1 Same as Fig. 2 but at different times and for different disc sizes.
In the text

	Fig. B.2 Same as Fig. B.1 but for a star with a mass of 1.5 M⊙.
In the text

	Fig. B.3 Same as Fig. B.1 but for a star with a mass of 0.5 M⊙.
In the text

Fig. C.1 Same as Fig. 4 but with αv = 10−3 and Ṁg,0 = 10−8 M⊙/yr. We find similar results in both the irradiated case (left) and the viscously heated case (right) compared to our nominal model. In the irradiated case, the fraction of super-Earths saturates at a lower stellar mass due to efficient inward migration. In the viscously heated case, we find that the fraction of close-in super-Earth still decreases with increasing stellar mass for large discs but the reduction is slightly lower compared to our nominal model.

In the text

Fig. C.2 Same as Fig. 4 but with pebble fragmentation velocities set to 1 m/s for dry pebble inside the water iceline and 10 m/s for icy pebbles outside the water iceline. The reduction of super-Earths for higher stellar masses is slightly stronger in the viscous case, as the dry region of the disc – where pebbles are small and planet formation is inefficient - extends farther away from the star for higher stellar masses. In the irradiated case, the results are relatively unchanged.

In the text