© ESO, 2013

1. Introduction

The standard scenario for the formation of planets in protoplanetary discs generally involves the following steps: i) coagulation and settling of dust in the disc midplane, followed by growth of km-sized planetesimals; ii) runaway growth of planetesimals (Greenberg et al. 1978; Wetherill & Stewart 1989) into ~ 10^-5  M_⊕ embryos; iii) oligarchic growth of these embryos (Kokubo & Ida 1998, 2000; Leinhardt & Richardson 2005) into planetary cores. Planetary cores forming oligarchically beyond the snow-line are expected to have masses ~10  M_⊕ (Thommes et al. 2003) and consequently are able to accrete a gaseous envelope to become a giant planet (Pollack et al. 1996) within the lifetime of protoplanetary discs. This however requires a relatively massive protoplanetary disc, equivalent to 10 times the mass of the minimum-mass solar nebula (MMSN; Hayashi 1981). Moreover, recent N-body simulations including the effects of fragmentation (Levison et al. 2010) indicate only a modest further growth of embryos once these have reached a mass of ~ 2  M_⊕. This occurs because as the embryos grow, they tend to scatter planetesimals outside of their feeding zone rather than accreting them. The action of gas drag then makes the orbits of these planetesimals quasi-circular, which prevents close encounters with embryos. These results emphasize the difficulty of forming giant planet cores within a few Myr in the context of the oligarchic growth scenario.

Recently, an alternative model for forming giant planet cores has been proposed by Lambrechts & Johansen (2012) and in which embryos grow by accretion of cm-sized pebbles. Compared with the classical scenario, the growth timescale to reach a critical core mass of 10  M_⊕ is typically reduced by a factor of 30–1000 at 5 AU in this model and naturally accounts for the preferential prograde spin of large asteroids (Johansen & Lacerda 2010). Using hydrodynamical simulations, Morbidelli & Nesvorny (2012) recently examined this process in more details and found that for a MMSN model, the mass doubling time of a 1  M_⊕ embryo accreting 20-cm pebbles at 1 AU is only ~5500 yr. They confirmed that the model of Lambrechts & Johansen (2012) is promising for forming embryos of a few Earth masses.

Subsequent giant impacts between embryos may produce giant planet cores. Giant impacts have been invoked during the late stages of terrestrial planet formation (Wetherill 1985) and may explain the origin of Earth’s Moon (e.g., Canup & Asphaug 2001). The fact that Uranus’ equatorial satellites are on prograde orbits despite the planet’s retrograde rotation can be explained by multiple giant impacts of roughly Earth-mass embryos during Uranus’ accretion (Morbidelli et al. 2012). Although Neptune’s rotation is prograde, its modest obliquity of 28.5° also appears to require a giant impact (Morbidelli et al. 2012). Finally, it is possible that giant impacts can stimulate runaway gas accretion (Broeg & Benz 2012).

Whether they form by accreting planetesimals or pebbles, embryos must undergo Type I migration due to their interaction with the gaseous disc (Ward 1997; Tanaka et al. 2002). The disc torque exerted on a low-mass planet and causing Type I migration consists of two components: i) the differential Lindblad torque due to the angular momentum exchange between the planet and the spiral density waves it generates inside the disc, which is generally negative and therefore responsible for inward migration; ii) the corotation torque exerted by the material located in the coorbital region of the planet, which scales with both the vortensity (i.e. the ratio between the vertical component of the disc vorticity and the disc surface density; Goldreich & Tremaine 1979) and the entropy gradients inside the disc (Baruteau & Masset 2008; Paardekooper & Papaloizou 2008). In particular, positive surface density gradients or negative entropy gradients give rise to a positive corotation torque which may eventually counteract the effect of the differential Lindblad torque and subsequently lead to outward migration. This arises provided that an amount of viscous/thermal diffusion is present inside the disc so that the corotation torque remains unsaturated, and that diffusion processes operate in such a way that the amplitude of the corotation torque is close to its fully unsaturated value. For non-isothermal, viscously heated protoplanetary discs, the torque experienced by a protoplanet is typically positive in the inner, optically thicks regions whereas the outer, optically thin regions give rise to a negative torque (Kretke & Lin 2012; Bitsch et al. 2013). In that case, one can expect protoplanetary discs to present locations where the Type I torque cancels and where protoplanets may converge. These are referred to as zero-migration lines or convergences zones and are generally considered as ideal sites for the growth of planetary embryos (Lyra et al. 2010; Hasegawa & Pudritz 2011).

A significant body of work has recently investigated the role of these zero-migration lines on the formation of giant planet cores. Sandor et al. (2011) studied this process in isothermal discs where zero-torque radii are located at dead-zone boundaries and found that 10  M_⊕ bodies can be formed in less than ~ 5 × 10⁵ yr through collisions of smaller embryos. The case of non-isothermal discs was investigated by Hellary & Nelson (2012) who performed N-body simulations of planetary growth in radiatively-inefficient protoplanetary discs. They showed that in non-isothermal discs, the convergent migration induced by corotation torques can indeed enhance the growth rate of planetary embryos. Similar results were obtained by Horn et al. (2012) who confirmed that giant planet cores can form at convergence zones from sub-Earth mass embryos in 2–3 Myr.

In this paper, we present the results of hydrodynamical simulations of the interaction of multiple protoplanets in non-isothermal disc models. Our simulations begin with 3  M_⊕ bodies with positions initially distributed around an opacity transition located just inside the snow-line. This opacity transition corresponds to a zero-torque radius for planets of masses 1 − 10  M_⊕, and inside (resp. outside) which Type I migration proceeds outward (resp. inward). The main aim of this work are: i) to determine the typical evolution outcome of a swarm of mutiple Earth-mass protoplanets which convergently migrate toward a zero-migration line; ii) to investigate whether giant planet cores can be formed at convergence zones from giant impacts between bodies of a few Earth masses. With respect to previous studies based on N-body simulations and which employ prescribed forces for migration, hydrodynamical simulations allow a self-consistent treatment of the interactions between the embryos and the gas disc. We performed simulations varying the initial number of embryos and tested the impact of including stochastic forces on the planets to mimic the effects resulting from disc turbulence. For laminar runs involving a modest initial number of embryos (N ≤ 5 − 7), we find that the system enters in a long resonant chain which remains stable for ~10⁴ yr whereas growth of embryos through collisions occurs when the initial number of objects is increased. As expected, including a small level or turbulence tends to break resonant configurations, which consequently enhances close encounters between embryos and promotes planetary growth.

In order to study the dynamical evolution on Myr timescales, we also present in this paper the results of N-body runs calibrated on our hydrodynamical simulations. These N-body runs confirm the results of hydrodynamical simulations and show that systems which are formed at convergence zones generally reach a quasi-stationary state with each body in resonance with its neighbours and evolving on a non-migrating orbit. In ~43% of the runs in which the initial mass of embryos is ~3  M_⊕, protoplanets of masses ≳10  M_⊕ were produced, suggesting thereby that giant planet cores can be formed at convergence zones through collisions between bodies of a few Earth masses.

This paper is organized as follows. In Sect. 2, we present the hydrodynamical model. In Sect. 3, we describe how our N-body runs are calibrated from hydrodynamical simulations. The results of hydrodynamical simulations are presented in Sect. 4. In Sect. 5, we discuss the results of the N-body runs. Finally, we draw our conclusions in Sect. 6.

2. The hydrodynamical model

2.1. Numerical method

Simulations were performed with the GENESIS numerical code (De Val-Borro et al. 2006) which solves the equations for the disc on a polar grid. This code employs an advection scheme based on the monotonic transport algorithm (Van Leer 1977) and includes the FARGO algorithm (Masset 2000) to avoid time step limitation due to the Keplerian orbital velocity at the inner edge of the grid. The energy equation that is implemented in the code reads: $\begin{array}{ccc}\frac{\mathit{\partial e}}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{e}{v}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{(}{\mathit{\gamma }}_{\mathrm{ad}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{\nabla }\mathrm{·}{v}\mathrm{)}\mathit{e}\mathrm{+}{\mathit{Q}}_{\mathrm{visc}}^{\mathrm{+}}\mathrm{-}{\mathit{Q}}_{\mathrm{rad}}^{\mathrm{-}}& & \end{array}$(1)where v is the gas velocity, e the thermal energy density and γ_ad the adiabatic index which is set to γ_ad = 1.4. Since we expect the effects resulting from stellar irradiation to be negligible in the inner parts of protoplanetary discs (Bitsch et al. 2013), we include only the contribution from viscous heating in the expression for the heating term ${\mathit{Q}}_{\mathrm{visc}}^{\mathrm{+}}$. In the previous equation, ${\mathit{Q}}_{\mathrm{rad}}^{\mathrm{-}}$ is the radiative cooling term which is given by: $\begin{array}{ccc}{\mathit{Q}}_{\mathrm{rad}}^{\mathrm{-}}\mathrm{=}\mathrm{2}{\mathit{\sigma }}_{\mathrm{B}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}& & \end{array}$(2)where σ_B is the Stephan-Boltzmann constant and T_eff is the effective temperature which is related to the central temperature T by: $\begin{array}{ccc}{\mathit{T}}^{\mathrm{4}}\mathrm{=}{\mathit{\tau }}_{\mathrm{eff}}{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}& & \end{array}$(3)where τ_eff is the effective optical depth given by (Hubeny 1990): $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{eff}}\mathrm{=}\frac{\mathrm{3}\mathit{\tau }}{\mathrm{8}}\mathrm{+}\frac{\sqrt{\mathrm{3}}}{\mathrm{4}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{4}\mathit{\tau }}\mathrm{·}& & \end{array}$(4)In the previous equation τ = κΣ/2 is the optical depth, where Σ is the disc surface density and κ the Rosseland mean opacity which was taken from Bell & Lin (1994).

