Formation of solar system analogs I: looking for initial conditions through a population synthesis analysis

TL;DR

This study develops a comprehensive planetary formation model to identify initial conditions that lead to solar system analogs, emphasizing the importance of disc mass, migration rates, and planetesimal size in forming habitable, water-rich planets.

Contribution

It introduces a detailed planet formation simulation focusing on initial conditions for solar system analogs, integrating physical processes like photoevaporation, migration, and water delivery.

Findings

Solar system analogs favor massive discs with low migration rates.

Small planetesimals increase the likelihood of forming analogs.

Habitable zone planets tend to be dry after disc dissipation.

Abstract

Population synthesis models of planetary systems developed during the last $\sim$15 years could reproduce several of the observables of the exoplanet population, and also allowed to constrain planetary formation models. We present our planet formation model, which calculates the evolution of a planetary system during the gaseous phase. The code incorporates relevant physical phenomena for the formation of a planetary system, like photoevaporation, planet migration, gas accretion, water delivery in embryos and planetesimals, a detailed study of the orbital evolution of the planetesimal population, and the treatment of the fusion between embryos, considering their atmospheres. The main goal of this work, unlike other works of planetary population synthesis, is to find suitable scenarios and physical parameters of the disc to form solar system analogs. We are specially interested in the final planet distributions, and in the final surface density, eccentricity and inclination profiles for the planetesimal population. These final distributions will be used as initial conditions for N-body simulations, to study the post-oligarchic formation in a second work. We then consider different formation scenarios, with different planetesimal sizes and different type I migration rates. We find that solar system analogs are favored in massive discs, with low type I migration rates, and small planetesimal sizes. Besides, those rocky planets within their habitables zones are dry when discs dissipate. At last, the final configurations of solar system analogs include information about the mass and semimajor-axis of the planets, water contents, and the properties of the planetesimal remnants.

Tables

Table 1. Table 1: Final percentages of planetary systems with respect to the total of the developed simulations.

Planetary system	Final percentage
Rocky planetary systems	64.46%
Icy giant planetary systems	13.4%
Gas giant planetary systems	22.14%
Solar system analogs	4.3%

Table 2. Table 2: Percentage of SSA found in each formation scenario.

Planetesimal size	$f_{\text{migI}}=0$	$f_{\text{migI}}=0.01$	$f_{\text{migI}}=0.1$	$f_{\text{migI}}=1$
100 m	26.90%	14.90%	1.10%	0%
1 km	16.40%	1.4%	0%	0%
10 km	5.10%	2.40%	0.3%	0%
100 km	0.2%	0.1%	0%	0%

Table 3. Table 3: Ranges, means and dispersion values of the disc parameters that formed SSA.

Disc	Range	Mean	$\sigma$
parameter
$M_{\text{d}}$	$0.04{\text{M}}_{\odot }$ - $1.15{\text{M}}_{\odot }$	$0.10{\text{M}}_{\odot }$	$0.027{\text{M}}_{\odot }$
$\gamma$	$0.5$ - $1.5$	0.95	0.27
$R_{\text{c}}$	20 au - 50 au	32.8 au	8.73 au
$\alpha$	${10}^{-4}$ - ${10}^{-2}$	$3.4\times {10}^{-3}$	$2.8\times {10}^{-3}$
$\tau$	1.65 Myr - 11.96 Myr	6.45 Myr	2.67 Myr

Equations

$${\Sigma }_{\text{g}}(R)={\Sigma }_{\text{g}}^0{\left(\frac{R}{R_{\text{c}}}\right)}^{-\gamma }{e}^{-{\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }},$$

$${\Sigma }_{\text{p}}(R)={\Sigma }_{\text{p}}^0{\eta }_{\text{ice}}{\left(\frac{R}{R_{\text{c}}}\right)}^{-\gamma }{e}^{-{\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }},$$

