© The Authors 2025

1 Introduction

Cold Jupiters are massive gaseous planets, usually categorized as those heavier than Neptune (≳30 M_⊕) that orbit their host star beyond the inferred ice lines (≳1 au). The occurrence rate of cold Jupiters is usually constrained with RV observations (Johnson et al. 2009; Montet et al. 2014; Rosenthal et al. 2021; Fulton et al. 2021; Hirsch et al. 2021; Ribas et al. 2023; Bonomo et al. 2023; Pass et al. 2023). The observational data suggest that cold Jupiters are more common around more massive stars with higher metallicity. While the occurrence rate of cold Jupiters is constrained around 20% around G-type stars, it decreases lower than 10% around M- and K-type stars (Fulton et al. 2021). The occurrence rate is rather low (less than 1.5%) around the low-mass (0.1-0.3 M_⊙) M dwarfs (Pass et al. 2023), but two confirmed planets, LHS 252 b (Morales et al. 2019) and GJ 83.1 b (Feng et al. 2020; Quirrenbach et al. 2022), and several candidates are detected there.

This result has been interpreted as a confirmation of the core accretion model (e.g. Mizuno 1980) since the solid material required for forming heavy-element cores is more abundant around those stars. In the core formation scenario, planetesimal and/or pebble accretion form a massive solid core (∼10 M_⊕). The solid core accretes surrounding disk gas and forms a gaseous envelope around the core. Once the protoplanet reaches the crossover mass, where the envelope mass equals the solid core mass, the protoplanet starts a rapid gas accretion (Pollack et al. 1996; Ikoma et al. 2000). The crossover mass is usually considered the starting point of runaway gas accretion. To form gas giant planets, protoplanets must reach the crossover mass before the disk dissipation. In the core accretion model, the occurrence rate of giant planets depends on the timescales of core formation and the subsequent gas accretion. Besides the core accretion model, the gravitational instability scenario (Boss 2006, 2019) is also suggested to form giant planets around low-mass M dwarfs (Morales et al. 2019). However, the relation between the occurrence rate and the stellar metallicity is not explained in this formation mechanism (Mercer & Stamatellos 2020). Therefore, the core accretion model would be the primary formation mechanism of cold Jupiters, while some of the giant planets around lower mass stars might be formed via gravitational instability.

Various formation models investigated how the occurrence rate of giant planets changes with the assumed stellar properties. The planetesimal accretion scenario is the classical formation model for a massive solid core (Inaba et al. 2003; Kobayashi & Tanaka 2021, 2023). Population synthesis models based on the planetesimal accretion scenario (Ida & Lin 2005; Alibert et al. 2011; Miguel et al. 2019; Burn et al. 2021; Guo et al. 2025) show that the formation rate of giant planets increases with stellar mass, which is consistent with the observations. However, the planetesimal accretion scenario has difficulty in the formation of giant planets around M dwarfs since the current planetesimal accretion model failed to form a massive solid core (Miguel et al. 2019; Schlecker et al. 2022).

In the pebble accretion scenario, protoplanets grow via accretion of small millimeter-to-centimeter-sized solids. Since the available pebble flux increases with the stellar metallicity and disk mass, Liu et al. (2019) found that giant planet formation is more efficient around higher mass stars with higher metal-licity. The solid core growth stops at the pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018), where pebbles are stalled from drifting inward at the pressure bump formed by the protoplanet. Since the pebble isolation mass is smaller around lower mass stars, the formation of solid cores that are massive enough to trigger runaway gas accretion is less efficient around lower mass stars (Liu et al. 2019, 2020). To form giant planets around lower mass stars, the pebble accretion scenario requires multiple planet formation and giant impacts between them to grow further than the pebble isolation mass (Pan et al. 2024, 2025).

The pebble accretion models successfully reproduced the observed trend in the occurrence rate of giant exoplanets. However, these studies mainly focus on the formation of solid cores, and many simplifications are used for the gas accretion process. Figure 1 shows our disk model at t = 0.1 Myrs. In the above pebble accretion models, protoplanets are assumed to start the Kelvin-Helmholtz contraction once the protoplanet reaches the pebble isolation mass. The Kelvin-Helmholtz contraction timescale depends on the boundary conditions of the planetary envelope (Piso & Youdin 2014) and the grain opacity (Ikoma et al. 2000), which strongly affects the giant planet formation (Bitsch & Savvidou 2021). Previous studies have used the fixed grain opacity and neglected the effects of boundary conditions. Therefore, it is still unclear how these parameters affect the formation of cold Jupiters, especially in the pebble accretion scenario.

Alibert et al. (2011) pointed out that the formation rate of giant planets is related to the disk’s lifetime since the protoplanets must undergo rapid gas accretion before disk dissipation. Despite the observations showing a large variety in the disk’s lifetime (Bayo et al. 2012; Ribas et al. 2015; Richert et al. 2018), how the occurrence rate of giant planets changes with the disk lifetime is not well discussed. While the detection of the inner dust ring in circumstellar disks shows that the disk’s lifetime rarely depends on the stellar mass (Richert et al. 2018), theoretical models for X-ray photoevaporation predict that the disk lifetime is longer around the lower mass stars because of the weaker X-ray irradiation (Picogna et al. 2021). Additionally, long-lived accretion disks are detected around M-type stars, known as Peter Pan disks (Silverberg et al. 2020). If the disk lifetime is longer around the lower mass stars, the giant planets may form in the core accretion scenario around M-dwarfs, even if the available solids are less abundant.

In this work, we model giant planet formation via pebble accretion while focusing on the crossover time, which is the time it takes for a protoplanet to reach the crossover mass. We investigate how the crossover time changes with the formation environment and how the occurrence rate changes with stellar type. In Sect. 2, we describe our formation model. In Sect. 3, we present the results of our numerical model and show how the crossover time changes with the formation environment. We discuss our results and model in Sect. 4. Our summary and conclusions are presented in Sect. 5.

Fig. 1 Disk model used in this work. Panels a–c show the surface density of the gaseous disk, disk temperature, and the disk’s aspect ratio around different stellar masses as shown in the color bar in panel a. The horizontal dashed line in panel b shows Tdisk = 170 K, where H2O condenses. Here, we show the case when t = 105 yr.

2 Methods

2.1 Disk model

Our disk model is based on the self-similar solution for the surface density profile of a gaseous disk (Lynden-Bell & Pringle 1974). We assume that the heating of the disk’s mid-plane gas is dominated by irradiation from the central star rather than viscous heating, where cold Jupiters form. Assuming a vertically optically thin and radially optically thick disk, the mid-plane temperature, T_disk, is given by (e.g., Chiang & Goldreich 1997; Ida et al. 2016) $$\begin{array}{r}{T}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}=150\mathrm{K}{\left(\frac{L_{\mathrm{s}}}{L_{\odot }}\right)}^{2/7}{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^{-1/7}{\left(\frac{r}{1\mathrm{a}\mathrm{u}}\right)}^{-3/7},\end{array}$$(1)

where L_s and M_s are the luminosity and mass of central star, and r is the radial distance from the central star. If the disk viscosity scales with the radial distance as ν ∝ r^γ, the surface density of disk gas, Σ_gas, is given as $$\begin{array}{r}{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{a}\mathrm{s}}=\frac{(2-\gamma )M_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k},0}}{2\pi {R_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}^2}{\left(\frac{r}{R_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}\right)}^{-\gamma }{T}^{-\frac{\frac{5}{2}-\gamma }{2-\gamma }}\exp \{-\frac{1}{T}{\left(\frac{r}{R_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}\right)}^{2-\gamma }\},\end{array}$$(2)

with $$\begin{array}{rl}\gamma & =\frac{3}{2}+\frac{\mathrm{d}\ln T_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}{\mathrm{d}\ln r}\end{array}$$(3) $$\begin{array}{rl}T& =1+\frac{t}{{\tau }_{\mathrm{v}\mathrm{i}\mathrm{s}}},\end{array}$$(4)

$$\begin{array}{rl}{\tau }_{\mathrm{v}\mathrm{i}\mathrm{s}}& =\frac{{R_{\mathrm{d}}}^2}{{\nu }_{\mathrm{d}}},\end{array}$$(5) where M_disk,0 is the disk total mass at t = 0, R_disk is a radial scaling length of the disk, τ_vis is the characteristic viscous timescale, and ν_d is a disk gas viscosity at r = R_d. We use the α-viscosity model (Shakura & Sunyaev 1973) and the disk gas viscosity is written as ν_d = α_acch²_gasΩ_K, where α_acc is the viscosity parameter, h_gas is the disk’s scale height, and Ω_K is the orbital angular velocity. The disk’s scale height is expressed as h_gas = c_s/Ω_K, where c_s is the sound speed of the disk gas. The sound speed is obtained using the ideal gas law, assuming a mean molecular weight of 2.34. In this work, we use the fixed α_acc = 10⁻³.

To investigate giant planet formation around different stellar masses, we introduce simple scaling laws following the numerical models available in the literature (Miguel et al. 2019; Burn et al. 2021; Chachan & Lee 2023; Venturini et al. 2024). Our scaling laws are given as $$L_{\mathrm{s}}=a{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^A,$$(6) $$M_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k},0}=b{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^B,$$(7) $$R_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}=c{\left(\frac{M_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k},0}}{0.1M_{\odot }}\right)}^C.$$(8)

