© ESO, 2017

1. Introduction

In the classical core-accretion scenario (Bodenheimer & Pollack 1986; Pollack et al. 1996), the formation of planets begins on the smallest scale, and lasts from sub-micron sized dust to Jupiter-like giant planets. This scenario successfully explains the formation of both terrestrial and giant planets. One key element of the theory is the formation of the meter-sized bodies via dust coagulation, since these objects are the building blocks of planetesimals. Terrestrial planets, and solid cores of giant planets form by subsequent collision and accretion of planetesimals. The giant planet formation is completed with a rapid gas accretion from the ambient disk in a process known as gaseous runaway, when the mass of the envelope equals the mass of the core when it has ~ 10 M_⊕. Despite its attractive completeness in predicting both the formation of terrestrial and giant planets, the classical core-accretion still has two groups of unresolved problems: the existence of the various barriers hindering the growth of dust particles, and the timescale problem of the giant planet formation.

In this paper we consider the latter problem, which is due to the combination of the relatively short lifetime of a protoplanetary disk with the fast type I inward migration of massive solid cores of giant planets. Observations indicate that the lifetime of gaseous protoplanetary disks is approximately a few million years (Haisch et al. 2001; Mamajek 2009; Pfalzner et al. 2014). Consequently, an approximately 10 M_⊕ solid core should form early enough to be able to rapidly accrete gas reaching the giant planet phase. However, during its growing phase the planetary core is also subject to the fast type I migration (Ward 1997), whose rate is linearly proportional to the core’s mass. According to theoretical models, the time needed to reach the critical mass for a planetary core may exceed its migration timescale. Recent new developments also show the possibility of outward type I migration (see Paardekooper et al. 2010, 2011). However, rapid inward or outward type I migration is still a danger to the full development of a Jupiter-like planet. In both cases the formation of a giant planet might be inhibited, because fast inward migration results in the quick loss of the embryo, while due to a rapid outward migration the embryo quickly reaches the outer part of the disk, where the surface density of planetesimals is very low, thus the timescale of the planetesimal accretion becomes very long, too.

In order to overcome the above problems, an extension to the classical core-accretion scenario has recently been suggested, namely, that there are particular places for planet growth. These places are the planet traps, where the torque responsible for the type I migration vanishes. If the location of a planet trap coincides with or is close to a local pressure maximum of gas where the inward radial drift of dust and planetesimals is stalled, the growing protoplanet may increase its mass quickly by accreting the solid material accumulated there. In the following, we deal with this situation by investigating the mass accretion process of a growing core. We intend to demonstrate with this investigation that the formation time of a giant planet at a density/pressure maximum is significantly reduced, being much shorter than the lifetime of the protoplanetary disk.

Our paper is organized as follows: first we describe the physical background behind the development of planet traps and planetesimal accumulation. In Sect. 3 we present our model of giant planet growth in a time-evolving gas and planetesimal disk having a dead zone. In Sect. 4 the paper continues with the description of our result, and finally closes with a conclusion and summary.

2. Physical background: the concept of the density/pressure maximum

The theoretical existence of a planet trap was first reported by Masset et al. (2006) demonstrating that a steep surface density jump (being a maximum in the gas surface density) in a protoplanetary disk can halt the type I migration of planetary cores of several Earth masses. Moreover, assuming the equation of state for gas in the form $\mathit{P}\mathrm{=}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, c_s being the local sound speed, a pressure maximum also develops at the surface density maximum. Since at a pressure maximum the gas orbital velocity becomes circular Keplerian, the drag force felt by a particle disappears, therefore dust and the most sensitive planetesimals to aerodynamical drag can accumulate there. A pressure maximum itself can be a place where planetesimals are born, since approaching the pressure maximum, the relative velocities between dust particles decrease, thus their collisions will no longer be destructive. Moreover, at a pressure maximum, several barriers of dust growth can be overcome via coagulation and sweep-up growth of small particles (Brauer et al. 2008; Windmark et al. 2012; Dra¸żkowska et al. 2013).

The ideal candidates for the developments of density/pressure maxima are the inner and the outer boundaries of the dead zone. It is widely accepted that gas accretion in a protoplanetary disk is driven by the turbulence caused by the magneto-rotational instability (MRI), which needs ionized plasma to be triggered (Balbus & Hawley 1991). However, only the very inner part of the disk is ionized (weakly) in its full column by thermal ionization of alkali metals, which cannot be sustained below temperatures of 900 K. In the absence of thermal ionization, the other ionization sources could be the X-rays from the stellar magnetosphere, or cosmic rays. At high gas-surface density the gas is self-shielded against these ionization sources, therefore gas accretion happens only in an upper layer (Gammie 1996). This part of the disk with reduced accretion is called the dead zone. When the surface density drops as the function of the distance from the star, the external ionizing radiation can penetrate through the disk, which is ionized again in its full column. Due to the reduced/enhanced accretion at the boundaries of the dead zone, density and consequently pressure maxima will develop. If certain physical conditions are fulfilled (Lyra et al. 2015), these density maxima can be manifested as large-scale vortices, being prone to the Rossby wave instability (Lovelace et al. 1999). With the current observation techniques, these vortices, if they are far enough from the star, can be observed, and indeed, the horseshoe-shaped patterns discovered in dust continuum by the Sub-millimeter Array (Brown et al. 2009) and recently by the Atacama Large Millimeter/Submillimeter Array (van der Marel et al. 2013; Casassus et al. 2013; Fukagawa et al. 2013; Isella et al. 2013; Pérez et al. 2014) might be attributed to the large dust collecting Rossby vortices (Regály et al. 2012) indifferently from whether they are assumed to be generated at the outer edge of the dead zone and the gap opened by a giant planet.

Another candidate for the development of a pressure maximum might be the water snowline (Kretke & Lin 2007; Brauer et al. 2008). The water snow line separates the regions of the disk in which the water is in vapour and solid form. As the temperature drops below a certain limit (being approximately 170 K for typical disk conditions) the water vapour condenses out, possibly to the surface of silicate dust grains. Thus, when crossing the snow line moving away form the star, the solid-to-gas ratio increases suddenly, affecting the strength of the MRI-driven turbulence as well as the number of free electrons, and thus the conductivity of the plasma also decreases suddenly (Sano et al. 2000; Ilgner & Nelson 2006). This change in the gas turbulence is also reflected in the gas accretion rate, being lower outside the snow line than inside of it. Therefore, similarly to the inner edge of the dead zone, a density and a corresponding pressure maximum may appear.

In our work, we consider two pressure maxima, one at the water snow line, and the other at the outer edge of the dead zone. We do not consider the inner edge of the dead zone as a possible location for giant planet formation, since the high temperature T ~ 1000 K would certainly inhibit the effective cooling and collapse of the gas envelope, which is needed to form a giant planet.

3. Our planet formation model

According to the classical core accretion scenario, the formation of giant planets begins with the sedimentation and coagulation of dust to the disk’s midplane, which is followed by the formation of planetesimals being larger objects (in the size regime ≲ 100 km) less sensitive to drag from the ambient gas disk. According to the most recent picture, planetesimals are the outcome of a process called as gravoturbulent formation. During this process, mm to cm sized dust grains are concentrated by transient high pressure regions in sufficiently high amount to trigger the streaming instability followed by gravitational collapse of the dust aggregates. Further collisions of planetesimals lead to the formation of planetary embryos (among them the larger ones can become the solid cores of giant planets), which continue growing by the accretion of planetesimals and the surrounding gas. As mentioned in Sect. 1, there is a timescale problem associated with the formation of giant planets, meaning that the building time of a giant planet is very close to the lifetime of protoplanetary disks. In recent years, however, several works, including different physical phenomena, demonstrated the possibility of the formation of giant planets before the disk dissipation. For example, Hubickyj et al. (2005) showed that a reduction of the grain opacity of the planet’s envelope significantly reduces the formation time of giant planets. On the other hand, Hori & Ikoma (2011) and Venturini et al. (2015, 2016) found that the pollution by icy planetesimals and the consequent enrichment of the envelope of the planet significantly reduces the critical core mass and speeds up the formation of giant planets. Moreover, models for planetary population synthesis based on the classical core accretion scenario can reproduce several of the main observed characteristics of the exoplanet population (Ida & Lin 2004; Ida et al. 2013; Alibert et al. 2013).

Concomitantly, in the past few years a new alternative model, which is also based on the accretion of solid material, has been proposed for the formation of giant planets. The basic assumption of this model is that the core of a giant planet can be formed rapidly as a few hundred km-sized body accretes cm-sized particles, known as pebbles. Ormel & Klahr (2010), Lambrechts & Johansen (2012) and Guillot et al. (2014) showed that pebbles, strongly coupled to the gas (with Stokes numbers lower than unity), can be accreted very quickly and efficiently to form massive cores. The main difference between pebble accretion and planetesimal accretion is that pebbles can be accreted by the full Hill sphere of the growing core while planetesimals can only be accreted by a fraction of the Hill sphere. Thus, pebble accretion rates can be significantly larger than planetesimal accretion rates. Lambrechts et al. (2014), Levison et al. (2015), and Bitsch et al. (2015) showed that solar system giant planets could be formed by the pebble accretion mechanism.