$${\eta }_{\text{ice}}=\{\begin{array}{ll}1& \text{ if R≥Rice},\\ & \\ \frac{1}{\beta }& \text{ if R<Rice},\end{array}$$

$$\frac{\partial {\Sigma }_{\text{g}}}{\partial t}=\frac{3}{R}\frac{\partial }{\partial R}\left[R^{1/2}\frac{\partial }{\partial R}\left(\nu {\Sigma }_{\text{g}}{R}^{1/2}\right)\right]+{\dot{\Sigma }}_{\text{w}}(R),$$

$$c_{\text{s}}=\sqrt{\frac{{\gamma }_{\text{g}}{k}_{\text{B}}\text{T}}{{\mu }_{{\text{H}}_2}{m}_{{\text{H}}_2}}},$$

$$\text{T}=280{\left(\frac{R}{1\text{au}}\right)}^{-1/2}\text{K},$$

$${\dot{\Sigma }}_{\text{w}}(R)=\{\begin{array}{ll}{\dot{\Sigma }}_{\text{w}}^g\exp \left[\frac{1}{2}\left(1-\frac{R_{\text{g}}}{R}\right)\right]{\left(\frac{R_{\text{g}}}{R}\right)}^2& \text{ if R≤Rg},\\ & \\ {\dot{\Sigma }}_{\text{w}}^g{\left(\frac{R_{\text{g}}}{R}\right)}^{5/2}& \text{ if R>Rg},\end{array}$$

$${\dot{\Sigma }}_{\text{w}}^g=1.16\times 10^{-11}{f}_{41}^{0.5}{\left(\frac{1\text{au}}{R_{\text{g}}}\right)}^{1.5}\frac{{\text{M}}_{\odot }}{{\text{au}}^2\text{yr}},$$

$$\frac{\partial {\Sigma }_{\text{p}}(R)}{\partial t}-\frac{1}{R}\frac{\partial }{\partial R}\left[R{v}_{\text{mig}}(R){\Sigma }_{\text{p}}(R)\right]={\dot{\Sigma }}_{\text{p}}^{\text{tot}}(R),$$

$${\dot{\Sigma }}_{\text{p}}^{\text{tot}}(R)={\dot{\Sigma }}_{\text{p}}^{\text{ac}}(R)+{\dot{\Sigma }}_{\text{p}}^{\text{sc}}(R)+{\dot{\Sigma }}_{\text{p}}^{\text{ej}}(R),$$

$${\dot{\Sigma }}_{\text{p}}^{\text{sc}}(R)=\frac{1}{4}{\left(\frac{a_{\text{P}}{M}_{\text{P}}}{M_{\star }{\tilde{R}}_{\text{C}}}\right)}^2{\dot{\Sigma }}_{\text{p}}^{\text{ac}}(R),$$

$$\frac{\text{d}{M}_{\text{C}}}{\text{d}t}=\frac{2\pi {\Sigma }_{\text{p}}(a_{\text{P}})R_{\text{H}}^2}{P}{P}_{\text{coll}}.$$

$${\tilde{R}}_{\text{C}}=\{\begin{array}{l}\frac{R_{\text{B}}}{1+\frac{2W({\sigma }_{\text{crit}}-1)+\ln ({\sigma }_{\text{crit}})}{{\gamma }_{\text{g}}}},\text{if}{\sigma }_{\text{crit}}<{\sigma }_t,\\ \\ \frac{R_{\text{B}}}{\frac{R_{\text{B}}}{R_t}+\frac{4}{{\gamma }_{\text{g}}}(4W{)}^{1/3}({\sigma }_{\text{crit}}^{1/3}-{\sigma }_t^{1/3})},\text{if}{\sigma }_{\text{crit}}>{\sigma }_t,\end{array}$$

$$\frac{R_{\text{B}}}{R_t}$$

$$W$$

$${\sigma }_t$$

$${\sigma }_{\text{crit}}$$

$$L$$

$$R_{\text{B}}$$

$$M_{\text{crit}}\sim 10{\left(\frac{{\dot{M}}_{\text{C}}}{10^{-6}{\text{M}}_{\oplus }{\text{yr}}^{-1}}\right)}^{0.25}{\text{M}}_{\oplus },$$

$${\dot{M}}_{\text{KH}}=\frac{\text{d}{M}_{\text{g}}}{\text{d}t}=\frac{M_{\text{P}}}{{\tau }_{\text{g}}}$$

$${\tau }_{\text{g}}=8.35\times 10^{10}(M_{\text{P}}/{\text{M}}_{\oplus }{)}^{-3.65}\text{yr}.$$

$${\dot{M}}_{\text{disc}}=\frac{\text{d}{M}_{\text{g}}}{\text{d}t}=3\pi \nu {\Sigma }_{\text{g}}(R).$$

$$\frac{\text{d}{M}_{\text{g}}}{\text{d}t}=\text{min}[{\dot{M}}_{\text{KH}},{\dot{M}}_{\text{disc}}].$$

$${\dot{M}}_{\text{GP}}=\frac{\text{d}{M}_{\text{g}}}{\text{d}t}=\dot{A}{\Sigma }_{\text{acc}},$$

$$\dot{A}=0.29{\left(\frac{H_{\text{g,P}}}{a_{\text{P}}}\right)}^{-2}{\left(\frac{M_{\text{P}}}{M_{\star }}\right)}^{4/3}{a}_{\text{P}}^2{\Omega }_{\text{k,P}}$$

$$\Sigma(x)=\begin{cases}\Sigma_{\text{g}}\text{exp}\left[-\left(\dfrac{x}{l}\right)^{-3}\right]&\text{ if $x>x_{\text{m}}$},\\ \\ \!\begin{aligned} &\Sigma_{\text{g}}\text{exp}\Bigg{[}-\dfrac{1}{2}\left(\dfrac{x}{H_{\text{g,P}}}-\dfrac{5}{4}\dfrac{x_{\text{m}}}{H_{\text{g,P}}}\right)^{2}\\ \\ &+\dfrac{1}{32}\left(\dfrac{x_{\text{m}}}{H_{\text{g,P}}}\right)^{2}-\left(\dfrac{x_{\text{m}}}{l}\right)^{-3}\Bigg{]}\\ \end{aligned}&\text{if $x\leq x_{\text{m}}$},\end{cases}$$

$$l=0.146{\left(\frac{\nu }{10^{-5}{a}_{\text{P}}^2{\Omega }_{\text{k,P}}}\right)}^{-1/3}{\left(\frac{M_{\text{P}}}{10^{-3}{M}_{\star }}\right)}^{2/3}{a}_{\text{P}},$$

$$x_{\text{m}}=0.207{\left(\frac{H_{\text{g}}}{0.1a_{\text{P}}}\frac{M_{\text{P}}}{10^{-3}{M}_{\star }}\right)}^{2/5}{\left(\frac{\nu }{10^{-5}{a}_{\text{P}}^2{\Omega }_{\text{k,P}}}\right)}^{-1/5}{a}_{\text{P}},$$

$$\frac{\text{d}{M}_{\text{g}}}{\text{d}t}=\text{min}[{\dot{M}}_{\text{KH}},{\dot{M}}_{\text{disc}},{\dot{M}}_{\text{GP}}].$$

$$M_{\text{C}}=M_{\text{C}}^i+M_{\text{C}}^j,$$

$$M_{\text{E}}=M_{\text{E}}^i(1-{\chi }_{\text{loss}}^i)+M_{\text{E}}^j(1-{\chi }_{\text{loss}}^j),$$

$$a={\left(\frac{a_i^{0.5}(M_{\text{C}}^i+M_{\text{E}}^i)+a_j^{0.5}(M_{\text{C}}^j+M_{\text{E}}^j)}{M_{\text{C}}+M_{\text{E}}}\right)}^2,$$

Full text

††pubyear: 2017††pagerange: Formation of solar system analogs I: looking for initial conditions through a population synthesis analysis– Formation of solar system analogs I: looking for initial conditions through a population synthesis analysis

1 Introduction

For a long time, the study of planetary systems was restricted to our own solar system. However, this situation had a drastic change since the discovery of the first exoplanet orbiting a sun-like star, 51 Peg b, in 1995 (Mayor & Queloz, 1995). 51 Peg b, named as Dimidium by the International Astronomical Union in December 2015, is what we know now as a hot-Jupiter, a gas giant planet orbiting very close to its parent star. Since this discovery, the observational field of extrasolar planet search had an explosion, leading to numerous additional discoveries of planets orbiting other stars. Up to date, 3610 exoplanets that show a wide range of masses, orbits and compositional properties have been discovered orbiting near stars (http://exoplanet.eu/). Some of them are hot and warm Jupiters, Jupiter-analogs, giant planets at wide orbits, hot-Neptunes, super-Earths, etc. Moreover, much of these planets form part of about 610 multiple planet systems. These systems represent the final stage of a series of complex processes, where a protoplanetary disc evolves into a planetary system, with a few planets and probably, reservoirs of small bodies, like in our solar system. These discoveries also triggered theoretical studies about the formation and evolution of planetary systems. During more than a decade, several models of planet formation have been developed to study the formation of planetary systems and population synthesis with the aim to reproduce the main properties of the observational sample of exoplanets, and to get a better understanding of the main processes of planetary formation (see Benz et al., 2014, for a detailed review).

The first models of planetary population synthesis were developed in the pioneer works of Ida & Lin (2004a, b, 2005, 2008a, 2008b). In these works, the authors could numerically reproduce several of the main observed properties of the populations of exoplanets, specially the mass versus semimajor-axis diagram (M-a diagram). Thommes et al. (2008) studied the formation of planetary systems with an hybrid model, including N-body interaction between embryos, with the aim to link a mature planetary system with the properties of the protoplanetary disc where it was born. They found that solar system analogs are not common. These systems can form when the timescale for the formation of the giant planets is similar to the timescale of the disc dissipation, undergoing only modest type I migration, and for massive discs. Otherwise, Miguel et al. (2011) found that solar system analogs might not be so rare in the solar neighbourhood, where the formation of such systems could be favoured by massive discs, and where there is not a large accumulation of solids in the inner region. In other case, very slow type I migration is needed to form giant planets. However, Miguel et al. (2011) did not analyze the post-oligarchic growth and long-term evolution of the planetary systems and adopted a different definition for solar system analogs. More recently, Alibert et al. (2013) introduced a distribution of embryos calculating the gravitational interactions between them in their planetary population synthesis model (Alibert et al., 2005; Mordasini et al., 2009; Fortier et al., 2013), which previously calculates planetary population synthesis but considering only one planet per disc. These authors also showed that the main observed properties of the exoplanets M-a diagram can be reproduced, focusing in low-mass planets. Then, Pfyffer et al. (2015) used the results of Alibert et al. (2013) as initial conditions to calculate the long-term evolution of the planetary population synthesis but only considering the distribution of embryos and without including the remnant of planetesimal distributions. They found that the M-a diagram does not experience significantly changes due to the long-term evolution of the systems, but that the eccentricity distributions of such systems do not match the observed ones in the known exoplanets. Moreover, Ida et al. (2013) also introduced gravitational interactions between embryos in their planetary population synthesis model using a Monte Carlo approach (not calculating N-body integrations) founding that the observed planetary mass-eccentricity and M-a diagrams can be reproduced. However, in most of the cases, the evolution of the planetesimal population is treated in a simplified way and is not taken into account in the calculations of the long-term evolution of the systems.

On the other hand, several works study the formation of terrestrial planets in the post-oligarchic growth regime and the water delivery via purely N-body calculations considering a distribution of embryos and planetesimals after the dissipation of the primordial nebula (Chambers, 2001; Raymond et al., 2004, 2005, 2006; Raymond et al., 2009; O’Brien et al., 2006, 2014; Walsh et al., 2011; Ronco & de Elía, 2014). In general, the orbital initial conditions for the distributions of embryos and planetesimals are typically selected ad hoc. However, Ronco et al. (2015) showed that more realistic initial conditions lead to different accretion histories of the planets that remain in the habitable zone compared to the case of ad hoc initial conditions.

The aim of this first work is to improve our model of planet formation during the gaseous phase to obtain accurate initial conditions for the distributions of planets and planetesimals, to perform future high resolution N-body calculations for the analysis of the formation of terrestrial planets and water delivery in solar system analogs. To do so, we first include the most relevant physical phenomena that take place in the formation of a planetary system during the gaseous phase. Then, we perform a planetary population synthesis calculation in order to obtain the main properties of the protoplanetary discs that lead to the formation of solar system analogs. This work is then organized as follow: in Sec. 2 we describe our planet formation model; then, in Sec. 3 we describe the initial conditions for the development of our planetary systems, in Sec. 4 we show the results, and in Sec. 5, we discuss the conclusions.

2 Description of our model of planet formation

We have improved a model of planet formation based on previous works (Guilera et al., 2010, 2011, 2014), which we now called PLANETALP. This code describes the evolution in time of a planetary system during the gaseous phase and incorporates many of the most important physical phenomena for the formation of planetary systems. PLANETALP is a code developed almost entirely in FORTRAN 95/2003 and has the advantage of being programmed in a modular way so as to be able to turn on or turn off the different physical phenomena according to the user’s choice.

Our model presents a protoplanetary disc, which is characterized by two components: a gaseous component, evolving due to an $\alpha$-viscosity driven accretion and photoevaporation processes, and a solid component represented by a planetesimal population being subject to accretion and ejection by the embryos, and radial drift due to gas drag. Besides, the model presents an embryo population which is also part of the solid component of the disc and which grows by accretion of planetesimals, gas, and due to the fusion between them taking into account their atmospheres. The new improvements made in PLANETALP with respect to previous versions of this code (Guilera et al., 2010, 2011, 2014) are described in detail in the next subsections.

2.1 Protoplanetary disc structure

The gaseous component is characterized by the gas surface density $\Sigma_{\text{g}}(R)$, where $R$ is the radial coordinate in the midplane of the disc, and the planetesimal population is characterized by the planetesimal surface density $\Sigma_{\text{p}}(R)$. Although our model allows to include a discrete planetesimal size distribution described by $\Sigma_{\text{p}}(R,r_{\text{p}})$, where $r_{\text{p}}$ represents the different sizes of the planetesimals, we will consider in this work only one size of planetesimals for each simulation. Both $\Sigma_{\text{g}}(R)$ and $\Sigma_{\text{p}}(R)$ determine the mass distribution along the protoplanetary disc. In this work, as in previous ones (de Elía et al., 2013; Ronco et al., 2015; Dugaro et al., 2016) we consider the surface density profile suggested by Andrews et al. (2010), who found that the discs observed in the Ophiuchus star-forming region can be represented by

[TABLE]

where $R_{\text{c}}$ is a characteristic radius of the disc, $\gamma$ is the exponent that represents the surface density gradient and $\Sigma_{\text{g}}^{0}$ is a normalization constant, which is a function of the mass of the disc $M_{\text{d}}$, $R_{\text{c}}$ and $\gamma$, that can be obtained by integrating Eq. 1 over the total area of the disc. Assuming that the metallicity along the disc is the same as that of the central star, and considering that the dust sediments and coagulates very quickly to form a mid-plane planetesimal disc, is that the planetesimal surface density profile is given by

[TABLE]

where $\Sigma^{0}_{\text{p}}=z_{0}\Sigma^{0}_{\text{g}}10^{[\text{Fe/H}]}$ is a normalization constant, $z_{0}=0.0153$ is the primordial abundance of heavy elements in the Sun (Lodders et al., 2009), $[\text{Fe/H}]$ is the stellar metallicity, and $\eta_{\text{ice}}$ is a parameter that represents an increase in the amount of solid material due to the condensation of water beyond the snow line, given by

[TABLE]

where $R_{\text{ice}}$ represents the position in the disc where water condensates, usually at $170^{\circ}$ K, and the factor $\beta$ represents the amount of material condensed. In our solar system, $\beta$ usually takes a value of approximately 2 (Lodders, 2003; Lodders et al., 2009).

The values of the free parameters involved in the structure of the protoplanetary disc, like the mass of the disc $M_{\text{d}}$, the characteristic radius $R_{\text{c}}$, the surface density gradient $\gamma$, and the factor $\beta$ will be discussed in the next section.

2.1.1 Evolution of the gaseous component

As it was mentioned before, we improved the evolution of the gaseous component considering a viscous accretion disc (Pringle, 1981) with photoevaporation.

The evolution in time of the gas surface density is then represented by a diffusion equation

[TABLE]

where $\nu=\alpha c_{\text{s}}H_{\text{g}}$ represents the Shakura and Sunyaev viscosity (Shakura & Sunyaev, 1973) with $\alpha$ a constant a-dimensional parameter along the disc, $c_{\text{s}}$ the sound speed, $H_{\text{g}}$ the height scale of the disc, and $\dot{\Sigma}_{\text{w}}(R)$ the photoevaporation sink term. The sound speed is given by

[TABLE]

where $\gamma_{\text{g}}=5/3$, $k_{\text{B}}$ is the Boltzmann constant, and $\mu_{\text{H}_{2}}$ and $m_{\text{H}_{2}}$ are the molecular weight and mass of molecular hydrogen, respectively. Following Hayashi (1981) and Ida & Lin (2004a), we consider a temperature profile given by

[TABLE]

where, according to this, $R_{\text{ice}}$ is defined as 2.7 au. Finally, the height scale of the disc is $H_{\text{g}}=\sqrt{2}c_{\text{s}}/\Omega_{\text{k}}$ where $\Omega_{\text{k}}$ is the keplerian frequency. It is important to note that we are considering a flare disc, where $H_{\text{g}}\propto R^{5/4}$.

Finally, $\dot{\Sigma}_{\text{w}}(R)$ represents the photoevaporation rate. There are different photoevaporation regimes, and they depend on which is the source that emits irradiation. This source could be the central star (Dullemond et al., 2007; D’Angelo & Marzari, 2012), or external stars that irradiate the environment in where the disc forms (Veras & Armitage, 2004). If our sun was formed in a densely populated environment, the high stellar density could lead to encounters with stars that could affect and modify the size of the disc by photoevaporation (see Moro-Martín et al., 2008, and references therein). However, since up to date the environment in which our solar system was formed is still under discussion, and due to the fact that photoevaporation by the central star is, probably, the most important photoevaporation regime that affects the planet formation region, we only consider this case. Thus, following the model of Dullemond et al. (2007), also adopted by D’Angelo & Marzari (2012), $\dot{\Sigma}_{\text{w}}(R)$ is given by

[TABLE]

where $R_{\text{g}}$ is the radius of the gas surface layer at which the sound speed equals the local keplerian velocity, and beyond which the gas is disengaged from the disc. The photoevaporation rate at $R_{\text{g}}$ is represented by $\dot{\Sigma}^{g}_{\text{w}}$ and is given by

[TABLE]

where $f_{41}$ is the rate of EUV ionizing photons emitted by the star in units of $10^{41}$s*-1*. It is worth noting that $R_{\text{g}}$ is usually taken to be 10 au for a solar-mass star.

Equation 4 is solved on a disc defined between 0.1 au and 1000 au, using 5000 radial bins logarithmically equally spaced. It is worth remarking that we consider that the disc is dissipated when the mass of gas becomes lower than $10^{-6}~{}M_{\odot}$.

2.1.2 Evolution of the planetesimal population

In our model, the planetesimal population evolves by the drift of planetesimals due to the nebular gas including the Epstein, Stokes and quadratic regimes, and the accretion and ejection by the embryos. We also consider that the evolution of the eccentricities and inclinations of the planetesimals are due to two principal processes: the embryo gravitational excitation (Ohtsuki et al., 2002), and the damping due to the nebular gas drag (Rafikov, 2004; Chambers, 2008). The treatment of these processes is described in detail in Guilera et al. (2014). Then, the evolution in time of the planetesimal surface density is modeled with a continuity equation given by

[TABLE]

[TABLE]

where $v_{\text{mig}}$ is the planetesimal migration velocity and $\dot{\Sigma}^{\text{ac}}_{\text{p}}(R)$, $\dot{\Sigma}^{\text{sc}}_{\text{p}}(R)$, and $\dot{\Sigma}^{\text{ej}}_{\text{p}}(R)$ represent the sink terms due to the accretion by the embryos, due to the probability of scattered planetesimals, and due to the ejection of planetesimals with high eccentricities, respectively. In this work, we incorporate two different ejection rates of planetesimals that play different roles. It is important to remark that the ejected planetesimals are directly removed from the planetary system and are not relocated in a different orbit. The first one, $\dot{\Sigma}^{\text{sc}}_{\text{p}}(R)$, considers the ejection rate of planetesimals proposed by Ida & Lin (2004a) and also adopted by Alibert et al. (2005), that takes into account the probability of planetesimal scattering rather than collisions during close encounters, and is

[TABLE]

where $a_{\text{P}}$ and $M_{\text{P}}$ are the semimajor-axis and the total mass of the planet, respectively, and $\tilde{R}_{\text{C}}$ is the planet’s enhanced radius (see next section). The second one, $\dot{\Sigma}^{\text{ej}}_{\text{p}}(R)$, considers the ejection of planetesimals that reach values of the eccentricity higher than $0.99$. As the gas disc evolves, the gas surface density decreases, and because of this, the eccentricities of the planetesimals that are in the vicinity of the massive planets are no longer efficiently damped by the gas, and therefore, they increase considerably. At this time, $\dot{\Sigma}^{\text{ej}}_{\text{p}}(R)$ becomes effective. As in Laakso et al. (2006), we consider that the probability of a particle to be ejected per unit time is constant. This implies that $\dot{\Sigma}^{\text{ej}}_{\text{p}}=-\rho\Sigma_{\text{p}}$, where we adopt ad hoc values of $\rho$ between 0.001 yr*-1* and 0.01 yr*-1*. These values represent the ejection of the 0.1% and 1%, respectively, of the available planetesimals in each radial bin around the planet at each time-step, finding no differences in the final results.

We use a full implicit upwind-downwind mix method to solve Eq. 9, considering density equal to zero as boundary conditions. We do not allow changes greater than 10% for the planetesimal surface density profile in each timestep and in each radial bin.

2.2 Growth and orbital evolution of the protoplanetary embryos

At the beginning of our model, an embryo population is immersed in the protoplanetary disc. Here, we describe the details on the embryos growth and orbital evolution.

2.2.1 Solid accretion by the embryos

The embryos grow within the oligarchic regime, and the accretion of planetesimals is well described by the particle in a box approximation (Inaba et al., 2001)

[TABLE]

In this equation $M_{\text{C}}$ is the mass of the core, $\Sigma_{\text{p}}({a_{\text{P}}})$ is the planetesimal surface density at the planet location, $R_{\text{H}}$ is the Hill radius, $P$ is the orbital period, and $P_{\text{coll}}$ is the collision probability. $P_{\text{coll}}$ is a function of the enhanced radius $\tilde{R}_{\text{C}}$, the Hill radius of the planet, and the relative velocity of planetesimals, $P_{\text{coll}}=P_{\text{coll}}(\tilde{R}_{\text{C}},R_{\text{H}},v_{\text{rel}})$. The fact that $P_{\text{coll}}$ is function of the enhanced radius $\tilde{R}_{\text{C}}$ instead of the core radius $R_{\text{C}}$ is due to the fact that we are considering the drag force experienced by planetesimals when entering the planetary envelope. Here, $v_{\text{rel}}=\sqrt{5e^{2}/8+i^{2}/2}~{}v_{\text{k}}$, where $e$ represents the eccentricity, $i$ the inclination, and $v_{\text{k}}$ the keplerian velocity of the planetesimal population (see Guilera et al. (2010) for further details). It is worth remarking that the embryos evolve in circular and coplanar orbits during all the simulation.

In contrast to our previous works, we do not solve the equations of transport and structure for the envelope, thus we calculate $\tilde{R}_{\text{C}}$ following the works of Ormel & Kobayashi (2012) and Chambers (2014) as

[TABLE]

where

[TABLE]

In the previous equations, $R_{\text{B}}$ is the Bondi radius of the planet, $L$ the luminosity generated due to the accretion of planetesimals, $G$ the universal gravitational constant, $R_{\text{C}}$ the radius of the core, $\rho_{\text{p}}$ the planetesimal density, $\sigma_{B}$ the Stefan-Boltzmann constant, $\kappa_{e}=0.01~{}\text{cm}^{2}\text{g}^{-1}$ the envelope opacity, and $P_{\text{neb}},T_{\text{neb}},\rho_{\text{neb}}$ the local pressure, temperature and density of the nebular gas, respectively.

2.2.2 Gas accretion by the embryos

As the embryos grow, they are able to retain a gaseous envelope that is capable to maintain hydrostatic equilibrium. Initially, the gas accretion rate is much lower than the planetesimal accretion rate, thus, the core of the embryo grows faster than the corresponding envelope. However, when the core of the planet reaches a critical mass, the gas accretion begins to be substantial. Following Ida & Lin (2004a) the critical mass is given by

[TABLE]

where, as we have already seen, $\dot{M}_{\text{C}}$ is the planetesimal accretion rate. The gas accretion rate can be estimated by

[TABLE]

where $M_{\text{P}}$ is the total mass of the planet and $\tau_{\text{g}}$ is the characteristic Kelvin-Helmholtz growth timescale of the envelope. According to the results of the giant planet formation model of Guilera et al. (2010, 2014), and following the prescription of Ida & Lin (2004a), $\tau_{\text{g}}$ is fitted by

[TABLE]

Note that, unlike Miguel et al. (2011), we found a different exponent in the dependence with the mass of the planet of $\tau_{\text{g}}$. These authors found an exponent of $-4.89$. This difference is due to the fact that the fit on the gas accretion has been done for slightly different giant planet formation models111Miguel et al. (2011) used the results of Fortier et al. (2009)..

It is important to remark that $\dot{M}_{\text{KH}}$ is valid as long as it is less than the maximum rate at which the gas can be delivered by the disc onto the planet, which, following (Mordasini et al., 2009), is

[TABLE]

Thus, we consider that the effective gas accretion rate is

[TABLE]

The process of gas accretion by the planet is also limited if the planet is able to open a gap. Through high resolution hidrodynamical simulations, Tanigawa & Ikoma (2007) found that the gas accretion rate decreases significantly when a planet opens a gap in what they called the “gap-limiting” case. Following this idea, Tanigawa & Ikoma (2007) developed an analytic formula for the accretion rate also adopted by Xiao & Jin (2015)

[TABLE]

where

[TABLE]

with $\Sigma_{\text{acc}}=\Sigma(x=2R_{\text{H}})$ and,

[TABLE]

and $l$ and $x_{\text{m}}$ are defined as

[TABLE]

[TABLE]

where $M_{\star}$ is the mass of the star, and the sub-index P represents at the location of the planet. Therefore, the gas accretion rate adopted to limit the gas accretion after the planet opens a gap is defined as

[TABLE]

As an example, in Fig. 1 we plot a comparison for the in situ formation of a giant planet at 5 au between the giant planet formation model (Guilera et al., 2010, 2014) and PLANETALP. For the first case, when the mass of the envelope reaches $\sim 40\text{M}_{\oplus}$, the giant planet formation model of Guilera et al. (2010, 2014) is not able to solve any more the equations of transport and structure for the envelope and the simulation stops. It is important to note that there is no limitation in the gas accretion for the first case. We can see that PLANETALP reproduces well the growth of the envelope, specially when the planet begins the gaseous runaway growth (see the zoom in Fig. 1). We also plot the growth of the envelope using the prescription of Miguel et al. (2011) for the gas accretion rate. Different gas accretion rates lead to different final masses of the planet. The small difference between the value of the mass of the cores for the cross-over mass (when the mass of the envelope equals the mass of the core) is because the enhanced radius is calculated in different ways in each model.

It is important to note that the gas accretion rate onto planets is not considered as a sink term in the time evolution of the gas surface density profile. Therefore our model controls that the mass of gas accreted by the planets of the system is not greater than 90% of the mass of gas available in the disc at each time step.

2.2.3 Fusion of embryos and planets with atmospheres

As we have already mentioned, embryos grow by accretion of gas and planetesimals, but they also grow due to collisions with each other. In the present study, we consider that when the distance between two embryos becomes smaller than 3.5 mutual Hill radii, they merge into one object. From this, it is important to specify the mass of the resulting object, which will depend on the physical properties of the interacting bodies. In fact, if the embryos that participate in a collision are not very massive, their individual gaseous envelopes are negligible. In this case, we assume perfectly inelastic collisions, from which the mass of the resulting embryo will be the sum of the individual masses of the interacting embryos. If one or both of the interacting bodies is sufficiently massive to have a significant gaseous envelope, it is necessary to specify the core and envelope masses of the resulting body. On the one hand, the individual cores are perfectly merged in order to define the final core. On the other hand, the evolution of the gaseous envelopes is computed from the studies developed by Inamdar & Schlichting (2015). These authors analyzed the formation of super-Earths and mini-Neptunes and examined how much gaseous envelope could have been accreted by planetary embryos before giant impacts and how much could be retained throughout the giant impact phase. Inamdar & Schlichting (2015) computed the global atmospheric mass loss fraction $\chi_{\text{loss}}$ for ratios of envelope mass to core mass $M_{\text{E}}/M_{\text{C}}$ between 10*-1* and 10*-6* as a function of $v_{\text{imp}}m_{\text{imp}}/(v_{\text{esc}}M_{\text{C}})$, being $v_{\text{imp}}m_{\text{imp}}$ the impactor momentum, and $v_{\text{esc}}$ the escape velocity of the target core of mass $M_{\text{C}}$. The values of $\chi_{\text{loss}}$ used in the present work are those exposed in Fig. 5 of Inamdar & Schlichting (2015). From these considerations, if bodies $i$ and $j$ collide, the core mass $M_{\text{C}}$ of the resulting object will be given by

[TABLE]

where $M^{i}_{\text{C}}$ and $M^{j}_{\text{C}}$ are the core mass of the body $i$ and $j$, respectively, while the final envelope mass $M_{\text{E}}$ will be

[TABLE]

where $M^{i}_{\text{E}}$ and $M^{j}_{\text{E}}$ are the envelope mass of the body $i$ and $j$, respectively, and $\chi_{\text{loss}}^{i}$ and $\chi_{\text{loss}}^{j}$ the atmospheric mass loss fraction of the body $i$ and $j$, respectively.

Finally, by conservation of angular momentum, we compute the semimajor-axis of the resulting body by

[TABLE]

being $a_{i}$ and $a_{j}$ the semimajor-axis of the body $i$ and $j$, respectively.

It is important to remark that the studies of Inamdar & Schlichting (2015) concerning the loss of atmospheric mass were developed for planets with masses in the range of the Neptunes and the so-called super-Earths and initial envelopes less than 10% of the core mass. Due to the absence of works about atmospheric mass loss by giant impacts in large atmospheres, we also used the results derived by Inamdar & Schlichting (2015) to treat collisions that involve giant planets. However, it is important to note that, in our study, the final mass of a giant planet is not sensitive to this consideration. This is due to the fact that a planet that is the final result of a merger between two previous planets, is still growing in the disc. Thus, if this final planet has a significant core mass, it quickly accretes large amounts of gas. In fact, Broeg & Benz (2012) studied in detail the effect of this situation on the gas accretion rate of a planet. They found that initially, most of the envelope can be ejected, but afterwards, the gas is re-accreted very fast.

2.2.4 Type I and type II migration

A planet immersed in a protoplanetary disc modifies the local gas surface density and the mutual gravitational interactions between them, which lead to an exchange of angular momentum between the gas and the planet. This produces torques that cause that the planet migrates towards or away from the central star. Different types of migration regimes have been studied by many authors during the last two decades (Ward, 1997; Tanaka et al., 2002; Masset et al., 2006; Paardekooper et al., 2010, 2011). These regimes modify the local density profile and the planets orbits in different ways depending on the mass of the planet and on the local properties of the disc.

Type I migration (Ward, 1997) affects those planets that are not massive enough to open a gap in the gaseous disc. In idealized vertically isothermal discs, type I migration produces always inward migration and very fast migration rates (Tanaka et al., 2002). Since type I migration is much faster than the gas dissipation timescales of the discs and therefore much faster than the formation timescales of the cores of giant planets, the survival of the planets, particularly of giant planets, results to be a difficult process to fulfill (Papaloizou et al., 2007). Although these migration rates have been corroborated by numerical simulations (D’Angelo et al., 2003; D’Angelo & Lubow, 2008), many authors had to reduce the migration rates by a constant factor in order to reproduce observations (Ida & Lin, 2004a; Alibert et al., 2005; Mordasini et al., 2009; Miguel et al., 2011).

It is important to note that more recent works which consider non isothermal discs (Paardekooper et al., 2010, 2011) have shown that type I migration could be outward. Even in isothermal discs, type I migration could produce outward migration if a full magneto-hydrodynamic (MHD) disc turbulence is considered (Guilet et al., 2013). Moreover, Benítez-Llambay et al. (2015) showed that the heating torques produced by a planet while the process of solid accretion is given, are capable of slowing, halting or even reversing the migration of the planet in the protoplanetary disc if the mass of the planet is lower than $5\text{M}_{\oplus}$.

In order to be consistent with our model of planet formation, PLANETALP, we use the migration rates for isothermal discs of Tanaka et al. (2002) but considering different reduction factors $f_{\text{migI}}$, given by

[TABLE]

where $L_{\text{P}}=M_{\text{P}}\sqrt{GM_{\star}a_{\text{P}}}$ is the angular momentum of the planet, and the total torque is

[TABLE]

Remember that the sub-index P represents the location of the planet.

When a planet is massive enough to open a clear gap in the gas surface density profile, it undergoes type II migration. Crida et al. (2006) showed that a planet of mass $M_{\text{P}}$ is able to clean more than 90$\%$ of the gas around its orbit if

[TABLE]

Once the gap is opened, the orbital evolution of the planet is completely tied to the viscous evolution of the gas disc and, as long as the local mass of the disc is greater than the mass of the planet (“disc-dominated type II migration” following Armitage, 2007), the planet migrates inward at a rate given by

[TABLE]

which is the same rate at which the gas moves in the disc (Ida & Lin, 2004a). When the mass of the planet is high enough to be comparable to the local mass of the disc (“planet-dominated type II migration” following Armitage, 2007), the migration rate is even lower and, as in Edgar (2007) and Mordasini et al. (2009), we calculate this migration rate as

[TABLE]

In this case, the migration rate not only depends on the viscosity of the disc but also on the gas surface density and on the mass of the planet.

2.3 Water distribution on embryos and planetesimals

We also incorporate a water radial distribution for the embryos and planetesimals in our model of planet formation. In order to be consistent with the initial surface density of planetesimals adopted previously (Eq. 2 and 3), the initial fraction per unit mass of water in embryos and planetesimals is given by

[TABLE]

The time evolution of the water radial distribution of the planetesimals is calculated as the Eq. 9 is solved. Moreover, the amount of water in the embryos is calculated self-consistently as they accrete planetesimals, and when two embryos merge between them we consider that the amount of water of the new embryo is the sum of the amounts of water of the previous ones. We remark that we do not consider loss of water either by the accretion of planetesimals or the merger between embryos. Thus, the final amounts of water represent upper limits.

3 Initial conditions for our planetary systems

The main goal of this work is to find favorable scenarios and disc properties for planetary systems similar to our own with PLANETALP, to obtain initial conditions for the post-oligarchic growth of these systems (Ronco $\&$ de Elía in prep., from now on, PII). To do so, PLANETALP has been automated to develop several thousand planetary systems and to make a statistical study of population synthesis. Figure 2 shows an schematic view of the operation of PLANETALP. In this section, we describe the initial conditions adopted to run PLANETALP taking into account that we are not interested in reproducing the actual exoplanet population, but we are interested in generating a great diversity of planetary systems without following any observable distribution in the disc parameters.

3.1 Scenarios and free parameters

As we have already mentioned, type I migration for idealized isothermal discs can be very fast, and thus, the timescales of migration are usually shorter than the disc lifetime. As in Ida & Lin (2004a) and Miguel et al. (2011), we then consider different scenarios for the type I migration introducing the reduction factor $f_{\text{migI}}$, which considers possible mechanisms that can slow down the planet migration 10 and 100 times ($f_{\text{migI}}=0.1$ and $f_{\text{migI}}=0.01$). We also consider the extreme cases in where type I migration is not slowed down at all ($f_{\text{migI}}=1$), and where it is not considered ($f_{\text{migI}}=0$). It is important to remark that we are not considering dynamical torques in our model (Paardekooper, 2014). These torques could slowed down significantly the migration rates estimated from the statical torques given by Tanaka et al. (2002). Sasaki & Ebisuzaki (2017) found that the M-a diagram calculated from a simple planet population synthesis model, based in the work of Ida & Lin (2004a) including statical and dynamical torques, can reproduce the observational exoplanet distribution around sun-like stars. It is important to note that this M-a diagram is similar to the M-a diagram obtained using statical torques and reduction factors in the type I migration rates.

Another important parameter for our model of planet formation is the size of the planetesimals. In the standard model of core accretion, both terrestrial planets and solid cores of gas giants are formed through the accretion of planetesimals (Safronov & Zvjagina, 1969; Wetherill, 1980; Hayashi, 1981). However, the mechanisms for the formation of planetesimals and their initial sizes are still under debate (Johansen et al., 2014; Testi et al., 2014; Birnstiel et al., 2016). While collisional growth models could lead to the formation of sub-kilometer/kilometer planetesimals (Okuzumi et al., 2012; Garaud et al., 2013; Kataoka et al., 2013; Dra̧żkowska et al., 2013; Kobayashi et al., 2016; Arakawa & Nakamoto, 2016), gravitational collapse models predict the formation of tens/hundreds kilometer planetesimals (Johansen et al., 2007, 2012; Cuzzi et al., 2008; Simon et al., 2016; Schäfer et al., 2017). Moreover, the size distribution of the asteroid belt can be reproduced by collisional evolution models using both, an initial population of big planetesimals (Morbidelli et al., 2009) or an initial population of small planetesimals (Weidenschilling, 2011).

In this work, as in Fortier et al. (2013), we consider scenarios with different sizes for the planetesimals ($r_{\text{p}}=100$ m, $1$ km, $10$ km and $100$ km), but we consider a single size per simulation. It is important to note that, while big planetesimals ($r_{\text{p}}\gtrsim 10$ km) are mainly governed by the quadratic regime for the determination of the radial drift and solid accretion rates, small planetesimals ($r_{\text{p}}\lesssim 1$ km) can be in different drag regimes along the disc.

We consider random values from uniform distributions for all the free parameters because we are interested in developing a great diversity of planetary systems without following any observable parameter distribution. However, we take into account that the bounds of the parameter ranges are obtained from previous works. Following Andrews et al. (2009, 2010), we consider the mass of the disc, $M_{\text{d}}$, between $0.01\text{M}_{\odot}$ and $0.15\text{M}_{\odot}$, the characteristic radius $R_{\text{c}}$ between 20 au and 50 au, and the density profile gradient $\gamma$ between 0.5 and 1.5. The factor $\beta$ is taken randomly between 1.1 and 3 given that, in our solar system, Lodders (2003) usually takes a value of approximately 2 and Hayashi (1981) adopts a value of 4. It is important to note that a value of $\beta=1.1$ represents an amount of 9% of water by mass, while a value of $\beta=3$ represents an amount of 66% of water by mass on embryos and planetesimals at the beginning of the simulations and beyond the snowline, following Eq. 39. The viscosity parameter $\alpha$ for the disc and the rate of EUV ionizing photons $f_{41}$ have a uniform distribution in log between $10^{-4}$ and $10^{-2}$, and between $10^{-1}$ and $10^{4}$, respectively (D’Angelo & Marzari, 2012).

Since the maximum value of the mass of the discs considered is relatively high and since this can lead to the occurrence of gravitational instabilities, our model checks, according to the Toomre criterion (Toomre, 1964) given by

[TABLE]

the stability of our discs. When $Q>1$ discs are considered unstable. It is important to remark that all the considered discs in our simulations result to be stable throughout their entire extent.

Other parameters are kept constant for all the simulations. As we are particularly interested in the formation of solar system analogs, the mass of the central star is always $1\text{M}_{\odot}$ and the metallicity is $[\text{Fe}/\text{H}]=0$. Also the planetesimal and embryo densities are considered constant and to be 1.5 g cm*-3* and 3 g cm*-3*, respectively.

Regarding the embryo population, we consider an embryo distribution between 0.1 au and 30 au. These embryos are separated by 10 mutual Hill radii assuming circular and coplanar orbits, and their initial masses, which are of the order of the mass of the Moon, are given by the transition mass between the runaway and the oligarchic regime (Ida & Makino, 1993) as

[TABLE]

where $m_{\text{p}}$ is the mass of the planetesimals.

Finally, with the aim of finding the most suitable scenarios and parameters that provide us with solar system analogs we perform 1000 simulations per each combination between $f_{\text{migI}}$ and $r_{\text{p}}$. This is, we form 16000 different planetary systems that evolve in time until the gas of the disc dissipates. From now on, when we talk about formation scenarios, we are referring to all the combinations between $f_{\text{migI}}$ and $r_{\text{p}}$.

3.2 The disc dissipation timescale

Another important parameter that determines the final configuration of a planetary system during the gaseous phase is the gas-disc dissipation timescale, $\tau$. Haisch et al. (2001) and Mamajek (2009) have observed young stars in clusters of different ages and they suggested that protoplanetary discs have characteristic lifetimes between 3 Myr and 10 Myr. More recently, Pfalzner et al. (2014) showed that, taking into account that the previous works present observational biases, gas discs could last much longer than 10 Myr.

In our model, unlike Miguel et al. (2011) for example, wherein the evolution of the gas is represented by an exponential decay law, the gas disc evolves due to the viscous accretion and photoevaporation processes. Both phenomena together determine the lifetime of the gas, $\tau$. Therefore, we need to form planetary systems whose combinations between the viscosity parameter $\alpha$, the EUV rate $f_{41}$, the mass of the disc $M_{\text{d}}$, the characteristic radius of the disc $R_{\text{c}}$, and the surface density gradient $\gamma$ make the gas discs of these systems to dissipate on timescales according to the results of the previous mentioned works. To do this, we run a version of PLANETALP that is limited to the study of the evolution in time of the gas disc. This is, we “turn off” the evolution of the planetesimal surface density and the evolution of the embryo population. We then generate, using a von Neumann acceptance-rejection method, as many gas discs with dissipation timescales between 1 and 12 Myr as we need, this is, 16000 combinations of all the parameters to develop complete planetary systems.

Fig 3 shows the timescale (in Myr) vs. mass of the disc (in solar masses) plane. Each point represents a set of parameters randomly taken. The 16000 red points represent sets whose combinations between all the parameters obtained a disc of gas that dissipated between 1 Myr to 12 Myr. We find a mean value of $\tau$ of 6.45 Myr with a dispersion of $\sigma=2.67$ Myr. The black points are those sets that resulted in discs dissipating in less than 1 Myr or in more than 12 Myr. These last sets are then discarded for our simulations.

4 Results

The main goal of this work is to find, through the population synthesis analysis developed with our model of planet formation, sets of appropriate initial conditions to study, in a future work (PII), the post-oligarchic evolution of planetary systems similar to our own with N-body simulations. However, through the development of the 16000 simulations, which were performed randomly and uniformly varying the parameters of the disc, we found a great diversity of planetary systems architectures, not only solar system analogs.

4.1 General results

To classify the diversity of planetary systems, we define general properties for the different types of planets. On the one hand, gas giant planets present an envelope which is greater than the mass of the core. Within this category we can find Jupiter-like or Saturn-like planets. The difference between them is that we consider that Jupiter-like planets manage to open a gap in the disc or their total mass is greater than $200\text{M}_{\oplus}$. On the other hand, icy giant planets present envelopes greater than 1% of the core mass but must have envelopes smaller than the core mass. At last, rocky planets present envelopes smaller than 1% of their core mass. Having classified the final planets, from all the diversity of planetary systems developed, we can distinguish three global kinds of planetary systems that we classify as:

• •

Rocky planetary systems: planetary systems which are formed only by rocky planets, without gas or icy giants.

• •

Icy giant planetary systems: these planetary systems can harbor rocky planets but they must present at least one icy giant planet analog to a Neptune-like planet, or a hot/warm Neptune-like planet.

• •

Gas giants planetary systems: planetary systems which must harbor at least one giant planet (in this global characterization we do not distinguish the location of the giant), but can also contain rocky and icy giant planets.

In this work we are particularly interested in those planetary systems which, at the end of the gas phase, present similarities to our own solar system, and so we classify them as:

• •

Solar system analogs (SSA): planetary systems that host only rocky planets in the inner zone of the disc and at least one giant planet beyond 1.5 au.

SSA represent a subkind of those planetary systems with gas giants. Table 1 shows the general percentages of each kind, including the percentage that SSA represent respect to the total of the performed simulations. Although the scope of this work is to focus only on SSA, we will mention some general results of the developed simulations that are in good agreement with previous works.

Following Table 1 it is clear that rocky planetary systems represent the vast majority, indeed, more than 60% of the total of the simulated planetary systems. Moreover, we can find this kind of systems in all the formation scenarios considered. In fact, rocky planetary systems represent more than 40% of the planetary systems in each formation scenario for all of them but one, wherein represent $\sim 20\%$. Finally, we can also find them in the whole range of disc masses considered, between $0.01\text{M}_{\odot}$ and $0.15\text{M}_{\odot}$. This result is in good agreement with the results obtained by Miguel et al. (2011) although they followed observable distributions for the most important disc parameters and although they only considered planetesimals of 10 km.

Fig. 4 shows maps of density of the different planetary systems for each formation scenario. Rocky planetary systems are significantly more favorable in formation scenarios with big planetesimals and low-moderate/null type I migration rates, although, as we mentioned before, we find them in each of all the formation scenarios. This is a natural consequence of the fact that planetesimal accretion rates are smaller for big planetesimals, and thus, the formation of massive cores generally requires formation timescales that exceed the characteristic dissipation timescales of the discs. Besides, the smaller the type I migration rate, the smaller the probability of merger between embryos and thus, the formation of massive cores by this mechanism. It is also important to note that a significant number of rocky planetary systems is formed by small planetesimals and low-moderate/null type I migration rates. This is due to the fact that Fig. 4 does not represents explicitly the dependence with the mass of the discs. As we adopted that the mass of the disc is uniformly distributed between $0.01\text{M}_{\odot}$ and $0.15\text{M}_{\odot}$, there is a significant percentage of systems in which the total mass of the disc is not enough to form massive cores, and consequently, icy giant or gas giant planets.

The opposite case is given by the icy giant planetary systems. These planetary systems are mainly formed in scenarios with small planetesimals and large type I migration rates. The planetesimal accretion rates are greater for smaller planetesimals, situation that lead to the quickly formation of several Earth-cores beyond the snowline. However, large type I migration rates produces that these cores quickly reach the inner edge of the disc, avoiding the formation of gas giants. Therefore, on the one hand, high type I migration rates could lead to the merge of embryos, an thus, the formation of massive cores as we mentioned before. But, on the other hand, high type I migration rates also cause the planets to quickly get into the inner zone, and fail in the gas giant formation.

Finally, gas giant planetary systems are majority for small planetesimals and low-moderate/null type I migration rates. As it was already mentioned, small planetesimals favor the formation of massive cores, and low-moderate/null type I migration rates avoid the planet to drop into the inner edge of the disc, allowing the planet to trigger the gaseous runaway growth. However it is important to note that we find gas giants in scenarios with high type I migration rates. About a $70\%$ of the planetary systems formed in scenarios with high type I migration rates, considering all the chosen planetesimal sizes, present Hot/Warm-Jupiters within 0.5 au. However, taking into account all the formation scenarios that form giant planets, not only those with high type I migration rates, the planetary systems with Hot/Warm-Jupiters represent the 43% of Gas Giant planetary systems. Moreover, those planetary systems with Hot/Warm-Jupiters represent only a 9.5% of the total of the developed simulations. Thus, although we find higher percentages of Hot-Jupiter occurences than those predicted observationally, which are $\sim 1\%$ (Marcy et al., 2005; Cumming et al., 2008; Mayor et al., 2011; Wright et al., 2012), Wang et al. (2015) concluded that the current knowledge of stellar properties and the stellar multiplicity rate is yet very limited to reach quantitative results for the Hot-Jupiter occurrence.

We remark again, that these density maps do not take into account the dependency with the properties of the disc, i.e. the mass of the discs, the exponent $\gamma$, the characteristic radius, and the factor $\beta$.

4.2 Solar system analogs

Of the performed simulations, we found that only 688 protoplanetary discs formed SSA, this represents only a 4,3% of the total. Table 2 shows the percentages of SSA for each formation scenario and Figure 5 shows $\gamma$ vs. $M_{\text{d}}$ planes for each formation scenario. In this figure, each point represents a formed planetary system with a particular set of disc parameters. The colored points represent SSA, the grey points are the rest of the planetary systems performed in which we are not going to focus in this work. The colorscale shows the dissipation timescale of the gas disc for each planetary system. As we can see in table 2 and figure 5, we do not find SSA in scenarios with non-reduced type I migration. This is due to the fact that, as we have already mentioned, type I migration without any reduction factor causes that giant-forming planets in the outermost zone of the disc migrate to within 1.5 au in timescales smaller than the dissipation of the gas disc, regardless of the size of the planetesimals considered. Clearly, formation scenarios of low (or none) type I migration rates and small planetesimals (100 m and 1 km) are the most favorable for the formation of SSA. Since the planetesimal accretion rates are higher for small planetesimals, the rapid formation of massive cores before the dissipation of the gas is favored. Otherwise, for large planetesimals, the formation timescales of massive nuclei are not high enough to trigger the growth of the gaseous envelope and to form giant planets in timescales according to the disc lifetime.

4.2.1 Solar system analogs architectures

As we already mentioned, SSA must harbor at least one giant planet beyond 1.5 au and rocky planets in the inner regions of the disc. However, the architectures of the SSA found are quite different, and not all of them present icy giants or only one gas giant planet. Taking into account all the SSA, without distinguishing between planetesimal sizes or migration rates, the most representative architectures found among the SSA are those that present between 1 to 3 gas giants, between 0 and 4 icy giants and between 100 and 200 rocky planets along the disc. In fact, the mean of the number of gas giants is 2, with a dispersion of $\sim 1$, the mean of the number of icy giants is also 2, with a dispersion of $\sim 2$, and the mean of the number of rocky planets is 150, with a dispersion of $\sim 47$. Particularly, Figure 6 displays histograms showing the percentage of SSA that present one, two, three, four or five gas giants, and also shows how many systems present different numbers of icy giants and rocky planets, but differentiating the systems according to the size of planetesimals and migration rates. For planetesimals of 100 km, only one giant planet is formed in the only three systems found for low or null type I migration. These systems also harbor only one, or none, icy giant, but a significant number of rocky embryos, some of them in the outer part of the disc. For these big planetesimals, we do not find SSA formed with moderate and non reduced type I migration rates. The formation of two giant planets is more likely for low or null type I migrations rates and scenarios with planetesimals of 10 km, while the formation of only one gas giant is more favorable for moderate type I migration rates. It is important to note that, in both cases, there is a not negligible probability to form systems with three or four giant planets. These systems also harbor preferentially a small number of icy giants. However, we found that for the case of null or low type I migration rates, it is possible to find SSA that harbor until 8 icy giants. These kind of systems also present preferably between 100 and 200 rocky embryos. For small planetesimals of 1 km and 100 m we find similar results. In both cases, planetary systems prefentially harbor one or two giant planets, with the difference that no giant planet is formed for planetesimals of 1 km of radius for moderate type I migration rates. In both cases, a small number of icy giants is also expected. Finally, a large number of rocky planets are expected in such systems when the gas of the disc is dissipated. It is important to remark that the stability and dynamical evolution of these kind of systems after the gas dissipation will be discussed in the forthcoming paper.

4.2.2 Favorable disc parameters

Taking into account all the SSA, table 3 shows the ranges, mean and dispersion values for the disc parameters. As we can see, we can find SSA for all the ranges considered of $\gamma$, $R_{\text{c}}$, $\alpha$ and $\tau$. However, we do not find SSA for low-mass discs ($M_{\text{d}}<0.04\text{M}_{\odot}$) because low-mass discs have not enough solid material to form massive cores before the disc dissipation timescales. Similar results were found by Thommes et al. (2008) and Miguel et al. (2011). Although we find SSA for almost all the values in the ranges of the free disc parameters, they are preferentially formed for massive disc ($<M_{\text{d}}>=0.1~{}\text{M}_{\odot}$), smooth surface densities profiles ($<\gamma>=0.95$), moderate values for the disc characteristic radius ($<R_{\text{c}}>\sim 33$ au), small values of the $\alpha$-viscosity parameter ($<\alpha>=3.4\times 10^{-4}$), and moderate timescales for the disc dissipation ($<\tau>\sim 6.5$ Myr). It is worth mentioning that only 20 SSA have all their disc parameters between the dispersion values, $\pm\sigma$. These planetary systems are represented with black diamonds in Figure 5. This represents only a $\sim 3~{}\%$ of the total SSA. All these systems are formed considering small planetesimals (100 m and 1 km) and adopting null or low type I migration rates.

4.2.3 Evolution and final configuration of solar system analogs

In this section we describe the general characteristics of the evolution in time of a SSA, and also show the final configurations of the most representative SSA of each formation scenario, which are composed of:

• •

An embryo distribution that provide us with information about the semimajor-axis, mass of the core, mass of the envelope, mass of water due to the accretion of planetesimals and mass of water due to the fusion with other embryos, of each final body.

• •

A planetesimal surface density, eccentricity and inclination profiles.

These final configurations are going to serve as initial conditions to analyze the post-oligarchic stage of planet formation, via N-body simulations (PII).

Figure 7 shows the evolution in time, represented by the colorscale, of a representative SSA for the formation scenario without type I migration and with planetesimals of 1 km. The parameters of this particular system are $M_{\text{d}}=0.13\text{M}_{\odot}$, $\gamma=0.92$, $R_{\text{c}}=34$ au, $\alpha=1.1\times 10^{-3}$, and a dissipation timescale of $\tau=6.65$ Myr. The factor that represents the amount of solid material condensed beyond the snowline, at 2.7 au, is $\beta=1.96$, in this example. Panel (a) represents the time evolution of the embryo distribution. The small red points linked by lines represent the initial embryo distribution. These embryos start to grow due to the accretion of planetesimals within their feeding zones. In time, as they grow, if the distance between them is less than 3.5 Hill radius they merge into one object. When the cores reach a critical mass, which is approximately $10\text{M}_{\oplus}$, they start to accrete significant amounts of gas until they achieve their final masses and the disc is dissipated. Panels (b) and (c) represent the evolution in time of the gaseous component of the disc. The gas surface density profile (b), which is plotted every 0.1 Myr, begins to decrease and to expand in the radial direction due to the disc angular momentum conservation. At approximately 2 Myr, a gap is opened in the disc due to the photoevaporation process. As a consequence, the gas in the inner zone of the disc, inside the gap, rapidly falls to the central star, while the gas in the outer zone continues its evolution for a few million years. The mass of the gas in the disc decreases as it is shown in panel (c). The planetesimal surface density (d) evolves in time due to the radial drift and due to the embryos, which accrete and excite them, increasing their eccentricities (panel e) and inclinations (panel f). As time advances, the planetesimal profile decreases, and the inner zone of the disc quickly runs out of planetesimals. This is mainly due to the fact that the planetesimal accretion rates are higher in the inner zone where the planetesimal surface density is initially higher. This favors the formation of many, low-mass rocky planets very close to each other, which accrete all the available solid material. These embryos do not merge between them due to the fact that their feeding zones are very narrow, an thus, they do not grow enough to be separated less than 3.5 mutual Hill radius (specially, very close to the star). Besides, the formation of two giant planets (see panel a) helps to remove planetesimals between 10 au and 15 au. At the end, when there is no more gas in the disc, there are only planetesimals behind 15 au. This remnant of planetesimals (represented with the blue line) present final eccentricity and inclination profiles according to the blue lines in panels (e) and (f), respectively. It is important to note that these eccentricity and inclination profiles along the disc have physical sense only in the region where there is still presence of planetesimals, for this particular case, beyond $\sim$ 15 au.

As we mentioned before, the aim of this work is to find suitable initial conditions for SSA to develop N-body simulations, thus, Figure 8 shows representative final configurations of SSA in all the formation scenarios in which they were found. Each panel shows the final embryo distribution with the final planetesimal density profile of each system at the top, and the final planetesimal eccentricity and inclination profiles in the bottom. The colorscale represents the percentage of water by mass.

Although all this planetary systems are SSA, they present global differences regarding on the location of the giant planets and on the final planetesimal surface density profiles. On the one hand, as the size of the planetesimals and the reduction factor of type I migration increase, the location of the giant planets changes, moving inward. On the other hand, also the planetesimal remnant until 30 au increases towards the central star. The final location of the giant planet depends on several phenomena. As it was shown by Guilera et al. (2011), the competition between the planetesimal accretion timescales and the planetesimal migration timescales, regulates the formation of a giant planet. These authors found that these two timescales depend on the planetesimal size and on the location of the planet with respect to the central star. Thus, the optimal location in the disc for the formation of a massive core, precursor of a giant planet, is different for different planetesimal sizes. Besides, the main properties of the protoplanetary disc, i.e. the slopes of the surface densities, the jump in the position of the snowline due to the condensation of volatile material, the characteristic radius, etc., can play an important role in the formation of a giant. In fact, Guilera et al. (2011) found that for smooth slopes, the formation of a giant planet is optimized in the outer part of the disc for scenarios with small planetesimals ($r_{\text{p}}<100$ m). Moreover, another important phenomenon that could change the optimal location of the formation of a giant planet, is the merger between several planets.

As an example, Figure 9 shows the time of cross-over as a function of the position of the planet in the disc. The time of cross-over is the time at when the mass of the envelope of a planet equals the mass of the core (mass of cross-over) and the gaseous runaway triggers. The black points and black diamonds in Figure 9 represent the in situ formation of an isolated planet, located at different positions between 2 au and 15 au, for the same disc parameters of the formation scenario of panels d) and i), of Figure 8, respectively. As we can see, for scenarios with planetesimals of 1 km, the optimal location for the formation of a giant planet is around 5 au, while for scenarios with planetesimals of 100 km, the optimal location is 3 au. It is worth mentioning that we find similar results between the scenarios of 1 km and 100 m, and between the scenarios of 10 km and 100 km. However, when we include the formation of the whole complete system, the optimal location for the formation of the giant can change, and can move away from the central star. As Figure 9 displays, for scenarios with planetesimals of 1 km the red point in Figure 9 shows the location at where the first planet of the system reached the mass of cross-over. This situation shows that, when the simultaneous formation is considered, the formation of giant planets can be optimized in outermost regions. This is due to the multiple mergers between embryos that have taken place during the formation of the system. For planetesimals of 100 km, the location of the giant planet is optimized a little bit away from the position found for the isolated formation. In this scenario, as the embryos do not grow enough to increase their Hill radius, the rate of mergers is lower.

Continuing with the description of Figure 8, if we only consider those planetary systems without type I migration, to distinguish the effect caused by the increase in the planetesimal size, (see panels a, d, f and i), we can see that for planetesimals of 100 m, the giant planets form in the outer regions of the disc. Since the accretion rates are higher than for big planetesimals, massive cores form in a shorter timescale, and thus, the planets in the outer region can merge between them before the dissipation of the gas disc, giving rise to giant planets in the outer zones. These giants clean their surroundings of planetesimals, leaving a very small remnant almost beyond 20 au. Despite both planets are Jupiter-analogs, the innermost giant opened a gap in the disc and moved inward through type II migration, while the outermost did not. A brief discussion about why a giant planet with $\sim 4.8\text{M}_{\text{J}}$ did not opened a gap in the disc is discussed in the next section.

Following with the scenario of planetesimals of 1 km (panel d), as the planetesimal accretion rates decrease with respect to the previous scenario, those planets in the outer zone grow achieving masses between $\sim 1.5\text{M}_{\oplus}$ and $\sim 6.5\text{M}_{\oplus}$. These planets do not grow enough to merge between them and thus, we do not find giants on the outer part of the disc. Moreover, they do not remove their feeding zones completely and do not empty them from planetesimals. Therefore, at the end of the gas phase, a population of planets of several Earth masses and a population of planetesimals coexist in the outermost regions of the disc. A similar situation occurs for scenarios with planetesimals of 10 km (panel f). But in this case, the optimal location for the formation of the giants moves inward a little bit. At last, for scenarios with planetesimals of 100 km, the formation of giant planets is optimized near the snowline, at $\sim 3$ au. In panel i), the only formed giant reached 2 au due to the fact that opens a gap in the disc and moves through type II migration.

Regarding the population of planetesimals, a remnant of planetesimals survived at the end of the gaseous phase in the majority of the performed simulations inside 30 au. However, it is worth mentioning that we found a few simulations of planetary systems without planetesimals inside 30 au, for scenarios with planetesimals of 100 m. How this result will affect, or not, the post-oligarchic growth of the planets in these system is an analysis that we are going to take into account in our next work.

Another important thing to note, with respect to the final planetesimal population, is that some planetary systems, present a small remnant of planetesimals inside the location of the giant planet. This can be seen, for example, in panels c) and h) of Figure 8. Particularly, the remnant of planetesimals in the formation scenario of panel h) is located between 2 au and 3 au. This seems to be of interest, since it could lead to the formation of an asteroid belt between the terrestrial-planet forming region and the giant planets. However, it will be important to analyze in a future if these planetesimals would survive in the disc, since after the dissipation of the gas they present eccentricities and inclinations greater than $\sim 0.4$ and $10^{\circ}$.