The coefficients we used in this work are summarized in the Table 1. Note that the coefficients are not strongly constrained by observations. For example, the mass-luminosity relation changes with the age of the stellar mass. Here, we begin with pre-main sequence stars and follow the numerical results obtained by Ramirez & Kaltenegger (2014); Choi et al. (2016). The relation between the disk’s mass and stellar mass is taken from Andrews et al. (2013), and that between the disk’s mass and size is taken from Andrews et al. (2010).

2.2 Pebble accretion model

Throughout this paper, we assume the protoplanet’s eccentricity and inclination are negligibly small. We adopt the pebble accretion model by Johansen & Lambrechts (2017) where the pebble accretion rate is given by $${\dot{M}}_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}=\pi {R_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}}^2{\rho }_{\mathrm{p},\mathrm{m}\mathrm{i}\mathrm{d}}\overline{S}\delta v,$$(9)

where R_peb,acc is the pebble accretion radius, ρ_p,mid is the pebbles’ midplane density, S¯ is the stratification integral of pebbles, and δv is the pebble’s approaching speed which is given by ∆v + Ω_KR_p,acc where ∆v is the sub-Keplerian speed. R_peb,acc is obtained by solving the equation $${\tau }_{\mathrm{f}}=\frac{{\xi }_{\mathrm{B}}\mathrm{\Delta }v+{\xi }_{\mathrm{H}}{\mathrm{\Omega }}_{\mathrm{K}}{R}_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}}{\mathcal{G}{M}_{\mathrm{p}}/{R_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}}^2},$$(10)

where τ_f is the Stokes number, G is the gravitational constant, and ξ_B and ξ_H are fitting parameters. Numerical results are well reproduced with ξ_B = ξ_H = 0.25 (Ormel & Klahr 2010; Lambrechts & Johansen 2012). The stratification integral S¯ is given by $$\begin{array}{r}\overline{S}=\frac{1}{\pi R_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}^2}{\int }_{z_{\mathrm{p}}-R_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}}^{z_{\mathrm{p}}+R_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}}2\exp (-\frac{z^2}{2{h_{\mathrm{p}\mathrm{e}\mathrm{b}}}^2})\sqrt{{R_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}}^2-{(z-z_{\mathrm{p}})}^2}\mathrm{d}z,\end{array}$$(11)

where z_p is the height of the protoplanet measured from the disk’s mid-plane, and h_peb is the scale height of the pebble layer. We set z_p = 0 in this study.

The mid-plane density of pebbles is ${\rho }_{\mathrm{p},\mathrm{m}\mathrm{i}\mathrm{d}}={\mathrm{\Sigma }}_{\mathrm{p}\mathrm{e}\mathrm{b}}/\sqrt{2}{h}_{\mathrm{p}\mathrm{e}\mathrm{b}}$ where Σ_peb is the surface density of pebbles, which is given by $${\mathrm{\Sigma }}_{\mathrm{p}\mathrm{e}\mathrm{b}}=\frac{{\dot{M}}_{\mathrm{p}\mathrm{e}\mathrm{b}}}{2\pi r{v}_{\mathrm{p}\mathrm{e}\mathrm{b}}},$$(12)

with the radial velocity of pebbles $$v_{\mathrm{p}\mathrm{e}\mathrm{b}}=\frac{2\eta v_{\mathrm{K}}}{{\tau }_{\mathrm{f}}+{{\tau }_{\mathrm{f}}}^{-1}}+\frac{\nu }{r},$$(13)

and $$\eta =-\frac{1}{2}{\left(\frac{h_{\mathrm{g}\mathrm{a}\mathrm{s}}}{r}\right)}^2\frac{\mathrm{\partial }\ln P_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}{\mathrm{\partial }\ln r}.$$(14)

Here, Ṁ_peb is the pebble flux, v_K is the Kepler velocity at the protoplanet’s orbit, and P_disk is the gas pressure of the disk, which is calculated assuming an ideal gas. The pebble layer’s scale height h_peb is $$\begin{array}{r}{h}_{\mathrm{p}\mathrm{e}\mathrm{b}}=h_{\mathrm{g}\mathrm{a}\mathrm{s}}\sqrt{\frac{{\alpha }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}{{\tau }_{\mathrm{f}}}},\end{array}$$(15)

where α_turb is the local turbulent viscosity. Note that this equation is valid in the regime α_turb ≪ τ_f . The local turbulent viscosity α_turb could be different from the accretion viscosity α_acc, which drives the global angular momentum transfer (e.g. Bai et al. 2016). Following Ida et al. (2018), we assume α_acc ≥ α_turb and set α_turb as an input parameter. The pebble’s size is obtained by equating the pebble growth timescale with the drift timescale (Lambrechts & Johansen 2014). The surface density of pebbles transforms into $$\begin{array}{r}{\mathrm{\Sigma }}_{\mathrm{p}\mathrm{e}\mathrm{b}}=\sqrt{\frac{2{\dot{M}}_{\mathrm{p}\mathrm{e}\mathrm{b}}{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{a}\mathrm{s}}}{\sqrt{3}\pi {ϵ}_{\mathrm{p}}r{v}_{\mathrm{K}}}},\end{array}$$(16)

where ϵ_p is the pebble’s sticking efficiency which we set to 0.5. The pebble flux M_peb is given by the model developed by Lambrechts & Johansen (2014): $${\dot{M}}_{\mathrm{p}\mathrm{e}\mathrm{b}}=2\pi r_{\mathrm{g}}\frac{\mathrm{d}{r}_{\mathrm{g}}}{\mathrm{d}t}{Z}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{a}\mathrm{s}}(r=r_{\mathrm{g}},t=0)$$(17)

with $$\begin{array}{r}{r}_{\mathrm{g}}={\left(\frac{3}{16}\right)}^{1/3}{\left(\mathcal{G}{M}_{\mathrm{s}}\right)}^{1/3}{\left({ϵ}_{\mathrm{D}}{Z}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}\right)}^{2/3}{t}^{2/3},\end{array}$$(18) $$\begin{array}{r}\frac{\mathrm{d}{r}_{\mathrm{g}}}{\mathrm{d}t}=\frac{2}{3}{\left(\frac{3}{16}\right)}^{1/3}{\left(\mathcal{G}{M}_{\mathrm{s}}\right)}^{1/3}{\left({ϵ}_{\mathrm{D}}{Z}_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}\right)}^{2/3}{t}^{-1/3},\end{array}$$(19)

where Z_disk is the disk’s metallicity, which is set to 0.01 throughout this paper, and ε_D = 0.05. Note that in Eq. (17) we use the surface density of solids at t = 0 and neglect the radial movement of small dust with the viscously evolving gas. We adopt this model for the mass conservation of solids because our model does not include dust crossing of the pebble condensation front by advection and diffusion. Therefore, our model overestimates the pebble flux if we calculate it using the time-evolving Σ_gas. More details on that point can be found in Appendix A.

Pebble accretion terminates when the growing planet reaches the pebble isolation mass (Lambrechts et al. 2014; Bitsch et al. 2018), which is given by $$\begin{array}{r}{M}_{\mathrm{i}\mathrm{s}\mathrm{o}}=25.0M_{\oplus }\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right){\left(\frac{h/r}{0.05}\right)}^3\\ & & \times (0.34{\left[\frac{-3}{{\log }_{10}({\alpha }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}})}\right]}^4+0.66)(1-\frac{\frac{\mathrm{\partial }\ln P_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}{\mathrm{\partial }\ln r}+2.5}{6}).\end{array}$$(20)

We introduce an exponential cutoff at the pebble isolation mass and a cap for the pebble flux, which limits the heavy-element accretion rate, and the solid accretion rate is finally given by $$\begin{array}{r}{\dot{M}}_{\mathrm{Z}}=min({\dot{M}}_{\mathrm{p}\mathrm{e}\mathrm{b}},{\dot{M}}_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{a}\mathrm{c}\mathrm{c}}\exp (-{\left[\frac{M_{\mathrm{p}}}{M_{\mathrm{i}\mathrm{s}\mathrm{o}}}\right]}^{10})),\end{array}$$(21)

where M_p is the protoplanet’s mass.

2.3 Gas accretion

To calculate the gas accretion rate onto the growing protoplanet, we use the MESA-extention code developed by Valletta & Helled (2020), which is based on the stellar evolution code MESA (Paxton et al. 2011, 2018). MESA’s version is r24.03.1. The planet is assumed to be in hydrostatic equilibrium and spherically symmetric. The equations of state are adapted from Müller et al. (2020). We assumed that the outer edge of the planetary envelope is connected to the local disk gas. The surface temperature T_surf and pressure P_surf were set to T_disk and P_disk. Before entering the runaway gas accretion (M_p ≲ 10 M_⊕), the protoplanet’s radius is comparable to or smaller than the disk’s scale height (see Eq. (22) below). In this regime, temperature and pressure gradients across the surface of the envelope are relatively small, justifying the assumption of spherical symmetry within the scope of this study.

The gas accretion rate Ṁ_env is calculated by the method developed in Pollack et al. (1996). In each timestep, we add a mass ∆m to fill the gap between the planetary radius, R_p, and the gas accretion radius, R_gas,acc, which is given by $$R_{\mathrm{g}\mathrm{a}\mathrm{s},\mathrm{a}\mathrm{c}\mathrm{c}}=\frac{\mathcal{G}{M}_{\mathrm{p}}}{{c_{\mathrm{s}}}^2+4\mathcal{G}{M}_{\mathrm{p}}/R_{\mathrm{H}}},$$(22)