According to (Chambers 2016), one important parameter of pebble accretion is the amount of the remaining pebbles after planetesimal formation has taken place. If at the onset of pebble accretion, the leftover pebble population is still significant, planetary systems with multiple gas giants beyond the snow line and small planets closer to the star can be formed. Otherwise, pebble accretion would not be effective enough, and no giant planets would be formed. In the latter case, the largest bodies have comparable sizes to Earth. Moreover, the outcome of planet formation is also sensitive to the sizes of planetesimals that form as a result of gravoturbulent collapse. If the largest planetesimals do not overgrow the critical size of 300 km before the depletion of the cm-sized pebble population, giant planet formation is also inhibited. Additionally, the formation size of planetesimals is still a debated issue. According to Morbidelli et al. (2009), large planetesimals form, having characteristic sizes of ~ 100 km. On the other hand, the accretion of sub-km sized planetesimals cannot be ruled out (Weidenschilling 2011).

In the present study we investigate giant planet formation, considering the growth of the protogiant cores by both accretion of planetesimals in the size interval between 0.1–100 km, and by cm-sized pebble accretion. As we show below, giant planet formation is possible in both cases as well as in cases where pebble accretion would be ineffective.

In a series of previous works (Guilera et al. 2010, 2011, 2014), we developed a model which calculates the formation of gaseous giant planets embedded in a time-evolving protoplanetary disk. In this work, we incorporate some modifications to our previous model with the aim of studying the formation of massive cores (which are the precursors of giant planets) at the pressure maxima of protoplanetary disks.

In our model, the protoplanetary disk is characterized by two components: a gaseous component, evolving due to an α-viscosity-driven accretion, and a solid component represented by a planetesimal population being subject to accretion by the planets and radial drift due to gas drag. The protogiant planets embedded in the disk grow by accretion of planetesimals and gas.

3.1. Initial radial profiles for gas and solid material

In our disk model, the computational domain is defined between 0.1 au and 1000 au, using 5000 radial bins equally logarithmically spaced, as using a classical 1D radial model. The gaseous component is characterized by the corresponding surface density Σ_g(R), where R is the radial coordinate, and the solid component is characterized by the planetesimal surface density Σ_p(R). Moreover, our model allows us to study a discrete planetesimal size distribution as well, thus in a more general view the planetesimal surface density can be characterized by Σ_p(R,r_p), where r_p represents the different sizes of the discrete distribution.

In order to define the initial surface density profiles, we follow the suggestions of Andrews et al. (2010), who while studying the Ophiuchus star-forming region found that the gas surface density of the disks observed can be represented by $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{g}}& \mathrm{=}& \end{array}$(1)where R_c is a characteristic radius, γ represents the surface density gradient, and ${\mathrm{\Sigma }}_{\mathrm{g}}^{\mathrm{0}}$ is a parameter function of the disk mass, $\begin{array}{ccc}{\mathit{M}}_{\mathrm{d}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{2}\mathit{\pi R}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{d}\mathit{R}\mathit{.}& & \end{array}$(2)Integrating Eq. (2), one can find that ${\mathrm{\Sigma }}_{\mathrm{g}}^{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{2}\mathrm{-}\mathit{\gamma }\mathrm{)}{\mathit{M}}_{\mathrm{d}}\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\right)}$.

A common assumption of planet formation models is that the metalicity along the disk is the same as that of the central star, and that dust sediments and coagulates very quickly to form a mid-plane planetesimal disk. Following this hypothesis, the initial planetesimal surface density is given by $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{p}}& \mathrm{=}& \end{array}$(3)where η_ice(R) takes into account the sublimation of water-ice given by $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{ice}}\mathrm{=}\{\begin{array}{c}\\ & \\ \\ & \end{array}& & \end{array}$(4)with R_ice = 2.7 au being the ice line (or snow line). For the solar system, the factor β can take values between ~ 2 and ~ 4 (Hayashi 1981; Lodders 2003). In this work, we adopt a value of β = 3. For numerical convenience, we smoothed the discontinuity at R = R_ice by $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{ice}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\beta }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\mathit{\beta }}\right)\left[\mathrm{1}\mathrm{+}\tanh \left(\frac{\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{ice}}}{{\mathrm{\Delta }}_{\mathrm{ice}}}\right)\right]\mathit{,}& & \end{array}$(5)where Δ_ice = H_g(R_ice), H_g being the scale height of the gas disk.

On the other hand, ${\mathrm{\Sigma }}_{\mathrm{p}}^{\mathrm{0}}$ is given by $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{p}}^{\mathrm{0}}\mathrm{=}{\mathit{z}}_{\mathrm{0}}{\mathrm{\Sigma }}_{\mathrm{g}}^{\mathrm{0}}\mathit{,}& & \end{array}$(6)with z₀ = 0.0153 being the initial abundance of heavy elements, referred to also as the dust-to-gas ratio (Lodders et al. 2009). Finally, we adopted some typical values for γ and R_c given by Andrews et al. (2010): γ = 1 and R_c = 25 au.

3.2. Evolution of the gas disk with a wide-range viscosity reduction

Fig. 1 Alpha-viscosity parameter as a function of the radial coordinate for a flat and a flared disk, with αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au. We consider a disk of 0.05 M⊙, and for both disk cin - dz = cout - dz = 1. For a flat disk, Hg(Rin - dz) = 0.135 au and Hg(Rout - dz) = 1 au, while for a flared disk Hg(Rin - dz) = 0.165 au and Hg(Rout - dz) = 2.013 au.

As we mentioned above, the gas surface density of the disk Σ_g evolves as an α accretion disk (Pringle 1981) $\begin{array}{ccc}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{=}\frac{\mathrm{3}}{\mathit{R}}\frac{\mathit{\partial }}{\mathit{\partial R}}\left[{\mathit{R}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial R}}\mathrm{(}\mathit{\nu }{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{R}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\right]\mathit{,}& & \end{array}$(7)where ν = αc_sH_g is the kinematic viscosity given by the dimensionless parameter α (Shakura & Sunyaev 1973). Usually, the parameter α ~ 10^-3...10^-2 is a constant value along the disk. In order to reproduce the effect of the dead zone, we have chosen α to be a particular function of R (see later). The sound speed is given by $\begin{array}{ccc}{\mathit{c}}_{\mathrm{s}}\mathrm{=}\sqrt{\frac{{\mathit{\gamma }}_{\mathrm{g}}{\mathit{k}}_{\mathit{B}}\mathrm{T}}{{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}}\mathit{,}& & \end{array}$(8)where γ_g = 5/3, k_B is the Boltzmann-constant, and μ_H2 and m_H2 are the molecular weight and mass of molecular hydrogen, respectively. The radial temperature profile is given by the following power-law function $\begin{array}{ccc}\mathit{T}\mathrm{=}\mathrm{280}{\left(\frac{\mathit{R}}{\mathrm{1}\mathrm{au}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{K}\mathit{.}& & \end{array}$(9)Regarding the geometry of the disk, we consider two cases, a flat and a flared disk. In the flat case, the aspect ratio is constant h = 0.05, so H_g = 0.05 R, while for the flared disk we assume that H_g = c_s/ Ω_k ∝ R^5/4, Ω_k being the Keplerian frequency.

In order to generate the inner and the outer pressure maximum at the water snow line and at the outer edge of the dead zone, respectively, we apply a reduction of the α parameter between them. Denoting by α_back the background α-viscosity parameter, and by α_dz, its reduced value, the functional form of α(R) is given as $\begin{array}{ccc}\mathit{\alpha }\mathrm{\left(}\mathit{R}\mathrm{\right)}& \mathrm{=}& [\left({\mathit{\alpha }}_{\mathrm{back}}\mathrm{-}{\mathit{\alpha }}_{\mathrm{d}\mathit{z}}\right)\{\mathrm{1}\mathrm{-}\mathrm{0.5}[\mathrm{1}\\ & & \mathrm{+}\tanh \left(\frac{\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{in}\mathrm{-}\mathrm{d}\mathit{z}}}{{\mathit{c}}_{\mathrm{in}\mathrm{-}\mathrm{d}\mathit{z}}{\mathit{H}}_{\mathrm{g}}\mathrm{\left(}{\mathit{R}}_{\mathrm{in}\mathrm{-}\mathrm{d}\mathit{z}}\mathrm{\right)}}\right)]\}\mathrm{+}\frac{{\mathit{\alpha }}_{\mathrm{d}\mathit{z}}}{\mathrm{2}}]\\ & & \mathrm{+}[\left({\mathit{\alpha }}_{\mathrm{back}}\mathrm{-}{\mathit{\alpha }}_{\mathrm{d}\mathit{z}}\right)\{\mathrm{1}\mathrm{-}\mathrm{0.5}[\mathrm{1}\\ & & \end{array}$(10)where R_{in - dz} and R_{out - dz} are the locations of the water snow line and outer edge of the dead zone, respectively, and c_{in - dz} and c_{out - dz} are constants that define the width of the transition in the viscosity profile. In the following, we refer to this region as the dead zone, but we should keep in mind that the real dead zone may extend much closer to the star, until the thermal ionization dominates accretion. On the other hand, the MRI driven turbulent viscosity qualitatively behaves similarly at the inner edge of the dead zone and at the water snow line, as the viscosity is suddenly reduced with increasing R.

Figure 1 represents the radial profiles of the alpha-viscosity parameter for a flat and a flared disk. We note that the differences in the widths of the transition regions are due to the fact that they are expressed in terms of the scale height of gas (see Eq. (10)), which is larger for a flared disk than for a flat disk. These differences play an important role in the evolution of the gas and planetesimal surface densities, as is discussed in following sections.

Finally, we mention that Eq. (7) is solved using a full implicit Crank-Nicholson method considering zero torques as boundary conditions. For each time-step, we do not allow changes greater than 10% for the gas surface density in each radial bin.

3.3. Evolution of the planetesimal population

The numerical treatment of the evolution of the planetesimal population is described in detail in Guilera et al. (2014). Here, we will only discuss the most relevant properties and some modifications of it with respect to our previous work. In our model, planetesimals are subject to radial drift due to the drag force arising from the ambient gaseous disk. Moreover, planetesimals are also accreted by the growing protoplanets, and affected by mutual collisions, though this latter effect is not included in our model.

The drag force between planetesimals and the gas depends on the relative velocities between gas and planetesimals and on the ratio between the planetesimal radii and the mean free path of the gas molecules. Similarly to our previous work, we consider three different regimes of the drag force: Epstein, Stokes and quadratic. The separate treatment of these regimes is important, because while big planetesimals are always in the quadratic regime, small planetesimals can change their regimes along the disk. The radial drift velocities of planetesimals are given by $\begin{array}{ccc}{\mathit{v}}_{\mathrm{mig}}\mathrm{=}\{\begin{array}{c}\\ & \\ \\ & \\ \\ & \end{array}& & \end{array}$(11)where t_stop and s are the stopping time and the Stokes number (for the corresponding regime), respectively (see Guilera et al. 2014), and η = (v_g−v_k) /v_k is the fraction by which the gas deviates from the Keplerian circular velocity given by $\begin{array}{ccc}\mathit{\eta }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{{\mathit{H}}_{\mathrm{g}}}{\mathit{R}}\right)}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathit{P}}{\mathrm{d}\ln \mathit{R}}\mathit{,}& & \end{array}$(12)where P is the gas pressure in the midplane of the disk. Considering a local isothermal equation state for the gas, $\mathit{P}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{g}}$, where ${\mathit{\rho }}_{\mathrm{g}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{/}\mathrm{\left(}\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{H}}_{\mathrm{g}}\mathrm{\right)}$ is the volumetric gas density at the midplane. Equation (12) can be expressed as $\begin{array}{ccc}\mathit{\eta }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{{\mathit{H}}_{\mathrm{g}}}{\mathit{R}}\right)}^{\mathrm{2}}\frac{\mathrm{d}\ln \mathrm{\left(}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\right)}}{\mathrm{d}\ln \mathit{R}}\mathit{,}& & \end{array}$(13)thus, if the gas density is a decreasing function of the distance R to the central star, the radial drift of the planetesimals is always inward (Eq. (11)). However, if there is some local maximum in the gas density, planetesimals can also drift outwards, when $\mathrm{d}\ln \mathrm{\left(}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{g}}\mathrm{\right)}\mathit{/}\mathrm{d}\ln \mathit{R}\mathit{>}\mathrm{0}$.

