© ESO, 2016

1. Introduction

The core accretion model (Perri & Cameron 1974; Mizuno 1980; Bodenheimer & Pollack 1986; Pollack et al. 1996) is the most widely accepted scenario to explain the formation of a vast diversity of planets (e.g. Alibert et al. 2005; Mordasini et al. 2009; Guilera et al. 2011). The central idea of the core accretion model can be summarised as follows. First, a solid core must be formed from the accretion of planetesimals or pebbles. Once this core reaches approximately a lunar mass, the core gravity is strong enough to start binding some gas from the protoplanetary disc. From this stage on, the protoplanet thus continues to grow by accreting both solids (planetesimals or pebbles) and gas (basically H-He). Planet formation models typically assume for simplicity that solids and gas do not mix: all the solids deposit their mass and energy at the top of the core, and the primordial H-He is collected above, building the atmosphere (or envelope). This is, of course, a very strong and unrealistic simplification: bolides that traverse a gaseous atmosphere undergo thermal ablation and mechanical breakup. Hence, volatile material can vaporise and mix with the primordial H-He, changing the composition of the envelope during the formation of a planet.

If planetesimal or pebble disruption did not occur during formation, then the envelope metallicities of planets should instead be sub-stellar because the gas accreted into the planets should in principle be metal-poor compared to the central star (the metals condense to form planetesimals or pebbles). This is not what is observed in the solar system, where the giant plants show some level of envelope enrichment (Irwin et al. 2014; Niemann et al. 1998; Guillot & Gautier 2014). The alternative hypothesis to planetesimal or pebble dissolution for explaining envelopes enriched in heavy elements is core erosion (Wilson & Militzer 2012). From an energetic point of view, it is possible to mix part of the core upwards (Guillot et al. 2004). In addition, core material is miscible in metallic hydrogen (Wilson & Militzer 2012), which allows for the heavy elements (if they can be lifted up) not to settle to re-form a core. For the mixing of an initial core within the gaseous envelope, Vazan et al. (2015) showed that this process is favoured when an initial compositional gradient exists in the interior of the planet. Thus, the mixing of heavy elements on the planetary envelope seems to be more likely if the heavy elements are not initially concentrated in well-defined core. The formation of such a diffuse core requires planetesimal dissolution in the deep envelope during the formation of the planet. Hence, even if core erosion could play a relevant role in mixing heavy elements in the planetary envelope, this process seems to demand as well that the envelope is initially enriched by planetesimal or pebble dissolution.

The problem of considering an envelope that is enriched with respect to stellar values during the formation of a giant planet has been raised since the very early studies of planet formation. Already in 1986, Bodenheimer & Pollack (1986) mentioned among another two problems that remained to be solved “the fact that the molecular weight of the envelope is expected to increase with time as some of the icy planetesimals dissolve in it”, and added that this problem “could significantly change the accretion scenario”. Stevenson (1982) showed that the critical core mass (the mass required to trigger rapid accretion of gas) is reduced when the envelope mean molecular weight increases. Moreover, Wuchterl (1993) showed a direct dependence of the critical core mass with the adiabatic gradient. The adiabatic gradient is expected to decrease when chemical reactions take place. This necessarily occurs when abundant elements such as H, C, and O are present in the envelope. Therefore, the effect of polluting the primordial envelopes reduces the critical core mass not only through increasing the mean molecular weight but also by reducing the adiabatic gradient (Hori & Ikoma 2011). Another effect that can reduce the adiabatic gradient is the condensation of species. We showed (Venturini et al. 2015) that if a planet forms in cold regions of the disc, such that the temperatures and pressures are low enough for certain species to condense in the atmosphere of the protoplanet, then this effect leads to an even larger reduction of the critical core mass.

Despite its expected importance for the formation of giant planets, the effect of envelope enrichment has so far never been implemented in a self-consistent way in any evolutionary calculation of planet formation. This is of course a consequence of the difficulty in modelling all the processes involved, which include planetesimal thermal ablation and dynamical breakup, mass and energy deposition in different envelope layers, and the inclusion of a self-consistent change in the envelope microphysics (opacities and equation of state) as the envelope metallicity evolves.

The magnitude of the enrichment depends on the accretion rate of solids and on their size and strength properties. The smaller and more porous the bolide, the easier it is to disrupt it and mix it with the gaseous envelope when crossing the atmosphere (Podolak et al. 1988). This tells us that the effect of envelope enrichment is more relevant for smaller bolides. Hence, when growing planets from cm-m size particles (the so-called pebbles, Lambrechts & Johansen 2012; Lambrechts et al. 2014), it is necessary to include this effect.

In this work we compute for the first time the in situ growth of a planet taking into account envelope enrichment by icy planetesimals or pebbles. Given the uncertainties in the initial size distribution of planetesimals, we test the effect of envelope enrichment corresponding to a wide range of particle sizes. We solve internal structure equations using global energy conservation arguments and taking into account the changes of envelope metallicity in the opacities and equation of state. In Sect. 2 we explain the numerical method, Sect. 3 is devoted to discussing our assumptions, in Sects. 4 and 5 we show and analyse our results, and in Sects. 6 and 7 we discuss our main results and predictions. We summarise our conclusions in Sect. 8.

2. Method

We used the numerical method developed in Venturini et al. (2015) to solve the equations of hydrostatic equilibrium, mass conservation, and heat transport. We applied the usual Schwarzschild criterion to distinguish between radiative and convective layers, adopting an adiabatic gradient for the latter case.

For the energy conservation, we used global energy conservation arguments (Mordasini et al. 2012; Fortier et al. 2013; Piso & Youdin 2014) to determine the total luminosity radiated away by the protoplanet. In the following subsection we show how we computed this total luminosity.

2.1. Computation of the total luminosity

We computed the total luminosity (L) by analysing the total energy of the system at time t and at time t + dt, as is sketched in Fig. 1.

The planetesimals are assumed to have started with zero velocity at infinity, therefore their total energy is zero. The release of their gravitational potential energy, which heats the envelope, appears as the term of the change in the gravitational potential energy of the core, as we show below. Hence, at a given time t, the total energy of the system is given by the sum of:

• the total energy of the envelope (gravitational plus internal): E_g,env(t) + E_i,env(t);

• the gravitational potential energy of the core (we neglected heat sources coming from the core ¹): E_g,core(t);

• the total energy of the gas layer that will be accreted, represented by the term e_gas,accdm, where e_gas,acc is the total energy of the layer per unit mass and dm the mass of gas that will be accreted;

• the energy that needs to be subtracted from the envelope to vaporise and thermalise the ices of the planetesimals that will be mixed in the envelope. The vaporisation is given by the term − l_vapdm_z, l_vap being the latent heat of vaporisation of the ice per unit mass and dm_z the mass of ice remaining in the envelope. We assumed the ices to be water ice. We explain this assumption in more detail in Sects. 3.1 and 3.2.2. The thermalisation term is given by c_pΔTdm_z, being c_p the specific heat capacity of the ice and ΔT the change of temperature from 0 °C to the mean envelope temperature.