where R_H is the Hill radius. By adding ∆m, R_p increases (or decreases if ∆m is negative), and the gap becomes narrower.

We iterate each step until δ = |(R_p - R_gas,acc)/R_gas,acc| becomes smaller than 10⁻³.

We used a dust grain opacity prescription of Valencia et al. (2013) where the grain opacity was scaled with a grain opacity factor, f_g. f_g = 1 corresponds to the interstellar medium’s grain opacity, while in planetary atmospheres, lower values are expected (Movshovitz et al. 2010; Mordasini 2014; Ormel 2014). In this study, we set f_g as a free parameter.

It is known that the heavy-element deposition in the planetary envelope accelerates the gas accretion process (Hori & Ikoma 2010, 2011; Venturini et al. 2015, 2016; Valletta & Helled 2020; Mol Lous et al. 2024). For simplicity, in this study, we assumed that the accreted heavy elements settle to the core and release the accretion energy there. Instead, we focus on the disk conditions and how the stellar mass affects the formation of giant planets. We discuss the importance of heavy-element deposition on giant planet formation in Sect. 4.1.

Convection was modeled using the mixing-length theory. The mixing length parameter α_MLT was set to 2, which is estimated from the stellar evolution models (Paxton et al. 2011). The value α_MLT in planetary condition is poorly constrained, and α_MLT could be significantely lower (Leconte & Chabrier 2012). However, gas accretion is rather independent of α_MLT if the envelope is not enriched with heavy elements (see Appendix B). Therefore, in our baseline simulations, we use α_MLT = 2.

2.4 Formation simulation setup

We start the simulation assuming a heavy-element core with a mass M₀ = 0.01 M_⊕ at t = t₀. The location of the protoplanet r was fixed at the same orbital period, which corresponds to 3 au around a solar mass star; namely, r = 3 au (M_s/M_⊙)^1/3. We neglected the planetary migration in our simulations. To form cold Jupiters, protoplanets need to enter the outward migration region (Coleman & Nelson 2016), or start their formation in the orbit farther than 10 au (Bitsch et al. 2015). Otherwise, rapid type-I migration allows the protoplanets to migrate to the disk’s inner edge and become warm/hot-Jupiters. The migration speed and the location of the outward migration region depend on various parameters such as the disk’s structure (Tanaka et al. 2002; Coleman & Nelson 2016), temperature profile, and the corotation torque (Paardekooper et al. 2010, 2011; Paardekooper 2014). We fixed the protoplanet’s orbit to reduce the number of free parameters and focus on the gas accretion process. The effects of planetary migration are discussed in Sect. 4.4. The heavyelement core increases its mass by accreting pebbles and disk gas. The simulation was terminated when the protoplanet reaches the crossover mass, when the planetary envelope mass, M_env, equals the core mass, M_core, or when the protoplanet’s growth timescale τ_grow = Ṁ_p/M_p exceeds 10⁹ yr.

We perform parameter studies when changing the stellar mass M_s, core formation time t₀, grain opacity factor f_g, and the disk’s turbulent viscosity α_turb. We summarise the parameter ranges used in this paper in Table 2.

3 Results

3.1 The role of grain opacity

Figure 2 shows the planetary formation when assuming different grain opacity factors f_g. Before reaching the pebble isolation mass (horizontal dashed line), the protoplanet’s mass increases mainly via pebble accretion (core growth phase). Since the rapid pebble accretion brings energy large enough to support the hydrostatic condition, the protoplanet’s envelope does not contract. The gas accretion rate is much lower than the solid accretion rate and is controlled by the expansion of the accretion radius R_gas,acc. Once the protoplanet reaches pebble isolation mass, the pebble accretion stops, and the envelope begins to contract. The cooling timescale controls the gas accretion rate (contraction phase). Since the cooling timescale depends on the atmospheric opacity, the gas accretion rate is shorter for the lower grain opacity factors.

Figure 3 shows the crossover time t_g,cross, defined as the time when the protoplanet reaches the crossover mass M_cross (solid lines). For comparison, we also plot the core growth time t_g,core, defined as the time required for the protoplanet to grow from the initial mass M₀ to the pebble isolation mass M_iso (dashed lines). Additionally, the envelope’s contraction time t_g,cont, representing the time between reaching M_iso and M_cross, is shown with dotted lines. Note that in Fig. 3 – as well as Figs. 5 and 6 – we plot the crossover time measured from the core formation time t_g,cross - t₀ to facilitate direct comparison with the core growth time t_g,core and the contraction time t_g,cont. Here, we show the cases with M_s = 1 M_⊕, t₀ = 10⁶ yr, and α_turb = 10⁻³. The crossover time increases with f_g because a larger f_g delays the contraction of the planetary envelope. Unlike in the contraction phase, the core growth phase is independent of f_g. If f_g is greater than 0.1, we find that the contraction phase becomes longer than the core growth phase; therefore, the crossover time depends on f_g. If f_g is smaller than 0.1, the protoplanet reaches the crossover mass shortly after the pebble isolation mass is reached. In this case, the formation of giant planets is regulated only by pebble accretion, and the crossover time rarely depends on f_g.

Fig. 2 Growth of protoplanets up to crossover mass around 1 M⊙ star with the core formation time t0 = 1 Myrs. Different colors show the cases with different grain opacity fg. The solid, dashed, and dotted lines show the total, envelope, and solid core mass, respectively. The horizontal gray dashed line shows the pebble isolation mass.

Fig. 3 Growth times obtained in the simulations with Ms = 1 M⊙ , t0 = 106 yr, and αturb = 10−3. The solid line shows the crossover time, measured from the core formation time. The dashed and dotted lines show the times required for the protoplanet to grow from M0 to Miso and from Miso to Mcross, respectively.

Fig. 4 Same as Fig. 2, but for different assumed turbulent viscosities αturb. The upper and lower panels show the cases with fg = 10−2 and fg = 100, respectively.

3.2 The role of turbulent viscosity

Next, we explore the importance of the assumed disk’s turbulent viscosity. The disk’s turbulent viscosity affects the pebble accretion rate and isolation mass. The pebble accretion rate is higher for lower α_turb because of the smaller pebble layer’s scale height (see Eq. (15)). On the other hand, the pebble isolation mass becomes smaller for lower α_turb as shown in Eq. (20). Figure 4 shows the time evolution of protoplanets with different α_turb. We also plot the crossover time presented in Fig. 5. Here, we show the cases with M_s = 1 M_⊕, t₀ = 10⁶ yr, and f_g = 0.01 (upper panels) and 1 (lower panels). Due to the higher pebble accretion rate, the protoplanet reaches the pebble isolation mass earlier for lower α_turb. The planetary core mass is smaller for lower α_turb because of the smaller pebble isolation mass. As shown in Fig. 5, the core growth time t_g,core is shorter for lower α_turb. On the other hand, the contraction time t_g,cont is longer for lower α_turb. This is because of the smaller pebble isolation mass. The envelope’s contraction timescale strongly depends on the mass of protoplanets (Ikoma et al. 2000; Piso & Youdin 2014).

When the grain opacity is small (for example, f_g = 0.01), the core growth time regulates the crossover time as discussed in Sect. 3.1, and the crossover time decreases for lower α_turb. On the other hand, when f_g = 1, the dominant growth phase is the envelope’s contraction phase, and the crossover time is longer for lower α_turb.

Fig. 5 Same as Fig. 3, but for different assumed turbulent viscosities αturb. The upper and lower panels show the cases with fg = 10−2 and fg = 100, respectively.

Fig. 6 Same as Fig. 3, but for different core formation times t0. The top and bottom panels show the cases with fg = 0.01 and 0.1, respectively. The different colors correspond to different stellar masses. The solid, dashed, and dotted lines show the crossover time measured from the core formation time tg,cross - t0, the core growth time tg,core, and the contraction time tgcont, respectively. Here, we show the cases with αturb = 10−3.

3.3 Core formation time and stellar mass

Figure 6 shows the crossover time as a function of the core formation time. We show the cases with low grain opacity f_g = 0.01 and high grain opacity f_g = 1 in the top and bottom panels, respectively. First, we focus on the cases around 1 M_⊙ stars (blue lines). We find that the protoplanets take longer to reach the pebble isolation mass (dashed lines) if the core formation time is longer. This is because both the total dust mass remaining in the protoplanetary disk and the pebble flux decrease with time. Unlike in the core growth phase, the growth time in the contraction phase t_g,cont (dotted lines) depends only weakly on t₀.

We also plot the cases around different stellar masses in Fig. 6. Because of the lower pebble flux, the protoplanet takes longer to reach the pebble isolation mass around lower-mass stars. When f_g = 0.01, the protoplanet reaches crossover mass soon after the pebble isolation mass is reached, and t_g,cont is rather independent of the mass of the host stars. However, t_g,cont depends on the stellar type when f_g = 1. The pebble isolation mass is larger for planets around more massive stars. Due to the larger core mass, the planetary envelope contracts quickly, and therefore, t_g,cont is shorter around more massive stars.

We find that protoplanets do not reach pebble isolation around very low mass stars M_s = 0.1 M_s. However, the protoplanets could reach the crossover mass if the core formation time is as short as 10⁵ yr and f_g = 0.01. In the case of t₀ = 10⁵ yr, the protoplanet’s core mass is 1.2 M_⊕. Despite the small core mass, the protoplanet can accrete disk gas due to the low atmospheric opacity. With higher grain opacity f_g = 1, protoplanets do not reach the crossover mass with M_s = 0.1 M_s.