Finally, as a consequence of the mass conservation, the evolution of the planetesimal surface density is described by the following advection equation $\begin{array}{ccc}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\left(}\mathit{R,}{\mathit{r}}_{{\mathit{p}}_{\mathit{j}}}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{-}\frac{\mathrm{1}}{\mathit{R}}\frac{\mathit{\partial }}{\mathit{\partial R}}\mathrm{[}\mathit{R}{\mathit{v}}_{\mathrm{mig}}\mathrm{\left(}\mathit{R,}{\mathit{r}}_{{\mathit{p}}_{\mathit{j}}}\mathrm{\right)}{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{(}\mathit{R,}{\mathit{r}}_{{\mathit{p}}_{\mathit{j}}}{\mathrm{)}}^{\mathrm{]}}\mathrm{=}\mathrm{ℱ}\mathrm{\left(}\mathit{R,}{\mathit{r}}_{{\mathit{p}}_{\mathit{j}}}\mathrm{\right)}\mathit{,}& & \end{array}$(14)where ℱ represents the sink terms due to the accretion by the growing embryos or cores, and r_pj emphasizes the fact that Eq. (14) is solved independently for each planetesimal size, when a planetesimal size distribution is considered. In this case, the total planetesimal surface density is given by $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}}{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\left(}\mathit{R,}{\mathit{r}}_{{\mathit{p}}_{\mathit{j}}}\mathrm{\right)}\mathit{.}& & \end{array}$(15)Equation (14) is solved using a full implicit upwind-downwind mix method considering zero density as boundary conditions. For each time-step, we do not allow changes greater than 10% for the planetesimal surface density in each radial bin.

3.4. Growth of the protoplanets

In our model, the planetary embryos embedded in the disk grow by the concurrent accretion of planetesimals and the surrounding gas (see Guilera et al. 2010, 2014, for a detailed explanation).

The mass increase of a core in the oligarchic growth regime due to the accretion of planetesimals (Inaba et al. 2001) is described by $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{M}}_{\mathit{C}}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\left(}{\mathit{a}}_{\mathrm{P}}\mathrm{\right)}{\mathit{R}}_{\mathrm{H}}^{\mathrm{2}}}{{\mathit{T}}_{\mathrm{orb}}}{\mathit{P}}_{\mathrm{coll}}\mathit{,}& & \end{array}$(16)where M_C is the core’s mass, Σ_p(a_p) is the surface density of solids at the location of the planet, R_H is the Hill radius, and T_orb is the orbital period. P_coll is a collision probability, which is a function of the core radius R_C, the Hill radius of the planet, and the relative velocity between the planetesimals and the planet v_rel, thus P_coll = P_coll(R_C,R_H,v_rel). In fact, as we also consider the drag force that planetesimals experience on entering the planetary envelope (following Inaba & Ikoma 2003), the collision probability is function of the enhanced radius $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{R}}_{\mathrm{C}}\end{array}$ instead of R_C.

The gas accretion rate and the thermodynamic state of the planet envelope are calculated by solving the standard equations of transport and structure (see Benvenuto & Brunini 2005; Fortier et al. 2007, 2009; Guilera et al. 2010, for details), $\begin{array}{ccc}\frac{\mathit{\partial r}}{\mathit{\partial }{\mathit{m}}_{\mathit{r}}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{\rho }}\mathit{,}\\ \frac{\mathit{\partial P}}{\mathit{\partial }{\mathit{m}}_{\mathit{r}}}& \mathrm{=}& \mathrm{-}\frac{\mathit{G}{\mathit{m}}_{\mathit{r}}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{4}}}\mathit{,}\\ & & \\ \frac{\mathit{\partial }{\mathit{L}}_{\mathit{r}}}{\mathit{\partial }{\mathit{m}}_{\mathit{r}}}& \mathrm{=}& {\mathit{ϵ}}_{\mathit{pl}}\mathrm{-}\mathit{T}\frac{\mathit{\partial S}}{\mathit{\partial t}}\mathit{,}\\ \frac{\mathit{\partial T}}{\mathit{\partial }{\mathit{m}}_{\mathit{r}}}& \mathrm{=}& \mathrm{-}\frac{\mathit{G}{\mathit{m}}_{\mathit{r}}\mathit{T}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{4}}\mathit{P}}\mathrm{\nabla }\mathit{,}\end{array}$(17)where ρ is the envelope density, G is the universal gravitational constant, ϵ_pl is the energy release rate by planetesimal accretion, S is the entropy per unit mass, and ∇ = dlnT/ dlnP is the dimensionless temperature gradient, which depends on the type of energy transport.

We also incorporated in our model the prescription of type I migration for a locally isothermal disk to calculate the change in semi-major axis of the planetary embryo $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{a}}_{\mathrm{P}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathrm{2}{\mathit{a}}_{\mathrm{P}}\frac{\mathrm{\Gamma }}{{\mathit{L}}_{\mathrm{P}}}\mathit{,}& & \end{array}$(18)where a_P represents the planetary embryo’s semi-major axis and ${\mathit{L}}_{\mathrm{P}}\mathrm{=}{\mathit{M}}_{\mathrm{P}}\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}{\mathit{a}}_{\mathrm{P}}}$ its angular momentum. Γ is the total torque, which is given by: $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{=}\mathrm{(}\mathrm{1.364}\mathrm{+}\mathrm{0.541}\mathit{\delta }\mathrm{)}{\left(\frac{{\mathit{M}}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}{\mathrm{\Omega }}_{\mathrm{P}}}{{\mathit{M}}_{\mathit{\star }}{\mathit{c}}_{{\mathit{s}}_{\mathrm{P}}}}\right)}^{\mathrm{2}}{\mathrm{\Sigma }}_{{\mathit{g}}_{\mathrm{P}}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{4}}{\mathrm{\Omega }}_{\mathrm{P}}^{\mathrm{2}}\mathit{,}& & \end{array}$(19)where Ω_P, c_sP, and Σ_gP are the values of the Keplerian frequency, the sound speed, and the gas surface density at the position of the planet, respectively (see Tanaka et al. 2002). The factor δ is defined by δ = dlog Σ_g/ dlog R evaluated at R = a_P. To follow the orbital migration and mass growth of a planetary embryo, Eqs. (7), (14), (16), (17), and (18) have to be numerically solved together self-consistently.

