© ESO 2020

1 Introduction

The formation of the solid core of gas giant planets is still an open question in regards to the core accretion model (Pollack et al. 1996). The timescale on which a core of critical mass must be formed in order for runaway gas accretion to set in is given by the disk lifetime. Depending on the stellar properties, this period has been found observationally to be between 3 and 10 Myr for Solar-type stars (Ribas et al. 2015). On this short timescale, interstellar medium μm dust grains must grow by almost 14 orders of magnitude in size and 40 orders of magnitude in mass. The critical core mass value at which the hydrostatic equilibrium breaks down is in the 5−20 M_⊕ range (Mizuno 1980; Béthune 2019), depending on how the envelope is able to cool (Rafikov 2006; Piso et al. 2015), and it is strongly dependent on the adopted opacity (Mordasini 2014; Ormel 2014; Brouwers & Ormel 2020). Furthermore, due to the interaction between the planetary core and its natal circumstellar disk, planets migrate and can be accreted onto the central star if no other physical processes stop them.

In order to match the observed population of gas giants, a faster growth timescale is needed in the core accretion paradigm. In recent years, two promising mechanisms have been studied. Under certain circumstances, the streaming instability might allow for the growth of planetesimal-sized objects directly from small dust particles (see e.g., Johansen et al. 2014, for a review), while pebble accretion might enhance the solid accretion on growing planetary cores (Ormel & Klahr 2010; Lambrechts & Johansen 2012).

Although these mechanisms have successfully reduced the time required to form giant planetary cores, a further barrier must be overcome in order to reach the core mass necessary to ignite the rapid gas accretion phase. A growing planet is able to open a partial gap in the surrounding gaseous disk by generating a pressure maximum, which effectively halts the inward drift of solids in the disk, in particular the pebble component. The mass at which the pebble flux is halted has been defined as the pebble isolation mass (Lambrechts et al. 2014), and it varies based on the local physical properties of the disk (Bitsch et al. 2018; Ataiee et al. 2018; Picogna et al. 2018).

Within this framework, we study the pebble accretion rate and isolation mass, in realistic 3D radiative hydrodynamical simulations, and how this affects the efficiency of gas giant planet formation. In Sect. 2 we describe the numerical set-up adopted to model the solid evolution and accretion onto a planetary core embedded in a disk with a realistic equation of state. We then discuss our results in Sect. 3 and draw the main conclusions in Sect. 4.

2 Numerical methods

We consider the same setup as in Picogna et al. (2018) to allow for a direct comparison; the only exception is that instead of using a locally isothermal equation of state for the gas, we adopted the ideal equation of state and included the radiative transfer in the flux-limited diffusion approximation (FLD; Kolb et al. 2013), and we consider laminar viscous disks.

2.1 Gas component

The disk was initially set in an axisymmetric configuration ranging from 2.08 to 13 au (0.4–2.5 in code units, where the unit of length is 5.2 au). In the vertical direction, the domain extends from θ_min = 83° to θ_max = 90° (half disk), corresponding to 5 disk scale heights H in the initial setup, and in the ϕ direction from ϕ_min = 0° to ϕ_max = 360° (full disk). In the three coordinate system’s directions (r, θ, ϕ), we used a 600 × 64 × 1024 grid resolution. The grid cells are spaced logarithmically in the radial, and linearly in the other directions. The main parameters of the simulations are summarized in Table 1.