One of the differences between our model model of planet formation and other population synthesis models is that we include a detailed treatment of the orbital evolution of the planetesimal population. This improvement allows us to obtain more realistic planetesimal surface density, eccentricity and inclination profiles at the end of the gaseous phase to use as initial conditions for the planetesimal population in the post-gas evolution of these systems.

4.2.4 Gap opening planets

In SSA the frequency of a gap opening planet in the disc is relatively low. For each forming-giant planet that opens a gap in the disc, there are 5.5 giant planets that did not. As we described in Sec. 2.2.4, the condition that a giant planet must fulfill to open a gap, is that proposed by Crida et al. (2006). This condition depends basically on the mass and the location of the planet, and on the viscosity and the height scale of the disc. More explicitly, this condition depends on the sum of two terms: the first term involves the ratio between the height scale of the disc and the Hill radius of the planet, this last is a function of the mass and position of the planet; the second term involves the ratio between the viscosity of the disc and a function of the location and mass of the planet. Thus, the opening of a gap is favored by massive planets, and low viscosity and flat/thin discs (i.e. discs with small aspect ratios $H_{\text{g}}/R$). However, our model considers a flare disc, hence, $H_{\text{g}}/R\propto R$ and it increases with the distance to the central star. This situation coupled with having moderate/high values of the viscosity of the discs, results in the low occurrence of gap opening for giant planets in our model. As an example, Figure 10 shows the minimum mass needed to open a gap, adopting the criteria developed by Crida et al. (2006), as a function of the planet semimajor-axis for discs with different viscosities. We show the limit cases for the values of the $\alpha$ parameter adopted in this work, and the mean value ($\pm$ the dispersion) of such parameter. We also plotted the curve corresponding to the planetary system described in panel (a) of Fig. 8 (grey line), and also the final location of both giant planets of that system. We can see that the inner Jupiter-analog has a final mass greater than the minimum needed to open a gap in such disc. However, the outer Jupiter-analog, despite of having a mass of $\sim 4.8~{}\text{M}_{\text{J}}$, did never achieved the minimum mass needed to open a gap.