Despite the developments of the last years aiming at improving analytic formulae for type I migration in more and more realistic disk models (e.g., Paardekooper et al. 2010, 2011; Bitsch et al. 2013, 2014a,b; Benítez-Llambay et al. 2015), also including the works Dittkrist et al. (2014) and Bitsch et al. (2015), in our study we use the torque prescription given by Eq. (19) to be consistent with our locally isotherm disk model. We also emphasize that in our formation scenario, the embryo mainly increases its mass being trapped in the zero torque location, and the mass growth during its migration is not significant, thus the outcome of simulations may be independent of the migration timescale.

4. Results

The aim of this work is to study the formation of massive planetary cores due to the accumulation of solids in the form of pebbles and planetesimals at pressure maxima developed in the disk. We recall that in this work we assume an inner pressure maximum that appears at the water snow line, and an outer pressure maximum that develops at the outer edge of the dead zone. We performed different sets of simulations varying the mass of the disk, the values of the α viscosity parameter inside and outside of the dead zone, and the size of the pebbles and planetesimals (considering a single-sized pebble/planetesimal population). The combinations of such parameters have also been considered for flared and flat disks. Table 1 summarizes the free parameters used in this work.

4.1. Disk evolution without planets

As a first step, we analyze the evolution of the disk without planets, namely, the evolution of Σ_g(R,t) and Σ_p(R,t). It is important to note that in contrast to previous works (e.g., Matsumura et al. 2009), for the sake of simplicity, we assume that neither the location of the snow line nor the outer edge of the dead zone evolves with time. Our choice for fixed pressure maxima/migration traps is discussed in more detail in Sect. 5. We intend to study the accumulation of planetesimals at the pressure maxima developed due to the viscosity reduction.

In Figs. 2 and 3 the time evolution of the surface density of gas Σ_g(R,t) and planetesimals Σ_p(R,t) are shown for a disk without a dead zone (pebble accumulation is analyzed in Sect. 4.3). We consider four different simulations in each using a single size distribution for planetesimal radii. While the time evolution of Σ_g(R,t) is the same for the four simulations (Fig. 2), there are significant differences in the time evolution of Σ_p(R,t) due to the different drift rates for the different planesimal sizes (Fig. 3). While big planetesimals, of 10 km and 100 km of radius, do not suffer a significant inward drift except in the inner part of the disk, small planetesimals, particularly planetesimals of 100 m of radius, undergo a significant radial drift all along the disk. We note that due to the inner boundary condition there is an accumulation of planetesimals in the inner edge of the disk due to the generation of a gas pressure maximum. We run all our simulations for 5 Myr, a characteristic protoplanetary disk lifetime (Mamajek 2009; Pfalzner et al. 2014). It is generally accepted that EUV/FUV/X-ray photo-evaporation plays an important role in disk dissipation. It has been showed that after a few Myr of viscous evolution, photo-evaporation becomes significant and the protoplanetary disk is dispersed in a time-scale of 10⁵ yr, (Alexander et al. 2006; Gorti et al. 2009; Owen et al. 2011). As we are interested in the first stages of planet formation, particularly in the formation of massive cores until the planet achieves the critical mass (when the mass of the envelope equals the mass of the core), we only considered the viscous evolution of the disk.

Fig. 2 Time evolution of the gas surface density radial profiles. The simulation corresponds to a disk of Md = 0.1 M⊙, using α = 10-3. The simulation is stopped after 5 Myr of viscous evolution.

Fig. 3 Time evolution of the planetesimal surface density radial profiles for planetesimal populations of different radii. Simulations correspond to a disk of Md = 0.1 M⊙, using α = 10-3. Simulations are stopped after 5 Myr of viscous evolution.

As we mentioned before, here we investigate the possibility of the rapid formation of massive planetary cores due to the accumulation of planetesimals around an inner and outer pressure maximum, respectively. To generate the inner and the outer pressure maximum, we implemented the functional form of the viscosity, given by Eq. (10) when numerically solving Eq. (7). Figure 4 shows the time evolution of the gas surface density for a flared disk (top panel) and a flat disk (bottom panel), considering that R_{in - dz} = 2.7 au, R_{out - dz} = 20 au, α_back = 10^-3, and α_dz = 10^-5. At first glance, the time evolution of the gas density profiles is very similar in both cases. The evolution around the inner edge of the viscosity reduction is practically the same for the flared and the flat disk. However, at the outer edge of the dead zone, the time evolution of the gas surface density turns out to be different. We can see that the gas density maximum (and the pressure maximum) disappear more quickly for the case of a flared disk. As we show in following sections, this plays an important role in the formation of massive cores at the outer edges of the dead zones. One reason of this difference is that we express the width of the transition region for the viscosity at the dead zone’s outer edge in terms of a disk’s scale height, which is larger for a flared disk (see Fig. 1). However, even adopting the same width, for example using c_{out - dz} = 0.5 for the flared disk, the pressure maximum also disappears faster than in the case of a flat disk.

The two important locations in the disk that might play an important role in the rapid formation of massive cores are the planet or migration traps and the pressure maxima. The growing cores are trapped at the zero-torque locations, the most drag-sensitive planetesimals are expected to be accumulated at pressure maxima. Figure 5 shows the time evolution of the locations of the migration traps and of the pressure maxima in the neighborhood of the inner (top panel) and outer (bottom panel) edge of the dead zone for the flared and the flat disk. We note that while the migration trap and pressure maximum at the inner edge of the dead zone (both for the flared and the flat disk) survives during the 5 Myr of viscous evolution, they disappear after some time at the outer edge of the dead zone. Usually, the migration trap disappears earlier than the pressure maximum. As we mentioned before, we can also see that the migration trap and the pressure maximum at the outer edge of the dead zone are vanishing first for the flared disk. It is also important to note that for the inner and outer edge of the dead zone, the locations of zero torque and maximum pressure do not exactly coincide. As we show in the following sections, this fact has important consequences for the presented scenario of planet formation.

Fig. 4 Time evolution of the gas surface density radial profiles assuming the existence of a dead zone for a flared disk (top panel) and a flat disk (bottom panel). Simulations corresponding to a disk of Md = 0.1 M⊙ using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5. Simulations are stopped after 5 Myr of viscous disk evolution.

Fig. 5 Time evolution of the zero torque locations (solid lines) and the locations of pressure maxima (dashed lines) at the inner egde (top panel) and outer edge (bottom panel) of the dead zone for a flared disk (red lines) and a flat disk (black lines). Simulations correspond to a disk of Md = 0.1 M⊙ using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5.

Fig. 6 Time evolution of the zero torque locations (solid lines) and pressure maximum locations (dashed lines) at the inner and outer edge of the dead zone for a flat disk using αback = 10-3 and αdz = 10-5 (black lines), and αback = 10-2 and αdz = 10-4 (blue lines).

A no less important parameter of our simulations is the value adopted for the α-viscosity inside and outside of the dead zone. Figure 6 shows the time evolution of the locations of the migration trap and pressure maximum at the inner and the outer edge of the dead zone for a flat disk of M_d = 0.01 M_⊙ using different values for the α-viscosity parameter inside and outside the dead zone. The evolution of the locations of the migration trap and the pressure maximum at the inner edge of the dead zone is practically the same for the different values of α_back and α_dz. However, the migration trap and the pressure maximum at the dead zone’s outer edge quickly disappear (in less than 0.5 Myr) for the case of α_back = 10^-2 and α_dz = 10^-4.

Finally, in Fig. 7 we show the time evolution of the planetesimal surface density radial profiles for different planetesimal sizes when a dead zone is considered in the disk. For smaller-sized planetesimals the accumulation of solids at the inner and outer edge of the dead zone is more significant due to the larger radial drift velocities. We also note that while the location of the planetesimal accumulation at the inner edge of the dead zone seems to be both fixed (or moving slightly outwards) and the same for the different planetesimal sizes, the location of the planetesimal accumulation at the outer edge moves inward at different rates. In fact, this is a consequence of the inward motion of the pressure maximum at the outer edge of the dead zone (see Fig. 5 bottom panel). However, it is interesting to analyze the relative migration between the location of the pressure maximum and the corresponding location of planetesimal accumulation. Planetesimals are trapped when dlnP/ dlnR = 0 (Eq. (12)). However, we can see in Fig. 8 that this situation only happens for planetesimals of 0.1 km, which have the larger radial drift velocities, and not for all the time that the maximum pressure exists in the disk. For larger planetesimals, the pressure maximum moves inwards faster than the average radial velocity of the bodies, thus they cannot be trapped there.