We employ N_R = 864 radial grid cells uniformly distributed between R_in = 0.6 and R_out = 3.2 and N_φ = 1800 azimuthal grid cells. For a 3  M_⊕ planet, this corresponds to the half-width of the horseshoe region being resolved by ~5 grid cells.

We adopt computational units such that the mass of the central star is M_⋆ = 1M_⊙, the gravitational constant G = 1 and the radius R₀ = 1 in the computational domain corresponds to 1 AU. We use closed boundary conditions at both the inner and outer edges of the computational domain and employ wave-killing zones for R < 0.75 and R > 2.7 to avoid wave reflections at the disc edges.

Evolution of planet orbits is computed using a fifth-order Runge-Kutta integrator (Press et al. 1992). Close encounter between the planets i and j is assumed to occur whenever their mutual distance is less than (3m_i/4πρ_s)^1/3 + (3m_j/4πρ_s)^1/3, where m_i,j is the mass of the planets and ρ_s the mass density which is set to ρ = 1 g cm^-3. In order to guarantee that two planets do not pass trough each other undetected, we set the time step size to Δt = min(Δt_h,Δt_n) where Δt_h is the hydrodynamical time step and Δt_n is a N-body timestep given by (Beaugé & Aarseth 1990): $\begin{array}{ccc}\mathrm{\Delta }{\mathit{t}}_{\mathit{n}}\mathrm{=}\mathrm{0.25}\underset{\mathit{i,j}}{\min }\left|\frac{{\mathit{r}}_{\mathit{ij}}^{\mathrm{2}}}{{{r}}_{\mathit{ij}}\mathrm{·}{{v}}_{\mathit{ij}}}\right|\mathrm{·}& & \end{array}$(5)In the previous expression, r_ij = R_i − R_j is the relative position vector between planets i and j and v_ij = v_i − v_j is their relative velocity.

Although a 2D disc model is adopted, we allow planets to evolve in the direction perpendicular to the disc midplane as well. With respect to coplanar orbits, this will reduce the collision rate between planets, increasing thereby the time during which planets can strongly interact. However, because of the 2D disc model used here, bending waves cannot be launched in the disc, and so there is no disc induced damping of inclination as it would be in a more realistic 3D disc model. To model the inclination damping due to the interaction with the disc we follow Pierens & Nelson (2008) and mimic the effect of bending waves by applying to each planet with mass m_p a vertical force F_z given by: $\begin{array}{ccc}{\mathit{F}}_{\mathit{z}}\mathrm{=}\mathit{\beta }\frac{{\mathit{m}}_{\mathrm{p}}{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\Omega }}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}\mathrm{(}\mathrm{2}{\mathit{A}}_{\mathit{z}}^{\mathit{c}}{\mathit{v}}_{\mathit{z}}\mathrm{+}{\mathit{A}}_{\mathit{z}}^{\mathit{s}}\mathit{z}{\mathrm{\Omega }}^{\mathrm{)}}& & \end{array}$(6)where c_s is the sound speed and where Ω and Σ_p are respectively the Keplerian angular velocity and the disc surface density at the position of the planet. ${\mathit{A}}_{\mathit{z}}^{\mathit{c}}$ and ${\mathit{A}}_{\mathit{z}}^{\mathit{s}}$ are dimensionless coefficients which are set to ${\mathit{A}}_{\mathit{z}}^{\mathit{c}}\mathrm{=}\mathrm{-}\mathrm{1.088}$ and ${\mathit{A}}_{\mathit{z}}^{\mathit{s}}\mathrm{=}\mathrm{-}\mathrm{0.871}$ (Tanaka & Ward 2004), and β is a free parameter which is chosen such that the inclination damping timescale t_i obtained in the simulations is approximately equal to the eccentricity damping timescale t_e. Test simulations show that choosing β = 0.33 give similar values for t_i and t_e.

2.2. Stochastic forces

The origin of turbulence is believed to be related to the magneto-rotational instability (MRI, Balbus & Hawley 1991). Here, turbulence effects are modelled as stochastic forcing using the turbulence model of Laughlin et al. (2004) and further modified by Baruteau & Lin (2010). This model employs a turbulent potential corresponding to the superposition of 50 wave-like modes and given by: $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{turb}}\mathrm{\left(}\mathit{R,\phi ,t}\mathrm{\right)}\mathrm{=}\mathit{\gamma }{\mathit{R}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\sum _{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathrm{50}}{\mathrm{\Lambda }}_{\mathit{k}}\mathrm{\left(}\mathit{R,\phi ,t}\mathrm{\right)}\mathit{,}& & \end{array}$(7)with $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathit{k}}\mathrm{=}{\mathit{\xi }}_{\mathit{k}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{(}\mathit{R}\mathrm{-}{\mathit{R}}_{\mathit{k}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}}}\cos \left({\mathit{m}}_{\mathit{k}}\mathit{\phi }\mathrm{-}{\mathit{\phi }}_{\mathit{k}}\mathrm{-}{\mathrm{\Omega }}_{\mathit{k}}\begin{array}{c}\mathrm{˜}\\ {\mathit{t}}_{\mathit{k}}\end{array}\right)\sin \mathrm{(}\mathit{\pi }\begin{array}{c}\mathrm{˜}\\ {\mathit{t}}_{\mathit{k}}\end{array}\mathit{/}\mathrm{\Delta }{\mathit{t}}_{\mathit{k}}\mathrm{)}\mathit{.}& & \end{array}$(8)In Eq. (8), ξ_k is a dimensionless constant parameter randomly sampled with a Gaussian distribution of unit width. R_k and φ_k are, respectively, the radial and azimuthal initial coordinates of the mode with wavenumber m_k, σ_k = πR_k/4m_k is the radial extent of that mode, and Ω_k denotes the Keplerian angular velocity at R = R_k. Both R_k and φ_k are randomly sampled with a uniform distribution, whereas m_k is randomly sampled with a logarithmic distribution between m_k = 1 and m_k = 150. Each mode of wavenumber m_k starts at time t = t_0,k and terminates when $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{t}}_{\mathit{k}}\end{array}\mathrm{=}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}\mathit{,k}}\mathit{>}\mathrm{\Delta }{\mathit{t}}_{\mathit{k}}$, where Δt_k = 0.2πR_k/m_kc_s denotes the lifetime of mode with wavenumber m_k. Such a value for Δt_k yields an autocorrelation time-scale τ_c ~ 0.5T_orb, where T_orb is the orbital period at R = 1 (Baruteau & Lin 2010).

Following Ogihara et al. (2007), we set Λ_k = 0 if m_k > 6 to save computing time. As noticed by Baruteau & Lin (2010), such an assumption is supported by the fact that a turbulent mode with wavenumber m has an amplitude decreasing as exp( − m²) and a lifetime ∝1/m, so that the contribution to the turbulent potential of a high wavenumber turbulent mode is relatively weak. In Eq. (7), γ denotes the value of the turbulent forcing parameter and is related to the value of the α viscous stress parameter (Shakura & Sunyaev 1973) by the relation (Baruteau & Lin 2010): $\begin{array}{ccc}\mathit{\gamma }\mathrm{~}\mathrm{8.5}\mathrm{\times }{\mathrm{10}}^{-2}{\mathit{h}}_{\mathrm{p}}\sqrt{\mathit{\alpha }}& & \end{array}$(9)where h_p is the aspect ratio. Since it is expected the typical amplitude of the surface density perturbations to scale with γ, the previous expression is consistent with the results of Yang et al. (2009) who found that these turbulent density pertubations scale with $\sqrt{\mathit{\alpha }}$. Although this parametrisation of turbulence does not capture all relevant physical effects like vortices (Fromang & Nelson 2006) or zonal flows (Lyra et al. 2008; Johansen et al. 2009), Baruteau & Lin (2010) have shown that applying the turbulent potential of Eq. (7) to the disc generates density perturbations that have similar statistical properties to those resulting from 3D MHD simulations.