While the crossover time t_g,cross ranges from 10⁵ to 10⁷ yr for lower grain opacity cases f_g = 0.01, t_g,cross exceeds 5 Myrs for the higher grain opacity cases f_g = 1.0. Since observations suggest that typical disk lifetimes are ∼1-3 Myrs (Richert et al. 2018), the frequency of giant planet formation is expected to be low if f_g = 1.0. In the next section, we further discuss the occurrence rate of giant planets.

3.4 Crossover time

Figure 7 shows the inferred crossover time t_g,cross as a function of stellar mass M_s and core formation time t₀ , for different f_g and α_turb. The crossover time shortens with faster core formation around more massive stars. When f_g = 1.0 and α_turb = 10⁻³, giant planet formation is rather slow, and crossover time exceeds 3 Myrs. In that case, only the long-lived protoplane-tary disks would enable the formation of giant planets. When the grain opacity factor is as small as f_g = 0.01, the crossover time becomes ∼1 Myr. Giant planets would form if the initial core formation ends within 1 Myrs, except around late M-dwarfs (M_s = 0.1 M_⊕).

As shown in Sect. 3.2, a smaller turbulent viscosity of α_turb = 10⁻⁵ shortens the crossover time if the grain opacity is small as f_g = 0.01, but lengthens if f_g = 1. A smaller turbulent viscosity α_turb = 10⁻⁵ and grain opacity f_g = 0.01 are required to form giant planets around late M-dwarfs. In this case, the crossover time weakly depends on the stellar type if the initial core forms within 1 Myr. Even if the initial core formation takes t₀ ∼ 3 Myrs (log₁₀t_g,cross = 6.5), giant planets could still form around high-mass stars with long-lived protoplanetary disks. On the other hand, no giant planets would form with a low turbulent viscosity α_turb = 10⁻⁵ and a high grain opacity f_g = 1.

Fig. 7 Crossover time obtained from simulations starting with an initial core of Mcore = 0.01 M⊕. We plot the obtained crossover time tg,cross as a color contour, and the color bar is shown on the right side of the panels. The horizontal and vertical axes are stellar mass Ms and the core formation time t0. The black lines show the contour lines and the number on each line shows log10 tg,cross. The left and right columns correspond to the low (fg = 0.01) and high (fg = 1) grain opacity cases. The top and bottom rows show the high (αturb = 10−3) and low (αturb = 10−5) turbulent viscosity cases.

3.5 From core formation to crossover

In the above sections, we investigate the formation of giant planets, setting the core formation time t₀ as an input parameter. Hereafter, we estimate t₀ by calculating the pebble accretion onto a protoplanet, which is set to be the largest planetesimal formed by streaming instability. Then, we estimate the crossover timescale using the numerical results obtained in Sect. 3.4.

3.5.1 Core formation time

Hydrodynamic simulations show that the maximum size of planetesimals formed by the streaming instability could be several hundred km in radius (Simon et al. 2017; Schäfer et al. 2017). We estimate the maximum size of planetesimals using two different models. The first model is the model used in Lau et al. (2022), where the characteristic planetesimal size is set by the typical clump size of pebbles determined by the local diffusivity δ, which is given by (Klahr & Schreiber 2020) $$\begin{array}{rl}{m}_{\mathrm{G}}& =\frac{1}{9}{\left(\frac{\delta }{{\tau }_{\mathrm{f}}}\right)}^{3/2}{\left(\frac{h_{\mathrm{g}\mathrm{a}\mathrm{s}}}{r}\right)}^3{M}_{\odot },\end{array}$$(23) $$\begin{array}{rl}& \sim 7.22\times {10}^{-3}{\left(\frac{\delta /{\tau }_{\mathrm{f}}}{{10}^{-4}/{10}^{-2}}\right)}^{3/2}{\left(\frac{h_{\mathrm{g}\mathrm{a}\mathrm{s}}/r}{0.058}\right)}^3{M}_{\oplus }.\end{array}$$(24)

The streaming instability simulations by Abod et al. (2019) show that the cumulative mass distribution of planetesimals exhibits an exponential cutoff ~0.3m_G. Although such massive planetesimals are relatively rare, planetesimals with masses of ~ m_G form via streaming instability. We adopt m_G as the maximum planetesimal size m_plt. Following Lau et al. (2022), we set δ as a free parameter independent of α_turb. We use log₁₀ δ = −4.0, −4.5, and −5.0. The second model is the model by Liu et al. (2020). In this model, the mass of the dust clump is determined by the self-gravity and the tidal shear, which is given by $$\begin{array}{r}{m}_{\mathrm{p}\mathrm{l}\mathrm{t}}=2\times {10}^{-3}{M}_{\oplus }\left(\frac{f_{\mathrm{p}\mathrm{l}\mathrm{t}}}{400}\right){\left(\frac{\gamma }{{\pi }^{-1}}\right)}^{1.5}{\left(\frac{h_{\mathrm{g}\mathrm{a}\mathrm{s}}/r}{0.05}\right)}^3\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)\end{array}$$(25)

with $$\begin{array}{r}\gamma =\frac{4\pi \mathcal{G}{\rho }_{\mathrm{g}}}{{\mathrm{\Omega }}_{\mathrm{K}}^2},\end{array}$$(26)

where f_plt is the ratio between the maximum mass and the characteristic mass of planetesimals formed by the streaming instability. Liu et al. (2020) estimates f_plt as 400 from the streaming instability simulations by Schäfer et al. (2017), but f_plt could change with the disk conditions (Simon et al. 2017). To account for this uncertainty, we also use a value of f_plt = 40, which is smaller by a factor of ten.

Figure 8 shows the inferred initial mass of the protoplanet. We assume that the first planetesimal forms at t = 10⁵ yr with M_p = m_plt. We calculate the planetary growth via pebble accretion without the gas accretion and estimate the core formation time t₀ by integrating it up to 10⁻²M_⊕. Then, we estimate the crossover time using the numerical results obtained in Sect. 3. We introduce a fitting function t_g,cross,fit(α_turb, f_g) to the obtained crossover time t_g,cross as a function of the core formation time t₀ and stellar mass M_s. We use scipy.interpolate.RBFInterpolator to determine the fitting function.

Fig. 8 Maximum mass of planetesimals mplt used in our model. The solid lines show the estimated mplt in the model by Lau et al. (2022). The different colors correspond to different δ. The dashed lines show the estimated mplt in the model by Liu et al. (2020). The different color corresponds to different fplt. The right y axis shows the radius of the protoplanet. We use the mean density of 3g/cm3 to estimate the radius.

3.5.2 Effects of initial protoplanet’s mass

Figure 9 shows the estimated crossover time. We plot the cases using different initial protoplanet masses m_plt . First, the protoplanet could not reach the crossover mass when we use Eq. (25) with δ = 10⁻⁵. In this case, the mass of the initial protoplanet is ≲10⁻⁵ M_⊕ and too small to initiate rapid pebble accretion. We present the mass-doubling timescale τ_peb in Fig. 10. As shown in this figure, the pebble accretion rate strongly depends on the mass of the growing protoplanet when ≲10⁻⁵ M_⊕. To initiate rapid pebble accretion, the initial protoplanet’s mass must be larger than ≳10⁻⁵ M_⊕.

The protoplanet can grow to the crossover mass with the other initial protoplanet mass. For f_g = 1, we find that the crossover time rarely depends on the initial protoplanet mass. This is because the growth time in the contraction phase is longer than in the core growth phase.

When the grain opacity factor is small as f_g = 0.01, the crossover time is shorter for the larger initial protoplanet mass. However, even if the initial mass was 10 times larger, the crossover time is shortened by 2–3. In the pebble accretion model, the mass-doubling timescale rarely changes with the protoplanet mass when M_p ≳ 10⁻⁵ M_⊕ as shown in Fig. 10. If the mass doubling timescale τ_peb is constant, the protoplanet mass is written as $$\begin{array}{rl}\frac{M_{\mathrm{p}}}{{\dot{M}}_{\mathrm{p}}}& ={\tau }_{\mathrm{p}\mathrm{e}\mathrm{b}},\\ \log \frac{M_{\mathrm{p}}}{m_{\mathrm{p}\mathrm{l}\mathrm{t}}}& =\frac{t}{{\tau }_{\mathrm{p}\mathrm{e}\mathrm{b}}}.\end{array}$$(27)

Since the contraction phase is negligibly shorter than the core growth phase, the time reaching the crossover is given by $$\begin{array}{r}{t}_{\mathrm{g},\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}}={\tau }_{\mathrm{p}\mathrm{e}\mathrm{b}}(\ln M_{\mathrm{i}\mathrm{s}\mathrm{o}}-\ln m_{\mathrm{p}\mathrm{l}\mathrm{t}}).\end{array}$$(28)

Therefore, the crossover time weakly depends on the initial protoplanet mass if the initial planetesimals are born big as ≳10⁻⁵ M_⊕.

From the above results, we conclude that the protoplanet’s initial size has a minor effect on the crossover time in the pebble accretion scenario, compared to the other parameters, like the grain opacity and stellar mass.