In the following sections, we study whether or not a planet is able to accrete the above mentioned accumulations of planetesimals at the pressure maxima generated by the dead zone.

4.2. Formation of massive cores at pressure maxima

As we have shown in the previous section, the inner pressure maximum survives during the whole evolution of the disk both for flared and flat disks, and, independently of the value of α_back and α_dz, the pressure maximum associated with the outer edge of the dead zone vanishes at some time, which depends on the width of the viscosity transition and on the magnitude of the viscosity reduction. Moreover, while planetesimals are accumulated around the inner pressure maximum, large planetesimals cannot accumulate at the outer edge of the dead zone, since they are drifted inwards slower than the pressure maximum moves.

4.2.1. In situ formation at pressure maxima

First we analyze the in situ formation of giant planets by fixing their positions to the locations of planetesimal accumulation. The inner pressure maximum efficiently accumulates planetesimals. However, as we have shown in the previous section, the pressure maximum generated at the outer edge of the dead zone migrates inward. Consequently, planetesimals with lower inward drift velocities could not accumulate there resulting in that the accumulation of planetesimals with larger sizes does not happen. Moreover, depending on their sizes, the generated planetesimal accumulation deviates differently from the location of the outer pressure maximum (see Fig. 8). Thus, for the outer pressure maximum, we adopted different planet locations for the different planetesimal sizes. In the following, we compare the in situ formation of a planet up until critical mass in a model with and without a dead zone.

As a demonstrative example, we calculate the formation of two planets embedded in a flat disk with mass M_d = 0.1 M_⊙. For a disk without a dead zone, we use α = 10^-3, and for the case with a dead zone we use α_back = 10^-3, α_dz = 10^-5, R_{in - dz} = 2.7 au and R_{out - dz} = 20 au. Initially, both embryos have a core mass of M_c = 0.01 M_⊕ and an envelope mass of ~ 10^-13M_⊕. For all sizes of the planetesimal population, the inner planet has been located at ~ 3.2 au (roughly corresponding to the planetesimal accumulation location), while the outer planet is located at 16 au, 16.5 au, 17 au, and 17.5 au, considering a planetesimal population size of 0.1 km, 1 km, 10 km, and 100 km, respectively. Figures 9 and 10 show the comparison of the core mass growth as a function of time for the inner and outer planet, respectively, between the cases with and without a dead zone for the different planetesimal radii. We run our simulations until the planets achieve the critical, or “cross-over” mass, which happens when the mass of the core is equal to the mass of the envelope, or for 5 Myr of viscous evolution. For the inner planet, except for the case of planetesimals of 100 km radius, the planet achieved the critical mass before 5 Myr when a disk without a dead zone was considered. When a dead zone is considered in the disk, the inner planet achieves the critical mass for all planetesimal sizes. In all of the disk models with a dead zone, the formation times (the time needed for the planet to achieve the critical mass) are significantly shorter than the time needed in a disk model without a dead zone. It is clear that the accumulation of planetesimals at the pressure maxima, and the consequent increment in the planetesimal surface density, significantly favors the formation of massive cores.

Similar results have been found for the outer planet achieving the critical mass in less than 5 Myr for planetesimals of 0.1 km, 1 km, and 10 km radius when a dead zone is considered in the disk. In all of the above cases the accumulation of planetesimals due to the pressure maxima generated by the dead zone clearly favors the formation of massive cores.

We should note however that the above presented in situ scenario does not take into account the planet-disk interaction, so the results obtained are idealizations, considering a maximal growing rate for the planetary core, and therefore the shortest formation times.

4.2.2. Formation with planet migration

The gravitational interaction between the gaseous protoplanetary disk and an embedded planet causes a migration of the planet due to exchange of angular momentum. Low-mass planets, up to a dozen Earth-masses, are subject to type I migration. In order to study the formation of massive cores in the presence of a dead zone, as we mentioned in Sect. 3.4, we incorporated into our model the prescription for the type I migration given by (Tanaka et al. 2002).

We analyze the formation of a giant planet whose core is located at different positions in the disk by the accretion of planetesimals of different sizes (Fig. 11). We initially start the planetary core at 2 au, 5 au, 10 au, 17.5 au, and 22 au. The first and the last position of the planet are outside the dead zone, while the other three locations are inside of it, where we choose two positions near the inner and outer edges of the dead zone. We run simulations for our fiducial model, that is, M_d = 0.1 M_⊙, α_back = 10^-3, α_dz = 10^-5, R_{in - dz} = 2.7 au and R_{out - dz} = 20 au, considering a flat disk.

If the starting position of the planetary core is at 2 au, for planetesimals of 0.1 km and 1 km radius, the planet becomes massive enough to migrate inwards approaching the inner edge of the disk. For planetesimals of 10 km and 100 km radius, the mass of the planet remains small enough to avoid a significant migration.

On the other hand, when the initial location of the core is at 5 au and 10 au, the formation history of the planet is similar. For planetesimals of 0.1 km of radius, the planet quickly becomes massive enough to obey a fast inward migration, and the planet quickly achieves the migration trap, namely the zero torque location, generated by the inner edge of the dead zone. Having been trapped there, the planet is able to accrete the already accumulated planetesimals at the pressure maximum. This phenomenon allows the planet to achieve the critical mass very quickly, in a time-scale of 10⁵ yr. For planetesimals of 1 km and 10 km, the planet grows and migrates inward until it achieves the inner zero torque location in less than 1 Myr. Then the planet continues growing being trapped at the zero torque location until it achieves the critical mass. For planetesimals of 1 km radius, the planet achieves the critical mass in less than 1 Myr, while for planetesimals of 10 km radius the planet achieves the critical mass in a few Myr. Finally, for planetesimals of 100 km radius the planet only achieves a few Earth masses.

When the initial location of the planet is near to the outer edge of the dead zone, in our case at 17.5 au, thus inside of it, the planet initially migrates outward until it becomes trapped at the outer migration trap. As we show in Figs. 5 and 6, the migration trap moves inward due to the diffusive evolution of the outer density/pressure maximum. Thus, during its slow inward motion with the migration trap, the planet gathers mass by accreting planetesimals. However, as shown in Figs. 5 and 6, the location of the pressure maximum deviates significantly from the migration trap’s position at the outer edge of the dead zone over a ~ 1 Myr timescale. Moreover, as shown in Fig. 7, the locations of planetesimal accumulation and the pressure maximum coincide only for the smaller planetesimals, which have the largest radial drift velocities. For these reasons, only for planetesimals of 0.1 km in diameter can the planet reach critical mass before the migration trap vanishes. For the rest of the planetesimal sizes, after the vanishing of the zero torque location, the planet migrates inward until it reaches the inner zero torque location, achieving the critical mass there.

Finally, for planetesimals of sufficiently small size, if the initial location of the planet is 22 au, beyond the outer edge of the dead zone, the planet becomes massive enough to undergo a significant inward migration. After the outer zero torque vanishes, the planet starts to migrate achieving the critical mass at the inner migration trap. We note that the planet significantly increases its mass when passing through the outer planetesimal accumulation location. However, for larger planetesimals, the planet does not become massive enough to perform a substantial migration.

Fig. 7 Time evolution of the planetesimal surface density radial profiles for planetesimal populations of different radii. Simulations correspond to a flat disk of Md = 0.1 M⊙, using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5. Simulations are stopped after 5 Myr of viscous evolution.

4.2.3. Oscillation of the planet’s semi-major axis around the zero torque location

As we have shown in the previous section, the migration of the planet towards the inner zero torque location, and its trapping there, favors the formation of a massive core. However, as we have also demonstrated in previous sections, the positions of the inner pressure maximum and of the migration trap do not coincide. Thus, the accumulation of planetesimals is at a different location than the migration trap.

Also, numerical magnetohydrodynamic simulations have revealed that at the inner edge of the dead zone the planet is not simply trapped but oscillates around the zero torque location (Faure & Nelson 2016). The oscillation of the semi-major axis of a massive core at the outer edge of the dead zone has also been found in the hydrodynamic simulations of Regály et al. (2013). In both cases, beside the turbulence, the oscillation of the planet’s semi-major axis might be the result of the planet-vortex interaction, since the sudden jumps in gas surface density are prone to the Rossby wave instability enabling the formation of a large vortex at the density jump’s position (Lovelace et al. 1999).

Thus, following the above-mentioned works, we mimicked the oscillation of the growing protoplanet around the migration trap. For this, we implemented a noisy oscillation of the semi-major axis of the protoplanet in the form $\begin{array}{ccc}{\mathit{a}}_{\mathrm{P}}\mathrm{=}{\mathit{a}}_{\mathrm{0}}\mathrm{+}\mathit{\psi }\mathrm{\Delta }\mathit{a}\sin \left(\frac{\mathit{t}\mathrm{-}{\mathit{t}}^{\mathrm{\ast }}}{\mathit{\xi }{\mathit{T}}_{\mathrm{orb}}}\right)\mathit{,}& & \end{array}$(20)where a₀ is the zero torque location, Δa is the oscillation amplitude, t^∗ is the time at which the planet achieves the zero torque location, and ψ and ξ are random numbers that take values between 1 ± f and 1 ± g, respectively. The constants f and g (0 <f,g< 1), are chosen in a way that enables the planet to accrete the accumulated planetesimals. To fulfill this criterion, the oscillation amplitude should be set large enough for the planet to repeatedly go through the region where planetesimals are accumulated.