4.2.5 Water accretion on planets of the habitable zone

A point of interest for our SSA is to analyze the final amounts of water by mass in those planets that remain, at the end of the gaseous phase, in the habitable zone. Taking into account the classical definition, the habitable zone is defined as the range of heliocentric distances at which a planet can retain liquid water on its surface (Kasting et al., 1993). Ten years later, Kopparapu et al. (2013a, b) proposed that our solar system presents a conservative habitable zone determined by loss of water and by the maximum greenhouse effect provided by a CO2 atmosphere, between 0.99 au and 1.67 au. Moreover, these authors proposed that there is also an optimistic habitable zone between 0.75 au and 1.77 au, limits that were estimated based on the belief that Venus has not had liquid water on its surface for the last 1 billion years, and that Mars was warm enough to maintain liquid water on its surface. Although the definition of the habitable zone is a topic under permanent debate, and although there is no guarantee that a planet inside this region will be potentially habitable, we consider this definition just as a first approximation. As it was mentioned before in Sec. 2.3, both embryos and planetesimals present, at the beginning of our simulations, a water radial distribution. As the embryos grow and the planetesimal surface density evolves, also the amount of water on each body changes. It is important to notice that, as we can see in Figure 8 those planets in SSA that remain within the habitable zone at the end of the gaseous phase do not present any water contents. This is due to the fact that the embryo distribution beyond the snowline act as a barrier preventing that water rich planetesimals from the outer zone are accreted by embryos of the inner zone. Moreover, the low/null type I migration rates also prevent that embryos from beyond 2.7 au reach the habitable zone. Although these planets are dry when the gas in the disc dissipates, this does not mean that they can not acquire significant amounts of water for the possible development of life during the post-oligarchic growth stage via, for example, accretion of wet embryos or planetesimals from beyond the snowline. This analysis is out of the scope of this work but is going to be taken into account in PII.