Fig. 9 The estimated crossover time in the pebble accretion scenario. The solid and dashed lines show the cases with the planetesimal formation model by Eqs. (24) and (25), respectively. The different colors show the cases with different parameters, as described in the legend box. Each panel shows the cases with different grain opacities. The gray area shows the typical disk lifetime inferred by the observations (Richert et al. 2018), which ranges from 1 to 3 Myrs. We plot the horizontal dotted gray lines at 5 and 10 Myrs for the eye guide.

Fig. 10 Mass doubling timescale in the pebble accretion model. The solid and dashed lines show the different times t = 105 and 106 yr, respectively. The color corresponds to the different masses of central stars, and the color bar is written on the panel’s upper right.

3.5.3 Crossover time across different stellar types

We find that the crossover time decreases with increasing stellar mass and changes the dependence on the stellar mass around M_s ~ 0.4 M_⊕ as shown in Fig. 9. For planets orbiting around stars of M_s ≲ 0.4 M_⊕, the crossover time strongly depends on M_s. These planets could not reach pebble isolation mass before the pebbles in the protoplanetary disk are depleted. The pebble flux is generated by the sweep of the pebble growth line r_g. When r_g < R_disk, the pebble flux slowly decreases with time (∝ t^−8/21). Once the pebble growth line reaches the disk’s typical radius R_disk, the pebble flux decreases exponentially because of the reduction in the dust surface density at the disk’s outer edge. By substituting r_g = 2 R_disk into Eq. (18), r = 2 R_disk gives a ten times reduction in Σ_gas because of the disk’s outer edge, we obtained the pebble depletion time as $$\begin{array}{r}{\tau }_{\mathrm{p}\mathrm{e}\mathrm{b},\mathrm{d}\mathrm{e}\mathrm{p}}=2.08\mathrm{M}\mathrm{y}\mathrm{r}{\left(\frac{M_{\mathrm{s}}}{1M_{\odot }}\right)}^{\frac{3}{2}BC-\frac{1}{2}}{\left(\frac{{ϵ}_{\mathrm{D}}}{0.05}\right)}^{-1}{\left(\frac{Z_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}}{0.01}\right)}^{-1}.\end{array}$$(29)

We plot Eq. (29) on the left panel in Fig. 9.

When f_g = 0.01, the envelope contraction phase is much shorter than the core growth phase, as we found in Sect. 3. t_g,cross < τ_peb,dep means that the protoplanet can reach the pebble isolation mass before the pebbles are depleted. The dependence of the crossover time on the stellar mass comes from the core growth time, which decreases with increasing stellar mass because the mass doubling timescale τ_peb is shorter around the higher mass stars. On the other hand, the protoplanets could not reach the pebble isolation mass in the cases with t_g,cross > τ_peb,dep. For these non-isolated protoplanets, the envelope’s contraction takes a longer time because of the smaller core mass. The protoplanet’s core mass is smaller around the smaller mass stars because of the shorter τ_peb,dep and longer τ_peb. Therefore, the crossover time steeply increases with decreasing stellar mass.

We plot the typical disk lifetime 1-3 Myrs with the gray area in Fig. 9. Using the dust observations in stellar clusters, Richert et al. (2018) found that the disk longevity rarely depends on the stellar mass, and the fraction of stars with disks is 60–70% at 1 Myrs and 30-40% at 3 Myrs. Currently, there is no evidence that the lifetimes of the dust disk and the gas disk are the same. However, it is reasonable to assume that both components could be dispersed by the same mechanisms, such as photoevaporation or disk winds. However, dust growth could lead to the selective removal of small dust grains while leaving the gas behind, and possibly leading to a gaseous disk that survives longer than dust disks inferred from observations. For simplicity, in this study the lifetimes of the dust and gaseous disks are assumed to be equal. Under this assumptions our results suggest that more than 30% of stars heavier than ~0.4 M_⊙ have cold Jupiters if f_g ≲ 0.1. This result is inconsistent with the observed occurrence rate of cold Jupiters (Fulton et al. 2021), where only ≲20% of G-type stars have cold Jupiters.

To be consistent with the low occurrence rate of cold Jupiters, a higher grain opacity of f_g > 0.1 is required. When f_g = 1, the crossover time ranges between 5 and 10 Myrs in M_s ≳ 0.4 M_⊙. Also, the crossover time decreases with the stellar mass. Due to the high grain opacity, the crossover time is regulated by the envelope’s contraction. The timescale in Kelvin-Helmholtz contraction is given by (Ikoma et al. 2000) $${\tau }_{\mathrm{K}\mathrm{H}}={10}^6\mathrm{y}\mathrm{r}{\left(\frac{M_{\mathrm{p}}}{10M_{\oplus }}\right)}^{-k_2}\left(\frac{{\kappa }_{\mathrm{e}\mathrm{n}\mathrm{v}}}{1\mathrm{c}{\mathrm{m}}^2/\mathrm{g}}\right),$$(30)

where κ_env is the effective opacity in the planetary atmosphere and k₂ is estimated 2.5-3.5. In our model, the pebble isolation mass scales as $$M_{\mathrm{i}\mathrm{s}\mathrm{o}}\sim 8M_{\oplus }{\left(\frac{M_{\mathrm{s}}}{M_{\odot }}\right)}^{\frac{3}{7}}.$$(31)

Using k₂ = 3.0 and substituting M_p = M_iso, we get τ_KH ∝ M_s^−9/7. The strong dependency between the crossover time and the stellar mass could explain the higher occurrence rate of cold Jupiter around more massive stars, as inferred from observations (Fulton et al. 2021).

The exact occurrence rate also depends on the distribution of the disk’s lifetime. However, the available observational data are limited to t ≤ 5 Myrs in Richert et al. (2018), and the fraction of stars with disks between 5 and 10 Myrs is unclear yet. Recent studies suggest that the characteristic timescale for the disk’s lifetime could be 5-10 Myr (Michel et al. 2021; Pfalzner et al. 2022), which prefers the scenario with the high grain opacity f_g = 1 to be consistent with the observed cold Jupiter’s occurrence rate. Since the formation of giant planets is controlled by the gaseous disk’s lifetime, a better understanding of the distribution of disk lifetimes is required.

Also, the gas distribution in the protoplanetary disks was assumed to be smooth. However, this is a clear simplification. Indeed, recent observations of protoplanetary disks show the common existence of substructures such as rings and spirals (e.g. Andrews 2020; Benisty et al. 2023). These substructures would affect the pebble flux and its accretion onto protoplanets. For example, Lau et al. (2022) shows that if planetesimals form and grow in a pressure bump, planetesimal accretion and rapid pebble accretion can occur, leading to the formation of a core of several M_⊕ within 10⁵ years. If such substructures shorten the core growth phase, it emphasizes the need for a slow contraction phase. It is therefore clear that a better understanding of the mechanisms by which substructures in protoplanetary disks accelerate or delay the formation of giant planets is required. This will allow us to better connect the properties of the planetary disk and the planetary formation history.

Fig. 11 The growth of the planetary core up to the crossover mass accounting for envelope enrichment caused by pebble evaporation. The blue lines show the case with Ms = 1 M⊙, t0 = 106 yr, αturb = 10−3 and fg = 1. The gray lines show the case without envelope enrichment. Panel a: time evolution of total mass (solid line), envelope mass (dashed line), and core mass (dotted line). Panel b: time evolution of the bulk metallicity of planetary envelope. Panel c: evolution of planetary masses as a function of the total mass. Panel b: evolution of the bulk metallic-ity of planetary envelope vs. the total mass of protoplanet. The vertical dashed lines in Panels c and d are pebble isolation mass.

4 Discussion and caveats

4.1 Heavy-element deposition in the envelope

In the simulations presented above, we neglect the ablation of the accreted pebbles and assume that all the heavy elements settle to the center. However, the enrichment of the envelope with heavy elements accelerates the contraction of the planetary envelope (Hori & Ikoma 2010, 2011; Venturini et al. 2015, 2016; Valletta & Helled 2020; Mol Lous et al. 2024). In this section, we investigate how envelope enrichment changes the crossover time.

We use the ablation model developed by Valletta & Helled (2020) and Mol Lous et al. (2024). We use the direct deposition model, where the ablated heavy elements from the accreted pebbles are deposited into the local envelope shell. We set the pebble composition to be 50% rock and 50% H₂O, and assume that only the ablated H₂O is deposited in the envelope. We include the condensation of H₂O using the saturation pressure and the maximum water enhancement of the supercritical water Z_sup,H2O = 0.9. We also introduce the pre-mixing algorithm for heavy elements if the Rayleigh-Taylor criterion is met. Further details of our heavy-element enrichment method are presented in Appendix D. To save computational time, we start with M₀ = 0.1 M_⊕, instead of M₀ = 0.01 M_⊕.

Figure 11 shows the planetary evolution including the envelope enrichment by the accreted pebbles. We find that the metallicity in the protoplanet envelope increases quickly and saturates to ~0.8. After reaching the pebble isolation mass, the atmospheric metallicity decreases because only gas accretes to the protoplanet in this phase. We also plot the cases without envelope enrichment with gray lines. Note that these cases differ from the simulations in Sect. 3 because of different M₀. Interestingly, we find that the crossover time is longer for the cases with envelope enrichment. The high metallicity of the envelope leads to higher gas accretion rates and larger envelope mass during the core growth phase. Since the total protoplanet mass regulates the pebble isolation mass, the core mass decreases at the pebble isolation compared to the cases without envelope enrichment. Piso & Youdin (2014) shows that the gas accretion rate in the following contraction phase strongly depends on the core mass. Due to the smaller core, the growth time in the contraction phase is longer when envelope enrichment is considered. This situation does not occur in the study of Mol Lous et al. (2024) because the pebble isolation mass was set to ∼25 M_⊕ and the protoplanet did not reach the pebble isolation in 10 Myrs.

