Issue 
A&A
Volume 558, October 2013



Article Number  A105  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201322123  
Published online  14 October 2013 
Making giant planet cores: convergent migration and growth of planetary embryos in nonisothermal discs
^{1} Université de Bordeaux, Observatoire Aquitain des Sciences de l’Univers, BP 89 33271 Floirac Cedex, France
email: arnaud.pierens@obs.ubordeaux1.fr
^{2} Laboratoire d’Astrophysique de Bordeaux, BP 89 33271 Floirac Cedex, France
Received: 24 June 2013
Accepted: 5 August 2013
Context. Rapid gas accretion onto gas giants requires the prior formation of ~ 10 M_{⊕} cores, and this presents a continuing challenge to planet formation models. Recent studies of oligarchic growth indicate that in the region around 5 AU growth stalls at ~ 2 M_{⊕}. Earthmass bodies are expected to undergo Type I migration directed either inward or outward depending on the thermodynamical state of the protoplanetary disc. Zones of convergent migration exist where the Type I torque cancels out. These “convergence zones” may represent ideal sites for the growth of giant planet cores by giant impacts between Earthmass embryos.
Aims. We study the evolution of multiple protoplanets of a few Earth masses embedded in a nonisothermal protoplanetary disc. The protoplanets are located in the vicinity of a convergence zone located at the transition between two different opacity regimes. Inside the convergence zone, Type I migration is directed outward and outside the zone migration is directed inward.
Methods. We used a gridbased hydrodynamical code that includes radiative effects. We performed simulations varying the initial number of embryos and tested the effect of including stochastic forces to mimic the effects resulting from disc turbulence. We also performed Nbody runs calibrated on hydrodynamical calculations to follow the evolution on Myr timescales.
Results. For a small number of initial embryos (N = 5–7) and in the absence of stochastic forcing, the population of protoplanets migrates convergently toward the zerotorque radius and forms a stable resonant chain that protects embryos from close encounters. In systems with a larger initial number of embryos, or in which stochastic forces were included, these resonant configurations are disrupted. This in turn leads to the growth of larger cores via a phase of giant impacts between protoplanets, after which the system settles to a new stable resonant configuration. Giant planets cores with masses ≥ 10 M_{⊕} formed in about half of the simulations with initial protoplanet masses of m_{p} = 3 M_{⊕} but in only 15% of simulations with m_{p} = 1 M_{⊕}, even with the same total solid mass.
Conclusions. If 2 − 3 M_{⊕} protoplanets can form in less than ~1 Myr, convergent migration and giant collisions can grow giant planet cores at Type I migration convergence zones. This process can happen fast enough to allow for a subsequent phase of rapid gas accretion during the disc’s lifetime.
Key words: accretion, accretion disks / planets and satellites: formation / hydrodynamics / methods: numerical
© 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 kmsized 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 snowline 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 minimummass solar nebula (MMSN; Hayashi 1981). Moreover, recent Nbody 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 quasicircular, 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 cmsized 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 20cm 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 Earthmass 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 lowmass 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 nonisothermal, 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 zeromigration 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 zeromigration lines on the formation of giant planet cores. Sandor et al. (2011) studied this process in isothermal discs where zerotorque radii are located at deadzone boundaries and found that 10 M_{⊕} bodies can be formed in less than ~ 5 × 10^{5} yr through collisions of smaller embryos. The case of nonisothermal discs was investigated by Hellary & Nelson (2012) who performed Nbody simulations of planetary growth in radiativelyinefficient protoplanetary discs. They showed that in nonisothermal 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 subEarth mass embryos in 2–3 Myr.
In this paper, we present the results of hydrodynamical simulations of the interaction of multiple protoplanets in nonisothermal disc models. Our simulations begin with 3 M_{⊕} bodies with positions initially distributed around an opacity transition located just inside the snowline. This opacity transition corresponds to a zerotorque 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 Earthmass protoplanets which convergently migrate toward a zeromigration 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 Nbody simulations and which employ prescribed forces for migration, hydrodynamical simulations allow a selfconsistent 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^{4} 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 Nbody runs calibrated on our hydrodynamical simulations. These Nbody runs confirm the results of hydrodynamical simulations and show that systems which are formed at convergence zones generally reach a quasistationary state with each body in resonance with its neighbours and evolving on a nonmigrating 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 Nbody 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 Nbody 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 ValBorro 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: (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 . In the previous equation, is the radiative cooling term which is given by: (2)where σ_{B} is the StephanBoltzmann constant and T_{eff} is the effective temperature which is related to the central temperature T by: (3)where τ_{eff} is the effective optical depth given by (Hubeny 1990): (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 halfwidth 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_{0} = 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 wavekilling 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 fifthorder RungeKutta 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 Nbody timestep given by (Beaugé & Aarseth 1990): (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: (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. and are dimensionless coefficients which are set to and (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 magnetorotational 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 wavelike modes and given by: (7)with (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 , where Δt_{k} = 0.2πR_{k}/m_{k}c_{s} denotes the lifetime of mode with wavenumber m_{k}. Such a value for Δt_{k} yields an autocorrelation timescale τ_{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^{2}) 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): (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 . 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 nonisothermal 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): (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). 

Open with DEXTER 
2.3. Initial conditions
The initial disc surface density profile is Σ = Σ_{p} × (R/R_{0})^{− σ} with σ = 0.5 and Σ_{p} = 850 g/cm^{2}.
The anomalous viscosity in the disc arising from MHD turbulence is modelled using a constant kinematic viscosity ν = 10^{14} g/cm^{2}, 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 zerotorque line remains approximately fixed in the course of the simulations.
The initial temperature profile is T = T_{0}(R/R_{0})^{− β} with β = 1 and T_{0} 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 steadystate 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 snowline is located at ~2.7 AU, which is consistent with estimations of the location of the snowline 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 steadystate. 

Open with DEXTER 
In the hydrodynamical simulations presented below, the initial mass of each protoplanet is assumed to be m_{p} = 3M_{⊕}. The motivation for choosing equalmass 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 equalmass bodies occur in the radiative case but not in the locally isothermal limit, which clearly demonstrates the role of the zeromigration 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 semimajor 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: (11)where a_{p,p′} denotes the semimajoraxes 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 Nbody 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 Nbody simulations
3.1. Prescription for Type I migration
The torque exerted on a protoplanet embedded in a nonisothermal 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 entropyrelated 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 timescales and the horseshoe libration timescale. Its optimal value, also referred to as the fully unsaturated corotation torque, is obtained when the diffusion timescales are approximately equal to half the horseshoe libration time (e.g. Baruteau & Masset 2013). In the limit where the diffusion timescales become shorter than the Uturn timescale, 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 timescales and the horseshoe libration timescale. Torque formulae as a function of viscosity and thermal diffusivity were proposed by Paardekooper et al. (2011).
In our Nbody runs, Type I migration is modelled by an extraforce a_{m} acting on each body and defined by: (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 m_{p} = 6 M_{⊕}. 

Open with DEXTER 
For a purely active disc, radiative diffusion and viscous timescales 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): (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 zerotorque 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 zeromigration line, the entropy gradient is weaker (S ∝ R^{0.14}) so that the (positive) entropyrelated corotation torque is not strong enough to counterbalance the (negative) differential Lindblad torque. We note that in Fig. 3, the zeromigration line at R ~ 1.7 AU exists provided that the corotation torque remains unsaturated. This arises when the diffusion timescale across the horseshoe region is shorter than the libration timescale t_{lib} = 8πa_{p}/3Ω_{p}x_{s} but longer than the Uturn timescale t_{U − turn} ~ h_{p}t_{lib}, where x_{s} is the halfwidth of the horseshoe region which is given by (Paardekooper et al. 2010): (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: (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 lowmass 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 halfwidth 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 Nbody 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 halfwidth 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: (16)In the previous equation, is defined by , where is the dimensionless halfwidth 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 .
Fig. 4 Corotation torque damping factor as a function of the ratio between the eccentricity and the dimensionless halfwidth of the horseshoe region for m_{p} = 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 dashdotted line corresponds to the prescription of Cossou et al. (2013). 

Open with DEXTER 
3.2. Eccentricity and inclination damping in Nbody simulations
Fig. 5 Orbital evolution of N = 5 embryos with mass m_{p} = 3 M_{⊕} and without stochastic forces included. The dashed line shows the location of the zeromigration line. 

Open with DEXTER 
Eccentricity and inclination damping resulting from the interaction with the disc are modelled by applying the following accelerations to each body: (17)and (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): (19)and (20)where: (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^{2} 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 halfwidth of the planet’s horseshoe region which is estimated to be . 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 outwardmigrating 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 threeplanet 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 outwardmigrating 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 zerotorque 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. 

Open with DEXTER 
Fig. 7 Left panel: orbital evolution of N = 7 embryos with mass m_{p} = 3 M_{⊕} and without stochastic forces included. Right panel: same but with N = 9. The dashed line shows the location of the zeromigration line. 

Open with DEXTER 
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 inwardmigrating bodies are captured in resonance later than outwardmigrating 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, outwardmigrating 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^{4} 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 fullyunsaturated corotation torque scales as (or equivalently as ) for lowmass 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. 

Open with DEXTER 
Not surprisingly, the significant resulting differential migration between the outermost 6M_{⊕} planet (blue) and its exterior inwardmigrating 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^{4} 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 longterm evolution outcomes using Nbody simulations.
In order to unveil the role of the zeromigration 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’ semimajoraxes and eccentricities as a function of time for the isothermal run. We note that although we consider here equalmass 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 semimajoraxes and eccentricities of protoplanets for an isothermal hydrodynamical simulation with N = 9. 

Open with DEXTER 
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 zeromigration line. 

Open with DEXTER 
Fig. 11 Top panel: time evolution of the semimajoraxes 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 zeromigration line. 

Open with DEXTER 
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’ semimajor 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 zeromigration radius, surface density perturbations inside the planets’ horseshoe regions and related to the coorbital dynamics are clearly visible. Due to the negative entropy gradient, coorbital 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 zeromigration 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 firstorder 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 snowline 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 zeromigration 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 zeromigration line, and which would promote close encounters between embryos, is couterbalanced by a stronger discinduced 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’ semimajoraxes for N = 5. Here, a 6 M_{⊕} planet (cyan) is formed at the outer edge of the swarm at t ~ 5 × 10^{3} yr and is likely to become trapped at the location of the zerotorque radius at later times. From t ~ 1.3 × 10^{4} 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^{3} 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 threeplanet system will cancel out (Cossou et al. 2013).
5. Results of Nbody runs
Fig. 12 Time evolution of the semimajoraxes 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 zeromigration line. 

Open with DEXTER 
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 Nbody 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 longterm 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 Nbody run with intially N = 9 embryos and without stochastic forces included. Right: same but with stochastic forces included. 

Open with DEXTER 
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 Nbody 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 zeromigration line by ~10^{4} years. It is interesting to notice that of these two 6 M_{⊕}, one is formed at ~8 × 10^{3} years through the merging of two coorbital planets that entered in a 1:1 resonance at ~ 2 × 10^{3} years. An additional 6 M_{⊕} body is produced at t ~ 1.8 × 10^{4} yr as a result of the collision between two 3 M_{⊕} objects. After ~ 3 × 10^{4} years, the systems attains a quasistationary state and consists of three 6 M_{⊕} objects plus two 3 M_{⊕} bodies trapped in a resonant chain and evolving on nonmigrating 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 Nbody 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 Nbody 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_{⊕}. 

Open with DEXTER 
The right panel of Fig. 13 illustrates the evolutionary outcome for a simulation in which two 12 M_{⊕} planets are produced. At t ~ 10^{5} yr, a nonmigrating compact system is formed with two coorbital 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^{5} 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 quasicircular 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^{5} yr. We note that although these planets evolve inside the zeromigration 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^{5} yr. This is in good agreement with the results of simulations with initially equalmass 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 zerotorque 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: (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)): (23)Given that τ = κΣ/2, Eq. (22) can be recast as: (24)Moreover, in the context of the standard α prescription of Shakura & Sunyaev (1973), the effective kinematic viscosity is given by ν = αc_{s}H, 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}: (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 Σ = Σ_{0}(R/R_{0})^{− σ} and Ω = Ω_{0}(R/R_{0})^{− 3/2}, where Ω_{0} is the angular velocity at R = R_{0}.
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): (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: (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 soliddashed 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 zerotorque radius in our hydrodynamical and Nbody 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 t_{dif} becomes equal to half of the horseshoe libration timescale (t_{lib}) whereas the soliddashed line shows the radius where t_{dif} equals the Uturn 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. 

Open with DEXTER 
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 JupiterSaturn 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 Earthmass 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 Earthmass 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 Nbody simulations of the evolution of a swarm of Earthmass protoplanets that gravitationally interact near a zeromigration 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 lowmass protoplanets convergently migrate toward a zeromigration 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 equalmass 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^{4} 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 Nbody 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^{5} yr. For a randomised mass distribution with μ_{m} = 1 M_{⊕}, however, only ~ 15% of the Nbody 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 cmsized pebbles (Morbidelli & Nesvorny 2012).
For a protoplanetary disc with mass typical of the MMSN, we find that the mechanism presented here may allow the insitu formation of giant planet cores from Earthmass 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 zeromigration 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 Nbody 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 zeromigration 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 Nbody simulations that account for the inward migration of the zerotorque 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 GENCIcines (c2012046957).
References
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Lin, D. N. C. 2010, ApJ, 709, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Masset, F. 2008, ApJ, 672, 1054 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Masset, F. 2013, Lect. Notes Phys., 861, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Beauge, C., & Aarseth, S. J. 1990, MNRAS, 245, 30 [NASA ADS] [Google Scholar]
 Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., & Kley, W. 2011, A&A, 536, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & DobbsDixon, I. 2013, A&A, 549, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Broeg, C. H., & Benz, W. 2012, A&A, 538, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canup, R. M., & Asphaug, E. 2001, Nature, 412, 708 [NASA ADS] [CrossRef] [Google Scholar]
 Cossou, C., Raymond, S. N., & Pierens, A. 2013, A&A, 553, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2008, A&A, 482, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Garaud, P., & Lin, D. N. C. 2007, ApJ, 654, 606 [NASA ADS] [CrossRef] [Google Scholar]
 Gladman, B. 1993, Icarus, 106, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Greenberg, R., Hartmann, W. K., Chapman, C. R., & Wacker, J. F. 1978, Icarus, 35, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hasegawa, Y., & Pudritz, R. E. 2011, MNRAS, 417, 1236 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Fundamental Problems in the Theory of Stellar Evolution (Dordrecht: R. Reidel), IAU Symp., 93, 113 [NASA ADS] [Google Scholar]
 Hellary, P., & Nelson, R. P. 2012, MNRAS, 419, 2737 [NASA ADS] [CrossRef] [Google Scholar]
 Horn, B., Lyra, W., Mac Low, M.M., & Sándor, Z. 2012, ApJ, 750, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I. 1990, ApJ, 351, 632 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475 [NASA ADS] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Kokubo, E., & Ida, S. 2000, Icarus, 143, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Kretke, K. A., & Lin, D. N. C. 2012, ApJ, 755, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489 [NASA ADS] [CrossRef] [Google Scholar]
 Lecar, M., Podolak, M., Sasselov, D., & Chiang, E. 2006, ApJ, 640, 1115 [NASA ADS] [CrossRef] [Google Scholar]
 Leinhardt, Z. M., & Richardson, D. C. 2005, ApJ, 625, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Thommes, E., & Duncan, M. J. 2010, AJ, 139, 1297 [NASA ADS] [CrossRef] [Google Scholar]
 Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 491, L41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lyra, W., Paardekooper, S.J., & Mac Low, M.M. 2010, ApJ, 715, L68 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icarus, 219, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Ogihara, M., Ida, S., & Morbidelli, A. 2007, Icarus, 188, 522 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper, S.J., & Papaloizou, J. C. B. 2008, A&A, 485, 877 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paardekooper, S.J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper, S.J., Baruteau, C., & Kley, W. 2012, MNRAS, 410, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2008, A&A, 478, 939 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., Baruteau, C., & Hersant, F. 2011, A&A, 531, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., Baruteau, C., & Hersant, F. 2012, MNRAS, 427, 1562 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, 2nd edn. (Cambridge: University Press), [Google Scholar]
 Sándor, Z., Lyra, W., & Dullemond, C. P. 2011, ApJ, 728, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] [Google Scholar]
 Thommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icarus, 161, 431 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 van Leer, B. 1977, J. Comput. Phys., 23, 276 [NASA ADS] [CrossRef] [Google Scholar]
 Ward, W. R. 1997, Icarus, 126, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Wetherill, G. W. 1985, Science, 228, 877 [NASA ADS] [CrossRef] [Google Scholar]
 Wetherill, G. W., & Stewart, G. R. 1989, Icarus, 77, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., Mac Low, M.M., & Menou, K. 2009, ApJ, 707, 1233 [NASA ADS] [CrossRef] [Google Scholar]
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). 

Open with DEXTER  
In the text 
Fig. 2 Surface density and temperature profiles at steadystate. 

Open with DEXTER  
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 m_{p} = 6 M_{⊕}. 

Open with DEXTER  
In the text 
Fig. 4 Corotation torque damping factor as a function of the ratio between the eccentricity and the dimensionless halfwidth of the horseshoe region for m_{p} = 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 dashdotted line corresponds to the prescription of Cossou et al. (2013). 

Open with DEXTER  
In the text 
Fig. 5 Orbital evolution of N = 5 embryos with mass m_{p} = 3 M_{⊕} and without stochastic forces included. The dashed line shows the location of the zeromigration line. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 7 Left panel: orbital evolution of N = 7 embryos with mass m_{p} = 3 M_{⊕} and without stochastic forces included. Right panel: same but with N = 9. The dashed line shows the location of the zeromigration line. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 9 Time evolution of the semimajoraxes and eccentricities of protoplanets for an isothermal hydrodynamical simulation with N = 9. 

Open with DEXTER  
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 zeromigration line. 

Open with DEXTER  
In the text 
Fig. 11 Top panel: time evolution of the semimajoraxes 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 zeromigration line. 

Open with DEXTER  
In the text 
Fig. 12 Time evolution of the semimajoraxes 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 zeromigration line. 

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

Open with DEXTER  
In the text 
Fig. 14 Left panel: distribution of mass of the most massive planet which is formed in Nbody 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_{⊕}. 

Open with DEXTER  
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 t_{dif} becomes equal to half of the horseshoe libration timescale (t_{lib}) whereas the soliddashed line shows the radius where t_{dif} equals the Uturn 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. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.