At time t + dt, all the mass of the planetesimals and gas that were accreted in the time elapsed (dt) are now part of the protoplanet, so that the total energy is the gravitational potential of the core and envelope, plus the internal energy of the envelope.

Our protoplanet is not a closed system, it is embedded in a protoplanetary disc, which pushes the outermost layer of the protoplanet, making work represented by − P_discdV. In addition, the protoplanet cools, radiating energy away by photons (− Ldt). Hence, ${\mathit{E}}_{\mathrm{tot}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{d}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{E}}_{\mathrm{tot}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{L}\mathrm{d}\mathit{t}\mathrm{-}{\mathit{P}}_{\mathrm{disc}}\mathrm{d}\mathit{V}\mathit{.}$(1)Regrouping terms, we find $\begin{array}{ccc}\mathit{L}& \mathrm{=}& \mathrm{-}\frac{\mathrm{d}{\mathit{E}}_{\mathrm{g,core}}}{\mathrm{d}\mathit{t}}\mathrm{-}\frac{\mathrm{d}{\mathit{E}}_{\mathrm{env}}}{\mathrm{d}\mathit{t}}\mathrm{+}{\mathit{e}}_{\mathrm{gas,acc}}\mathit{Ṁ}\mathrm{gas}\mathrm{-}{\mathit{P}}_{\mathrm{disc}}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{t}}\\ & & \mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\mathrm{(}{\mathit{l}}_{\mathrm{vap}}\mathrm{+}{\mathit{c}}_{\mathrm{p}}\mathrm{\Delta }\mathit{T}\mathrm{)}\mathit{Ṁ}\mathit{z}\mathit{.}\end{array}$(2)The first term is the change in the gravitational potential energy of the core, the second is the change in the total energy of the envelope, the next two are the surface terms explained above, and the last are the “ice heating” terms, also explained above. We call β the mass fraction of planetesimals that sinks to the core, that is why a factor 1−β (what remains in the envelope) is needed in front of the “ice heating” terms.

We develop below the first two terms, which are always at least two orders of magnitude larger than the others. To compute the first term (which we call hereafter L_core), we made use of the fact that the gravitational potential energy of a sphere of uniform density is ${\mathit{E}}_{\mathrm{core}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{5}}\frac{\mathit{G}{\mathit{M}}_{\mathrm{core}}^{\mathrm{2}}}{{\mathit{R}}_{\mathrm{core}}}\mathrm{·}$(3)Differentiating with respect to time and using the fact that the core density is constant, we can readily find ${\mathit{L}}_{\mathrm{core}}\mathrm{=}\mathrm{-}\frac{\mathrm{d}{\mathit{E}}_{\mathrm{g,core}}}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{core}}\mathit{Ṁ}\mathrm{core}}{{\mathit{R}}_{\mathrm{core}}}\mathrm{·}$(4)We assume throughout this work a constant density for the core, an usual practice in most planet formation works. If the core were compressible, an additional term in L_core should be included, reflecting the change in core gravitational energy associated with the change in density. Nevertheless, since the core is in a non-gaseous phase, that term is not relevant.

If we call Ṁ_Z the total accretion rate of solids, then Ṁ_core = βṀ_Z since β is the mass fraction of planetesimals/pebbles that sinks to the core.

We express the second term as L_gas: ${\mathit{L}}_{\mathrm{gas}}\mathrm{=}\mathrm{-}\frac{{\mathit{E}}_{\mathrm{env}}\mathrm{(}\mathit{t}\mathrm{+}\mathrm{d}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{E}}_{\mathrm{env}}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{·}$(5)The problem with computing L_gas is that the energy of the envelope at time t + dt is not known before obtaining the structure at t + dt. Therefore, we must make a guess of E_env(t + dt) to be able to have L(t + dt) for computing the structure at t + dt. This guess was performed following Mordasini et al. (2012) and Fortier et al. (2013). These authors assumed a similar functional form for the total energy of the envelope as for the gravitational potential energy of the core of uniform density. Hence, ${\mathit{E}}_{\mathrm{env}}\mathrm{=}\mathrm{-}{\mathit{k}}_{\mathrm{env}}{\mathit{M}}_{\mathrm{env}}\mathit{g,}$(6)where g is a mean gravitational potential, defined by $\mathit{g}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{G}\mathrm{(}{\mathit{M}}_{\mathrm{core}}\mathit{/}{\mathit{R}}_{\mathrm{core}}\mathrm{+}{\mathit{M}}_{\mathrm{P}}\mathit{/}{\mathit{R}}_{\mathrm{P}}\mathrm{)}$.

For a given converged structure at time t, we know the total energy of the envelope at this time. Hence, Eq. (6) defines k_env at time t, which is then used for the total energy of the envelope at time t + dt. This assumption is quite good during most of the growth of the planet, since k_env practically does not vary from one timestep to another. To test this, we ran simulations using a predictor-corrector scheme: we first computed the structure of the envelope using the k_env from the previous time step, k_env(t−dt), then computed the resulting k_env of the actual time step, k_env,corr, and finally recomputed the structure of the planet at time t using k_env = 0.5 (k_env(t−dt) + k_env,corr).

2.2. Envelope enrichment and iteration scheme

For the initial conditions, we assumed that the gas in the disc, and therefore the planetary envelope, is made of hydrogen and helium. This is justified by the fact that we are forming planets beyond the ice line, which means that most of the metals have been condensed into solids², which exist in the disc in the form of embryos, planetesimals, or pebbles.

We started all simulations with M_core = 0.01M_⊕. The core grows at a given accretion rate Ṁ_Z (see Sect. 3.3 for the different schemes adopted to compute Ṁ_Z). When the core reached a threshold value of M_thresh, we assumed that the core continues to grow by an amount ΔM_core = βṀ_Zdt, while the amount ΔM_Z,env = (1−β)Ṁ_Zdt remains uniformly mixed in the envelope. The dependence of M_thresh on the planetesimal properties is discussed in Sect. 3.5.

At the beginning of each time step, an initial guess of envelope metallicity was assumed to compute the corresponding equilibrium envelope mass (the envelope mass depends on the core mass and radius, total luminosity, and envelope metallicity). To compute M_env we used the second numerical scheme described in Venturini et al. (2015): iteration on M_env for a given core mass and radius, and Z_env. The total luminosity was recalculated at each iteration on M_env because L_gas depends on M_env, as explained in the above subsection. The envelope metallicity was taken to be uniform throughout the envelope and was defined as ${\mathit{Z}}_{\mathrm{env}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{Z,env}}}{{\mathit{M}}_{\mathrm{env}}}\mathit{,}$(7)M_Z,env being the mass of volatiles that remains mixed in the envelope.

Once convergence in M_env was achieved, the envelope metallicity was updated to the new M_env. Afterwards, an additional iteration on Z_env was required to define the metallicity self-consistently (and thus guarantee mass conservation of metals).