5 Summary and discussion

During the last years, advances in observational astronomy have allowed us to increase the sample of confirmed exoplanets very quickly. This sample, which continues to grow day by day, has motivated the theoretical astronomers to improve their models of planet formation to reproduce the most important properties of this population.

In this first work, we improved our model of planet formation, PLANETALP, incorporating relevant physical phenomena to the formation during the gaseous phase. In our model, the gaseous component of the disc evolves as a viscous accretion disc with photoevaporation due to the central star. We adopted an EUV-photoevaporation model based in the works of Dullemond et al. (2007) and D’Angelo & Marzari (2012). Many theoretical studies link the photoevaporation process with the observed inner holes in transition discs (Owen et al., 2012; Owen, 2016; Terquem, 2017). Particularly, Ercolano et al. (2017), were able to show that the gap in the disc around TW-Hya, which is the closest protoplanetary disc to the Earth, is consistent with their photoevaporation models. However it is important to note that these models are a combination of EUV and X-rays (or just X-rays) as photoevaporation sources. The percentage of observed transitional discs, which are discs with an inner hole of several au, respect to the total number of observed discs, is relatively small, of about $10\%-20\%$. The low percentage of observed transition discs suggests that once the gap is opened in the disc by photoevaporation, the disc should dissipate in a $10\%-20\%$ of the dissipation timescale. In this sence, it is important to note that our model of purely EUV-photoevaporation source could be overestimating the dissipation time after the gap is opened.