Fig. 8 Time evolution of the location of the pressure maximum planetesimal accumulations for the different planetesimal sizes at the outer edge of the dead zone for a flat disk of Md = 0.1 M⊙ using αback = 10-3 and αdz = 10-5.

Fig. 9 Time evolution of the core’s mass in the in situ formation of a planet located at ~ 3.2 au for different planetesimal sizes. Simulations correspond to a flat disk of Md = 0.1 M⊙ considering a disk with (black lines) and without a dead zone (red lines). When a dead zone is considered we use αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au, while if the dead zone is not taken account we use α = 10-3 along the disk. Simulations end when the planet achieves the critical mass or after 5 Myr of viscous disk evolution.

Fig. 10 Time evolution of the core’s mass in the in situ formation of a planet located near the outer edge of the dead zone. The planet is located at 16 au, 16.5 au, 17 au, and 17.5 au for planetesimals of 0.1 km, 1 km, 10 km, and 100 km of radius, respectively. Simulations correspond to a flat disk of Md = 0.1 M⊙ considering a disk with (black lines) and without a dead zone (red lines). When a dead zone is considered we use αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au, while if the dead zone is not taken into account we use α = 10-3 along the disk. Simulations end when the planet achieves the critical mass or after 5 Myr of viscous disk evolution.

Fig. 11 Time evolution of the planet’s semi-major axis (left y-axis) and planet’s core mass (right y-axis) for different initial planet locations: 2 au, 5 au, 10 au, 17.5 au, and 22 au. Black lines represent the position of the zero torque. Simulations correspond to a flat disk with a mass of Md = 0.1 M⊙ using αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au. Simulations end when the planet achieves critical mass or after 5 Myr of viscous disk evolution.

Fig. 12 Top: time evolution of the planet’s semi-major axis (left y-axis) and planet’s core mass (right y-axis) being initially located at 5 au. The red dashed line represents the location of the pressure maximum, and the black dashed line represents the zero torque position. Bottom: time evolution (color palette) of the planetesimal surface density radial profiles. The left panel represents the case where the planet is trapped at the zero torque location, while the right panel represents the case where the semi-major axis of the planet oscillates around the zero torque location after being trapped. Simulations correspond to a flat disk with mass Md = 0.1 M⊙ using αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au, Rout - dz = 20 au, and planetesimals of 1 km radius. Simulations end when the planet achieves critical mass.

Figure 12 (top-left panel) shows the time evolution of the semi-major axis of the planet’s core and mass, started initially from 5 au, considering a population of 1 km sized planetesimals. The core grows and migrates quickly until reaching the migration trap. Having been trapped there, the protoplanet continues growing until achieving critical mass. However, as we can see, the pressure maximum (red dashed line) is not at the same location as the migration trap. Analyzing the time evolution of the planetesimal surface density profiles for this simulation (bottom-left panel) one can see that while the planet is able to accrete a large quantity of the accumulated planetesimals, it is not able to accrete all the available mass and, a significant number of planetesimals remain in the disk. When applying the random oscillation of the growing protoplanet around the zero torque location, the formation time of the critical mass core is reduced by more than 50% with respect to the previous case. One can see that the random oscillation allows the growing planet to accrete all the mass accumulated at pressure maximum (bottom-right panel).

This phenomenon is more significant for less massive disks. Figure 13 shows the results of the same simulations as before, but considering a disk with mass M_d = 0.03 M_⊙ (top panel) and with M_d = 0.05 M_⊙ (middle panel). Left and right panels show the behavior of the semi-major axis and the mass of the growing protoplanet as functions of time with and without the random oscillation around the zero torque location, respectively. For a disk of M_d = 0.03 M_⊙, the planet does not achieve the critical mass in less than 5 Myr unless the random oscillation around the zero torque location is considered. For a disk of M_d = 0.05 M_⊙, the formation time is practically halved in the case of random oscillation.

Finally, the bottom panel of Fig. 13 shows again the same results but now for a disk of M_d = 0.05 M_⊙ considering α_back = 10^-2 and α_dz = 10^-4. In this case, the formation times are ~ 1 Myr greater than the case of α_back = 10^-3 and α_dz = 10^-5 (middle panel). This is due to the fact that the accumulation of planetesimals around the pressure maximum location is more effective for smaller values of the alpha parameter.

Fig. 13 Same as top panel of Fig. 12 for a disk with Md = 0.03 M⊙ (top panel) and a disk with Md = 0.05 M⊙ (middle panel). The bottom panel represents the case of a disk with mass Md = 0.05 M⊙ but considering that αback = 10-2, αdz = 10-4.

Regarding the outer pressure maximum, the planetesimal accumulation location coincides with the pressure maximum location only for planetesimals with a 100 m radius. For the other planetesimal sizes, the radial drift velocities are lower than the shifting velocity of the outer pressure maximum, so the initial planetesimal accumulations remain far away from the pressure maximum and zero torque location (Fig. 5 bottom panel). Thus, except for planetesimals of 100 m radius we need a large oscillation amplitude in order to accrete the planetesimal accumulation. For this reason, we only analyze the case of planetesimals of 100 m radius. Figure 14 shows both the time evolution of the planet’s semi-major axis (left y-axis) and the mass of the planetary core (right y-axis). When we do not consider a dead zone in the disk, the planet achieves the critical mass at the inner part of the disk due to type I migration for both disk masses (0.1 M_⊙ and 0.05 M_⊙). When a dead zone is considered in the disk, but the planet does not oscillate around the zero torque location, and, for the massive disk, the planet is able to achieve critical mass before the disappearance of the zero torque location. However, for the disk of 0.05 M_⊙, the planet is not able to achieve critical mass at the outer region until the zero torque disappears, and it quickly migrates inward reaching the inner migration trap where it achieves critical mass. When the noisy oscillation of the planet’s semi-major axis is considered, the planet achieves critical mass before the zero torque disappears. Thus, the planet becomes a giant planet at wide orbit, even for a moderate-mass disk. However, for the disk of 0.1 M_⊙, the formation history of the planet is practically the same independently of whether or not the noisy oscillation of its semi-major axis has been considered around the zero torque location.

Fig. 14 Time evolution of the planet’s semi-major axis (left y-axis ) and planet’s core mass (right y-axis) being initially located at 17.5 au for a disk of 0.1 M⊙ (top panel) and a disk of 0.05 M⊙ (bottom panel). Red dashed line represents the location of the pressure maximum, and black dashed line represents the zero torque location. The left panel represents the case where the planet is trapped at the zero torque location (blue lines) and the case where a dead zone is not present in the disk (pink lines), while the right panel represents the case where the semi-major axis of the planet randomly oscillates around the zero torque location after becoming trapped. Simulations correspond to a flat disk with αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au for the disk with the dead zone and α = 10-3 for the disk without the dead zone, and using planetesimals of 100 m radius. Simulations end when the planet achieves critical mass.

4.3. Pebble accretion at pressure maxima