3. Assumptions

3.1. Choice of β

We recall that β is the mass fraction of the incoming planetesimals that is assumed to reach the core (hence, the fraction 1−β remains homogeneously mixed in the envelope). For our nominal model, we took β = 0.5. This choice is based on the refractory to volatile content of comets (Lambrechts et al. 2014; Thiabaud et al. 2015). Of the material falling in, we assumed that only ices mix in the envelope and that the refractory sinks into the core. We assumed for simplicity that the ices are water ice because it is the main volatile species in comets (Kofman et al. 2015) and because of self-consistency with our equation of state (see Sect. 3.2.2). The choice of β = 0.5 is a strong assumption, but we expect the ices to be the first component of the incoming planetesimals to be vaporised in the envelope. Moreover, it has been shown that water and hydrogen remain well mixed in the interior of giant planets for pressures higher that 1 GPa and temperatures higher than 2000 K (Soubiran & Militzer 2015). At least for the mentioned ranges, the assumption of water being homogeneously mixed in the envelopes is therefore fair. However, these ranges do not necessary cover all the pressure-temperature values during formation (see Fig. 5 of Venturini et al. 2015), which can extend from P ~ 10^-2 Pa to P ~ 10 GPa and from T ~ 100 K to T ~ 10 000 K. The miscibility of water into hydrogen for low pressure and temperature needs to be tested in the future.

3.2. Envelope microphysics

3.2.1. Opacities

Venturini et al. (2015) assumed a total opacity resulting from contributions from dust and gas. We used dust opacities from Semenov et al. (2003) and a gas opacity from tables computed by Ferguson for arbitrary metallicity, based on the original calculations of Alexander & Ferguson (1994). The dust opacities from Semenov et al. (2003) are suited for protoplanetary discs and therefore do not consider grain growth and settling in a planetary envelope. This is taken into account in the analytical opacities calculated by Mordasini (2014), which we adopted in this work. The dust opacities of Mordasini (2014) have the additional advantage of being independent of the accretion rate of metals (an equilibrium among the deposition of grains, grain growth, and settling was found in those calculations), making these opacities independent of the envelope metallicity and therefore suited for any value of Z_env.

For the gas opacities, new calculations by Freedman et al. (2014) are suited for planetary envelopes and for a wide range of metallicities (0 ≲ Z_env ≲ 0.5). Compared to the opacities computed by Ferguson, they have the advantage of not needing to be extrapolated for T ≤ 1000 K, and an analytical fit is also available, which is very useful to reduce computational time. The domain of validity in pressure and temperature of these opacities is 1 μbar to 300 bars and 75 K to 4000 K. These ranges do not reach the high temperatures and pressures that are typical at the surface of the core, but in practice, the envelopes are always convective in this domain of high pressure and temperature. Therefore, even taking a constant extrapolation in these regions does not affect the computation of the internal structure.

It is well known that opacities drastically affect the timescale on which a giant planet forms (Pollack et al. 1996; Hubickyj et al. 2005; Ikoma et al. 2000). We therefore ran comparison tests using different extreme combinations of dust and gas opacities. On the one hand, we took the combination of Semenov and Ferguson, and on the other that of Mordasini and Freedman. The results for the crossover mass are summarised in Table 1 for an envelope of solar composition and a constant accretion rate of metals of 10^-6M_⊕/yr. The notorious decrease in the crossover mass when using the new opacities mainly arises because the dust opacities from Mordasini (2014) are much lower than those from Semenov et al. (2003; whose values are very similar to the interstellar ones). We repeated this comparison test in Sect. 4.4 taking envelope enrichment into account. For all results shown before, we adopted the new opacities of Mordasini for the dust and the opacities reported by Freedman for the gas.

3.2.2. Equation of state

We assumed in this work that the volatile content of the planetesimals (or pebbles) that are destroyed while traversing the envelope is what remains well mixed in the envelope, thanks to ice sublimation. We also assumed that this volatile component is entirely made of water. This is justified mainly by two reasons. First, as we mentioned before, water is thought to be the main volatile molecule in planetesimals. Second, there is no accurate equation of state (EOS) available in the literature for an arbitrary mixture of volatiles that covers the wide ranges of temperature and pressure in planetary interiors. The only EOS suited for this purpose is that of water.

We therefore adopted an EOS for a mixture of H, He, and H₂O that takes degeneracy due to free electrons into account. For the H-He component we implemented the Saumon et al. (1995) EOS, and for the H₂O component we used an improved version of ANEOS (Thompson 1990). An important drawback in the standard version of ANEOS is the assumption that the substance is monoatomic in the gas phase. We corrected this for water (Benitez et al., in prep.) to include the proper degrees of freedom, following the approach of Melosh (2007). We implemented this new version of ANEOS for our water component. The EOS of the mixture of H-He with water was obtained, as in Baraffe et al. (2008), by means of the additive volume rule, which has been proven to yield adequate results for mixtures of H, He, and H₂O (Soubiran & Militzer 2015).

3.3. Accretion rate of solids

We adopted different models for the accretion rate of solids.

• 1.

  Constant accretion rate. This basic scheme allowed us to analyse the effect of envelope enrichment in its simplest form. For nominal results, we adopted a Ṁ_Z = 10^-6M_⊕/yr. Other values of accretion rates are tested in Sect. 4.5.

• Accretion rate following Pollack. We implemented our enrichment code in a Pollack-scheme (Pollack et al. 1996, hereafter P96) planetesimal accretion rate (see Sect. 5 for an explanation). The purpose of this was to show the effect of envelope enrichment in a more realistic formation scenario (although it is still not very accurate, since the proper excitation of the planetesimals in the oligarchic growth regime is not included as it is in Fortier et al. 2007, 2013). The crucial parameters assumed for this scenario are summarised in Table 2. The different choices of Σ₀ for the different sets of opacities were set to obtain similar formation timescales, as is explained later in Sect. 5. For the capture radius we used the prescription of Inaba & Ikoma (2003; see also Guilera et al. 2011).

  Table 2

  Parameters used for results where the Pollack-scheme is implemented and nominal opacities (Mordasini and Freedman) are adopted.

3.4. Boundary conditions

To define the radius of the protoplanet, we used the definition proposed by Lissauer et al. (2009), which takes flow circulation from 3D hydrodynamical simulations into account, ${\mathit{R}}_{\mathrm{P}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{P}}}{\mathrm{(}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathit{G}{\mathit{M}}_{\mathrm{P}}\mathit{/}{\mathit{R}}_{\mathrm{H}}\mathrm{)}}\mathit{,}$(8)R_P being the planet radius, M_P the planet mass, c_s the sound speed of the gas, and R_H the Hill radius.

For nominal results, we adopted boundary conditions suited for the present position of Jupiter. The values of these boundary conditions are given in Table 3. The results are very insensitive to the boundary conditions (as long as the planet is beyond the ice line and the accretion rate of solids is constant), therefore we do not show results for other boundary conditions.