The model also presents a solid component, composed by a planetesimal and embryo population. The code calculates in detail the orbital migration of the planetesimal population due to the gas drag, and the evolution of their eccentricities and inclinations due to the embryo gravitational excitation, and damping by the gas. The embryos grow in the disc due to the accretion of planetesimals and the surrounding gas. PLANETALP also considers the fusion between embryos, and planets with significant atmospheres. It is worth mentioning that the gas accretion of the planets is not calculated by solving the equations of transport and structure for the envelope. Instead, we fit the results of the giant planet formation model of Guilera et al. (2010, 2014). Type I and type II migration are also included for the embryos. We also incorporate a water radial distribution for both populations.

The principal aim of this work, is to find scenarios and disc parameters that favor the formation of solar system analogs. Obtaining final embryo distributions and final planetesimal density, eccentricity and inclination profiles of these solar system analogs at the end of the gas phase will allow us to study and analyze in a future work (PII), the post-oligarchic growth of these systems, using N-body simulations.

First, we define initial conditions for our simulations, considering that we are interested in generating a great diversity of planetary systems without following observable distributions for the disc parameters. However, the bounds of the disc parameter ranges considered, come from previous observational works (Andrews et al., 2010). We also define different scenarios for the type I migration, introducing a reduction factor $f_{\text{migI}}$, which takes into account that the planet migration could be not considered ($f_{\text{migI}}=0$), slowed down 10 times ($f_{\text{migI}}=0.1$), 100 times ($f_{\text{migI}}=0.01$), or not slowed down at all ($f_{\text{migI}}=1$). Besides, we consider different planetesimal sizes of 100 m, 1 km, 10 km and 100 km. In order to find gas discs that dissipate in timescales according to the observed ones, we run a version of PLANETALP limited to the study of the evolution of the gas disc, without considering the evolution of the embryo and planetesimal population. We then generate as many gas discs with dissipation timescales between 1 and 12 Myr as we need. Finally, we automate our planet formation model to generate a population synthesis of 16000 planetary system simulations, separated in 16 blocks of different planetesimal sizes and type I migration rates.