In this section we analyze the formation of massive cores by pebble accretion at the pressure maxima of the disk. To do so, we consider a population of pebbles of 1 cm in size incorporating pebble accretion in our model of planet formation. As in Guilera (2016), we adopt the pebble accretion rates given by Lambrechts et al. (2014), $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{M}}_{\mathit{C}}}{\mathrm{d}\mathit{t}}{|}_{\mathrm{pebbles}}\mathrm{=}\{\begin{array}{c}\\ \\ \\ \end{array}& & \end{array}$(21)where St = t_stopΩ_k is the Stokes number, t_stop being the stopping time, which depends on the drag regime (Rafikov 2004; Chambers 2008), and where we introduce the factor β = min(1,R_H/H_p) in order to take into account a reduction in the pebble accretion rates if the scale height of the pebbles, H_p, becomes greater than the Hill radius of the planet. The scale height of the solids at a given distance R from the central star is given by (Youdin & Lithwick 2007) $\begin{array}{ccc}{\mathit{H}}_{\mathrm{p}}\mathrm{=}{\mathit{H}}_{\mathrm{g}}\sqrt{\frac{\mathit{\alpha }}{\mathit{\alpha }\mathrm{+}\mathrm{St}}}\mathrm{·}& & \end{array}$(22)First, we calculate the evolution of a disk without any planet immersed in it. Figure 15 shows the time evolution of the surface density of the pebble population for a low-mass flat disk of 0.03 M_⊙. The top panel represents the case where the dead zone is not considered. Initially, the migration of the pebbles from the outer part of the disk increases the pebble surface density in the planet formation region (R ≲ 10 au). However, as time advances, all the solid material is deposited in the pressure maximum generated by the inner boundary condition (see Fig. 2). The bottom panel shows the case where the dead zone is included in the model. We can see that the solid material is quickly accumulated at the pressure maxima of the disk. The pebbles between the star and the inner edge of the dead zone are concentrated at the pressure maximum generated by the inner boundary condition. The pebbles in the dead zone quickly migrate to the inner pressure maximum generated by the snow line (at ~ 3 au), increasing the surface density of solids there by approximately three orders of magnitude. The outer pressure maximum being developed at the outer edge of the dead zone concentrates all the inward drifting solid material of the outer disk. When this pressure maximum disappears, solids are drifted again toward the star until reaching the inner pressure maximum. In Fig. 16, we show that particles of 1 cm behave as pebbles (with St ≤ 1) for distances R ≲ 20 au, including the region of the disk where pressure maxima are developed.

Finally, in Fig. 17 we show the time evolution of the semi-major axis and the mass of the giant planet’s core for different starting positions. When the starting location of the planet is at 2 au, the core does not grow or migrate. This is due to the fact that pebbles at approximately 2 au are drifted very quickly to the inner radius of the disk meaning that all the outer solid material moving inward is collected at the pressure maximum at the inner edge of the dead zone depleting the feeding zone of the core very quickly. The situation for the other cases is very different. When the initial position of the core is at 5 au and 10 au, the core grows rapidly reaching a few Earth masses due to the high pebble accretion rates. Thus the growing core also migrates quickly achieving the inner zero torque location in ~ 10⁵ yr. When reaching the zero torque location, a rapid growth of the core begins due to accretion of the pebbles accumulated at the pressure maximum being close to the zero torque’s place. A similar result is obtained when the core is initially placed at 17.5 au. In this case, the planet slightly migrates outward until reaching the outer zero torque location, then a rapid growth of the core begins due to accretion of the pebbles accumulated at the outer pressure maximum. Finally, when the initial location of the core is at 22 au it reaches the cross-over mass in less than 10⁵ yr before it achieves the outer zero torque location. This is due to the high pebble accretion rates and the inward migration of all the solid material of the outer part of the disk.

Fig. 15 Time evolution of the pebble surface density radial profiles. The top panel corresponds to a disk without a dead zone while the bottom panel corresponds to the case of a disk with a dead zone. Simulations correspond to a low-mass flat disk of Md = 0.03 M⊙ and using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5, when the dead zone is considered.

Fig. 16 Time evolution of the Stokes number as a function of the distance to the central star for particles of 1 cm radius. Simulation corresponds to a low-mass flat disk of Md = 0.03 M⊙ with a dead zone using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5.

Fig. 17 Time evolution of the planet’s semi-major axis (solid lines) and planet’s core mass (dashed lines) for different initial locations. Black dashed line represents the zero torque location. Simulations correspond to a low-mass flat disk of Md = 0.03 M⊙ with a dead zone using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5. Simulations stop after 5 Myr or when the planet reaches the cross-over mass.

5. Summary and conclusions

In the present study we investigate the efficiency of giant planet formation near the density/pressure maxima that may exist in protoplanetary disks. We assumed two locations in a disk, where density/pressure maxima can be developed, these are the water snowline and the outer edge of the dead zone. The attractive features of these locations are twofold: at density maxima the torque responsible for the type 1 migration of planetary cores can vanish, thus they can act as migration traps. Moreover, near to a density maximum, a pressure maximum can also be developed. Approaching a pressure maximum, the gas orbital velocity tends to the circular Keplerian one implying that the gas drag on the solid particles becomes zero. Thus a pressure maximum is a preferential place for dust coagulation and concentration of planetesimals. Although the pressure maximum and migration trap do not exactly coincide, their proximity makes these places favorable for the oligarchic growth of embryos, which can lead to rapid formation of giant planets in a considerably shorter time than the disk’s lifetime.

Our physical model incorporates (i) the time evolution of a viscosity-driven axissymmetric gaseous disk; (ii) the radial drift of planetesimals due to the aerodynamic drag force arising from the ambient sub/super-Keplerian gas disk; (iii) type 1 migration of a growing planetary embryo; and finally (iv) its hydrostatic growth by planetesimal and gas accretion until the gaseous runaway phase, which happens when the mass of the gaseous envelope becomes larger than the mass of the solid core. The differential equations describing the above processes have been numerically solved simultaneously and self-consistently. The two pressure maxima were developed by applying a strong reduction of the α viscosity between the radial distances of the water snow line and the outer edge of the dead zone, mimicking the effect of the almost accretionally inactive dead zone. In some of our simulations, a random oscillation of the growing embryo’s semi-major axis around the density maximum has also been incorporated. It is noteworthy to mention that such oscillations appear in HD/MHD simulations due to the gravitational interaction between the planetary core and the large-scale vortex, the latter being the 2D/3D manifestation of a density/pressure maximum.

First we studied the development of the density/pressure maxima in flared and flat disk models. We find that at the onset of simulations, the density and pressure maxima develop in both cases. In the case of a flared disk, both the zero torque location and the pressure maximum exist throughout the numerical simulations (being 5 × 10⁶ yr) at the water snow line. It is also clearly seen, however, that these locations do not coincide, but the distance between them remains almost constant. Regarding the outer edge of the dead zone, the zero torque location vanishes in a million years timescale, while the pressure maximum lasts for a longer time, but their positions deviate relatively quickly. In the flat disk models, the inner pressure maximum and the migration trap also survive the whole simulation time. The pressure maximum and the zero torque location also disappear, but in a longer timescale than in the case of a flared disk. Their positions also deviate much slower than in the case of a flared disk. Thus as a conclusion, we can state that both in flared and flat disks, at the inner density/pressure maximum, the formation of a giant planet is supported during the disk’s lifetime. Another conclusion is that in flat disks, giant planet formation might be more likely at the outer edge of the dead zone than in flared disks, since in flat disks the migration trap and the pressure maximum exists for a longer time.

Simulations with planetesimals of different sizes have revealed that their accumulation at the pressure maxima is size dependent. We have performed simulations with planetesimal sizes of 0.1, 1, 10, and 100 km. Since the 0.1 km sized planetesimals are the most sensitive to the drag force, their accumulation is the most effective, being followed by the other planetesimal sizes in increasing order. Since the inner pressure maximum moves slightly outwards, all of the inward drifting planetesimals are trapped there, irrespective of their sizes. The density/pressure maximum at the outer edge of the dead zone moves slightly inwards, thus only those planetesimals can be trapped there which are drifted faster than the inward shift of the outer pressure maximum. We have found that practically only the planetesimals with sizes of 0.1 km can be trapped efficiently in the outer pressure maximum.

Having investigated planetesimal accumulation at the pressure maxima, we have studied the in situ formation of giant planets meaning that we have not considered planet-disk interaction and placed embryos to fixed locations, exactly where planetesimal accumulation occurs. We have found in all of our simulations (except the outer dead zone edge and 100 km sized planetesimals) that the formation time of a giant planet is much shorter at pressure maxima than in disk models without dead zone. (We recall that by formation time we mean the time needed until the solid core’s mass is equal to the gaseous envelope’s mass). In situ giant planet formation is however not physical, since due to the planet-disk interaction, the growing core either migrates or is trapped at a zero torque location, which for the outer edge of the dead zone deviates from the accumulation of planetesimals. Therefore as a next step, we have taken into account type 1 migration of the growing core. According to our results, except of the case of 100 km size planetesimals, the formation of giant planets at the water snow line happens well before the disk’s lifetime. Giant planet formation at the inner density/pressure maximum is even more effective, if a random oscillation of the semi-major axis of the planetary core is assumed. At the outer edge of the dead zone, the giant planet formation is only effective for 0.1 km sized planetesimals. If the core’s semi-major axis randomly oscillates around the outer edge of the dead zone, the formation time is also shorter than in the case without oscillation.

Finally, we also studied the formation of giant planets at the pressure maxima of the disk by the accretion of pebbles of 1 cm sizes. As for the case of planetesimals, pressure maxima act as locations of accumulation of solid material significantly increasing the pebble surface density. We found that due to the high pebble accretion rates, and the accumulated pebbles at the pressure maxima locations, massive cores are formed in a timescale of 10⁵ yr in the inner and outer edges of the dead zone, even for low-mass disks, meaning that giant planet formation via pebble accretion at a pressure maximum is the fastest and most efficient formation scenario in the core accretion model.

In our study we find that the inner pressure maximum is always a favorable place for giant planet formation for a wide range of physical disk parameters, meaning that the formation time is much shorter that the disk’s lifetime. The outer edge of the dead zone can also promote giant planet formation but only for pebbles or smaller, sub-kilometer sized planetesimals, and in disk models when the lifetime of the migration trap is long enough to enable the trapped core to accrete enough material for the onset of the runaway gas accretion.

It is important to note that during our investigations the positions of the inner and outer pressure maxima are kept fixed in time. We are aware of the fact that this might be a simplification of a more complex problem not addressed in this work. As we mentioned in Sect. 2, we locate the inner pressure maximum to the position of the water snow line that develops due to the condensation of water resulting in a sudden increase of the solid-to-gas ratio. The condensation of water reduces the number of free electrons and thus increases the resistivity of the gas suppressing the MRI driven turbulence. Water condensation happens when the disk’s midplane temperature drops below 170 K. During the disk’s lifetime its temperature profile evolves as a function of the distance from the star, therefore the position where water condensation takes place evolves also (Garaud & Lin 2007; Oka et al. 2011). On the other hand, the time-evolution of snow line position might be a more complex issue than simply monitoring the place where T(R_drop) ~ 170 K. For example, Ciesla & Cuzzi (2006) also took into account the sublimation/condensation of the appropriate amount of ice or vapor to maintain the equilibrium of the vapor pressure considering the radial drift of ice-rich dust and planetesimals. More recently, Morbidelli et al. (2016) found that at a certain epoch, t_crit of disk evolution, the radial inward velocity of gas is larger than the speed at which the condensation front moves inward. Thus at larger times t>t_crit, the radius of the temperature drop R_drop moves in water-poor gas, and no water condensation can take place. In that case, the pressure maximum might be developed at the interface of the water-poor and water-rich gas.

We note, however, that neither of the above scenarios are directly applicable to our case, since the pressure maximum, that we assume to develop, traps icy dust grains, and also reduces the gas radial inward velocity. Therefore we intend to investigate the above complex issue in a separate study. Regarding the outer edge of the dead zone, its location depends on the penetrating depth of the X-rays and cosmic rays. Therefore, when the gas surface density becomes less than some critical value, the disk can be ionized in its full vertical extent. Thus, it is expected that the outer edge of the dead zone moves in time. Matsumura et al. (2009) showed that the outer edge of the dead zone can reach a few au in a few Myr.

Nevertheless the motion of the snow line and the outer edge of the dead zone may not influence the growth of the embryo itself, since the planetesimals and pebbles will be trapped in the pressure maxima developed at those locations. Being trapped, the growing cores are moving presumably together with the snow line/dead zone’s outer edge as they are changing their place. Therefore an interesting issue that merits further investigation is the distance of the snow line/dead zone’s outer edge to the star when the planet traps are no longer able to hold in place the growing cores. We assume that the moving snow line and the outer dead zone edge would certainly affect the final formation place of the giant planet.

The efficiency of giant planet formation at the pressure maxima in more realistic disk models will be the subject of a forthcoming study.

Acknowledgments

We thank W. Lyra, the referee of this work, for his suggestions and comments helping us to improve the work. This work was initiated during the visit of Zs. Sándor at the Faculty of Astronomical and Geophysical Sciences of the National University of La Plata, therefore Zs. Sándor thanks Prof. P. Cincotta and Dr. N. Maffione for their kind invitation, arrangements, and hospitality. O. M. Guilera is supported by grants from the National Scientific and Technical Research Council and National University of La Plata, Argentina. Zs. Sándor also thanks the support of the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.

References

All Tables

All Figures

	Fig. 1 Alpha-viscosity parameter as a function of the radial coordinate for a flat and a flared disk, with αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au. We consider a disk of 0.05 M⊙, and for both disk cin - dz = cout - dz = 1. For a flat disk, Hg(Rin - dz) = 0.135 au and Hg(Rout - dz) = 1 au, while for a flared disk Hg(Rin - dz) = 0.165 au and Hg(Rout - dz) = 2.013 au.
In the text

	Fig. 2 Time evolution of the gas surface density radial profiles. The simulation corresponds to a disk of Md = 0.1 M⊙, using α = 10-3. The simulation is stopped after 5 Myr of viscous evolution.
In the text

	Fig. 3 Time evolution of the planetesimal surface density radial profiles for planetesimal populations of different radii. Simulations correspond to a disk of Md = 0.1 M⊙, using α = 10-3. Simulations are stopped after 5 Myr of viscous evolution.
In the text

	Fig. 4 Time evolution of the gas surface density radial profiles assuming the existence of a dead zone for a flared disk (top panel) and a flat disk (bottom panel). Simulations corresponding to a disk of Md = 0.1 M⊙ using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5. Simulations are stopped after 5 Myr of viscous disk evolution.
In the text

	Fig. 5 Time evolution of the zero torque locations (solid lines) and the locations of pressure maxima (dashed lines) at the inner egde (top panel) and outer edge (bottom panel) of the dead zone for a flared disk (red lines) and a flat disk (black lines). Simulations correspond to a disk of Md = 0.1 M⊙ using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5.
In the text

	Fig. 6 Time evolution of the zero torque locations (solid lines) and pressure maximum locations (dashed lines) at the inner and outer edge of the dead zone for a flat disk using αback = 10-3 and αdz = 10-5 (black lines), and αback = 10-2 and αdz = 10-4 (blue lines).
In the text

	Fig. 7 Time evolution of the planetesimal surface density radial profiles for planetesimal populations of different radii. Simulations correspond to a flat disk of Md = 0.1 M⊙, using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5. Simulations are stopped after 5 Myr of viscous evolution.
In the text

	Fig. 8 Time evolution of the location of the pressure maximum planetesimal accumulations for the different planetesimal sizes at the outer edge of the dead zone for a flat disk of Md = 0.1 M⊙ using αback = 10-3 and αdz = 10-5.
In the text

Fig. 9 Time evolution of the core’s mass in the in situ formation of a planet located at ~ 3.2 au for different planetesimal sizes. Simulations correspond to a flat disk of Md = 0.1 M⊙ considering a disk with (black lines) and without a dead zone (red lines). When a dead zone is considered we use αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au, while if the dead zone is not taken account we use α = 10-3 along the disk. Simulations end when the planet achieves the critical mass or after 5 Myr of viscous disk evolution.

In the text

Fig. 10 Time evolution of the core’s mass in the in situ formation of a planet located near the outer edge of the dead zone. The planet is located at 16 au, 16.5 au, 17 au, and 17.5 au for planetesimals of 0.1 km, 1 km, 10 km, and 100 km of radius, respectively. Simulations correspond to a flat disk of Md = 0.1 M⊙ considering a disk with (black lines) and without a dead zone (red lines). When a dead zone is considered we use αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au, while if the dead zone is not taken into account we use α = 10-3 along the disk. Simulations end when the planet achieves the critical mass or after 5 Myr of viscous disk evolution.

In the text

Fig. 11 Time evolution of the planet’s semi-major axis (left y-axis) and planet’s core mass (right y-axis) for different initial planet locations: 2 au, 5 au, 10 au, 17.5 au, and 22 au. Black lines represent the position of the zero torque. Simulations correspond to a flat disk with a mass of Md = 0.1 M⊙ using αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au. Simulations end when the planet achieves critical mass or after 5 Myr of viscous disk evolution.

In the text

Fig. 12 Top: time evolution of the planet’s semi-major axis (left y-axis) and planet’s core mass (right y-axis) being initially located at 5 au. The red dashed line represents the location of the pressure maximum, and the black dashed line represents the zero torque position. Bottom: time evolution (color palette) of the planetesimal surface density radial profiles. The left panel represents the case where the planet is trapped at the zero torque location, while the right panel represents the case where the semi-major axis of the planet oscillates around the zero torque location after being trapped. Simulations correspond to a flat disk with mass Md = 0.1 M⊙ using αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au, Rout - dz = 20 au, and planetesimals of 1 km radius. Simulations end when the planet achieves critical mass.

In the text

	Fig. 13 Same as top panel of Fig. 12 for a disk with Md = 0.03 M⊙ (top panel) and a disk with Md = 0.05 M⊙ (middle panel). The bottom panel represents the case of a disk with mass Md = 0.05 M⊙ but considering that αback = 10-2, αdz = 10-4.
In the text

Fig. 14 Time evolution of the planet’s semi-major axis (left y-axis ) and planet’s core mass (right y-axis) being initially located at 17.5 au for a disk of 0.1 M⊙ (top panel) and a disk of 0.05 M⊙ (bottom panel). Red dashed line represents the location of the pressure maximum, and black dashed line represents the zero torque location. The left panel represents the case where the planet is trapped at the zero torque location (blue lines) and the case where a dead zone is not present in the disk (pink lines), while the right panel represents the case where the semi-major axis of the planet randomly oscillates around the zero torque location after becoming trapped. Simulations correspond to a flat disk with αback = 10-3, αdz = 10-5, Rin - dz = 2.7 au and Rout - dz = 20 au for the disk with the dead zone and α = 10-3 for the disk without the dead zone, and using planetesimals of 100 m radius. Simulations end when the planet achieves critical mass.

In the text

	Fig. 15 Time evolution of the pebble surface density radial profiles. The top panel corresponds to a disk without a dead zone while the bottom panel corresponds to the case of a disk with a dead zone. Simulations correspond to a low-mass flat disk of Md = 0.03 M⊙ and using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5, when the dead zone is considered.
In the text

	Fig. 16 Time evolution of the Stokes number as a function of the distance to the central star for particles of 1 cm radius. Simulation corresponds to a low-mass flat disk of Md = 0.03 M⊙ with a dead zone using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5.
In the text

Fig. 17 Time evolution of the planet’s semi-major axis (solid lines) and planet’s core mass (dashed lines) for different initial locations. Black dashed line represents the zero torque location. Simulations correspond to a low-mass flat disk of Md = 0.03 M⊙ with a dead zone using Rin - dz = 2.7 au, Rout - dz = 20 au, αback = 10-3, and αdz = 10-5. Simulations stop after 5 Myr or when the planet reaches the cross-over mass.

In the text