We set the conservative number of 3.2 g/cm³ for the core density. This number should be larger when we consider that only rocky material sinks to the core, but we tested that increasing the core density has a very weak effect on the overall evolution.

3.5. Choice of M_thresh

The amount of material that is deposited in the envelope when a planetesimal is disrupted is a complicated function that depends upon many planetesimal properties (mass, density, material strength, etc.). Therefore an accurate knowledge of planetesimal properties is needed when the mass deposition is to be computed self-consistently. Since this is still unknown, we encoded this lack of information in the parameter M_thresh, which represents the lowest core mass that can bind an envelope massive enough to destroy the incoming planetesimals completely. For nominal results, we took M_thresh = 2M_⊕. We also test the effect of considering M_thresh = 0.5,1, and 4 M_⊕ in Sect. 4.2. To have an order of magnitude in mind, a core mass of 1 M_⊕ would have an envelope massive enough to disrupt completely icy planetesimals of ~100 m (Podolak et al. 1988), a core of 2 M_⊕ would disrupt icy planetesimals of ~1–10 km, and a core of 4 M_⊕ would disrupt icy planetesimals of ~10–100 km before they reach the core (Mordasini, priv. comm.; based on calculations presented in Mordasini et al. 2006). These values would be upper limits, planetesimals smaller that these would be fully disrupted in the envelope as well. Moreover, a core of 0.5 M_⊕  has an envelope mass equivalent to two Earth atmospheres in our models, which means that pebbles would surely sublimate but km size planetesimals would not.

4. Results for a constant accretion rate of solids

4.1. Differences in growth between standard H-He envelopes and enriched ones

Figure 2 shows the growth of a planet assuming Ṁ_Z = 10^-6M_⊕/ yr for a standard case where all the solids go to the core (envelope composed of H and He), and for a case where envelope enrichment is taken into account. Throughout this work, we refer to the first case as the non-enriched case and to the latter as the enriched case. The first remarkable fact we observe is that the time needed to form a giant planet is reduced when envelope enrichment is taken into account. If we use the classical definition of crossover time (time when the core mass equals the envelope mass, t_cross, Pollack et al. 1996) as the characteristic time to form a giant planet, then we find t_cross = 3.8 Myr for the enriched case against a t_cross = 7.0 Myr for the non-enriched case.

We analysed the enriched case in more detail and show the evolution of the envelope metallicity in Fig. 3. When the core mass reaches a threshold value M_thresh = 2M_⊕, half of the metals (water in our assumptions) remains uniformly mixed in the envelope and the other half is deposited in the core (see changing slope in M_core in Fig. 2). Since the envelope is still quite thin at this M_core (see Table 4), the envelope metallicity rises very rapidly at the beginning of the enrichment. In 3.5 × 10⁵ yr, Z_env reaches its maximum and then dilution slowly starts. There is an important fact to remark: rapid gas accretion is not triggered immediately as a consequence of envelope enrichment. This is well illustrated in Fig. 4, where we observe that the accretion rate of H-He remains lower than that of metals until t ~ 3 Myr. In the bottom panel of the same figure, we show the accretion timescale of H-He (defined as τ_HHe = M_HHe/Ṁ_HHe) also as a function of time for M_thresh = 2M_⊕. Immediately after the onset of enrichment, the accretion timescale of H-He is very short, but it starts to increase until it reaches a maximum at t ≈ 3.5 Myr.

Analysing the timescale to accrete H-He allows us to infer when runaway of gas starts. When this timescale increases, it means that gas accretion is slowing down. Conversely, when τ_HHe starts to decrease, HHe-accretion accelerates. We therefore define the onset of the runaway of gas as the time when τ_HHe reaches its maximum³, and we denote this time as t_run.

The time elapsed between the onset of enrichment (t ≈ 2 Myr) and the onset of runaway of gas is 1.5 Myr. During that time, the envelope metallicity is always higher than 33% (Fig. 3). This scenario implies that a very high gas-to-core ratio (GCR = M_env/M_core) is achievable before runaway gas accretion is triggered. Figure 5 depicts this more clearly. It shows the accretion timescale of H-He as a function of the GCR. For the non-enriched case, the maximum of τ_HHe occurs at a GCR ~ 0.1. For GCR higher than this, runaway of gas sets in. Therefore, protoplanets with a GCR higher than 0.1 are expected to become gas giants. In other words, low- and intermediate-mass planets with a high GCR are difficult to explain in the framework of the standard core accretion model, where envelopes are assumed to remain non-enriched (Lee & Chiang 2016). Nevertheless, when envelope enrichment is taken into account, a much higher CGR (~0.8) can be achieved before the onset of runaway (Fig. 5).

In this context we may wonder why giant planets are formed faster when envelope enrichment is taken into account. When envelope enrichment sets in, the structure of the envelope changes radically: the mean molecular weight increases very rapidly as a result of the increase in Z_env. Since the boundary conditions are the same, an increase in the mean molecular weight translates into an increase in the density profile. Therefore, the self-gravity of the envelope is stronger, and the planet becomes more prone to contracting and accreting more gas.

4.2. Dependence on M_thresh

In this section we test the effect of changing M_thresh. As we stated in Sect. 3.5, this would correspond to considering different sizes (and/or compositions, but mainly sizes) for the disrupted planetesimals or pebbles. Table 4 summarises some aspects of the simulations for values of M_thresh = 0.5, 1, 2, and 4 M_⊕. The first row shows the mass of the envelope when envelope enrichment starts. This information is important when trying to link the value of M_thresh with the corresponding maximum size of planetesimals to be fully disrupted.

In Fig. 6 (top) we show the evolution of the GCR and of the total mass fraction of H-He for the different M_thresh. The black dots indicate the time of the onset of the runaway of gas (maximum in τ_HHe) and the numbers above are the corresponding value of GCR at this time.

It is interesting to note that the lower M_thresh, the higher the GCR that can be achieved before runaway of gas. Likewisse, the period during which enrichment takes place but fast accretion of gas still does not occur, tends to be longer the lower M_thresh (last row of Table 4). This means that the smaller the planetesimals, the more likely it is that at the time the disc disappears, a small planet with high GCR remains.

Another surprising fact is that the total mass fraction of H-He (which we denote f_HHe) at the onset of runaway of gas is always ~30% for the enriched cases (see Table 5), regardless of the value of M_thresh. This would be the highest attainable value for planets that do not become gas giants, lower values of f_HHe are also possible, as Fig. 6 depicts. This means that envelope enrichment seems to be a natural way to explain the formation of low-mass, low-density planets (mini-Neptunes) and also the recently reported “super-puffs”: low-mass planets that have a bulk composition of H-He of presumably ≳20% of H-He by mass (Lee & Chiang 2016; Lopez & Fortney 2014). We emphasise that the claim that envelope enrichment is a more natural scenario to explain the formation of objects with ~20% of H-He is not just due to the fact that these values are reached before the onset of runaway of gas, but also because these high amounts of H-He last longer when envelope enrichment is included, as Fig. 6 shows.