We also performed simulations with a lower grain opacity factor f_g = 0.01 (Fig. 12). Due to the rapid envelope contraction in the contraction phase, the effect of envelope enrichment is weaker. The inferred difference in the crossover time is only <0.3 Myrs.

We find that heavy-element enrichment in the planetary envelopes does not necessarily shorten the formation timescale of giant planets. However, we investigated only a few parameter studies using the specific enrichment parameters (50%-rock-50%-water and Z_sup,H2O = 0.9). As shown in Mol Lous et al. (2024), the effect of envelope enrichment is more profound if the accreted pebbles are composed of pure-water because the energy deposited by the accreted pebbles is small. Even with 50%-rock-50%-water, the ablation and deposition of rock, which we assume non-missible to the H-He envelope, would accelerate gas accretion and would change the mass of the planetary core. Future studies should investigate in further detail how the planetary formation timescale changes with envelope enrichment when considering various compositions.

Fig. 12 Same as Fig. 11, but with fg = 0.01.

4.2 Extending the envelope contraction phase

We find that a high grain opacity factor of f_g > 0.1 is required to explain the low occurrence rate of cold Jupiters. However, dust growth models in planetary envelopes suggest that small dust grains grow to larger sizes and quickly settle to deeper regions, reducing grain opacity (Movshovitz et al. 2010; Mordasini 2014; Ormel 2014). The estimated grain opacity factor is f_g ≤ 0.01 (Mordasini 2014). Additional mechanisms that delay the envelope’s contraction, such as gas recycling or a continuous accretion of dust and planetesimals after the protoplanet reached pebble isolation mass, are required.

3D hydrodynamical simulations show that the accreted gas enters the Hill sphere through the poles and exits to the disk’s midplane (Tanigawa et al. 2012; Ormel et al. 2015; Kurokawa & Tanigawa 2018). The recycled gas keeps bringing small dust to the radiative-advective boundary region of the planets, which would lead to a dust opacity comparable to the (local) dust opacity in the disk (Lambrechts & Lega 2017). After pebble isolation mass is reached, the drifting pebbles pile up at the outer pressure bump. This can trigger fragmentation, and dust grains that are small enough to be coupled to the disk gas can flow into the planetary envelope (Stammler et al. 2023). The accretion of small dust could increase the grain opacity in the planetary envelope and may delay the envelope’s contraction (Chen et al. 2020). Further investigation of the grain opacity in protoplan-etary disks, including the grain growth around a protoplanet, should be conducted in future studies.

Another possible mechanism is planetesimal accretion. In this study, we neglected planetesimal accretion after the formation of initial planetesimals. No protoplanets grow via pebble accretion if the initial protoplanet mass m_plt is smaller than 10⁻⁵ M_⊕. However, planetesimal accretion could grow the small planetesimals to a point where they trigger rapid pebble accretion (e.g., Johansen & Lambrechts 2017). The timing of when rapid pebble accretion starts could be later than the time of planetesimal formation, and the crossover time could be more delayed than that obtained in our simulations. Also, we neglected the process of planetesimal accretion, which could take place in parallel to pebble accretion. Unlike pebble accretion, planetesimal accretion could continue even after the gap formation, and accreted planetesimals could provide additional energy to the planetary envelope, delaying the contraction of the planetary envelope (Alibert et al. 2018; Helled 2023). Since planetesimals need to form to start rapid pebble accretion, planetesimal accretion seems to be inevitable during giant planet formation, although the accretion rate could be smaller than that of pebble accretion. In fact, it is possible that pebble accretion and planetesimal accretion dominate at different stages during planet formation as well as in different disk environments. Understanding the interplay between planetesimal accretion and pebble accretion is desirable.

4.3 Effects of disk scaling models

The pebble accretion process strongly depends on the assumed disk model. In our disk model, protoplanets could not reach pebble isolation around stars of ≲0.4 M_⊕, and the break point of the crossover time where the dependency on stellar mass changes appears around 0.4 M₀. However, as shown in Eq. (29), the pebble depletion time decreases around low-mass stars, and the break point would shift to larger mass if the exponent BC gets larger.

In our model, the disk accretion rate M_d,acc scales as $$\begin{array}{r}\frac{\mathrm{d}\log {\dot{M}}_{\mathrm{d},\mathrm{a}\mathrm{c}\mathrm{c}}}{\mathrm{d}\log M_{\mathrm{s}}}\sim \frac{2}{7}A+B\{1-C(2-\gamma )\}=\frac{39}{112}.\end{array}$$(32)

Here, we neglect the viscous evolution and radial decay terms in Σ_gas. The relation between the disk accretion rate and the stellar mass can be estimated from observations (e.g. Manara et al. 2022), and our scaling law seems to be shallower than that in observations. Single power-law fitting suggests a slope of 1.6-2.. Such a steep slope requires higher A, B, or smaller C than those in our model. However, the exact relation between Ṁ_d,acc and M_⊙ is still actively discussed. Several studies show that a double power-law fit for Ṁ_d,acc - M₀ is a more appropriate representation than a simple power-law fit (Alcalá et al. 2017; Manara et al. 2017). In this case, the slope is flatter for stars more massive than 0.2-0.3 M₀. Note that the observed disk accretion rate is time-dependent and that the exact scaling law at time zero remains unknown. A more robust determination of the properties of disks around different stellar types is important for constraining the origin of cold Jupiters and predicting their occurrence rate.

Fig. 13 A summary of our results and the inferred giant planet formation regimes.

4.4 Planetary migration

For simplicity, planetary migration is not included in this study. However, it is clear that planet migration plays an important role in planet formation. The migration regime typically changes from type I to type II around the pebble isolation mass (∼10 M_⊕), and the migration speed peaks around this transition (Kanagawa et al. 2018; Ida et al. 2018). When the grain opacity is as high as f_g = 1, the protoplanet enters a long contraction phase after reaching pebble isolation mass. During this phase, the protoplanet can undergo significant radial migration, potentially becoming a warm or hot Jupiter. In contrast, when the grain opacity is low, for example, f_g = 0.01, the protoplanet grows rapidly after pebble isolation mass is reached (see Fig. 3) and quickly shifts to the slower type II migration. Figure 13 is the summary of our findings. This quick transition to slower migration could help prevent significant inward migration. Note, however, that it remains unclear whether such low grain opacities are common.

While inward migration is expected to occur, outward migration (Paardekooper et al. 2010, 2011; Paardekooper 2014), and migration trap (Masset et al. 2006; Coleman & Nelson 2016) could also take place. Outward migration occurs near the transition radius, where the dominant disk heating mechanism shifts from viscous accretion to stellar irradiation (Liu et al. 2019, 2020). However, this mechanism weakens as the planet grows in mass and the disk evolves. Alternative mechanisms are disk substructures that can act as migration traps, enabling the survival of cold Jupiters (Coleman & Nelson 2016). Both outward migration and migration traps are sensitive to the evolving disk structure. A more comprehensive understanding of the formation of “cold” giant planets will require future investigations that incorporate planetary migration self-consistently, ideally through disk evolution models and hydrodynamical simulations.

If significant planetary migration occurs, a quick transition to slower type II migration might be required. To be consistent with the lower occurrence rate of cold Jupiters, planetary cores should form less efficiently. Our model assumes that icy pebbles are sticky and continue to grow until the gas drag from the disk causes the growing pebbles to drift inward. However, recent observations of dust in protoplanetary disks (Jiang et al. 2024) and laboratory experiments (Musiolik et al. 2016) suggest that icy pebbles are more fragile and pebbles initiate fragmentation even outside the water ice line. If this is the case, the pebble accretion efficiency is reduced due to the lower Stokes number. Since the planetary core growth takes longer, giant planets would form less efficiently, even if the grain opacity is lower than f_g = 1. Fragile icy pebbles may support the formation of cold Jupiters.

It remains unclear how planet migration takes place in disks around different stellar types and we are still lacking a good understanding of how the disk properties (for example, density, lifetime) change with the mass of the central star. As a result, it is still unclear whether different migration histories are expected. Future studies should investigate how planetary migration affects planet formation across different stellar masses.

5 Summary and conclusions

We investigated the formation of cold Jupiters around different mass stars. We explored how the crossover time changes with stellar mass M_s, core formation time t₀, grain opacity factor f_g, and disk turbulent viscosity α_turb, assuming pebble accretion. Our key conclusions can be summarized as follows:

• The grain opacity strongly affects the planetary growth. If the grain opacity is lower than f_g ≲ 0.1, the crossover time is comparable to the growth timescale of the heavy-element core. Contrary, if f_g ≳ 0.1, the envelope contraction dominates the formation process, and the protoplanets experience a long envelope contraction phase after reaching the pebble isolation mass.

• Low disk turbulent viscosity α_turb accelerates (delays) reaching crossover mass when the grain opacity is lower (higher) than f_g ~ 0.1. When f_g < 0.1, the core growth time decreases for lower α_turb because of the smaller pebble layer’s scale height and higher pebble accretion rate. On the other hand, the envelope contraction timescale increases because of the lower pebble isolation mass if f_g > 0.1.