Our general results show that the most common planetary systems are those with only rocky planets. Moreover, these systems represent more than $60\%$ of the total of the performed simulations. This general result is in agreement with previous works of population synthesis analysis (Thommes et al., 2008; Mordasini et al., 2009; Miguel et al., 2011). Distinguishing the formation scenarios, rocky planetary systems are majority in scenarios with big planetesimals and low-moderate/null type I migration rates. On the contrary, icy giants are mostly in scenarios with small planetesimals and high type I migration rates. Finally, to most suitable formation scenarios to form giant planets are those with small planetesimals and low-moderate/null type I migration rates. Those giants formed with high migration rates end up as hot-Jupiters, very near the central star.

Regarding solar system analogs (SSA), classified as those planetary systems with only rocky planets in the inner zone of the disc and with at least one giant planet beyond 1.5 au, they represent only a $4.3\%$ of the total of performed simulations. The most representative SSA architectures found are those that present between 1 to 3 gas giants, between 0 to 4 icy giants, and between 100 to 200 rocky planets along the disc. The most favorable formation scenarios for this kind of systems are those with small planetesimals and low/null type I migration rates. SSA were preferentially formed in discs with smooth surface densities profiles ($<\gamma>=0.95$), moderate values for the characteristic radius ($<R_{\text{c}}>\sim 33$ au), small values of the $\alpha$-viscosity parameter ($<\alpha>=3.4\times 10^{-4}$), moderate timescales for the disc dissipation ($<\tau>\sim 6.5$ Myr) and massive discs ($<M_{\text{d}}>\sim 0.1\text{M}_{\odot}$). Although we found SSA in almost all the disc range parameters considered, we did not find them in low-mass discs, with $M_{\text{d}}<0.04\text{M}_{\odot}$. Regarding water delivery, we found that those embryos that remained within the habitable zone, were totally dry when the gas disc had dissipated.