In the context of non-isothermal disc models, an important effect resulting from these turbulent fluctuations is that the initial temperature is progressively alterated by turbulent heating (Pierens et al. 2012). In order for the temperature profile to remain fixed in the course of simulations, we follow Ogihara et al. (2007) and Horn & Lyra (2012) and rather apply the turbulent potential of Eq. (7) on the planets. In that case, the turbulent force F_turb acting on each body is related to Φ_turb by (Ogihara et al. 2007): $\begin{array}{ccc}{\mathit{F}}_{\mathrm{turb}}\mathrm{=}\mathrm{-}\mathit{C}\frac{\mathrm{64}{\mathrm{\Sigma }}_{\mathrm{p}}{\mathit{R}}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}{\mathit{M}}_{\mathit{\star }}}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{turb}}\mathrm{·}& & \end{array}$(10)In the previous equation, C is a constant which is set to a value such that the mean deviation of the turbulent torque distribution coincides with that obtained in the case where the turbulent potential is applied directly to the disc. In order to estimate C, we have measured the torque experienced by a 3  M_⊕ planet for i) an isothermal turbulent simulation in which the turbulent potential of Eq. (7) with γ = 10^-4 is applied to the disc; and ii) a series of laminar isothermal simulations in which the turbulent force of Eq. (10) is applied to the planet, and which differ by the value of C. In the case with C = 17, Fig. 1 shows that the distribution of the specific torque acting on the planet is in good agreement with that obtained in the turbulent simulation.

Fig. 1 Distribution of the specific turbulent torque exerted on a 3  M⊕ planet obtained from a turbulent simulation in which the turbulent potential of Eq. (7) with γ = 10-4 is applied to the disc (solid line), and obtained from a laminar simulation in which the turbulent force of Eq. (10) with C = 17 is applied to the planet (dashed line).

2.3. Initial conditions

The initial disc surface density profile is Σ = Σ_p × (R/R₀)^− σ with σ = 0.5 and Σ_p = 850 g/cm².

The anomalous viscosity in the disc arising from MHD turbulence is modelled using a constant kinematic viscosity ν = 10¹⁴ g/cm², which corresponds to a value for the α viscous stress parameter of α ~ 2 × 10^-3 at 1 AU. The choice of a constant viscosity is justified by the fact that for σ = 0.5, there is no viscous evolution of the disc surface density profile so that the zero-torque line remains approximately fixed in the course of the simulations.

The initial temperature profile is T = T₀(R/R₀)^− β with β = 1 and T₀ is the initial temperature at 1 AU. Under the action of the source terms in Eq. (1), the temperature profile evolves and reaches an equilibrium state once viscous heating equilibrates radiative cooling. The surface density and temperature profiles at steady-state are plotted in Fig. 2. The change in the temperature structure at R ~ 1.7 AU is related to a change in the opacity regime. For R ≤ 1.7 AU, the opacity is dominated by metal grains and varies as κ ∝ T^1/2 whereas for R ≥ 1.7 AU, melting of ice grains causes the opacity to drop with temperature and to vary as κ ∝ T^-7 (Bell & Lin 1994). By equating the viscous heating and radiative cooling terms and assuming an optically thick disc, it is straightforward to show that T ∝ R^− 8/7 inside AU whereas T ∝ R^− 4/11 for AU. This corresponds to an entropy S = p/Σ^γ_ad, where p is the pressure, decreasing as R^-0.92 inside the opacity transition and as R^-0.14 for AU. For this disc model, we notice that the snow-line is located at ~2.7 AU, which is consistent with estimations of the location of the snow-line at the epoch of planetesimal formation, although large excursions from this value are expected due to disc evolution (Lecar 2006; Garaud & Lin 2007).

Fig. 2 Surface density and temperature profiles at steady-state.

In the hydrodynamical simulations presented below, the initial mass of each protoplanet is assumed to be m_p = 3M_⊕. The motivation for choosing equal-mass embryos is based on the fact that this minimizes the rate of convergent migration and therefore the probability of close encounters between embryos when using an isothermal equation of state. As we see in Sect. 4, collisions between equal-mass bodies occur in the radiative case but not in the locally isothermal limit, which clearly demonstrates the role of the zero-migration line in forming more massive objects through collisions. Isothermal simulations with embryos of initially different masses would probably lead to close encouters between embryos (Cresswell & Nelson 2006).