An interesting aspect in Table 5 is that the total mass of the planet at the onset of runaway of gas grows with increasing M_thresh but decreases for M_thresh = ∞ (the non-enriched case). This is related to the fast accretion of H-He that is triggered when enrichment sets in. When we analyse the core masses of the different M_thresh at the onset of runaway (Table 5), we note that for M_thresh = ∞ the core is larger than for M_thresh = 4 M_⊕. However, because the mass fraction of H-He is much greater in the enriched cases that in the non-enriched one, the total planetary mass of the M_thresh = 4 M_⊕ case is higher than the M_thresh = ∞ one.

Finally, in Fig. 7 we show the total mass of solids of the planet as a function of the planetary mass acquired during growth. With the low opacities of our nominal model (Mordasini and Freedman), the highest mass of solids that can be attained during growth is M_Z ≈ 7M_⊕ even in the extreme case of non-enrichment. Higher accretion rates can lead to higher core masses (or total mass of solids, in general). We show this in Sect. 4.5.

4.3. Smoothing the transition in β

It might be thought that since in reality there should be a distribution of planetesimal sizes, there is not one exact value of M_thresh for which the envelope starts to become enriched, but a range of M_thresh. We have tested the effect of adopting different M_thresh in the previous section, but considering a distribution of planetesimal sizes also means that β would not change abruptly at a given M_thresh, but would transition for a range of M_thresh.

We therefore performed the following test. We assumed β = 1 for M_core ≤ 1 M_⊕, β = 0.5 for M_core ≥ 2M_⊕, and a linear variation of β in between. The growth of the planet is shown in Fig. 8 together with the evolution of envelope metallicity. Even though the transition in β is smoother than in the previous cases, the evolution is very similar to the case M_thresh = 2M_⊕. The envelope metallicity still grows quite abruptly because the envelope at M_core = 1M_⊕ is quite thin. Since the envelope at M_core = 1M_⊕ is thinner than at M_core = 2M_⊕, in this case of smoother transition of β, Z_env reaches a higher maximum than for M_thresh = 2M_⊕. The overall formation time is practically the same as for M_thresh = 2M_⊕, however. We therefore conclude that assuming a sharp transition of β at a given M_core is not a relevant simplification.

4.4. Other opacities

As we mentioned in Sect. 3.2.1, the time needed to form a giant planet strongly depends on the choice of opacities. Even though in our nominal results we used the latest opacities published, the opacities might not exactly follow these low values. For instance, the recondensation of upstreaming gas into grains would increase the opacities, and this is not taken into account for the computation of Mordasini (2014) nor Freedman et al. (2014).

To test how choosing higher opacities would affect our results, we ran a simulation with the other extreme set of opacities that we introduced in Sect. 3.2.1: dust opacities from Semenov et al. (2003) and gas opacities from Ferguson (based on Alexander & Ferguson 1994, but including arbitrary metallicities). We considered here the enriched case with Ṁ_Z = 10^-6M_⊕/yr and M_thresh = 2M_⊕.

The results of the planet growth are shown in Fig. 9. The total mass of the planet is plotted with a colour bar that indicates the timescale to accrete hydrogen and helium (τ_HHe). We note that the maximum of this timescale in this case is reached after crossover mass. It occurs for a total mass of ~12 M_⊕, of which M_core ~ 5.4M_⊕ and M_env ~ 6.6M_⊕. At that time, half of the mass of the envelope is H-He and the other half is water. This means that when τ_HHe reaches it maximum (t ≈ 9.5 Myr), we have a planet with a total mass and total metallicity (Z_tot ~ 72%) similar to Uranus.

The time shown here to form an ice giant seems quite long for a typical disc lifetime, but this is a consequence of the choice of the planetesimal accretion rate. For instance, if instead of using Ṁ_Z = 10^-6M_⊕/yr we were to choose Ṁ_Z = 3 × 10^-6M_⊕/yr, then t_run ≈ 4 Myr, and the overall evolution in terms of mass and metallicities is quantitatively similar to the results shown in Fig. 9. For a higher accretion rate of solids, the total planetary mass is that of Neptune at t = t_run, with a total metallicity of Z_tot ~ 70%. Hence, for the high opacities used in this section, higher accretion rates lead to shorter formation times, yielding still, as a typical output, a Neptune- or Uranus-like planet.

4.5. Other accretion rates of solids

For completeness, we test in this section the effect of considering other accretion rates of solids. We ran simulations with envelope enrichment and the microphysics of our nominal model (low opacities) but using Ṁ_Z = 10^-5M_⊕/yr and 5 × 10^-5M_⊕/yr instead of Ṁ_Z × 10^-6M_⊕/yr. We summarise the results in Tables 6 and 7.

Increasing the accretion rate of solids increases the core luminosity. Thus, the higher the accretion rate of solids, the hotter the envelope, which prevents gas accretion. Therefore, by increasing Ṁ_Z, we expect to obtain planets with higher planetary masses before the onset of runaway of gas. We note that with an accretion rate as high as 5 × 10^-5M_⊕/yr (more typical of pebble accretion than of planetesimal accretion), when τ_HHe reaches its maximum, the total mass of the planet is M_P ≈ 11M_⊕, with a core of M_core = 5M_⊕ and total fraction of H-He of ~30%. Clearly, this might be another mechanism to form intermediate-mass planets. This is exactly the mechanism proposed by Lambrechts et al. (2014) to form Uranus and Neptune through pebble accretion. The problem with this scenario is that the formation timescale is very short (~0.2 Myr), which makes it hard to justify why a planet that reached the onset of runaway of gas in 0.2 Myr would not continue accreting gas to become a gas giant. We can only assume that the embryo was formed late in the disc, but then this turns into the classical time fine-tuning problem: the embryo has to form when the disc is fading, but before the gas is entirely depleted. Therefore, this is an unlikely scenario to explain the formation of Neptunes. We discuss the formation of ice giants in the framework of pebble accretion in more detail in Sect. 6.3.

4.6. Water condensation

In all the simulations presented in this work we have included the effect of water condensation (as it should be). Venturini et al. (2015) showed that water condensation decreases the adiabatic gradient, and this plays a relevant role in diminishing the critical core masses to very low values, especially for very high envelope metallicities (Z_env ≳ 0.6).

We have repeated some of the simulations presented before, but excluding the effect of water condensation in the computation of the adiabatic gradient. To do this, we ran simulations using CEA for the equation of state (Gordon et al. 1994) as in Venturini et al. (2015), because with this package it is possible to switch the effect of water condensation off when the adiabatic gradient is computed. The results presented in this section were therefore all performed with CEA (those where water condensation is included, and those where it is not).

In the nominal case where we used Mordasini and Freedman opacities, the difference between including and excluding water condensation does not affect the evolution much because with these low opacities the outer layers of the envelope are radiative, so that the structure in these layers is independent of the value of the adiabatic gradient.