The initial density profile created by force equilibrium is given by $$\rho (R,Z)={\rho }_0{\left(\frac{R}{R_{\text{p}}}\right)}^p\exp \left[\frac{G{M}_{\text{s}}}{c_{\text{s}}^2}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$$(1)

where R = r sin θ describes the cylindrical radius, ρ₀ is the gas midplane density at the planet location R = R_p = 1, p = −1.5 the density exponent, and c_s is the isothermal sound speed. The disk’s initial aspect ratio was set to h = H∕R = 0.05, which corresponds to a temperature profile of $$T(R)=T_0{\left(\frac{R}{R_{\text{p}}}\right)}^q,$$(2)

with q = −1 and T₀ = 121 K.

The gas moves with an azimuthal velocity given by the Keplerian speed around a 1 M_⊙ star, corrected by the pressure support (Nelson et al. 2013): $$\Omega (R,Z)={\Omega }_{\text{K}}{\left[(p+q)h^2+(1+q)-\frac{qR}{\sqrt{R^2+Z^2}}\right]}^{\frac{1}{2}},$$(3)

where Ω_K is the Keplerian orbital frequency. At the inner and outer boundary, we adopted reflective conditions, and we damped the density as well as the radial and vertical velocity to the initial values at the timescale of a fraction of a local orbit in order to prevent reflection of the spiral wave caused by the planet-disk interaction onto the boundary. The damping was applied in the intervals [0.4, 0.5] r_p and [2.3, 2.5] r_p. For the vertical boundaries, a mirror condition was implemented at the disk midplane, while an open boundary was applied in the disk atmosphere, and in the azimuthal direction a periodic condition was applied. We adopted a constant of α = 5 × 10⁻⁴, which is consistent with values generated by hydrodynamical turbulences (see e.g., Stoll et al. 2017).

We employed the radiation hydrodynamics module developed by Kolb et al. (2013) for the PLUTO code (Mignone et al. 2007). The solver is based on the FLD approximation in the two-temperature approach. The equations were solved in the comoving frame in the frequency-independent (gray) approximation.

The motion of the gas is described by the Navier-Stokes equations (Eqs. (4)–(6)) that are coupled with radiation transport (Eq. (7)) $$\begin{array}{l}\frac{\partial }{\partial t}\rho +\nabla \cdot (\rho \upsilon )=0,\\ \frac{\partial }{\partial t}\rho \upsilon +\nabla \cdot (\rho \upsilon \otimes \upsilon -\sigma )+\nabla p=-\rho \nabla \Phi ,\\ \frac{\partial }{\partial t}e+\nabla \cdot [(e+P)\upsilon ]=-\rho \upsilon \cdot \nabla \Phi +(\sigma \cdot \nabla )\upsilon -{\kappa }_{\text{p}}\rho c(a_R{T}^4-E),\\ \frac{\partial }{\partial t}E+\nabla \cdot F={\kappa }_{\text{p}}\rho c(a_{\text{R}}{T}^4-E),\end{array}$$

where the first three equations describe the evolution of the gas motion where ρ is the gas density, P corresponds to the thermal pressure, υ is the velocity, e = ρɛ + 1∕2ρυ² the total energy density (i.e., the sum of internal and kinetic energy) of the gas without radiation, ɛ = C_VT is the specific internal energy in which C_V is the specific heat capacity (assumed constant here), σ is the viscous stress tensor, and Φ is the gravitational potential. This system of equations is closed by the ideal gas equation of state $$P=(\gamma -1)\rho \varepsilon =\rho \frac{k_{\text{B}}T}{\mu {\text{m}}_{\text{H}}},$$(8)

where γ is the ratio of specific heats, T represents the gas temperature, k_rmB is the Boltzmann constant, μ = 2.35 the mean molecular weight, and m_H is the hydrogen mass.

The evolution of the radiation energy density E is given by Eq. (7), where F denotes the radiative flux, which was computed in the FLD approximation as in Kolb et al. (2013), κ_P is the Planck mean opacity, c corresponds to the speed of light, and a_R is the radiative constant. This implementation does not include the advective transport terms for the radiation energy and radiative pressure work in Eqs. (6) and (7), since these terms are of minor importance for low-temperature disks. For the computational boundaries, we adopted, for the radiative part, reflective conditions in the radial direction, a symmetric condition at the disk midplane, and a fixed temperature of 5 K in the upper boundary. Period boundaries were applied in the azimuthal direction.

2.2 Dust component

The solid fraction of the disk is modeled with 10⁶ Lagrangian particles divided into ten size bins as reported in Table 1. This approach has the great advantage of modeling a broad range of dynamical behaviors self-consistently, using the same modeled particles. The trade-off is that in the regions of low density, the resolution of the dust population is lower. However, for our study, this is not a problem since we are mainly interested in the dynamical evolution of dust particles; thus we do not take into account collisions between particles or the back-reaction of the dust onto the gas. We study particles with sizes s from 0.1 mm up to 1 km and internal density ρ_d = 1 g cm⁻³. The particlesizes were chosen to cover a wide range of different dynamical behaviors. The initial surface density profile of the dust particles is $${\Sigma }_{\text{d}}(r)\propto R^{-1}.$$(9)

This particle distribution leads to an equal number of particles in each radial ring as the grid is spaced logarithmically in the radial direction. The resulting dust profile is steeper than the gas profile, and it allows for a better sampling of the dust dynamics in the vicinity of the planet location. The dust particles were initially placed with a vertical distribution given by the local disk scale height and the dust diffusion coefficient (see e.g., Youdin & Lithwick 2007). The evolution of dust particles is given by the gravitational interaction with the planet and central star, turbulent kicks from the gas that resemble a realistic turbulent behavior, and the drag force from the interaction with the gaseous disk $$F_{\text{drag}}=-\frac{m_{\text{d}}}{t_{\text{s}}}{v}_{\text{r}},$$(10)

where m_d is the dust particle mass, and t_s is the stopping time which represents the timescale on which the embedded dust particle approaches the gas velocity and, for well-coupled particles, it is given by $$t_{\text{s}}=\frac{s{\rho }_{\text{d}}}{\rho {\overline{v}}_{\text{th}}},$$(11)

or, in its dimensionless form (hereafter, Stokes number), as $${\tau }_{\text{s}}=t_{\text{s}}{\Omega }_{\text{K}}(r),$$(12)

which describes the effect of a drag force acting on a particle independently of its location within the disk. Here, ${\overline{v}}_{\text{th}}$ represents the mean gas thermal velocity. Dust particles were introduced at the beginning of the full 3D simulation in thermal equilibrium, and they evolve with two different integrators depending on their τ_s (see Picogna et al. 2018, for further details). We do not consider the effect of the disk self-gravity on the particle evolution. Particles that leave the computational domain at the inner boundary re-enter at the outer boundary. This solution allowed us to keep a constant number of particles, which is beneficial from a numerical point of view, and it does not affect the load balance between the computational cores. Adding particles at a constant rate from the outer boundary would only make sense if we were considering the evolution of small planets for longer times. As shown in Fig. 5, the disk close to the planet location does not run out of particles within the simulated time even for the fast evolving pebbles. Accreted particles are flagged but were otherwise kept in the simulations (see Sect. 3 for further details).

2.3 Planets

We embed a planet, with a mass in the range of [5, 10, 100] M_⊕, orbiting a solar mass star on a circular orbit with semi-major axes R_p = 1 (5.2 au) in code units. The planet does not migrate and its mass is kept fixed. Its gravitational potential is smoothed with a cubic expansion inside its Hill sphere $$r_{\text{H}}=r{\left(\frac{M_{\text{p}}}{3M_{\star }}\right)}^{1/3}$$(13)

(Klahr & Kley 2006). In order to obtain the initial conditions, we evolved the disk for ~ 200 orbits, with lower resolution in the azimuthal direction, until its scale height was not changing considerably. This allowed us to start our simulation with the disk in thermal equilibrium. Subsequently we increased the resolution to the one reported in Table 1, and added the dust component. After the dust component has been evolved for six orbits in the computational domain, the planetary mass slowly increases over an additional 20 orbits to allow for a smooth initial phase. Each simulation was evolved for 100 orbital periods of the planet when the disk structure is close to a stationary state and further changes are small, so we do not expect them to affect our results.

Fig. 1 Azimuthal gas velocity in units of the Keplerian speed as a function of radius for the different planetary masses and models. When the gas speed becomes super-Keplerian outside the planet location, the dust-filtration process occurs and the pebble isolation mass is reached. The purple line shows the corresponding simulation from Picogna et al. (2018) at the same time span. The 100 M⊕ planet corresponds to the right axis and it is shown with a green dotted line.

3 Results

3.1 Gas evolution

The presence of a planet significantly alters the gas structure close to its location even if it is not massive enough to open up a gap (see e.g., Lin & Papaloizou 1993; Lambrechts et al. 2014). In Fig. 1 we show the azimuthal gas velocity in units of the Keplerian speed as a function of radius. In this case, we compared our results with those of Picogna et al. (2018) for a viscous locally isothermal disk with a disk scale height of H∕R = 0.05 and with an isothermal case with H∕R = 0.02. This helps us to understand the influence of radiative cooling. The first thing to notice is that changing the disk aspect ratio changes the pressure support of the gas as well, so that the disk becomes more Keplerian on average, going from H∕R = 0.05 (purple line) to H∕R = 0.02 (black line). This plays an important role in the dust evolution because whenever the gas becomes super-Keplerian, the head-wind felt by the dust particles is reduced, meaning that the particles slow their inward migration down or stop it altogether. Thus, for a colder (thinner) disk, smaller planets can halt the dust evolution more easily. As seen from the parameter study probed, a 5 M_⊕ planet can create a super-Keplerian flow outside its location for an H∕R = 0.02, while the flow remains sub-Keplerian everywhere for the same planet in a hotter disk (H∕R = 0.05).

The radiative case lies in the middle of the two cases presented above even though the initial aspect ratio of the disk was 0.05 in this case and the 5 M_⊕ planet almost pushes the gas to the Keplerian speed (blue line). This points to the fact that the disk was able to effectively cool and adapt to a new equilibrium aspect ratio. This effect is highlighted in Fig. 2 where the two isothermal cases are plotted with purple dotted lines for reference. The radiative case settles to a disk scale height that is close to 0.03 at the planet location, with small perturbances in the planet vicinity depending on its mass. The features in the disk scale height can be directly related to the temperature structure (shown in Fig. 2, lower panel). There, we see that the region close to the location of the planet is significantly hotter than the surrounding disk. We show this effect in more detail in Fig. 3 where we plotted the temperature map at the disk midplane for the 10 M_⊕ planet. For a planet at a larger separation, this effect could be even more pronounced and indirect evidence for the planet (see e.g., Tsukagoshi et al. 2019).

Fig. 2 Upper panel: disk scale height (H∕R, solid lines) as a function of radius for the different planetary masses, 5 M⊕ (blue) and 10 M⊕ (red), for themodels with radiative transfer. The purple dotted lines at 0.05 and 0.02 represent the H∕R of the isothermal simulations performed, respectively. Lower panel: midplane gas temperature (T, solid lines) for the corresponding models as on the top panel. The black dashed line represents the initial temperature.

Fig. 3 Midplane gas temperature (T, solid lines) for the 10 M⊕, close to the planet location. The temperature close to the planet location increases significantly.

Fig. 4 Spatial distribution of the dust particles after 100 planetary orbits for the different planetary masses and for three representative particle sizes. The Stokes numbers at the planet location from left to right are, 0.08, 1.23, and 67.2, respectively.

3.2 Dust distribution

The different particle sizes modeled encompass a wide range of stopping times. In Fig. 4 we show the dust surface number density after 100 planetary orbits for three representative sizes and the three planetary masses. In particular, the middle column shows the distribution of pebble-sized particles, which are defined here as particles with a Stokes number of order unity and correspond to meter-sized objects in our setup.

The 5 M_⊕ planet (first row) is not massive enough to open up a gap in the gas and particle disk as seen for the small (coupled) dust in the left column. The only exception is for the planetesimal-sized objects (right column), which do not feel a strong gas drag, and the planetary core can gravitationally perturb their orbits, depleting the co-orbital region (Dipierro & Laibe 2017). In the middle row, on the other hand, the 10M_⊕ planet has stopped the inflow of pebble-sized particles (middle column), creating an overdensity at the location of the outer pressure bump and a depleted inner disk (Paardekooper & Mellema 2006; Rice et al. 2006; Pinilla et al. 2012). The 100 M_⊕ planet is able to carve a deep gap in both the gas and the small particles (last row, first column), where the spiral arms launched by the planet are visible as well. The gap edges also become Rossby wave unstable, generating vortices (Lovelace et al. 1999), as is visible in the larger particle distribution. In the pebble-sized objects (last row, middle column), the distribution of solids shows two concentric rings in the inner side of the planet caused by the perturbance in the gas Keplerian speed by the massive planets, which generates two confined regions where the gas becomes super-Keplerian (see Fig. 1).

Fig. 5 Radial distribution (averaged over azimuth) of the dust particles as a function of time for the different planetary masses and for three representative particle sizes. The 10 M⊕ planet is able to stop the pebbles efficiently, thus creating a transition disk, suggesting that the pebble isolation mass can be less than 10 M⊕. The Stokes numbers at the planet location from left to right are 0.08, 1.23, and 67.2 after 100 planetary orbits, respectively.

3.3 Dust evolution

In Fig. 5 we show the temporal evolution of the radial distribution for the same dust particles presented earlier as a function of planetary mass. The plot shows that a quasi-equilibrium state has been obtained for the dust distribution at the end of the simulated period. The 5 M_⊕ planet (first row) is not able to stop the radial inflow of the well-coupled dust particles (left column), while a small overdensity is building up outside its location for the pebble-sized particles (middle column). On the other hand, a 10 M_⊕ planet is able to stop the flux of pebble-sized objects efficiently (second row, middle column), thus creating a transition disk where the inner disk is depleted of dust as they are not able to cross the pressure bump created by the planet. The mass at which a planet is able to effectively stop the influx of pebble-sized objects is called pebble isolation mass (see e.g., Morbidelli & Nesvorny 2012). This effect becomes even stronger for the 100 M_⊕ planet (last row) where already after 60 planetary orbits, the pebble particles (middle column) inside the planet location are only present in two confined regions (as described earlier). In looking at the right column, one can see the strong mass dependence of the gap width for planetesimal-sized objects as well (Ayliffe et al. 2012; Weber et al. 2018).

3.4 Solid accretion

The solid accretion rate onto the planetary cores is computed in a post-processing phase using two different criteria, depending on the ratio between their Stokes number and the time they need to cross the planet’s Hill sphere (see Picogna et al. 2018, for a detailed description). A sample of particles for each size was chosen outside the planet location and followed for 50 planetary orbits, after which we checked the fraction of accreted and nonaccreted particles (see Fig. 10 of Picogna et al. 2018). Mass is not added to the planet, and even if accreted particles can remain close to the planet or reach a high speed, they do not influence the simulation since the interaction between dust particles and back reaction onto the gas are neglected. We verified that the total solid accretion rate reached a quasi-equilibrium state after 100 planetary orbits when the simulation was stopped.

We can then compute an effective accretion efficiency, P_eff, that is, the number of accreted particles onto the planet divided by the number of particles that would otherwise drift across the location of the planet in an unperturbed disk (Ormel & Klahr 2010) $$P_{\text{eff}}=\frac{{\dot{M}}_{\text{acc}}}{{\dot{M}}_{\text{drift}}},$$(14)

where Ṁ_acc is the measured accretion rate through the criteria explained above, while Ṁ_drift is the particle drift through the disk Ṁ_drift = 2πrΣ_pv_drift, where Σ_p is the particle surface density and v_drift is the unperturbed dust radial drift speed (see e.g., Nakagawa et al. 1986).

In Fig. 6, we plotted the accretion efficiency for two different planetary masses, 5 M_⊕ and 10 M_⊕, and we excluded the 100 M_⊕ planet because it is capable of carving a gap very quickly, thus strongly depleting the solid accretion for all the studied particle sizes. Again,we compare the radiative case (blue) with two locally isothermal ones with different aspect ratios (red & black). As already pointed out by Picogna et al. (2018), we can only trust the middle part of this plot (in the Stokes number range from 0.1 to 10) because for larger or smaller particles, the planet moves faster than the dust and the accretion efficiency cannot be computed correctly with our method, since the sample of particles outside the planet location chosen for the efficiency calculation is not able to cross its location. Therefore, we highlight the relevant region with a gray color.

For the 5 M_⊕ planet (Fig. 6, left panel), the accretion efficiency in the radiative case (blue line) has a similar steepness but it has increased with respect to the locally isothermal case with H∕R = 0.05 (red line), by a factor of three to four; whereas, for the case with a lower aspect ratio, the “pebble isolation mass” has already been reached, effectively stopping the inflow of pebble-sized particles. The difference becomes striking for the 10 M_⊕ planet (Fig. 6, right panel). The pebble isolation mass is reached for the radiative case (blue line), and the overall accretion efficiency is reduced. While in the isothermal case with a higher aspect ratio (red line), the planet is still accreting pebble-sizedobjects efficiently. The reason for this difference can be attributed to the change in the disk aspect ratio at equilibrium for the radiative disk compared to the locally isothermal ones (see Fig. 2). This effect has been predicted byBitsch et al. (2018), who found that the pebble isolation mass scales strongly with the disk aspect ratio for a locally isothermal disk $$M_{\text{iso}}\propto {\left(\frac{H/R}{0.05}\right)}^3.$$(15)

In Fig. 7 we compare the pebble isolation mass found for our simulations (with error bars given by the bin size in planetary mass used) with the analytical prescription from Bitsch et al. (2018), and we find good agreement when taking the reduced disk scale height at equilibrium for the radiative case into account.

We fit the accretion efficiency with an analytical prescription in the Stokes number range considered, (see green line in Fig. 6) as $$P_{\text{eff}}=\frac{3}{\pi \eta }{\left(\frac{{\tau }_s}{0.1}\right)}^{-2/3}{\left(\frac{r_{\text{H}}}{r}\right)}^2$$(16)

where r_H is the Hill radius, and $$\eta =-\frac{1}{2}{\left(\frac{H}{r}\right)}^2\frac{\partial \ln P}{\partial \ln r}=8.56\times {10}^{-4}$$(17)

is the value at 5.2 au in our simulation, just before the planet was included.

This relation is similar to the one obtained by Lambrechts & Johansen (2014, see their Eq. (33)), but it shows a steeper dependence on the Stokes number τ_s. Their analytical derivation gives a correlation to the Stokes number as $P_{\text{eff}}\propto {\tau }_{\text{s}}{}^{-1/3}$. There are two main reasons to explain why our study yields a steeper relation. First of all, their relation only holds up to τ_s = 0.1. This is because the ratio of the pebble accretion rate (see Eq. (28) of Lambrechts & Johansen 2014), given by $\stackrel{.}{M_{\text{c}}}\propto {\tau }_{\text{s}}^{2/3}$, and the pebble flux through the planet $\stackrel{.}{M_{\text{f}}}\propto {\upsilon }_{\text{r}}\propto \frac{{\tau }_{\text{s}}}{{\tau }_{\text{s}}^2+1}\simeq {\tau }_{\text{s}}$ gives this proportionality. However, this last relation breaks down for τ_s > 0.1. In our case, the proportionality of − 2∕3 holds for all of the values of τ_s that are in agreement with our method (τ_s from 0.1–1.0, gray area in Fig. 6). Furthermore, we are treating a full 3D disk with turbulent kicks on the dust particles, which is difficult to model using an analytic prescription, as in Lambrechts & Johansen (2014). In our 3D simulations, the small particles are lifted to a higher disk scale height compared to the larger particles that settle in the midplane. Hence, the accretion rate increases for smaller dust particles, and a steeper gradient of the efficiency as a function of the particle Stokes number develops.

Fig. 6 Efficiency of accreted particles as a function of Stokes number for different planetary masses. The blue line (rad) refers to the new simulations performed with radiative transfer, the black line corresponds to the isothermal case with H∕R =0.02, and the red line (H∕R = 0.05) represents the isothermal case from Picogna et al. (2018). The fit from Eq. (16) is overplotted with a green line. We shaded the intermediate τs values in gray for which the method adopted to calculate the efficiency holds (see text).

Fig. 7 Pebble-isolation mass as a function of α and for different aspect ratios H∕R. The lines show the pebble isolation mass for different α values and they follow the fit from Bitsch et al. (2018). Squares and error bars correspond to the data from Picogna et al. (2018) and the new models with radiative transfer.

4 Conclusions

We can summarize the main results of this study with the following points.

• The main effect of radiative transfer within the bulk of the disk is to produce a cooler and thinner disk (lower H∕R, see also Kley et al. 2009). This allows the growing protoplanet to enhance its efficiency at halting the flux of pebbles and decrease the pebble isolation mass. Since this limit becomes lower than 10 M_⊕, which corresponds to the mass at which the core can start its runaway gas accretion, this result can simply explain the observed, relatively low gas giant planet frequency.

• We found a new relation for the pebble efficiency as a function of the particle Stokes number for a 3D disk (see Eq. (16)), where the dependence on the Stokes number is steeper than the one obtained by the analytical derivation of Lambrechts et al. (2014).

• The planetary cores at the end of the disk lifetime would have had time to grow until this limit, which is consistent with thefinding that the majority of exoplanets are in the mass range between Earth and Neptune.

• Decreasing the limit of planetary masses that are able to stop the inflow of pebbles from the outer to the inner disk, with respect to the planet location, can help to explain the population of transition disks observed around young stars (with high accretion rates), in a natural way.

The main limitations of this study are the lack of dust back-reaction, particle growth, and irradiation from the central star. Dust coagulation and fragmentation can increase the fraction of solids crossing the planet location by grinding them down at the pressure bump and regrowing them in the inner disk (Dra̧żkowska et al. 2019). Along the same lines, dust back-reaction would modify the pressure gradient outside the planet location, reducing the planet efficiency in stopping the pebble flux (Kanagawa et al. 2018). In this respect, our work, which does not take these two effects into account, underestimates the pebble isolation mass. However, Dra̧żkowska et al. (2019) found that the effect of back-reaction can be strongly reduced when taking a full coagulation model into account. Stellar irradiation is expected to play an important role in modifying the disk aspect ratio and potentially the pebble isolation mass. This effect is, nevertheless, only stronger in the outer, irradiation dominated, disk. We looked in detail to understand the impact of stellar irradiation in Appendix A, confirming that its contribution at the location studied is negligible. We did focus on a planet at 5.2 au where viscous heating is still the main heating term (see e.g., Ziampras et al. 2020), and so our results can be safely generalized to a more complex model.

The fact that our conclusions do not hold for larger radii could imply that giant planets only form in the irradiation dominated region, and the signature of this effect can still be present in the bulk composition of the observed giant planets. Finally, we modeled a laminar viscous disk. It has been shown that, even in magnetically inactive regions, hydrodynamical instabilities can develop (see e.g., Stoll et al. 2017, for the vertical shear instability, VSI). As seen by Stoll & Kley (2016), VSI has a significant effect on reducing the effective viscosity and the disk scale height while increasing the dust vertical spreading. Its effect on solid accretion rates and isolation masses has been studied in the context of locally isothermal disks in Picogna et al. (2018) who found a negligible difference. However, for radiative disks, the damping effect of VSI on the disk scale height can further strengthen our main conclusions.

Acknowledgements

We acknowledge funding by the DFG Research Unit “Transition Disks” (FOR 2634/1, ER 685/8-1, KL 650/29-1), and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 714769), and by the German Research Foundation cluster of excellence ORIGINS (EXC 2094, www.origins-cluster.de). The authors gratefully acknowledge the compute and data resources provided by the Leibniz Supercomputing Centre (www.lrz.de). We thank Til Birnstiel and Joanna Dra̧żkowska for proofreading the paper, and the anonymous referee for the constructive feedback.