Despite our model of planet formation includes important physical phenomena for the formation of a planetary system during the gas phase, it yet does not include many fundamental interactions that may affect our final results. Perhaps, the most important are:

Type I migration rates for non ideal isothermal discs. Paardekooper et al. (2010, 2011) found that type I migration rates for non isothermal discs can be different with respect to those calcuted previously for isothermal discs. Type I migration rates for isothermal discs only depend on the disc surface density profile, while for more realistic discs, these migration rates also depend on the temperature profile of the disc. Moreover, while ideal isothermal discs in general lead to a rapid inward migration, type I migration for non isothermal discs can be outward, depending on the temperature structure of the disc, and on the mass and semimajor-axis of the planet. Dittkrist et al. (2014) studied the impact of planet migration models on planetary population synthesis incorporating type I migration rates for non isothermal discs. Comparing population synthesis adopting migration rates for isothermal (without any reduction factor) and non isothermal discs, they found that despite the M-a diagrams are similar for both cases, the percentage of planets that remain outside 0.1 au in the first case is $\sim$ the half, compared to the second case. Thus, type I migration rates inferred for more realistic discs could change the percentage of our SSA. Otherwise, Coleman & Nelson (2014, 2016) also incoporated type I migration rates for non isothermal discs in a planet formation model. These authors found that the formation of giant planets is favored by the accretion of small planetesimals ($r_{\text{p}}<100$ m) and that the existence of radial disc structures is necessary for the survival of such planets in the outer parts of the discs. However, we note that the planet formation model of Coleman & Nelson (2014) is quite different with respect to our own, or the one developed by Dittkrist et al. (2014), specially in the treatment of the planetesimal accretion and growth of the planets. Finally, it is important for us to remark again that, to be consistent with our assumption of isothermal disc, we followed Tanaka et al. (2002) prescriptions for type I migration rates. The impact of considering more realistic discs and type I migration rates for non isothermal discs in the formation of SSA will be subject of a future work. 2. 2.

Gravitational interactions between planets and mean-motion resonances. A key effect, that might significantly alter the results of our simulations, is the fact of considering gravitational interactions between the planets. These interactions could cause two different effects. On the one hand, when the gas is almost removed in the inner regions of the disc due to photoevaporation, the dispersion between gaseous giant planets can lead to the ejection of one or two planets exalting the eccentricities of the remaining ones, although there is still gas in the outer zones. Moreover, Matsumura et al. (2013) showed that even the terrestrial planets of the inner zone of the disc can be affected by the giants, despite being far from them, in discs without gas. These effects can alter the final configurations of our SSA. On the other hand, gravitational interactions could allowed the planets to be trapped in mean motion resonances (MMR). The migration of planets trapped in MMR along the gaseous disc is a complex phenomenon. Towards this mechanism, planets trapped in MMR could be able to avoid a fast orbital decay into very inner zones of the disc (Masset & Snellgrove, 2001). Moreover, regarding the formation of our solar system, if Jupiter was able to open a gap in the disc and migrate inward through type II migration, and at the same time Saturn was able to migrate faster towards Jupiter, they could have been locked in a 2:3 (MMR). This effect could have stopped or even reversed the migration of both planets toghether (Morbidelli & Crida, 2007). Then, Morbidelli et al. (2007) showed that if the gas giant planets were trapped in a MMR, the icy giant planets could also be trapped in MMR and all the system could evolve locked in a MMR chain. These results are of great importance, since they represent the initial orbital configuration of the outer solar system after the gas disc is dissipated, assumed by the Nice model (Tsiganis et al., 2005; Gomes et al., 2005; Morbidelli et al., 2005). However, it is important to note that several phenomena like, for example, disc turbulences, could break the MMR configurations (Adams et al., 2008). The inclusion of the treatment of these interactions in PLANETALP is one of our future goals. 3. 3.

Planetesimal fragmentation. Planetesimal collisional evolution could have and important impact in the population synthesis results. In fact, Guilera et al. (2014) developed a planetesimal fragmentation model to study the role of planetesimal fragmentation in giant planet formation. In line with the results of the pionner work of Inaba et al. (2003), and Ormel & Kobayashi (2012), the authors also found that substantial amounts of mass could be lost in the planetesimal collisional evolution process reducing the efficiency of planet formation, specially for small planetesimals, which have a lower specific impact energy. However, it is important to remark that these works considered that all the mass generated by the planetesimal fragmentation process distributed below some minimum particle size, is lost. However, Chambers (2014) found that relaxing this hypothesis, considering that the very small particles below such minimum particle size can quickly coagulate avoiding the loss of material by the collisional process, planetesimal fragmentation could favour the planet formation process. It is important to remark that the inclussion of planetesimal fragmentation in a population synthesis analysis is very costly numerically speaking and imposes a very important limitation.

Finally, it is important to remark that our results provide us with planet distributions and planetesimal density, eccentricity and inclination profiles at the end of the gaseous phase. These distributions will be use as initial conditions to develop N-body simulations, with the aim of studying the post-oligarchic stage of formation, focusing on the formation of terrestrial planets and their final amounts of water. The planet distributions provide us with information about the location, core and envelope mass of the final planets, as well as water contents due to the planets primordial contents and due to the accretion of embryos and planetesimals. We also obtain detailed characteristics of the planetesimal population regarding their final eccentricities, inclinations and also, water contents. The role and importance that these details may have in the final results and the stability and dynamical evolution of these systems in the post-gas phase are subjects that will be developed in the next work (PII).

Acknowledgements

We thank the referee, John Chambers, for constructive comments that helped us to improved the manuscript. This work is supported by grants from the National Scientific and Technical Research Council through the PIP 0436/13, and National University of La Plata, Argentina, through the PIDT 11/G144.