The distribution of semi-major axes a_i is such that about half of initial population of planetary embryos is located on each side of the convergence zone, with an initial orbital separation between two adjacent bodies p and p′ of 4.5  R_mH, where R_mH is the mutual Hill radius: $\begin{array}{ccc}{\mathit{R}}_{\mathrm{mH}}\mathrm{=}\frac{{\mathit{a}}_{\mathrm{p}}\mathrm{+}{\mathit{a}}_{{\mathrm{p}}^{\mathrm{\prime }}}}{\mathrm{2}}{\left(\frac{{\mathit{m}}_{\mathrm{p}}\mathrm{+}{\mathit{m}}_{{\mathrm{p}}^{\mathrm{\prime }}}}{\mathrm{3}{\mathit{M}}_{\mathit{\star }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}& & \end{array}$(11)where a_p,p′ denotes the semimajor-axes of planets p and p′. Although the initial separation between bodies is greater than the critical value of 3.46  R_mH below which rapid instability occurs for two planets on initially circular orbits (Gladman 1993), it is smaller to what is expected from the oligarchic growth scenario (Kokubo & Ida 1998). The adopted value for the initial embryo separation is chosen to make the hydrodynamical simulations computationally tractable, but N-body runs performed with an initial separation of 6  R_mH show consistent results in comparison with those obtained using the fiducial value of 4.5  R_mH (see Sect. 5).

Planetary embryos initially evolve on circular orbits with inclinations randomly sampled according to a Gaussian distribution with mean μ_i = 0° and standard deviation σ_i = 1°.

3. Calibration of N-body simulations

3.1. Prescription for Type I migration

The torque exerted on a protoplanet embedded in a non-isothermal disc can be decomposed into two components: the differential Lindblad torque which results from angular momentum exchange between the planet and the spiral waves it generates inside the disc plus the corotation torque which is due to the torque exerted by the material located in the coorbital region of the planet. A linear analysis reveals that the corotation torque consists of a barotropic part which scales with the vortensity plus an entropy-related part which scales with the entropy gradient. In the absence of any diffusion processes inside the disc, however, vortensity and entropy gradients tend to flatten through phase mixing, which causes the two components of the corotation torque to saturate. Consequently, desaturating the corotation torque requires that some amount of viscous and thermal diffusions are operating inside the disc. In that case, the amplitude of the corotation torque depends on the ratio between the diffusion time-scales and the horseshoe libration time-scale. Its optimal value, also referred to as the fully unsaturated corotation torque, is obtained when the diffusion time-scales are approximately equal to half the horseshoe libration time (e.g. Baruteau & Masset 2013). In the limit where the diffusion time-scales become shorter than the U-turn time-scale, the corotation torque decreases and approaches the value predicted by linear theory. Therefore, the corotation torque can be considered as a linear combination of the fully unsaturated corotation torque and the linear corotation torque, with coefficients depending on the ratio between the diffusion time-scales and the horseshoe libration time-scale. Torque formulae as a function of viscosity and thermal diffusivity were proposed by Paardekooper et al. (2011).

In our N-body runs, Type I migration is modelled by an extra-force a_m acting on each body and defined by: $\begin{array}{ccc}{{a}}_{\mathrm{m}}\mathrm{=}\mathrm{-}\frac{{{v}}_{\mathrm{p}}}{{\mathit{t}}_{\mathrm{m}}}\mathrm{·}& & \end{array}$(12)In the previous expression, v_p is the planet velocity and t_m   =   J/2Γ is the migration timescale, where J is the specific planet angular momentum and Γ the specific disc torque. To estimate Γ, we use the analytical prescription of Paardekooper et al. (2011) but multiplied by a factor of 0.7. As we will see shortly, very good agreement with hydrodynamical simulations is obtained in that case. We notice that this is consistent with the results of Cresswell & Nelson (2006) who found that for isothermal disc models, analytical torque formulae given in Tanaka et al. (2002) predict migration times that are faster than those observed in the simulations by about a factor of three.

Fig. 3 Upper panel: torque exerted on a 3  M⊕ planet as a function of its orbital distance. The solid line corresponds to the analytical prescription of Paardekooper et al. (2011). Triangles and diamonds are the results from hydrodynamical simulations with/without the effects due to the waves generated by the other planets included. Lower panel: same but for mp   =   6  M⊕.

For a purely active disc, radiative diffusion and viscous time-scales are expected to be equal (Bitsch & Kley 2011) so when computing the saturation parameters in the analytical formulae, we set ν = χ, where χ is the thermal diffusion coefficient. In the case where thermal diffusion is only due to radiative effects, χ is given by (e.g. Bitsch & Kley 2011): $\begin{array}{ccc}\mathit{\chi }\mathrm{=}\frac{\mathrm{16}{\mathit{\gamma }}_{\mathrm{ad}}\mathrm{(}{\mathit{\gamma }}_{\mathrm{ad}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{\sigma }{\mathit{T}}^{\mathrm{4}}}{\mathrm{3}\mathit{\kappa }{\mathit{\rho }}^{\mathrm{2}}{\mathit{H}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}& & \end{array}$(13)where H is the disc scale height and ρ = Σ/2H the gas density. In Fig. 3, we compare for the disc model described in Sect. 2.3 and for m_p = 3,6  M_⊕ the analytical torque of Paardekooper et al. (2011), which we multiplied by a factor of 0.7, with the numerical torque obtained using GENESIS. Clearly, very good agreement is obtained between the analytical prediction and the torques derived from numerical simulations. Inside the zero-torque radius located at R ~ 1.7 AU, the torque is positive due to the strong (negative) entropy gradient there (S ∝ R^-0.92, see Sect. 2.3) whereas outside the zero-migration line, the entropy gradient is weaker (S ∝ R^-0.14) so that the (positive) entropy-related corotation torque is not strong enough to counterbalance the (negative) differential Lindblad torque. We note that in Fig. 3, the zero-migration line at R ~ 1.7 AU exists provided that the corotation torque remains unsaturated. This arises when the diffusion timescale across the horseshoe region ${\mathit{t}}_{\mathrm{dif}}\mathrm{=}{\mathit{x}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathit{\chi }$ is shorter than the libration timescale t_lib = 8πa_p/3Ω_px_s but longer than the U-turn timescale t_U − turn ~ h_pt_lib, where x_s is the half-width of the horseshoe region which is given by (Paardekooper et al. 2010): $\begin{array}{ccc}{\mathit{x}}_{\mathrm{s}}\mathrm{=}\frac{\mathrm{1.1}{\mathit{a}}_{\mathrm{p}}}{{\mathit{\gamma }}_{\mathrm{ad}}^{\mathrm{1}\mathit{/}\mathrm{4}}}\sqrt{\frac{\mathit{q}}{{\mathit{h}}_{\mathrm{p}}}}& & \end{array}$(14)where q = m_p/M_⋆ the planet mass ratio. This condition yields an estimation of the planet mass range for which the corotation torque remains unsaturated. We find: $\begin{array}{ccc}\mathrm{2.8}{\left(\frac{\mathit{\chi }}{{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{p}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{h}}_{\mathrm{p}}^{\mathrm{5}\mathit{/}\mathrm{3}}\lesssim \mathit{q}\lesssim \mathrm{2.8}{\left(\frac{\mathit{\chi }}{{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{p}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{h}}_{\mathrm{p}}& & \end{array}$(15)which gives 3.5 × 10^-6 ≲ q ≲ 3 × 10^-5 for our disc model. This implies that convergent migration toward the opacity transition is expected for planets with masses in the range 1 ≲ m_p ≲ 10  M_⊕. Bodies more massive than ~10  M_⊕ or less massive than ~1  M_⊕ will rather experience inward migration.

We also examined the issue of whether the torque experienced by a protoplanet can be altered by the close proximity of other bodies. Horn & Lyra (2012) indeed speculated that effects resulting from planet wakes may lead to a more rapid saturation of the corotation torque. To achieve this aim, we performed i) one simulation in which we measured the torques experienced by N = 9 protoplanets separated by 4.5  R_mH and held on a fixed circular orbit; ii) an additional set of calculations with N = 1 protoplanet and which differ in the value for the planet radial position. From Fig. 3 we see that the torques derived from these two series of runs differ only slightly, which suggests that the wakes generated by other low-mass planets have only a marginal effect on the saturation of the corotation torque.

When the planet eccentrity reaches a value such that its radial excursion bebomes larger than the half-width of the horseshoe region, we expect the corotation torque to be strongly attenuated. To model this effect, we follow Hellary & Nelson and multiply in the N-body runs the analytical corotation torque by a damping factor E. In order to estimate how E depends on e_p, we have performed a subset of calculations with a m_p = 3  M_⊕ protoplanet evolving on a fixed orbit with e_p in the range 0 ≤ e_p ≤ h_p. Since the half-width of the horseshoe region is a fraction of the disc scaleheight, we expect a null corotation torque when e_p ~ h_p, leaving only the differential Lindblad torque. Given that the differential Lindblad torque depends weakly on the eccentricity for e_p < h_p, the corotation torque can be determined in each simulation by simply substracting the differential Lindblad torque to the total torque. We plot the results of these simulations in Fig. 4. Superimposed is the function that is found to best reproduce the simulation results and which is given by: $\begin{array}{ccc}\mathit{E}\mathrm{=}\exp \left(\mathrm{-}\frac{{\mathit{e}}_{\mathrm{p}}^{\mathrm{2}}}{{\mathit{X}}_{\mathrm{s}}^{\mathrm{\prime }\mathrm{2}}}\right)\mathrm{·}& & \end{array}$(16)In the previous equation, ${\mathit{X}}_{\mathrm{s}}^{\mathrm{\prime }}$ is defined by ${\mathit{X}}_{\mathrm{s}}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}{\mathit{x}}_{\mathrm{s}}^{\mathrm{\prime }}$, where ${\mathit{x}}_{\mathrm{s}}^{\mathrm{\prime }}\mathrm{=}{\mathit{x}}_{\mathrm{s}}\mathit{/}{\mathit{a}}_{\mathrm{p}}$ is the dimensionless half-width of the horseshoe region. Compared with the prescriptions of Hellary & Nelson (2012) and Cossou et al. (2013) which are also plotted in Fig. 4, our formula tends to enhance the corotation torque by a factor of ~3 for ${\mathit{e}}_{\mathrm{p}}\mathrm{~}{\mathit{x}}_{\mathrm{s}}^{\mathrm{\prime }}$.

Fig. 4 Corotation torque damping factor as a function of the ratio between the eccentricity and the dimensionless half-width of the horseshoe region for mp = 3  M⊕, and deduced from hydrodynamical simulations (triangles). The solid line corresponds to the function that best fits the simulation results. The dashed line is the prescription of Hellary & Nelson (2012) and the dash-dotted line corresponds to the prescription of Cossou et al. (2013).

3.2. Eccentricity and inclination damping in N-body simulations

Fig. 5 Orbital evolution of N = 5 embryos with mass mp = 3  M⊕ and without stochastic forces included. The dashed line shows the location of the zero-migration line.

Eccentricity and inclination damping resulting from the interaction with the disc are modelled by applying the following accelerations to each body: $\begin{array}{ccc}{{a}}_{\mathrm{e}}\mathrm{=}\mathrm{-}\mathrm{2}\frac{\mathrm{(}{{v}}_{\mathrm{p}}\mathrm{·}{{R}}_{\mathrm{p}}\mathrm{)}{{R}}_{\mathrm{p}}}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{e}}}& & \end{array}$(17)and $\begin{array}{ccc}{{a}}_{\mathrm{i}}\mathrm{=}\mathrm{-}\frac{\mathrm{(}{{v}}_{\mathrm{p}}\mathrm{·}{ẑ}\mathrm{)}}{{\mathit{t}}_{\mathrm{i}}}{ẑ}& & \end{array}$(18)where ẑ is a unit vector in the vertical direction. t_e and t_i are the eccentricity and inclination damping timescales for which we use the prescriptions of Cresswell & Nelson (2008): $\begin{array}{ccc}{\mathit{t}}_{\mathrm{e}}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{wave}}}{\mathrm{0.78}}\left[\mathrm{1}\mathrm{-}\mathrm{0.14}{\left(\frac{{\mathit{e}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{0.06}{\left(\frac{{\mathit{e}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)}^{\mathrm{3}}\mathrm{+}\mathrm{0.18}\left(\frac{{\mathit{e}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right){\left(\frac{{\mathit{i}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\right]& & \end{array}$(19)and $\begin{array}{ccc}{\mathit{t}}_{\mathrm{i}}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{wave}}}{\mathrm{0.544}}\left[\mathrm{1}\mathrm{-}\mathrm{0.3}{\left(\frac{{\mathit{i}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{0.24}{\left(\frac{{\mathit{i}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)}^{\mathrm{3}}\mathrm{+}\mathrm{0.18}{\left(\frac{{\mathit{e}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)}^{\mathrm{2}}\left(\frac{{\mathit{i}}_{\mathrm{p}}}{{\mathit{h}}_{\mathrm{p}}}\right)\right]& & \end{array}$(20)where: $\begin{array}{ccc}{\mathit{t}}_{\mathrm{wave}}\mathrm{=}\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathit{m}}_{\mathrm{p}}}\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathrm{\Sigma }}_{\mathrm{p}}{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}}{\mathit{h}}_{\mathrm{p}}^{\mathrm{4}}{\mathrm{\Omega }}_{\mathrm{p}}^{-1}\mathit{.}& & \end{array}$(21)

4. Results of hydrodynamical simulations

4.1. Laminar viscous simulations

The orbital evolution of N = 5 embryos with mass m_p = 3M_⊕ and initially separated by 4.5R_mH is displayed in Fig. 5. At early times, the two innermost (resp. outermost) bodies located inside (resp. outside) the opacity transition tend to undergo outward (resp. inward) migration. Regarding the third body (cyan) initially located at R_p = 1.7, it tends to experience only a weak positive torque due to its close proximity to the opacity transition. Nevertheless, the strong differential migration between this body and the fourth one (green) quickly leads to the formation of a 9:8 mean motion resonance (MMR) between these two protoplanets at t ~ 3 × 10² yr. As revealed by Fig. 6 which displays the temporal evolution of the torques experienced by each planet, eccentricity growth due to resonant trapping makes the torque experienced by the third planet become negative, in such a way that the third and fourth planets subsequently migrate inward together while maintaining their 9:8 resonance. This arises because, just after resonant trapping, the value reached by the eccentricity of the third planet e_p ~ 0.015 is comparable to the dimensionless half-width of the planet’s horseshoe region which is estimated to be ${\mathit{x}}_{\mathrm{s}}^{\mathrm{\prime }}\mathrm{~}\mathrm{1.1}{\mathit{a}}_{\mathrm{p}}\sqrt{\mathit{q}\mathit{/}\mathit{h}}\mathrm{~}\mathrm{0.02}$. Consequently, the radial excursion that the planet undergoes is eventually larger than the horseshoe region, causing thereby the (positive) corotation torque to be significantly attenuated (Bitsch & Kley 2010; Cossou et al. 2013).

These two bodies then catch up with the outward-migrating second protoplanet (blue) and enter in a 9:8 MMR with it at t ~ 600 yr. Again, eccentricity growth due to resonant capture causes the positive torque exerted on the second body to weaken. Although it remains positive, its amplitude is not sufficient to counterbalance the negative torques exerted on the third and fourth planets, and this three-planet system consequently suffers a slow, inward resonant migration. This proceeds until the fifth body (orange) catches up with the fourth protoplanet (green) and enters a 9:8 resonance with it at t ~ 1000 yr. At t ~ 2000 yr, the outward-migrating innermost planet (black) becomes trapped in a 9:8 MMR with the second protoplanet (blue) which, from this time onward, undergoes a marginally positive torque due to the high value reached by its eccentricity. However, as can be seen in the lower panel of Fig. 5, eccentricity pumping due to resonant interaction is relatively modest for the innermost core, with an equilibrium value of e_p ~ 5 × 10^-3. Consequently, the fraction of the corotation torque operating on the innermost planet is large enough for the total torque exerted on this body to remain positive, which is confirmed by looking at the evolution of the torque exerted on that planet in Fig. 6. This effect, however, is clearly not sufficient to couterbalance the negative torques experienced by the outer bodies so that at late times, the five bodies tend migrate inward in lockstep with each neighbouring pair of planets forming a 9:8 resonance. As the outer planets pass through the zero-torque radius, however, we expect the disc torques exerted on these planets to become positive so that it is likely that the swarm will ultimately stop migrating once the net torque acting on the whole system cancels (Cossou et al. 2013).

Fig. 6 Time evolution of the torque experienced by the protoplanets in the case with N = 5 and without stochastic forces included.

Fig. 7 Left panel: orbital evolution of N = 7 embryos with mass mp = 3  M⊕ and without stochastic forces included. Right panel: same but with N = 9. The dashed line shows the location of the zero-migration line.

In Fig. 7 we present, the orbital evolution of simulations with N = 7 (left panel) and N = 9 (right panel) embryos embedded in the same disc model. We remind the reader that the embryos are initially located in such a way that bodies of the inner half migrate outward whereas bodies of the outer half migrate inward. Because they are initially located on either side of the convergence line, the fourth (cyan) and fifth (green) bodies are the first to become trapped in a 9:8 MMR. Again, eccentricity growth due to resonant trapping causes the corotation torque operating on the fourth body to be significantly attenuated so that these two planets tend to migrate inward at later times. This is exemplified in the upper panel of Fig. 8 which shows the evolution of the torque exerted on each planet as a function of time. At t ~ 600 yr, the third body (blue) catches up with the fourth body (cyan) and enters a 9:8 MMR with it. A general trend is that inward-migrating bodies are captured in resonance later than outward-migrating bodies. This is a direct consequence from the corotation damping effect, which tends to make the resonant swarm migrate inward, strengthening thereby the differential migration rate with the inner, outward-migrating bodies.

Once again, the evolution outcome consists of inward migration of a group of N = 7 members which are in mutual mean motion resonances, with the resonance being 9:8 except for the two innermost bodies which are in 8:7 resonance. This arises because prior to resonant capture of the second body (purple), the torque exerted on this planet is slightly stronger in comparison with that experienced by the innermost one (black), resulting in divergent migration between these two cores. Comparing Figs. 6 and 8, we see that after ~10⁴ yr, only one planet feels a positive torque in the case where N = 5 while two cores are subject to a positive torque in the simulation with N = 7. This confirms the expectation that the resulting resonant system becomes more compressed as N increases, and consequently more prone to dynamical instability.

Indeed, increasing the number of initial embryos to N = 9 resulted in a more chaotic behaviour where protoplanets suffered close encounters and collisions, as illustrated in the right panel of Fig. 7 which displays the planets’ orbital evolution for that case. At early times, evolution proceeds similarly to that corresponding to N = 7, with a system of 9 protoplanets with each body locked in a 9:8 MMR with its neighbours being formed at t ~ 3000 yr. From that time onward, it can be seen in the lower panel of Fig. 8, which displays the temporal evolution of the torque exerted on each planet, that the five innermost embryos located inside the convergence line undergo a positive torque whereas the others undergo a negative torque, implying a significantly compressed resonant system. At t ~ 3500 yr, this leads to a physical collision between the sixth (green) and seventh (light green) embryos, forming thereby a 6  M_⊕ planet which subsequently undergoes inward migration since it is located outside the convergence line. The resulting perturbation then propagates to the inner system and causes two additional collisions at later times, between the third (dark blue) and fourth (blue) planets at t ~ 4000 yr and between the two innermost cores (black+purple) at t ~ 4100 yr. As revealed by the lower right panel of Fig. 7, these two newly formed ~6  M_⊕ planets are located inside the convergence line and have relatively modest eccentricity. This, combined with the fact that the fully-unsaturated corotation torque scales as ${\mathit{m}}_{\mathrm{p}}^{\mathrm{2}}$ (or equivalently as ${\mathit{x}}_{\mathrm{s}}^{\mathrm{4}}$) for low-mass planets, implies that these planets are subject to a stronger positive torque.

Fig. 8 Upper panel: time evolution of the torque experienced by the protoplanets in the case with N = 7 and without stochastic forces included. Lower panel: same but for N = 9.

Not surprisingly, the significant resulting differential migration between the outermost 6M_⊕ planet (blue) and its exterior inward-migrating 3  M_⊕ neighbour (cyan) leads to capture in a 6:5 resonance at t ~ 7000 yr. Eccentricity growth due to resonant trapping causes the corotation torque operating on the 6M_⊕ to be partially attenuated, causing the two planets migrate inward together at later times while maintaining their 6:5 resonance. This process enables the innermost 6  M_⊕ planet (black) to catch up with the outermost 6  M_⊕ body (blue) and to enter in a 5:4 resonance with it at t ~ 8500 yr. Protoplanets located outside the convergence line then become trapped in MMRs with the inner ones at later times so that the evolution outcome for that case corresponds again to the formation of a stable resonant system where the two inner planets are the more massive. Since the two inner planets are trapped in resonance with sustained eccentricities, they only feel a weak positive torque despite being located interior to the convergence zone. With no strong outward push, the whole system of embryos migrates inward. The high computational cost makes it impossible to perform hydrodynamical simulations for more than ~10⁴ yr so that the final fate of the system remains uncertain. Altough it can not be excluded that additional collisions will arise at later times, a possible issue is that the planets will reach stationary orbits once the net torque acting on the resonant system, which consists of the sum of the attenuated corotation torques and unattenuated differential Lindlbad torques exerted on each planet, cancels (Cossou et al. 2013). In Sect. 5, we examine in more details the possible long-term evolution outcomes using N-body simulations.

In order to unveil the role of the zero-migration line on the collisions events that were observed in this simulation, we have performed a similar simulation with initially N = 9 embryos but using an isothermal equation of state. In that case, all protoplanets are expected to experience inward migration due to their interaction with the gas disc. Figure 9 shows the evolution of the planets’ semimajor-axes and eccentricities as a function of time for the isothermal run. We note that although we consider here equal-mass planets, these tend to undergo convergent migration because the surface density profile has index σ < 3/2 (e.g. Pierens et al. 2011). Compared with the radiative simulation, however, the convergent migration rate is sufficiently weak in the isothermal case to prevent close encounters between embryos. This clearly demonstates that near the transition between the outward and inward migration regimes, close encounters and collisions can be stimulated due to the strong differential migration experienced by the embryos there.

Fig. 9 Time evolution of the semimajor-axes and eccentricities of protoplanets for an isothermal hydrodynamical simulation with N = 9.

4.2. Effects of disc turbulence

Fig. 10 Snapshot of the perturbed surface density profile in the hydrodynamical simulation with initially N = 5 embryos and with stochastic forces included. The dashed circle shows the location of the zero-migration line.

Fig. 11 Top panel: time evolution of the semimajor-axes of protoplanets in a hydrodynamical simulation with N = 5 and with stochastic forces included. Middle panel: same but for N = 7. Bottom panel: same but for N = 9. The dashed line shows the location of the zero-migration line.

The results presented above indicate that, in the limit of a moderate initial number of objects, the formation of mean motion resonances prevents embryos from undergoing close encouters with other bodies. In this section, we examine how the stability of these resonant configurations is affected when a level of disc turbulence is accounted for. Previous work (e.g. Pierens et al. 2011) suggested that mean motion resonances are likely to be disrupted by stochastic torques in the active regions of protoplanetary discs and within disc lifetimes.

To examine this issue, we have perfomed a series of simulations for N = 5,7,9 in which each body is subject to a an additional stochastic force F_turb =  −∇Φ_turb where Φ_turb is given by Eq. (7). The temporal evolution of the planets’ semi-major axis for these three simulations is displayed in Fig. 11. For N = 5, Fig. 10 presents a contour plot of the perturbed surface density distribution at the beginning of the simulation. For protoplanets initially located inside the zero-migration radius, surface density perturbations inside the planets’ horseshoe regions and related to the co-orbital dynamics are clearly visible. Due to the negative entropy gradient, co-orbital dynamics leads to a positive (resp. negative) surface density pertubation ahead (resp. behind) of the planet, giving rise to a positive corotation torque. For protoplanets located outside the zero-migration line, these additional surface density perturbations are weaker, indicating a relatively faint corotation torque in that case. For this run, a sequence of 9:8 MMRs is ultimately formed so that the final fate of the system is similar to that obtained in the laminar simulation. Although these resonances are stable on average, the corresponding resonant angles are observed to oscillate between periods of circulation and libration, implying a weaker resonant locking in presence of turbulence. We notice that a similar behaviour was observed in previous studies on the capture in resonance of pairs of planets embedded in a turbulent isothermal disc (Pierens et al. 2011). This stengthens the expectation that turbulence can break resonant configurations when the typical amplitude of the stochastic density fluctuations is large enough, and consequently that the resonant systems obtained using viscous laminar disc models are more prone to destabilization in presence of turbulence. This is therefore not too surprising that collisions arise in the turbulent run with N = 7, as can be seen in the middle panel of Fig. 11. Here, the resonant system that is formed at the end of the simulations is composed of two 6  M_⊕ bodies (black+blue lines) located at the inner edge of the swarm and trapped in a 5:4 resonance plus three exterior 3  M_⊕ protoplanets trapped in a first-order p+1:p resonance with its neighbours, where we observe a clear tendency for p to increase as one moves out through the swarm.

The lower panel of Fig. 11 shows the evolution of the system for the run with N = 9. Again, two 6  M_⊕ planets (purple+blue) are formed in the course of the evolution and the final fate corresponds again to the formation of a resonant chain with the two more massive planets located in the inner half of the swarm. Although three collisions ocurred in the laminar run with N = 9, it is clear, when comparing Fig. 11 with Fig. 7, that the evolution of the system is much more chaotic in the turbulent case.

4.3. Effect of the disc model

To investigate how the disc model affects the results presented above, we have performed a series of simulations using a disc model with mass twice that of the fiducial disc model, which corresponds to a the snow-line located at ~4.2 AU. Again, the initial mass of embryos is 3  M_⊕ and their initial positions are chosen such that half of the initial population of embryos is located inside the zero-migration line, which lies at ~2.4 AU for this disc model, whereas the other half is located outside. Although not shown here, runs with stochastic forces not included show evolution outcomes consistent with those of the fiducial disc model, with growth of embryos and formation of 6  M_⊕ protoplanets observed when the initial number of bodies is N ≥ 7. This is not surprising since as the disc mass increases, the effect of a faster differential migration toward the zero-migration line, and which would promote close encounters between embryos, is couterbalanced by a stronger disc-induced eccentricity damping.

For runs with stochastic forces included, however, the outcomes are found to differ noticeably from those of the fiducial disc model due to more vigorous turbulent fluctuations. In that case, collisions are found to arise even for low values of the initial number of embryos, as illustrated in Fig. 12 where is displayed, for this disc model, the evolution of the planets’ semimajor-axes for N = 5. Here, a 6  M_⊕ planet (cyan) is formed at the outer edge of the swarm at t ~ 5 × 10³ yr and is likely to become trapped at the location of the zero-torque radius at later times. From t ~ 1.3 × 10⁴ yr, this planet tends to separate from the three innermost 3  M_⊕ bodies which experience resonant inward migration with each planet forming a 8:7 resonance with its neighbours. Prior to the formation of this inner resonant system, it is interesting to note that the two innermost bodies (black+blue) entered in a coorbital 1:1 resonance and which remained stable for ~ 5 × 10³ yr.

Again, the final fate of the system remains uncertain but it is likely that the group composed of the three innermost planets will end up on a stationary orbit once the net torque acting on this three-planet system will cancel out (Cossou et al. 2013).

5. Results of N-body runs

Fig. 12 Time evolution of the semimajor-axes of protoplanets in a turbulent hydrodynamical simulation with initially N = 5 embryos and for a disc mass equivalent to two times that of the fiducial disc model. The dashed line shows the location of the zero-migration line.

In Sect. 3, we have presented how analytical formulae for the Type I torques can be calibrated using the results of hydrodynamical simulations. Here, we present the main results that emerge from ~600 N-body runs of protoplanets embedded in a radiative protoplanetary disc that we have performed using the prescriptions we obtained for Type I migration. The aims of this alternative approach are to i) study the long-term evolution of a swarm of protoplanets migrating toward the convergence line; and ii) examine the statistical properties of the planetary systems which are formed through this process.

Fig. 13 Left: time evolution of protoplanets’ orbital positions (top panel) and masses (lower panel) for a N-body run with intially N = 9 embryos and without stochastic forces included. Right: same but with stochastic forces included.

5.1. Comparison with hydrodynamical simulations

5.1.1. Runs without stochastic forces included

We first describe the range of outcomes that are observed in simulations with initial conditions chosen as close as possible to those for the hydrodynamical simulations presented in Sect. 4. We performed 100 N-body simulations for situations with N = 7 and N = 9 initial embryos. These simulations did not include stochastic forcing and differed only in the initial azimuthal positions of the embryos. For N = 7, collisions occured in 53% of the simulations. Objects as large as 9 M_⊕ formed in 6% of the simulations. In the runs with N = 9, virtually all (95%) of the simulations included collisions, as expected from the hydrodynamical simulations presented in Sect. 4. Most (77%) simulations only produced 6 M_⊕ planets, but a significant number (17%) formed 9 M_⊕ planets and one simulation formed a 12 M_⊕ core.

The left panel of Fig. 13 displays the evolution as a function of time of the planets’ orbital positions and masses for a run representative of the range of the observed outcomes. Overall, the early stages of evolution are consistent with what is seen in the corresponding hydrodynamical simulation, involving outward (resp. inward) migration of the innermost (resp. outermost) bodies plus formation of two 6  M_⊕ protoplanets inside the zero-migration line by ~10⁴ years. It is interesting to notice that of these two 6  M_⊕, one is formed at ~8 × 10³ years through the merging of two co-orbital planets that entered in a 1:1 resonance at ~ 2 × 10³ years. An additional 6  M_⊕ body is produced at t ~ 1.8 × 10⁴ yr as a result of the collision between two 3  M_⊕ objects. After ~ 3 × 10⁴ years, the systems attains a quasi-stationary state and consists of three 6  M_⊕ objects plus two 3  M_⊕ bodies trapped in a resonant chain and evolving on non-migrating orbits. Moving from inward to outward, the resonances which are formed are 5:4, 7:6, 7:6, 8:7 and 6:5. This stationary configuration is reached when the positive torque exerted on the two innermost 6  M_⊕ bodies (blue and black lines in Fig. 13) counterbalances the torques experienced by the three other bodies, leading to a zero net torque acting on the whole system (Cossou et al. 2013).

5.1.2. Effect of stochastic forces

To illustrate the dependence of the results presented above on the presence of disc turbulence, we have performed a set of simulations but with stochastic forces acting on the protoplanets included. We remind the reader that hydrodynamical simulations in which stochastic forces are included (see Sect. 4.2) indicated that collisions are more likely to occur in that case due to the general tendency for turbulent fluctuations to break resonant configurations. In agreement with this result, we find that collisions occured in 97% of the N-body runs performed with N = 7 and in 100% of the simulations performed with N = 9. Therefore, it is not surprising to observe a clear trend for forming more massive embryos when stochastic forces are included. For example, 36% of the runs with N = 7 resulted in the production of 9  M_⊕ bodies whereas this number increases to 50% for N = 9. In this latter case, two simulations resulted in the formation of 12  M_⊕ protoplanets. Moreover, we note that we performed a series of turbulent runs with initial separations of 6  R_mH, and which resulted in the production of 9  M_⊕ cores in 30 of a total of 100 simulations, which suggests that these statistics are relatively robust with regards to the value for the initial separation of embryos.

Fig. 14 Left panel: distribution of mass of the most massive planet which is formed in N-body runs with stochastic forces included and in which the mass of each embryo is sampled from a Gaussian distribution with μm = 3  M⊕ and σm = 1  M⊕. Right panel: same but for μm = 1  M⊕ and σm = 0.5  M⊕.

The right panel of Fig. 13 illustrates the evolutionary outcome for a simulation in which two 12  M_⊕ planets are produced. At t ~ 10⁵ yr, a non-migrating compact system is formed with two co-orbital planets of 3  M_⊕ and 6  M_⊕ located at the inner edge of the swarm and in a 7:6 resonance with an exterior 3  M_⊕ protoplanet. At t ~ 1.2 × 10⁵ years, the formation of a 12  M_⊕ body through the collision of two exterior 6  M_⊕ bodies (yellow+red) destabilizes the inner part of the swarm, which subsequently leads to the merging of the three inner planets into an additional 12  M_⊕ embryo. From this time, the system is composed of two 12  M_⊕ bodies plus an exterior 3  M_⊕ embryo which evolves on a quasi-circular orbit at the nominal convergence line. The two 12M_⊕ planets enter in a 6:5 resonance at later times and resonantly migrate inward until they reach the inner edge of the disc. This configuration remained stable until the end of the simulation which was evolved for 5 × 10⁵ yr. We note that although these planets evolve inside the zero-migration line, they experience a negative torque from the disc because the (positive) corotation torque is saturated for such a planet mass (see Sect. 3).

5.2. Effect of changing the initial mass distribution and statistical overview

In order to examine the dependence of our results on the initial mass distribution, we performed two additional sets of 100 simulations using a randomised mass distribution. In the first series of runs, the mass of each embryo is sampled from a Gaussian distribution with mean μ_m = 3  M_⊕ and standard deviation σ_m = 1  M_⊕, while in the second set of simulations, we set μ_m = 1  M_⊕ and σ_m = 0.5  M_⊕. The total mass of embryos is chosen to be 30  M_⊕ in both cases. Moreover, we choose the orbital position of the inner planet to be uniformely distributed between 1.1 ≤ a_p ≤ 1.2 AU whereas the initial separations of other planets are set to nR_mH, where 4 ≤ n ≤ 5 is randomly chosen according to a uniform distribution.

Figure 14 shows the distribution of mass of the most massive planet which is formed in these two series of simulations. In the case where μ_m = 3  M_⊕ and σ_m = 1  M_⊕ (left panel), we see that ~7  M_⊕ embryos are most common but of the 100 simulations, a significant number of them (43%) resulted in the formation of protoplanets with masses in the range 8 ≤ m_p ≤ 10  M_⊕ within 5    × 10⁵ yr. This is in good agreement with the results of simulations with initially equal-mass embryos presented in Sect. 5.1.2. Giant planet cores with masses ≥10  M_⊕ were produced in 8% of the simulations.

In the case where μ_m = 1  M_⊕ and σ_m = 0.5  M_⊕ (right panel), 14% of runs resulted in the formation of bodies with masses in the range 8 ≤ m_p ≤ 10  M_⊕, with only one run leading to the formation of a core with mass ≥10  M_⊕. We believe this is related to the fact that the initial population of embryos with initial masses ≲1  M_⊕ experience corotation torque saturation (see Sect. 3). Consequently, they undergo inward Type I migration even these are located inside the opacity transition, which tends to substantially reduce the process of convergent migration at the planet trap.

These results indicate that forming giant planet cores at the zero-torque radius is likely to occur provided it involves massive impacts between bodies of a few Earth masses which do not experience corotation torque saturation and which formed earlier through an alternative process.

6. Discussion

In this section, we discuss how the location of the opacity transition we considered here depends on the disc model that is adopted and examine under which conditions collisional planetary growth at this opacity transition is expected to arise.

In order to determine the location of the opacity transition, we balance radiative losses with viscous heating, which gives: $\begin{array}{ccc}\mathrm{2}\mathit{\sigma }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{9}}{\mathrm{4}}\mathit{\nu }\mathrm{\Sigma }{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{·}& & \end{array}$(22)For an optically thick disc, we remind that the midplane temperature T is related to the effective temperature T_eff by (see Eq. (4)): $\begin{array}{ccc}{\mathit{T}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{8}}\mathit{\tau }{\mathit{T}}_{\mathrm{eff}}^{\mathrm{4}}\mathit{.}& & \end{array}$(23)Given that τ = κΣ/2, Eq. (22) can be recast as: $\begin{array}{ccc}\mathit{\sigma }{\mathit{T}}^{\mathrm{4}}\mathrm{=}\frac{\mathrm{27}}{\mathrm{128}}\mathit{\nu \kappa }{\mathrm{\Sigma }}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}\mathit{.}& & \end{array}$(24)Moreover, in the context of the standard α prescription of Shakura & Sunyaev (1973), the effective kinematic viscosity is given by ν = αc_sH, where c_s = (γ_adℛT/μ)^1/2 (ℛ and μ are the gas constant and the mean molecular weight respectively) is the sound speed and H = c_s/Ω the disc scale height. Evaluating the previous equation at the opacity transition gives the following expression for the location of the transition radius R_t: $\begin{array}{ccc}\frac{{\mathit{R}}_{\mathrm{t}}}{{\mathit{R}}_{\mathrm{0}}}\mathrm{=}{\left(\frac{\mathrm{27}}{\mathrm{128}}\frac{{\mathit{\gamma }}_{\mathrm{ad}}\mathrm{ℛ}}{\mathit{\mu \sigma }}\frac{{\mathit{\kappa }}_{\mathrm{t}}{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{0}}}{{\mathit{T}}_{\mathrm{t}}^{\mathrm{3}}}\mathit{\alpha }\right)}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{2}\mathit{\sigma }\mathrm{+}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}}& & \end{array}$(25)where T_t and κ_t denote the values for the temperature and opacity at the transition respectively. Here, we have used the fact that Σ = Σ₀(R/R₀)^− σ and Ω = Ω₀(R/R₀)^− 3/2, where Ω₀ is the angular velocity at R = R₀.

Furthermore, convergent migration at the opacity transition for bodies of a few Earth masses arises provided that the corotation torque remains unsaturated. We expect the corotation torque to remain close to its fully unsaturated value provided that (e.g. Baruteau & Masset 2013): $\begin{array}{ccc}{\mathit{t}}_{\mathrm{U}\mathrm{-}\mathrm{turn}}\mathrm{\le }{\mathit{t}}_{\mathrm{dif}}\mathrm{\le }{\mathit{t}}_{\mathrm{lib}}\mathit{/}\mathrm{2.}& & \end{array}$(26)For given values of α and the mass ratio q = m_p/M_⋆, this condition provides an estimation of the range of radii for which the corotation torque is fully unsaturated. We find: $\begin{array}{ccc}\mathrm{0.52}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{R}}_{\mathrm{0}}}{{\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{t}}}\right)}^{\mathrm{2}}{\mathit{q}}^{\mathrm{6}\mathit{/}\mathrm{7}}{\mathit{\alpha }}^{\mathrm{-}\mathrm{4}\mathit{/}\mathrm{7}}\mathrm{\le }\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{0}}}\mathrm{\le }\mathrm{0.025}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{R}}_{\mathrm{0}}}{{\mathit{c}}_{\mathrm{s}\mathit{,}\mathrm{t}}}\right)}^{\mathrm{2}}{\mathit{q}}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{\alpha }}^{\mathrm{-}\mathrm{4}\mathit{/}\mathrm{9}}& & \end{array}$(27)where c_s,t is the value of the sound speed at the opacity transition. In the top panel of Fig. 15 we plot the location of R_t as a function of the α viscous stress parameter and for the disc model that we employed in the simulations. The region located in between the dashed and solid-dashed lines in Fig. 15 corresponds to that where the corotation torque remains unsaturated for a 3  M_⊕ planet, and which is defined by Eq. (27). Therefore, the intersection between the solid line and the shaded area in Fig. 15 represents the range of radii where convergent migration of 3  M_⊕ protoplanets can occur. For this disc model, this process is expected to arise in the region from ~1.5 to ~2.2 AU, in good agreement with the location of the zero-torque radius in our hydrodynamical and N-body simulations.

Fig. 15 Top panel: location of the opacity transition (solid line) as a function of the value for the α viscous stress parameter and for our fiducial disc model. The dashed line shows the location where the viscous timescale tdif becomes equal to half of the horseshoe libration timescale (tlib) whereas the solid-dashed line shows the radius where tdif equals the U-turn timescale. The intersection between the shaded area and the solid line corresponds to the range of radii where convergent migration for 3  M⊕ bodies is expected to occur. Lower panel: same but for a disc model ten times more massive than the MMSN.

For a given value of α, Eq. (25) predicts that the location of the opacity transition will increase with Σ. However, looking at the upper panel of Fig. 15, it is clear that as Σ increases, the range of α values for which the corotation torque remains unsaturated is shifted toward lower values. In the context of the formation of the giant planets in the solar system, this implies that both a significantly massive disc and a relatively low value for α will be needed for the convergent migration mechanism to operate in the Jupiter-Saturn region. This is illustrated in the lower panel of Fig. 15 which shows that for a disc model with σ = 1.5, convergent migration at ~4 AU arises provided that the disc mass corresponds to ten times the MMSN. Of course, these results will strongly depend on the opacity transition that is considered. We nevertheless expect a similar mechanism to occur at any opacity transition interior to which the temperature gradient is steep enough to allow for outward migration of Earth-mass bodies. Testing other opacity transitions is beyond the scope of that paper but we will discuss in a future paper the influence of changing the opacity table. Finally, we notice that formation of giant planet cores through collisions of Earth-mass bodies may also be possible at other kinds of planet traps like those arising at a dead zone or at the radius where stellar heating begins to take over viscous heating (Hasegawa & Pudritz 2011).

7. Conclusion

We have presented the results of both hydrodynamical and N-body simulations of the evolution of a swarm of Earth-mass protoplanets that gravitationally interact near a zero-migration line. The embryos are initially located around an opacity transition which plays the role of a planet trap where the Type I torque cancels. For bodies with masses in the range 1 − 10  M_⊕ and evolving inside the opacity transition, Type I migration proceeds outward whereas it is directed inward if these are located further out. The main aim of this work is to examine the possible outcomes that arise when low-mass protoplanets convergently migrate toward a zero-migration radius and in particular whether giant planet cores can be formed at such places through giant impacts between embryos of a few Earth masses.

Hydrodynamical simulations show that equal-mass embryos with mass of 3  M_⊕ located on both sides of a convergence zone tend to enter in a resonant chain whose stability depends on the initial number of objects and whether or not planets experience stochastic forces due to turbulence. For a limited number of bodies and in the absence of stochastic forcing, a sequence of resonances appears to be stable so that close encounters betweeen embryos are prevented. Increasing the initial number of protoplanets however leads to a significant compression of the system and eventually to the destabilization of these resonant chains. Formation of 6  M_⊕ protoplanets from 3  M_⊕ embryos is observed to occur in that case on a timescale ≤10⁴ yr.

Not surprisingly, including a moderate level of turbulence corresponding to a value for the αviscous stress parameter of α ~ 2 × 10^-3 clearly enhances this process of collisional planetary growth. Interestingly, we find that a significant fraction (~ 50%) of the N-body runs with stochastic forces included and performed with masses randomly sampled from a Gaussian distribution with μ_m = 3  M_⊕ resulted in the formation of giant planet cores with mass ≥8 − 10  M_⊕ in ~ 5 × 10⁵ yr. For a randomised mass distribution with μ_m = 1  M_⊕, however, only ~ 15% of the N-body runs produced giant planet cores. We conclude that forming giant planet cores at convergence zones is efficient provided that it involves collisions between embryos with mass ≳2  M_⊕ and which formed earlier according to the classical runaway/oligarchic growth scenario (Levison et al. 2010) or through accretion of cm-sized pebbles (Morbidelli & Nesvorny 2012).

For a protoplanetary disc with mass typical of the MMSN, we find that the mechanism presented here may allow the in-situ formation of giant planet cores from Earth-mass bodies at ~2 AU whereas very massive discs (≥10 times the MMSN) are required to form ~ 10  M_⊕ bodies at 4 − 5 AU. We note however that we considered here a zero-migration line corresponding to a particular change in the opacity regime. Other planet traps located further than ~4 AU are expected to arise in typical protoplanetary disc models. For example, the transition where stellar irradiation begins to provide most of the heating of the disc gives rise to an additional planet trap located at 2–40 AU, depending on the mass accretion rate (Hasegawa & Pudritz 2011). Planet traps can also exist at locations where the thermal diffusion timescale becomes longer than the horseshoe libration timescale, leading to a saturation of the corotation torque in the outer regions. Since the horseshoe libration timescale depends on the planet mass, planets with different masses tend to converge toward different radii. We will focus on the role of these additional planet traps on the formation of giant planet cores in a future publication.

The N-body runs that we have presented are the simplest we can perform. One limitation of our work is that the embryos do not accrete planetesimals or pebbles as they migrate. If there is a sufficient supply of pebbles, accretion may be very rapid for embryos of ~1  M_⊕ (Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012), with the consequence that embryo masses may change significantly over the timescales of the simulations. Morevover, the position of the zero-migration line remains fixed in the course of the simulations. In a more realistic scenario, we expect the convergence zone to move inward as the disc disperses (Lyra et al. 2010; Horn et al. 2012). We will present in a forthcoming paper the results of both hydrodynamical and N-body simulations that account for the inward migration of the zero-torque radius due to photoevaporation and irradiation from the central star effects.

Acknowledgments

Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and by HPC resources from GENCI-cines (c2012046957).

References

All Figures

	Fig. 1 Distribution of the specific turbulent torque exerted on a 3  M⊕ planet obtained from a turbulent simulation in which the turbulent potential of Eq. (7) with γ = 10-4 is applied to the disc (solid line), and obtained from a laminar simulation in which the turbulent force of Eq. (10) with C = 17 is applied to the planet (dashed line).
In the text

	Fig. 2 Surface density and temperature profiles at steady-state.
In the text

	Fig. 3 Upper panel: torque exerted on a 3  M⊕ planet as a function of its orbital distance. The solid line corresponds to the analytical prescription of Paardekooper et al. (2011). Triangles and diamonds are the results from hydrodynamical simulations with/without the effects due to the waves generated by the other planets included. Lower panel: same but for mp   =   6  M⊕.
In the text

Fig. 4 Corotation torque damping factor as a function of the ratio between the eccentricity and the dimensionless half-width of the horseshoe region for mp = 3  M⊕, and deduced from hydrodynamical simulations (triangles). The solid line corresponds to the function that best fits the simulation results. The dashed line is the prescription of Hellary & Nelson (2012) and the dash-dotted line corresponds to the prescription of Cossou et al. (2013).

In the text

	Fig. 5 Orbital evolution of N = 5 embryos with mass mp = 3  M⊕ and without stochastic forces included. The dashed line shows the location of the zero-migration line.
In the text

	Fig. 6 Time evolution of the torque experienced by the protoplanets in the case with N = 5 and without stochastic forces included.
In the text

	Fig. 7 Left panel: orbital evolution of N = 7 embryos with mass mp = 3  M⊕ and without stochastic forces included. Right panel: same but with N = 9. The dashed line shows the location of the zero-migration line.
In the text

	Fig. 8 Upper panel: time evolution of the torque experienced by the protoplanets in the case with N = 7 and without stochastic forces included. Lower panel: same but for N = 9.
In the text

	Fig. 9 Time evolution of the semimajor-axes and eccentricities of protoplanets for an isothermal hydrodynamical simulation with N = 9.
In the text

	Fig. 10 Snapshot of the perturbed surface density profile in the hydrodynamical simulation with initially N = 5 embryos and with stochastic forces included. The dashed circle shows the location of the zero-migration line.
In the text

	Fig. 11 Top panel: time evolution of the semimajor-axes of protoplanets in a hydrodynamical simulation with N = 5 and with stochastic forces included. Middle panel: same but for N = 7. Bottom panel: same but for N = 9. The dashed line shows the location of the zero-migration line.
In the text

	Fig. 12 Time evolution of the semimajor-axes of protoplanets in a turbulent hydrodynamical simulation with initially N = 5 embryos and for a disc mass equivalent to two times that of the fiducial disc model. The dashed line shows the location of the zero-migration line.
In the text

	Fig. 13 Left: time evolution of protoplanets’ orbital positions (top panel) and masses (lower panel) for a N-body run with intially N = 9 embryos and without stochastic forces included. Right: same but with stochastic forces included.
In the text

	Fig. 14 Left panel: distribution of mass of the most massive planet which is formed in N-body runs with stochastic forces included and in which the mass of each embryo is sampled from a Gaussian distribution with μm = 3  M⊕ and σm = 1  M⊕. Right panel: same but for μm = 1  M⊕ and σm = 0.5  M⊕.
In the text

Fig. 15 Top panel: location of the opacity transition (solid line) as a function of the value for the α viscous stress parameter and for our fiducial disc model. The dashed line shows the location where the viscous timescale tdif becomes equal to half of the horseshoe libration timescale (tlib) whereas the solid-dashed line shows the radius where tdif equals the U-turn timescale. The intersection between the shaded area and the solid line corresponds to the range of radii where convergent migration for 3  M⊕ bodies is expected to occur. Lower panel: same but for a disc model ten times more massive than the MMSN.

In the text