Appendix A Importance of stellar irradiation

In this work, we ignore the contribution of the stellar irradiation to the disk temperature structure. Since one of our main resultsis the reduction of the pebble isolation mass due to the decreased disk aspect ratio, we tested our assumption. We performed an extra run (without planet and dust particles) including stellar irradiation (with T_⋆ = 5772 K), until thermal equilibrium was reached at the planet location. As shown in Fig. A.1, the disk scale height at the planet location is unaffected by the contribution of the stellar irradiation (red line) when comparing it to the run for the 5 M_⊕ planet (and no stellar irradiation, blue line). On the other hand, this effect plays a major role in an extended disk (see Fig. A.2).

Fig. A.1 Disk scale height for a stellar irradiated disk (red) compared with the disk scale height of the run with a 5 M⊕ planet (blue).

Fig. A.2 Disk scale height for a stellar irradiated disk with a 100 au outer radius (red) compared with the disk scale height of the run with a 5 M⊕ planet (blue).

References

All Tables

All Figures

Fig. 1 Azimuthal gas velocity in units of the Keplerian speed as a function of radius for the different planetary masses and models. When the gas speed becomes super-Keplerian outside the planet location, the dust-filtration process occurs and the pebble isolation mass is reached. The purple line shows the corresponding simulation from Picogna et al. (2018) at the same time span. The 100 M⊕ planet corresponds to the right axis and it is shown with a green dotted line.