• The crossover time depends on the initial size of protoplanets. However, this dependency is weaker than linear. The grain opacity plays a bigger role than the initial size of protoplanets in determining the crossover time.

• Around stars of M_s ≲ 0.4 M₀, the crossover time rapidly increases with decreasing stellar mass. Around these stars, pebbles in the protoplanetary disk deplete before the protoplanets reach pebble isolation mass. Because of the small mass of the forming solid core, the envelope’s contraction timescale (before reaching the crossover mass) is long. A very long disk lifetime ≳10⁷ yr, as found in Peter Pan disks, is required to form giant planets around such stars.

• Around stars with masses M_s ≳ 0.4 M₀, protoplanets reach crossover mass within 3 Myrs if the grain opacity is f_g ≲ 0.1, which is comparable or shorter than the typical disk’s lifetime. To be consistent with the low occurrence rate of cold Jupiters around G-type stars (≲20%), a high grain opacity (f_g > 0.1) is required. With a high grain opacity f_g = 1, the crossover timescale increases around lower-mass stars, which could explain the inferred low-occurrence rate of cold Jupiters around small stars.

We find that the crossover time can be divided into two regimes by comparing the core formation time t_g,core and the pebble depletion timescale τ_peb,dep. In the fast pebble depletion regime where τ_peb,dep < t_g,core, the protoplanet does not reach the pebble isolation mass, and the crossover time exceeds 10 Myrs because of the small core mass. Giant planets may form in a long-lived disk, such as the Peter-Pan disk, but this formation path would not be the primary formation path of cold Jupiters around FGK-type of stars.

If the planetary core forms before the pebbles are depleted (slow pebble depletion regime where τ_peb,dep > t_g,core), the protoplanet can grow to the pebble isolation mass. To account for the observed low occurrence rate of cold Jupiters around G-type stars, either slow envelope contraction or slow core formation – lasting longer than the typical disk lifetim – is required. In our pebble accretion model, protoplanets quickly reach the pebble isolation mass and therefore, a slow contraction of the planetary envelope is required to suppress giant planet formation and explain the low occurrence rate of cold Jupiters. This result is consistent with recent formation models that suggest that the formation of Jupiter and Saturn could have taken a few million years (Alibert et al. 2018; Helled 2023). However, rapid contraction may be required to form giant planets around lower-mass stars. This result may indicate that the metallicity and dust content of protoplanetary disks increase with the mass of the central stars. Indeed, it was recently shown that the bulk metallicity of giant planets is lower for giant planets orbiting small stars in comparison to giant planets around sun-like stars (Müller & Helled 2025). The relationship between the heavy-element mass in the disk and stellar mass is yet to be determined, and is likely to control giant planet formation around different stellar masses.

Interestingly, we found that envelope enrichment with heavy elements released from accreted pebbles could delay envelope contraction, a different result from that obtained in previous studies. Because of the higher mean molecular weight, the protoplanets could have a more massive envelope at the pebble isolation mass. Consequently, the core mass at the pebble isolation is smaller when envelope enrichment is considered. Due to the smaller core mass, the envelope contracts more slowly, resulting in a delayed crossover time. Our results also do not rule out the possibility of slow core formation. The timescale for core growth depends sensitively on disk structure, metallicity, and the dust fragmentation velocity.

We note that the durations of core formation and envelope contraction depend on the formation history and, therefore, on various processes and parameters. The core formation timescale strongly depends on the local conditions in the disk and the types of solids that are accreted. Similarly, while in this study, a slow contraction phase is achieved with a high grain opacity, delaying the envelope contraction could be the result of other mechanisms, such as atmospheric recycling or planetesimal accretion. It is therefore clear that a better understanding of the conditions and processes taking place during planet formation is crucial. Overall, our findings highlight the complex interplay between core accretion, envelope physics, and disk evolution in shaping planetary outcomes. We suggest that future studies should include these processes self-consistently and investigate a wide parameter space. Only then can we reveal the origin of cold Jupiters in our galaxy.

Acknowledgements

The numerical computations were carried out on PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. R.H. acknowledges support from SNSF under grant 200020_215634.

Appendix A Solid mass conservation

Fig. A.1 Time evolution of the cumulative pebble flux. The solid lines show the cumulative pebble flux calculated with Eq. 17. The dashed lines are those obtained with Eq. 17 using the time evolving Σgas(t = t′) instead of Σgas(t = 0). The horizontal dotted lines show the initial total mass of solids in the protoplanetary disk. The different colors correspond to disks around different stellar masses.

In Eq. 17 we use the solid surface density at t = 0 instead of that at t = t′, where t′ is the disk’s age. Using Σ_gas(t = t′) in Eq. 17 would account for the radial transport of small dust grains due to gas advection and diffusion, and thus seems more realistic than using Σ_gas(t = 0). However, this approach overestimates the pebble flux. This is because the gas that diffuses outward from r < R_disk to r > R_disk may be depleted of solids. For simplicity, we do not consider this effect in our model.

To illustrate this point, we plot the cumulative pebble flux calculated with Eq. 17 using Σ_gas(t = 0) (solid lines) and Σ_gas(t = t′) (dashed lines). We find that the two models yield similar evolution up to t′ ~ τ_vis, but the model using Σ_gas(t = t′) eventually overestimates the pebble flux and violates mass conservation by exceeding the initial total solid mass. To preserve mass conservation, we adopt Σ_gas(t = 0) in Eq. 17.

Our pebble flux model is valid only if the pebble condensation front reaches the outer edge of the disk before the gaseous disk starts viscous evolution, when τ_peb,dep ≲ τ_vis. For disks with higher viscosity or larger radial extent, a more sophisticated model of dust evolution is required to compute the pebble flux accurately.

Appendix B Mixing length parameter

To check the effect of different mixing length parameter α_MLT, we compare the simulations with α_MLT = 2 (solid lines) and α_MLT = 0.1 (dashed lines) in Fig. B.1. We found that the difference in the obtained crossover time is less than 5%. We conclude that the effect of α_MLT is negligibly smaller than that of the other parameters, such as stellar mass M_s and the grain opacity factor f_g.

Fig. B.1 The evolution of protoplanets around different mass stars using different αMLT = 2 (solid lines) and αMLT = 0.1 (dashed lines). The top and bottom panels show the cases with different grain opacity factors fg = 0.01 and 1, respectively. The color corresponds to the stellar mass, as shown in the color bar.

Appendix C Supply-limited gas accretion

In our simulations, the gas accretion rate is obtained by iterating the planetary radius and the gas accretion radius. The energy balance in the planetary envelope determines the gas accretion rate in this phase. However, the gas accretion rate onto the protoplanet could be limited by the supply of local disk gas to the planetary orbital region. Once the envelope’s contraction and the following gas accretion become fast enough to deplete local disk gas, the planetary envelope detaches from the protoplane-tary disk, and the accretion regime shifts to supply-limited gas accretion. Here, we compare our numerical results to the maximum supply of disk gas and discuss how supply-limited gas accretion affects our results.

To estimate the maximum supply of disk gas, we follow the approach developed by Tanaka et al. (2020). The flux of disk gas entering the protoplanet’s Hill sphere is given by (Tanigawa & Watanabe 2002) $$\begin{array}{r}{\dot{M}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},\mathrm{a}\mathrm{c}\mathrm{c}}=D{\mathrm{\Sigma }}_{\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}},\end{array}$$(C.1)

with: $$\begin{array}{r}D=0.29{\left(\frac{M_{\mathrm{p}}}{M_s}\right)}^{4/3}{\left(\frac{h_{\mathrm{g}\mathrm{a}\mathrm{s}}}{r}\right)}^{-2}{r}^2{\mathrm{\Omega }}_{\mathrm{K}}.\end{array}$$(C.2)

Σ_band is the gas surface density at the gas accretion band. Protoplanets massive enough to open a gap in the gaseous disk reduce the local gas surface density and Σ_band. Using the gap structure model obtained by Kanagawa et al. (2017), we get $$\begin{array}{r}\frac{{\mathrm{\Sigma }}_{\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}}{{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{a}\mathrm{s}}}=\frac{1}{1+0.04K},\end{array}$$(C.3)

with $$\begin{array}{r}K={\left(\frac{M_{\mathrm{p}}}{M_s}\right)}^2{\left(\frac{h_{\mathrm{g}\mathrm{a}\mathrm{s}}}{r}\right)}^{-5}{{\alpha }_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}^{-1}.\end{array}$$(C.4)

Equation C.1 shows the limit of gas disk supply through the local gas flow. In addition, the global disk accretion M_d,acc limits the supply of disk gas onto protoplanets. Finally, the maximum supply of disk gas is $$\begin{array}{r}{\dot{M}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{s}\mathrm{u}\mathrm{p}}=min({\dot{M}}_{\mathrm{d},\mathrm{a}\mathrm{c}\mathrm{c}},{\dot{M}}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},\mathrm{a}\mathrm{c}\mathrm{c}}).\end{array}$$(C.5)

Fig. C.1 Gas accretion rate obtained in our numerical simulations of Ms = 1M⊙ and log10t0/yr = 6. The blue and green lines show the cases with fg = 1 and 0.01, respectively. The solid and dashed lines correspond to the cases with αturb = 10−3 and 10−5, respectively. We also plot the maximum gas accretion rate in the supply-limited case with black (t = 106 yr) and gray (t = 107 yr) lines. Also, the solid and dashed lines correspond to cases with αturb = 10−3 and αturb = 10−5, respectively. The thin gray dotted lines