For higher opacities and/or accretion rates of solids, the envelopes are more prone to be convective, so that for these cases water condensation plays an non-negligible role. For instance, considering the nominal case but increasing the dust opacities by a factor of 300, and taking an accretion rate of solids of Ṁ_Z = 3 × 10^-6M_⊕/yr, we find that the total formation time (when the planet reaches a mass of ~40 M_⊕) is 5 Myr for the case where water condensation is taken into account, and 6 Myr when it is not. The maximum Z_env attained for the former case is 0.75, and for the latter it is 0.85. This is consistent with what we found in Venturini et al. (2015): water condensation makes envelopes thicker for the same core mass, which explains that a lower envelope metallicity is reached and gas is accreted faster than when water condensation is neglected.

5. Results with an accretion rate of metals as in P96

In this section we study the effect of envelope enrichment in the framework of a Pollack et al. (1996) accretion scheme. The main difference with fixing a given accretion rate of planetesimals is that now the accretion rate of planetesimals depends on the total mass of the protoplanet and that the availability of planetesimals to be accreted depends on the initial amount in the neighbourhood of the embryo (the initial surface density of solids).

We first implemented an accretion rate of solids following Pollack with initial surface solids of Σ₀ = 4 g/cm² and nominal opacities (Mordasini and Freedman). It is important to note that the value of Σ₀ is arbitrary and was chosen so as to obtain a formation timescale for the non-enriched case of 8 Myr, as in Pollack et al. (1996). However, in our calculations this value has no influence on the opacities even though it is likely that dust-to-gas ratio and total disc mass would change opacities. Recent results by Mordasini (2014) and Ormel (2014) tend to show that a least regarding the dust, an increase in the accretion rate of solids (or Σ₀ in the context of the P96 accretion scheme) does not increase the dust opacity because the higher the flux of dust received by the protoplanet, the faster the coagulation and settling within the envelope.

The results for the opacities and surface density of solids are illustrated in the upper panel of Fig. 10. The enriched case assumes M_thresh = 3 M_⊕ and β changing at this core mass from 0 to 0.5. The choice of M_thresh was made to ensure that the envelope mass is appropriate to dissolving icy planetesimals of 100 km radius (which is the choice of planetesimal size for the P96 scheme, see Table 2). In this case, when phase 1 ends, the core mass is 2.9 M_⊕, which means that M_thresh is reached during phase 2 of P96. We note that the decrease in formation time, compared to the non-enriched case, is of a factor ~2.

We also tested the effect of using the combination of high opacities (Semenov and Ferguson) in this P96 accretion scheme (Fig. 10, bottom). To reach the same overall formation time as before, we chose an initial surface density of solids of Σ = 10 g/cm² for the non-enriched case (as in the baseline model of P96). Because of the high opacities and the high accretion rate of solids achieved in phase 1, the envelope required to destroy icy planetesimals of 100 km size now corresponds to M_thresh = 11M_⊕. We note that in this case, the formation time is reduced by a factor ~6 when envelope enrichment is included. The difference with the low-opacity case is that now M_thresh is achieved just before phase 1 ends. This means that when envelope enrichment takes place during phase 1, there are feedback mechanisms acting that favour more rapid gas accretion. In order to understand these mechanisms, we proceed to analyse the behaviour of the solid accretion rate for the high-opacity scenario.

The upper panel of Fig. 11 shows the change in the accretion rate of solids as a function of time for the enriched and non-enriched cases of the high-opacitiy simulation (Fig. 10, bottom). A bump in the accretion rate of solids occurs when the envelope begins to be enriched (at M_core = 11M_⊕). The evolution of the capture radius (Fig. 11, bottom) shows that at this moment, it also grows considerably. While the outer boundary conditions are the same as in the non-enriched case, the mean molecular weight of the envelope increases as accretion proceeds. Therefore, the density increases, so that the planetesimals are more efficiently slowed down when crossing the atmosphere. This increases the capture radius, which translates into a higher accretion rate.

In principle, the same effect occurs when enrichment sets in during phase 2. However, in that case, the accretion rate of solids is much lower and the availability of the protoplanet to increase the accretion rate of planetesimals is limited to the amount of planetesimals that can enter the feeding zone at each time step. Hence, it makes sense that the increase in the accretion rate of solids is lower than when enrichment starts in phase 1.

6. Implications on the formation of different types of planets

6.1. Formation of gas giants

We have shown that by including the effect of envelope enrichment during the growth of a planet, the time needed to form a gas giant is shorter than in the standard case of no envelope enrichment. In the classical context of planetesimal accretion, the formation of gas giants is challenging when the oligarchic growth regime is taken into account (Fortier et al. 2007, 2013) and even more challenging when planetesimal fragmentation is included (Guilera et al. 2014).

Fortier et al. (2013) showed (for the classical scenario of H-He envelopes) that the formation of giant planets required small planetesimals (~100 m). They found that km size planetesimals (and larger) were not sufficiently damped by gas drag, and therefore, because of their high relative velocities, were not efficiently accreted by the protoplanet. We have shown that envelope enrichment reduces the formation timescales. Therefore, including this effect in simulations that consider oligarchic growth might help the formation of gas giants from km size planetesimals. The reason for this is not related to envelope enrichment playing any role in the damping of planetesimals eccentricity, but to the fact that it accelerates the growth of the planetary envelope, producing a positive feedback in the increase of the capture radius. We showed in Sect. 5 that envelope enrichment can reduce the time needed to form a giant planet by a factor of ~5. The same effect should reduce the formation time of giants when oligarchic growth is taken into account. This scenario needs to be tested in the future.

In the framework of pebble accretion, the problem with gas giants is the opposite of the problem with planetesimals: giant planets are extremely easy to form. In this context, envelope enrichment would exacerbate the situation. It is important to stress, however, that the current status of pebble accretion (in giant planet formation) is that the embryos must be arbitrarily allowed to start to grow after the disc is significantly evolved (Bitsch et al. 2015). Otherwise, most of the embryos growing beyond the ice line end up as gas giants. The feasibility of pebble accretion to reproduce (statistically) the giant planets observed still needs to be proven, and envelope enrichment should be included for those models to be physically acceptable.

6.2. Formation of mini-Neptunes and Neptunes

Another important result we obtained is that if the disc disappears before the onset of runaway of gas, the type of planets we can form with envelope enrichment is quite different than in the non-enriched counterpart. With envelope enrichment we expect failed giants to have high GCRs (~1), high mass fractions of H-He (up to ~30%), and relatively high envelope metallicities (≳40%). This result is very interesting because one ongoing problem with the classical picture of the core accretion model (without envelope enrichment) is the formation of low- (~1–10 M_⊕) and intermediate-mass (~10–20 M_⊕) objects with high contents of H-He and high total metallicity. This has been a recurrent problem in the formation of Uranus and Neptune in the solar system, for instance, because even if the disc disappears at the moment when the mass of the planet is in the appropriate range (15–20 M_⊕), a core-dominated planet with ~20% in mass of H-He is difficult to obtain (Helled & Bodenheimer 2014). The same regarding mini-Neptunes: the highest GCRs expected from simulations of H-He envelopes typically is ≲10% (Lee & Chiang 2016).