In the text

Fig. 2 Upper panel: disk scale height (H∕R, solid lines) as a function of radius for the different planetary masses, 5 M⊕ (blue) and 10 M⊕ (red), for themodels with radiative transfer. The purple dotted lines at 0.05 and 0.02 represent the H∕R of the isothermal simulations performed, respectively. Lower panel: midplane gas temperature (T, solid lines) for the corresponding models as on the top panel. The black dashed line represents the initial temperature.

In the text

	Fig. 3 Midplane gas temperature (T, solid lines) for the 10 M⊕, close to the planet location. The temperature close to the planet location increases significantly.
In the text

	Fig. 4 Spatial distribution of the dust particles after 100 planetary orbits for the different planetary masses and for three representative particle sizes. The Stokes numbers at the planet location from left to right are, 0.08, 1.23, and 67.2, respectively.
In the text

Fig. 5 Radial distribution (averaged over azimuth) of the dust particles as a function of time for the different planetary masses and for three representative particle sizes. The 10 M⊕ planet is able to stop the pebbles efficiently, thus creating a transition disk, suggesting that the pebble isolation mass can be less than 10 M⊕. The Stokes numbers at the planet location from left to right are 0.08, 1.23, and 67.2 after 100 planetary orbits, respectively.

In the text

Fig. 6 Efficiency of accreted particles as a function of Stokes number for different planetary masses. The blue line (rad) refers to the new simulations performed with radiative transfer, the black line corresponds to the isothermal case with H∕R =0.02, and the red line (H∕R = 0.05) represents the isothermal case from Picogna et al. (2018). The fit from Eq. (16) is overplotted with a green line. We shaded the intermediate τs values in gray for which the method adopted to calculate the efficiency holds (see text).

In the text

	Fig. 7 Pebble-isolation mass as a function of α and for different aspect ratios H∕R. The lines show the pebble isolation mass for different α values and they follow the fit from Bitsch et al. (2018). Squares and error bars correspond to the data from Picogna et al. (2018) and the new models with radiative transfer.
In the text

	Fig. A.1 Disk scale height for a stellar irradiated disk (red) compared with the disk scale height of the run with a 5 M⊕ planet (blue).
In the text

	Fig. A.2 Disk scale height for a stellar irradiated disk with a 100 au outer radius (red) compared with the disk scale height of the run with a 5 M⊕ planet (blue).
In the text