Fig. C.1 shows the gas accretion rate obtained in our simulations and the maximum supply of disk gas M_max,sup. In most of the simulation phase, the gas accretion rate is slower than M_max,sup. In simulations with f_g = 0.01, the gas accretion rate reaches M_max,sup just before the crossover mass is reached. Around the crossover mass, the gas accretion timescale is much shorter than the disk’s lifetime. Therefore, even when we include the effects of supply-limited gas accretion, the crossover time obtained in our simulations remains similar.

Appendix D Heavy-element deposition in the gaseous atmosphere

We use the ablation model presented by Valletta & Helled (2020) and Mol Lous et al. (2024). We use the direct deposition model, where the ablated heavy elements (from the pebbles) are deposited into the local envelope shell. We set the composition of pebbles as 50 % rock and 50 % H₂O, and assume that only the ablated H₂O is deposited in the envelope.

Valletta & Helled (2020) and Mol Lous et al. (2024) introduced the condensation of H₂O by putting the maximum metal-licity of a local layer Z_max,i (index-i increases toward the center of the protoplanet). The first limit is the saturation of H₂O Z_sat,H2O. Using the saturation pressure P_sat,H2O, Z_sat,H2O is defined as P_sat,H2O/P_i where P_i is the local total pressure. The second limit is the maximum water enhancement of the supercritical water Z_H2O, which is set to 0.9, as used in Mol Lous et al. (2024).

In addition, Mol Lous et al. (2024) introduced a smoothing factor f_smooth to limit the jump in the local metallicity Z_i. Pebbles are usually evaporated in the upper atmosphere and sink to the lower envelope. A large fraction of the heavy elements are deposited in the layer where Z_sat,H2O becomes ≳ 0.1. This sink of heavy elements leads to a rather large jump in Z_i. As a result, it is not easy for MESA to find a hydrostatic solution. To avoid a large jump in Z_i, Mol Lous et al. (2024) introduced a smoothing factor Z_smooth,i = Z_pre,i + f_smooth, where Z_pre,i is the local metallicity in the previous step.

The maximum metallicity of the local layer i is given by $$\begin{array}{r}{Z}_{\mathrm{m}\mathrm{a}\mathrm{x},i}=min(Z_{\mathrm{s}\mathrm{a}\mathrm{t},{\mathrm{H}}_2\mathrm{O}},Z_{\mathrm{s}\mathrm{u}\mathrm{p},{\mathrm{H}}_2\mathrm{O}},Z_{\mathrm{s}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h},i}).\end{array}$$(D.1)

If the deposition of the accreted heavy elements leads to a local metallicity Z_i that is higher than Z_max,i, the overabundant heavy elements settle to the lower layer. The smoothing factor coefficient f_smooth is set to 0.2. The smoothing factor makes the composition gradient slightly shallower than that without it. This difference in the composition gradient has a weak effect on the crossover time. We found that with the smoothing factor, the crossover times could be 5% shorter at most than those obtained without it. This is because the solid core mass at the pebble isolation becomes slightly larger due to the shallower compositional gradient and the smaller envelope metallicity. This difference does not alter our main finding that the crossover time could be longer if the heavy element deposition is included in the planetary envelope.

Pebbles evaporated in the upper atmosphere create an increasing composition gradient toward the outer envelope. MESA tries to erase this inverted composition gradient by triggering the convection under the Ledoux criterion. However, it makes the timestep too small to proceed, and numerical convergence can be difficult. To avoid this, we introduce a pre-mixing algorithm that assumes Rayleigh-Taylor instability. We cannot use the actual criterion for Rayleigh-Taylor, which compares the density of adjacent layers, because our enrichment model increases the metallicity of local shells without changing the mass there (see Valletta & Helled 2020). Instead, we use this simplified criterion that checks the metallicity of adjacent layers. We mix heavy elements in and distribute them equally to the adjacent two shells if Z_i > Z_{i + 1}.

References

All Tables

All Figures

	Fig. 1 Disk model used in this work. Panels a–c show the surface density of the gaseous disk, disk temperature, and the disk’s aspect ratio around different stellar masses as shown in the color bar in panel a. The horizontal dashed line in panel b shows Tdisk = 170 K, where H2O condenses. Here, we show the case when t = 105 yr.
In the text

	Fig. 2 Growth of protoplanets up to crossover mass around 1 M⊙ star with the core formation time t0 = 1 Myrs. Different colors show the cases with different grain opacity fg. The solid, dashed, and dotted lines show the total, envelope, and solid core mass, respectively. The horizontal gray dashed line shows the pebble isolation mass.
In the text

	Fig. 3 Growth times obtained in the simulations with Ms = 1 M⊙ , t0 = 106 yr, and αturb = 10−3. The solid line shows the crossover time, measured from the core formation time. The dashed and dotted lines show the times required for the protoplanet to grow from M0 to Miso and from Miso to Mcross, respectively.
In the text

	Fig. 4 Same as Fig. 2, but for different assumed turbulent viscosities αturb. The upper and lower panels show the cases with fg = 10−2 and fg = 100, respectively.
In the text

	Fig. 5 Same as Fig. 3, but for different assumed turbulent viscosities αturb. The upper and lower panels show the cases with fg = 10−2 and fg = 100, respectively.
In the text

Fig. 6 Same as Fig. 3, but for different core formation times t0. The top and bottom panels show the cases with fg = 0.01 and 0.1, respectively. The different colors correspond to different stellar masses. The solid, dashed, and dotted lines show the crossover time measured from the core formation time tg,cross - t0, the core growth time tg,core, and the contraction time tgcont, respectively. Here, we show the cases with αturb = 10−3.

In the text

Fig. 7 Crossover time obtained from simulations starting with an initial core of Mcore = 0.01 M⊕. We plot the obtained crossover time tg,cross as a color contour, and the color bar is shown on the right side of the panels. The horizontal and vertical axes are stellar mass Ms and the core formation time t0. The black lines show the contour lines and the number on each line shows log10 tg,cross. The left and right columns correspond to the low (fg = 0.01) and high (fg = 1) grain opacity cases. The top and bottom rows show the high (αturb = 10−3) and low (αturb = 10−5) turbulent viscosity cases.

In the text

Fig. 8 Maximum mass of planetesimals mplt used in our model. The solid lines show the estimated mplt in the model by Lau et al. (2022). The different colors correspond to different δ. The dashed lines show the estimated mplt in the model by Liu et al. (2020). The different color corresponds to different fplt. The right y axis shows the radius of the protoplanet. We use the mean density of 3g/cm3 to estimate the radius.

In the text

Fig. 9 The estimated crossover time in the pebble accretion scenario. The solid and dashed lines show the cases with the planetesimal formation model by Eqs. (24) and (25), respectively. The different colors show the cases with different parameters, as described in the legend box. Each panel shows the cases with different grain opacities. The gray area shows the typical disk lifetime inferred by the observations (Richert et al. 2018), which ranges from 1 to 3 Myrs. We plot the horizontal dotted gray lines at 5 and 10 Myrs for the eye guide.

In the text

	Fig. 10 Mass doubling timescale in the pebble accretion model. The solid and dashed lines show the different times t = 105 and 106 yr, respectively. The color corresponds to the different masses of central stars, and the color bar is written on the panel’s upper right.
In the text

Fig. 11 The growth of the planetary core up to the crossover mass accounting for envelope enrichment caused by pebble evaporation. The blue lines show the case with Ms = 1 M⊙, t0 = 106 yr, αturb = 10−3 and fg = 1. The gray lines show the case without envelope enrichment. Panel a: time evolution of total mass (solid line), envelope mass (dashed line), and core mass (dotted line). Panel b: time evolution of the bulk metallicity of planetary envelope. Panel c: evolution of planetary masses as a function of the total mass. Panel b: evolution of the bulk metallic-ity of planetary envelope vs. the total mass of protoplanet. The vertical dashed lines in Panels c and d are pebble isolation mass.

In the text

	Fig. 12 Same as Fig. 11, but with fg = 0.01.
In the text

	Fig. 13 A summary of our results and the inferred giant planet formation regimes.
In the text

Fig. A.1 Time evolution of the cumulative pebble flux. The solid lines show the cumulative pebble flux calculated with Eq. 17. The dashed lines are those obtained with Eq. 17 using the time evolving Σgas(t = t′) instead of Σgas(t = 0). The horizontal dotted lines show the initial total mass of solids in the protoplanetary disk. The different colors correspond to disks around different stellar masses.

In the text

	Fig. B.1 The evolution of protoplanets around different mass stars using different αMLT = 2 (solid lines) and αMLT = 0.1 (dashed lines). The top and bottom panels show the cases with different grain opacity factors fg = 0.01 and 1, respectively. The color corresponds to the stellar mass, as shown in the color bar.
In the text

Fig. C.1 Gas accretion rate obtained in our numerical simulations of Ms = 1M⊙ and log10t0/yr = 6. The blue and green lines show the cases with fg = 1 and 0.01, respectively. The solid and dashed lines correspond to the cases with αturb = 10−3 and 10−5, respectively. We also plot the maximum gas accretion rate in the supply-limited case with black (t = 106 yr) and gray (t = 107 yr) lines. Also, the solid and dashed lines correspond to cases with αturb = 10−3 and αturb = 10−5, respectively. The thin gray dotted lines

In the text