Our nominal model is able to produce mini-Neptunes, but not Neptunes. This is because of the low opacities. Gas opacities from Freedman et al. (2014) include the effect of a changing metallicity (and therefore, increase with increasing Z_env), but the dust opacities of Mordasini (2014) are still very low. Therefore, with these opacities, crossover masses are low (even without envelope enrichment, see Table 1). It is not possible with these low opacities to obtain the static structure of a ~15 M_⊕ planet with a core of ~10 M_⊕ (regardless of the envelope metallicity). Nevertheless, a Neptune-type planet can be formed with envelope enrichment if higher opacities are invoked, as we showed in Sect. 4.4. This shows the importance of continuously improving the opacity models. The formation of clouds could be one mechanism for higher opacities in the envelope of protoplanets, and we have shown that water clouds can be present in the atmospheres of protoplanets formed beyond the ice line (Venturini et al. 2015). Efforts to include these physics into the computation of opacities need to continue.

6.3. Formation of ice giants through pebble accretion

In Sect. 4.5 we showed that when high accretion rates are considered, such as are invoked in the context of pebble accretion, the formation of ice giants was difficult because formation times were too short. In principle, the situation could be improved considering higher opacities (as we showed in Sect. 4.4 for planetesimal-like accretion rates). However, we note that when we ran the non-enriched case with the high opacities of Semenov and Ferguson and a solid accretion rate of 5 × 10^-5M_⊕/yr, runaway of gas is triggered at t_run ≈ 0.45 Myr. This means that if we were to consider envelope enrichment with these high opacities, the formation times of giant planets would still be too short (of ~0.1 Myr).

Lambrechts et al. (2014) claimed that they were able to explain the formation of ice giants in the context of pebble accretion. Their argument was that pebble accretion rates are so high that they provide a means for increasing the critical core mass to very high values (Stevenson 1982; Papaloizou & Terquem 1999) so that a Neptune-mass planet is in this context still sub-critical, and therefore will not accrete large amounts of gas.

Moreover, they determined an isolation mass for pebbles, which has values of M_iso ≳ 20M_⊕ for distances from the star a ≳ 5 AU. The isolation mass is lower the closer the planet is to the star, therefore they claimed that Jupiter and Saturn reached isolation mass, which caused the onset of runaway of gas (halting solid accretion promotes gas accretion, Ikoma et al. 2000); whereas Uranus and Neptune did not.

Lambrechts et al. (2014) showed that they were able to compute the correct structure of the four giant planets of the solar system (in terms of total mass and total metallicity). However, they did not show any time evolution. Despite the increase in critical core mass, if this mass is reached in 0.1 Myr, then fast gas accretion from this stage on will be inevitable. Perhaps a way to avoid runaway of gas is to consider that planets accrete both pebbles and planetesimals: when pebble isolation mass is reached, the planet could continue to accrete planetesimals that provide the necessary luminosity to delay runaway of gas (as in phase 2 of Pollack et al. 1996). This will be the subject of a future work.

7. Predictions and comparison with observations

7.1. Solar system

We considered in this work enrichment by icy planetesimals. Therefore our predictions are relevant for planets formed beyond the ice line. Planetesimal disruption inside the ice line could take place as well, but given that the composition would be more silicate-iron rich, the thermal ablation of planetesimals would either require very small planetesimal sizes or thick envelopes. The last option would mean that enrichment would probably take place after the planet has entered the phase of runaway of gas. In addition, even if thermal ablation and mechanical disruption occur, the fate of a mixture of silicates with hydrogen is expected to be miscible for T ≳ 10 000K (Wilson & Militzer 2012). Water has been proven theoretically to remain homogeneously mixed throughout the envelope of giant planets (Soubiran & Militzer 2015), but for refractories, the fate of the fragments of planetesimals that are composed of them is less certain, given the immiscibility just mentioned. If they tend to sink to the core, then the effect of envelope enrichment would not operate during formation, meaning that critical core masses would have the usual values of ~10–20 M_⊕. This would imply that forming giant planets from rocky planetesimals (i.e. inside the ice line) would be difficult.

For the formation beyond the ice line, our results show that we can expect both gas giants and Neptune-like planets to form through envelope enrichment, as explained above.

We have shown that by including envelope enrichment by icy planetesimals during formation, we expect the planets formed to present a decrease in envelope (and total) metallicity with increasing total mass. The decrease in envelope metallicity with planetary mass was pointed out empirically by Kreidberg et al. (2014). The decrease in total metallicity with planetary mass is a natural outcome of core accretion (e.g. Alibert et al. 2005; Mordasini et al. 2009). To illustrate this behaviour better, we show in Fig. 12 the result of planet formation for different choices of opacities and accretion rates of planetesimals. In all cases, the envelope starts to become enriched when it reaches a mass of M_env = 5 × 10^-3M_⊕, which corresponds to the full disruption of icy planetesimals of approximately 1–10 km-size (Sect. 3.5). We have overplotted the total metallicities of the giant planets of the solar system for reference. This figure suggests that rather high opacities are preferable to explain the giants of the solar system. Of course, the behaviour of the curves depends somewhat on the choice of parameters. Still, when we analysed the nominal model (Mordasini and Freedman opacities) for different M_thresh, for instance, we saw that even the extreme case of non-enrichment could not lead to total high-Z content of more than ~7 M_⊕ (Sect. 4). It is clear therefore that with these low opacities we cannot explain Jupiter, which has at least 15 M_⊕ of heavy elements (Baraffe et al. 2014). This situation can of course be solved, at least for Jupiter, by increasing the accretion rate of solids (Sect. 4.5). However, for intermediate-mass planets, an increase in the solid accretion rate would cause formation times to be so short that they would likely become gas giants. The combination of low opacities and high accretion rates of solids could work for forming Jupiters, but forming Neptunes would be more challenging.

We note that the dust opacities from Mordasini (2014) are analytical and therefore many simplifications were considered to compute them. Perhaps the most important one is the assumption of a predominant grain size. Mordasini (2014) mentioned that if planetesimal ablation is important in all layers of the envelope, then the constant supply of small grains could increase the dust opacities. It is likely that this plays a relevant role in our enriched scenario.

7.2. Exoplanets

Our formation model assumes that envelopes become enriched in water during formation because icy planetesimals are expected to be water rich. Water has indeed been detected in all the atmospheres of the giant planets of the solar system, although in very small amounts, even less than what is expected based on the detection of other volatile molecules, such as CH₄ (Niemann et al. 1998; Irwin et al. 2014). The problem with measuring water in the outermost layers of the giant planets is that because of the low temperatures, this species is expected to have condensed deeper in, and therefore be present in the form of clouds.

In this sense, transiting exoplanets offer a better opportunity to detect water in their atmospheres: since the equilibrium temperatures are much higher (these planets are much closer to their central star than the giants of the solar system as a result of observational bias), water is expected to be present in the form of vapour. Of the 19 transiting planets whose spectra have been measured, ten present signatures of water vapour in their atmospheres (the remaining nine are thought to posses water as well, but their signature is probably obscured by clouds, Iyer et al. 2016). Of these ten exoplanets, nine have masses in the range of 0.5–2 Jovian masses, which means that these planets are expected to be hydrogen-helium rich. In other words, the amount of water present is expected to be low, as our results suggest (envelope metallicity should decrease with increasing planet mass).

Two works have reported precise water abundances of hot-Jupiter atmospheres: Kreidberg et al. (2014) for WASP-43b and Madhusudhan et al. (2014) for HD 189733b, HD 209458b, and WASP-12b. We have converted the abundances they give of water mixing ratios into mass fraction of water (Z_env in our models, see Appendix A for an explanation on this conversion) and plotted the predictions of our models together with these estimates in Fig. 13. We note that for WASP-43 b the predictions of both high- and low-opacity models work surprisingly well to explain the estimated abundance of water vapour. For the hot Jupiters reported by Madhusudhan et al. (2014), we find that all our models predict a higher mass fraction of atmospheric water, at least by a factor of 2. Madhusudhan et al. (2014) also claimed that the water abundances they found seemed too low and suggested that clouds might be obscuring some of the molecular features of the spectra (and clouds were not included in the retrieval model of that work). Madhusudhan et al. (2014) suggested alternatively that the atmospheres of these hot Jupiters might bear more carbon-rich species than oxygen-rich. It is important to remark that in our formation model, we assumed that all the volatiles deposited in the planet envelope were only water (because of the equation of state). In this sense, it might perfectly well be that we overestimate the amount of water because our model neglects the possibility of initial amounts of CH₄, CO₂, and CO, for example. Nevertheless, more recent works (Sing et al. 2016; Benneke 2015) suggest that indeed, the non-inclusion of clouds in the retrieval modelling makes the water abundances of Madhusudhan et al. (2014) to be considerably underestimated.

Another interesting case to link our results with observations is the detection of water in HAT-P-11b (Fraine et al. 2014), the only Neptune-mass exoplanet where water vapour has been inferred from transmission spectroscopy. Fraine et al. (2014) reported from retrieval models that the atmosphere of HAT-P-11b is expected to have a metallicity of 1 to 700 times solar, with no real preference for a specific value within this range (Benneke, priv. comm.). Unfortunately, this wide range could imply an envelope metallicity from solar up to Z_env ~ 0.9. HAT-P-11b has a total mass of M_P ≈ 26M_⊕, which means that our models would predict a Z_env in the range of, approximately, 0.05–0.4 (see Fig. 12). This would correspond to values of “metallicity times solar” between 1 and 100. More precise measurements are needed to determine whether our models are able to predict the metal atmospheric content of this planet.

For the detection of water in exoplanets, it is finally important to remark that even though the sample of transmission spectra measured is still small, water is shown to be common in the atmospheres of planets. This reinforces a formation scenario beyond the ice line, with volatile-rich planetesimals being dissolved in the atmospheres during formation. Even when it might be found in the future that statistically the atmospheres are water poor, this could be a hint of an initial planetesimal composition that is more carbon rich (Madhusudhan et al. 2014). Therefore it is important to determine precisely other volatile abundances in addition to water.

8. Conclusions

We have performed the first self-consistent calculation of the growth of a planet including the effect of envelope enrichment due to the dissolution of icy planetesimals or pebbles. We implemented suitable equations of state and opacities taking different metallicities into account. Moreover, we considered two different sets of opacities to test the effect on our results, which we summarise below.

• Envelope enrichment notably accelerates the formation of gas giants. This is mainly a consequence of the increase in mean molecular weight of the envelope. The thinner the envelope (i.e. the smaller the planetesimal), the sooner envelope enrichment sets in and the shorter the time in which a giant planet is formed.

• When envelope enrichment is taken into account, low- and intermediate-mass planets (mini-Neptunes to Neptunes) can be formed, with total mass fractions of H-He up to 30%, this number is independent of the choice of opacities.

• Low opacities allow for the formation of mini-Neptunes, while high-opacities lead to the formation of Neptunes.

• High opacities are preferable for explaining the total mass and metallicity of the giant planets of the solar system.

• We were able to quantify the amount of volatile material remaining in the primordial atmospheres as a result of formation. This allowed us to compare our results with water abundances inferred from the transmission spectra of transiting exoplanets. We find good agreement for WASP-43b, but for the other three cases, the measured water abundance is lower than what is predicted by our models.

• We predict that envelope metallicity and total metallicity decrease with planetary mass.

---

Acknowledgments

We thank M. Ikoma for fruitful discussions and F. Benitez for providing the EOS tables and revising this manuscript. J.V. acknowledges financial support from the Swiss-based MERAC Foundation. This work has, in part, been carried out within the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge financial support from the SNSF.

References

Appendix A: Conversion of water mixing ratios to water mass fractions

It is a common practice in papers that estimate abundances of different compounds to give them in terms of volume mixing ratios. In particular, for estimates on water abundances derived from transmission spectroscopy measurements, the volume mixing ratio of water (m_r) is defined as (Madhusudhan et al. 2014) ${\mathit{m}}_{\mathit{r}}\mathrm{=}\frac{\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathrm{\right)}}{\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}}\mathit{,}$(A.1)where n(H₂O) and n(H₂) are the number density of water molecules and hydrogen, respectively.

Because of the molar masses of H₂O and H₂, in terms of mass, the ratio of water to hydrogen reads $\frac{{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}\mathrm{=}\frac{\mathrm{18}}{\mathrm{2}}\frac{\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathrm{\right)}}{\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}}\mathrm{=}\mathrm{9}{\mathit{m}}_{\mathit{r}}\mathit{.}$(A.2)In the formation models presented here, we defined the envelope metallicity as the mass fraction of water in the envelope. Thus, $\begin{array}{ccc}{\mathit{Z}}_{\mathrm{env}}& \mathrm{=}& \frac{{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{+}{\mathit{m}}_{\mathrm{He}}\mathrm{+}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}\\ & \mathrm{=}& \frac{{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\mathit{/}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}{\mathrm{1}\mathrm{+}\mathit{Y}\mathit{/}\mathit{X}\mathrm{+}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\mathit{/}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}\\ & \mathrm{=}& \frac{\mathrm{9}{\mathit{m}}_{\mathit{r}}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{Y}\mathit{/}\mathit{X}{\mathrm{)}}_{\mathrm{\odot }}\mathrm{+}\mathrm{9}{\mathit{m}}_{\mathit{r}}}\mathit{,}\end{array}$(A.3)where we assumed the ratio of helium to hydrogen to be solar. This last equation is the one we used to estimate the values of Z_env shown in Fig. 13 for the different hot Jupiters whose water mixing ratios has been determined (Madhusudhan et al. 2014; Kreidberg et al. 2014).

All Tables

All Figures