In situ accretion of gaseous envelopes on to planetary cores embedded in evolving protoplanetary discs

TL;DR

This study models the in situ formation and evolution of gaseous envelopes on planetary cores within protoplanetary discs, revealing conditions for runaway accretion and implications for the origins of Hot Jupiters and smaller exoplanets.

Contribution

It provides new insights into how planetary cores accrete gaseous envelopes in situ and their subsequent evolution, challenging the necessity of migration for certain exoplanet types.

Findings

Cores with 5 M⊕ do not undergo runaway accretion.

Runaway accretion occurs for cores ≥10 M⊕ at certain orbital distances.

Post-disc dispersal, planets cool and contract over Gyr timescales.

Abstract

The core accretion hypothesis posits that planets with significant gaseous envelopes accreted them from their protoplanetary discs after the formation of rocky/icy cores. Observations indicate that such exoplanets exist at a broad range of orbital radii, but it is not known whether they accreted their envelopes in situ, or originated elsewhere and migrated to their current locations. We consider the evolution of solid cores embedded in evolving viscous discs that undergo gaseous envelope accretion in situ with orbital radii in the range $0.1-10\rm au$. Additionally, we determine the long-term evolution of the planets that had no runaway gas accretion phase after disc dispersal. We find: (i) Planets with $5 \rm M_{\oplus}$ cores never undergo runaway accretion. The most massive envelope contained $2.8 \rm M_{\oplus}$ with the planet orbiting at $10 \rm au$. (ii) Accretion is more efficient onto $10 \rm M_{\oplus}$ and $15 \rm M_{\oplus}$ cores. For orbital radii $a_{\rm p} \ge 0.5 \rm au$, $15 \rm M_{\oplus}$ cores always experienced runaway gas accretion. For $a_{\rm p} \ge 5 \rm au$, all but one of the $10 \rm M_{\oplus}$ cores experienced runaway gas accretion. No planets experienced runaway growth at $a_{\rm p} = 0.1 \rm au$. (iii) We find that, after disc dispersal, planets with significant gaseous envelopes cool and contract on Gyr time-scales, the contraction time being sensitive to the opacity assumed. Our results indicate that Hot Jupiters with core masses $\lesssim 15 \rm M_{\oplus}$ at $\lesssim 0.1 \rm au$ likely accreted their gaseous envelopes at larger distances and migrated inwards. Consistently with the known exoplanet population, Super-Earths and mini-Neptunes at small radii during the disc lifetime, accrete only modest gaseous envelopes.

Tables

Table 1. Table 1: Disc and stellar model parameters. Note that the value for Σ g subscript Σ g \Sigma_{\rm g} applies only to the 1 × MMSN 1 MMSN 1\times{\rm MMSN} model and should be scaled appropriately for the 3 and 5 × MMSN 5 MMSN 5\times{\rm MMSN} models.

Parameter	Value
Disc inner boundary	0.05 $\mathrm{au}$
Disc outer boundary	40 $\mathrm{au}$
${\Sigma }_g$(1 $\mathrm{au}$)	$1731g/{\mathrm{cm}}^{-2}$
Metallicity	Solar
$\alpha$	$2\times {10}^{-3}$
${\alpha }_{\mathrm{active}}$	$1\times {10}^{-2}$
Stellar Mass	$1M_{⨀}$
$R_S$	$2R_{⨀}$
$T_S$	4280 K
$f_{41}$	10

Table 2. Table 2: Initial parameters for models evolved after disc dispersal

Model	Total mass	Core mass	Orbital radius
	$M_p/M_{\oplus }$	$M_{\mathrm{core}}/M_{\oplus }$	$a_p/au$
A1	5.216	5.00	0.2
A2	5.398	5.00	1.0
A3	5.776	5.00	5.0
B1	11.821	10.00	0.2
B2	13.380	10.00	1.0
B3	16.578	10.00	5.0

Table 3. Table 3: The first column gives the mean density in the final state with zero internal luminosity for each of the models listed in table ( 2 ). The second column gives the ratio of the radius after 2 Gyr and the end state radius.

Model	Final mean density	Radius at 2Gyr/ (Final radius)
	$\overline{\rho }$	$R_{\mathrm{surf}}/R_{\mathrm{surf},f}$
A1	0.674	1.37
A2	1.10	1.63
A3	1.13	2.02
B1	0.68	1.58
B2	0.83	1.83
B3	0.81	1.52

Equations

$$\frac{dP}{dR}=-\rho g$$

$$L_{\mathrm{rad}}=-\frac{64\pi \sigma R^2{T}^3}{3\kappa \rho }\frac{dT}{dR}$$

$${\nabla }_{\mathrm{rad}}=\frac{3\kappa L_{\mathrm{tot}}P}{64\pi \sigma GM{T}^4}$$

$${\nabla }_{\mathrm{ad}}={\left(\frac{\partial \ln T}{\partial \ln P}\right)}_{\mathrm{S}},$$

$$L_{\mathrm{conv}}=\pi R^2{C}_p{\Lambda }^2{\left[{\left(\frac{dT}{dR}\right)}_{\mathrm{S}}-\left(\frac{dT}{dR}\right)\right]}^{3/2}\times \sqrt{\frac{1}{2}\rho g\left|{\left(\frac{\partial \rho }{\partial T}\right)}_{\mathrm{P}}\right|}$$

$$R=R_{\mathrm{core}}={\left(\frac{3M_{\mathrm{core}}}{4\pi {\rho }_{\mathrm{core}}}\right)}^{1/3}.$$

$$R_{\mathrm{surf}}=\frac{2}{3}{\left(\frac{M_{\mathrm{p}}}{3M_{\ast }}\right)}^{1/3}{a}_{\mathrm{p}}$$

$$T_{\mathrm{surf}}=(T_{\mathrm{disc}}^4+T_{\mathrm{effp}}^4{)}^{1/4}$$

$${\tau }_{\mathrm{surf}}=\kappa ({\rho }_{\mathrm{disc}}{T}_{\mathrm{disc}}){\rho }_{\mathrm{disc}}{R}_{\mathrm{surf}}.$$

$$E={\int }_0^{M_{\mathrm{p}}}\left(U-\frac{GM}{R}\right)dM.$$

$$\frac{dE}{d{M}_{\mathrm{p}}}\frac{d{M}_{\mathrm{p}}}{dt}=-L,$$

$$P\kappa =g.$$

$$T_{\mathrm{surf}}=(T_{\ast }^4+T_{\mathrm{effp}}^4{)}^{1/4},$$

$$\left(dE/d{R}_{\mathrm{surf}}\right)\left(d{R}_{\mathrm{surf}}/dt\right)=-L_{\mathrm{tot}}.$$

$$\frac{P}{T}\frac{dT}{dP}=\frac{3\kappa P{L}_{tot}}{64\pi G{M}_p\sigma T^4}={\nabla }_{\mathrm{ad}}\equiv \frac{{\Gamma }_2-1}{{\Gamma }_2}.$$

$$\frac{1}{4-\eta }>\frac{{\Gamma }_2-1}{{\Gamma }_2}.$$

Full text

1 Introduction

Two distinct paradigms have been proposed to explain the origin of gas giant planets: core nucleated instability (Perri & Cameron, 1974; Mizuno et al., 1978; Stevenson, 1982) and gravitational instability (e.g. Helled et al., 2014, and references therein). In principle, given the presence of a solid core of sufficient mass, core nucleated instability can lead to the formation of gas giant planets at almost any orbital location within a protoplanetary disc during its lifetime.

For gravitational instability to operate it is required that the disc gas can cool efficiently on orbital time-scales. Accordingly this mode of planet formation is most likely to occur at large stellocentric distances ($\gtrsim 50\,{\rm au}$) where cooling times are short (e.g. Gammie, 2001). On the other hand, formation of more massive objects such as brown dwarfs may be preferred (see e.g. Kratter & Lodato, 2016).

In this paper, we examine the accretion of gas on to already formed cores located at a wide variety of stellocentric distances, ranging between $0.1\,{\rm au}$ and $10\,{\rm au}$ within evolving protoplanetary disc models, to determine where gas giant planets can form via in situ gas accretion, and the likely envelope masses of planets that do not achieve runaway growth. Most previous studies of gas accretion on to planetary cores have used either static models of planetary envelopes, where there is a thermal equilibrium for which planetesimal accretion on to the core provides the energy source that balances radiative losses from the envelope (Perri & Cameron, 1974; Mizuno et al., 1978; Rafikov, 2006), or have used quasi-static models where accretion and Kelvin-Helmholtz contraction of a gaseous envelope proceeds through a sequence of phases for which the model is assumed to be in hydrostatic equilibrium (Bodenheimer & Pollack, 1986; Pollack et al., 1996). In the former case, for a given planetesimal accretion luminosity, there is a core mass above which an envelope in hydrostatic and thermal equilibrium cannot exist, and it is normally assumed that once this critical core mass has been exceeded then rapid gas accretion will ensue. This provides one definition of the critical core mass that is often cited in the literature. However in the latter case, the accretion and settling of gas provides a source of gravitational energy that replaces the energy that is radiated away by the planet when there is no planetesimal accretion. Continuing radiative losses then enable further settling and gas accretion to occur. The quasi-static models show that giant planet formation goes through a number of phases that are characterised by different time-scales. The core is normally built by rapid planetesimal accretion in phase 1. In phase 2, the planetesimal accretion rate drops and the atmosphere grows slowly, where the rate of gas accretion is controlled by the rate of the envelope’s Kelvin-Helmholtz contraction. As the envelope mass approaches that contained in the core (the so-called cross-over mass), phase 3 begins and gas accretion increases quasi-exponentially because of the increasing importance of the envelope self-gravity, leading to rapid growth of the planet.

In this paper we adopt the approach developed in Papaloizou & Nelson (2005) to examine the accretion of gaseous envelopes onto planet cores using a sequence of quasi-static models. When discussing the critical core mass in this work, we refer to the fact that the sequence of hydrostatic models that we compute eventually comes to an end when a hydrostatic model cannot be obtained for a given core and envelope mass that has a luminosity that is sourced by the contraction and continuing accretion of the envelope. At this stage the core is said to be of a critical mass. We assume that at this point the planet enters runaway gas accretion and will eventually detach from the disc as the envelope contracts to become smaller than the planet Hill sphere. We also consider the long term evolution of the models that do not undergo runaway gas accretion after the gas disc has dispersed by evolving them for up to Gyr time-scales in the absence of gas accretion and under the assumption that the planet now has a free surface through which it radiates rather than being embedded in the protoplanetary disc.

Our motivation for this work is manifold. First, we have incorporated the gas accretion models of Papaloizou & Nelson (2005) into the combined N-body and protoplanetary disc evolution code presented in Coleman & Nelson (2014, 2016a), in order to increase the realism of the planetary system formation simulations that can be computed with this numerical tool. The purpose of this paper is to provide a simplified examination of gas accretion onto isolated, non-interacting planetary cores at a variety of orbital locations using the evolving disc models from Coleman & Nelson (2016a), as a precursor to presenting simulations that consider the contemporaneous accretion of solids and gas onto protoplanets under the influence of migration within the evolving disc models. A second motivation is to examine how accreting planets with different initial core masses at different locations in the disc evolve within disc lifetimes. Observations indicate the presence of exoplanets (super-Earths, Neptunes and Jupiters) that are likely to be gas bearing with a broad range of orbital radii, covering the interval $0.04\,{\rm au}\lesssim a_{\rm p}\lesssim 100\,{\rm au}$ (e.g. Borucki et al., 2011; Marois et al., 2008). Analysis of the occurrence rates as a function of planet radius in the Kepler data suggests that planets with radii in the range $0.85\leq R_{\rm p}\leq 4\,{\rm R}_{\oplus}$ are present around $\sim 50\%$ of Sun-like stars (e.g. Fressin et al., 2013), which is consonant with the results from earlier radial velocity surveys that focussed on low mass planets (Mayor et al., 2011). Transit timing variations (Wu & Lithwick, 2013; Hadden & Lithwick, 2016) and radial velocity measurements (Marcy et al., 2014; Weiss et al., 2016) indicate that the super-Earths and Neptunes in the Kepler data with radii $R_{\rm p}<4\,{\rm R}_{\oplus}$ have masses in the range $2\lesssim M_{\rm p}\lesssim 15\,{\rm M}_{\oplus}$, and show a tendency for the bulk density to increase for $R_{\rm p}<1.5\,{\rm R}_{\oplus}$ and decrease for $1.5<R_{\rm p}<4\,{\rm R}_{\oplus}$, indicating that the smaller planets are composed almost entirely of solids whereas the larger planets have larger radii due to the presence of significant gaseous envelopes that have been accreted during formation. Comparing theoretical models with the data suggests that the gaseous envelopes contribute little to the total masses of the super-Earths and Neptunes (Lopez & Fortney, 2014), with recent estimates indicating that planets with $R_{\rm p}\leq 1.2\,{\rm R}_{\oplus}$ contain $\sim 1\%$ gas by mass, rising to $\sim 5\%$ for planets with $R_{\rm p}\sim 3\,{\rm R}_{\oplus}$ (Wolfgang & Lopez, 2015). Considering the gas giant exoplanets, numerous attempts have been made to constrain the core masses and total heavy element abundances of those transiting planets that have mass estimates from radial velocity measurements using evolutionary models of the internal structure (Guillot et al., 2006; Burrows et al., 2007; Miller & Fortney, 2011; Thorngren et al., 2016). These studies provide a strong indication that the majority of gas giant exoplanets are substantially enriched with heavy elements relative to their host stars, with heavy element masses typically being $\sim 20\,{\rm M}_{\oplus}$ or larger. These results are consistent with the idea that giant planets form when a sufficiently massive (i.e. in the range 10 – $15\,{\rm M}_{\oplus}$) core forms that is able to accrete gas and undergo runaway gas accretion sufficiently early in the disc lifetime, in basic agreement with the theoretical core accretion models. In the absence of a sufficiently massive core that forms early enough, a more common outcome is the formation of super-Earths or Neptunes that fail to become gas giants because of the relatively slow rate of gas accretion until the envelope mass becomes comparable to the core mass. It is clearly of interest to address the above observational results through theoretical modelling of in situ planet formation, and to compare the outcomes of our calculations with recent work that has considered gas accretion onto planets with short (Lee et al., 2014; Batygin et al., 2016) and long orbital periods (Piso & Youdin, 2014; Piso et al., 2015). Our third motivation is to consider the longer term cooling and contraction of the gaseous envelopes for those planets that remain in a state of quasi-static equilibrium throughout the period of gas accretion when the disc is present. Since the observations are typically relevant to exoplanets that have existed around their stars for Gyrs, they will have evolved for time periods of this duration after the epoch of formation. This is the first time that our models have been used to examine this longer term evolution, and our approach is similar to that used in recent population synthesis studies (see the recent review by Mordasini et al., 2015).

This paper is organised as follows. We present the basic equations and their methods of solution in Sections 2 and 3. In Section 4 we present the evolution of the disc models, and in Section 5 we present calculations of the accretion of gas onto cores of varying mass located at different disc locations. In Section 6 we present calculations of the long term evolution of planetary envelopes, and in Section 7 we discuss our results and draw conclusions.

2 Calculating gas accretion rates

The temporal evolution of the accreting gaseous envelopes in our models are computed using the methods described in Papaloizou & Nelson (2005). This adopts the approach presented in previous works (e.g. Papaloizou & Terquem, 1999; Pollack et al., 1996; Bodenheimer & Pollack, 1986) in approximating the gaseous envelope structure as being non-rotating, spherically symmetric and in hydrostatic equilibrium, similar to the approach used in modelling stellar evolution. The equation of hydrostatic balance is given by

[TABLE]

where $P$ is the pressure, $\rho$ is the density, and $g$ is the acceleration due to gravity, $g=GM(R)/R^{2}$ with $M(R)$ being the mass, including that of the solid core, interior to the radius $R$, and $G$ is is the gravitational constant.

We adopt the equation of state for a hydrogen and helium mixture given by Saumon et al. (1995). The mass fractions of hydrogen and helium are 0.7 and 0.28 respectively. The luminosity, $L_{\rm rad}$, transported by radiation is given by

[TABLE]

where $\sigma$ is the Stefan-Boltzmann constant, $T$ is the temperature, and $\kappa$ is the opacity. We assume that the opacity is given by the Rosseland mean opacity, with the temperature and density dependencies being calculated using the formulae in Bell et al. (1997) for temperatures below 3730 K, and by Bell & Lin (1994) above 3730 K.

2.1 Inner convection zones

When the planet is radiating on account of the envelope settling onto the core, most of the mass in the envelope is unstable to convection so a theory of energy transport by convection is required. We adopt the conventional mixing length theory (e.g. Cox & Giuli, 1968).

The radiative and adiabatic temperature gradients $\nabla_{\rm rad}$ and $\nabla_{\rm ad}$ are defined as

[TABLE]

and

[TABLE]

with the subscript S denoting evaluation at constant entropy. The total luminosity is $L_{\rm tot}$. Within our models this is equal to the rate of release of gravitational energy of gas accreted from the surrounding protoplanetary disc, as well as that liberated from the gravitational settling of the gaseous envelope (see sect. 2.3). We include the binding energy of the planet core in the initial calculation of the luminosity, and in future work we will include the accretion luminosity obtained from the release of gravitational energy of accreted planetesimals. In this paper, however, the accretion of solids in the form of planetesimals, boulders or pebbles is neglected.

When the gas is convectively stable, such that $\nabla_{\rm rad}<\nabla_{\rm ad}$, the energy is transported entirely through radiation. However when $\nabla_{\rm rad}>\nabla_{\rm ad}$, the gas is convectively unstable and so some of the energy is transported through convection. The total luminosity passing through a sphere of radius $r$ is equal to $L_{\rm tot}(r)=L_{\rm rad}+L_{\rm conv}$, where $L_{\rm conv}$ is the luminosity associated with convection. Adopting the mixing length theory (Cox & Giuli, 1968), this is equal to

[TABLE]

where $\Lambda=\left|\alpha P/(dP/dR)\right|$ is the mixing length, $\alpha$ being a constant parameter expected to be of order unity, $(dT/dR)_{\rm S}=\nabla_{\rm ad}T(d{\rm ln}P/dR)$, and the subscript P denotes evaluation at constant pressure. We adopt the mixing length parameter $\alpha=1$. As for our models, convection is efficient near the core where the gravitational energy is liberated, the temperature gradient will be close to the adiabatic gradient with the luminosity dropping out of the equations. For this reason we take $L_{\rm tot}$ to be constant which is the approximation usually made when considering stellar envelopes in corresponding circumstances.

2.2 Boundary conditions

2.2.1 Inner boundary

We assume that there is a solid core of mass $M_{\rm core}$ with a uniform density, $\rho_{\rm core}=3.2$ gcm*-3* (Bodenheimer & Pollack, 1986). The models that calculate the structure of the gaseous envelope must therefore satisfy that the total mass $M(R_{\rm core})=M_{\rm core}$ when

[TABLE]

2.2.2 Outer boundary

We assume that the envelope structure extends to the Roche lobe, or the boundary of the Hill sphere beyond which material must be gravitationally unbound from the planet. The surface radius of the core-envelope structure, $R_{\rm surf}$, is then given by

[TABLE]

where $M_{\rm p}$ is the total planet mass (gaseous envelope and solid core), and $a_{\rm p}$ is the orbital radius of the planet. We assume that the outer layers of the gaseous envelope eventually join smoothly to the surrounding protoplanetary disc, where $T=T_{\rm disc}$, $P=P_{\rm disc}$ and $\rho=\rho_{\rm disc}$ respectively. Therefore the boundary conditions are that at $R=R_{\rm surf}$, $M(R_{\rm surf})=M_{\rm p}$, $P_{\rm surf}=P_{\rm disc}$ and the temperature is given by

[TABLE]

where $T^{4}_{\rm effp}=3\tau_{\rm surf}L_{\rm tot}/(16\pi\sigma R^{2}_{\rm surf})$ with $L_{\rm tot}$ denoting the total luminosity escaping from the surface. We approximate the additional optical depth above the planet’s atmosphere, through which radiation passes, by (Papaloizou & Terquem, 1999)

[TABLE]

Equation 8 expresses the fact that $T_{\rm surf}$ must exceed the disc temperature $T_{\rm disc}$ at $R=R_{\rm surf}$ in order that the luminosity be radiated away from the planet and into the surrounding protoplanetary disc.

2.3 Accretion, settling and envelope evolution sequences

The models we present here provide an evolutionary sequence where the mass $M_{\rm p}$ increases through the accretion of gas from the surrounding protoplanetary disc. As the gas settles, it also liberates gravitational energy, increasing the planet’s luminosity. To calculate this energy attributed to the planet’s evolution, we consider the total energy of the planet within the Roche lobe

[TABLE]

Here $U$ is the internal energy per unit mass and we neglect the energy involved in bringing material from infinity to the Roche lobe. We justify this as the majority of the mass is concentrated well inside the Roche lobe, close to the planet’s core where the specific energies are significantly higher.

Since the planet radiates energy into the surrounding disc, this energy must come from either the settling of the atmosphere described above, or through the release of gravitational energy from gas that is accreted by the planet. This change in the energy of the planet as its mass increases is $dE=(dE/dM_{\rm p})dM_{\rm p}$. As this energy change balances the loss of energy through radiation, then the conservation of energy requires that

[TABLE]

with $L$ being the luminosity at the surface. This then determines the evolution of the atmosphere as it settles and continues to accrete gas from the surrounding disc. We note that the accretion of planetesimals, neglected in this work, can also provide an additional energy source. If included in the models, then this would reduce the gas accretion rate as it would supply some fraction of the energy required to balance that radiated by the planet. The final masses of the gas envelopes would then be reduced accordingly.

To summarise the procedure used to generate a sequence of models and accretion rates, we remark that by analogy with stellar structure, for given core and envelope masses we have three state variables, $(P,T,R)$ to be determined as a function of $M(R)$ through equations 1–5. We have boundary conditions for each of these at the surface and one for $R,$ at the surface of the core. The problem thus appears to have one too many boundary conditions and thus be overdetermined. This is resolved by allowing $L_{\rm tot}$ to be determined so that all the conditions can be satisfied. Thus specification of the total mass together with the local disc properties specifies $L_{\rm tot}$. This is then used to determine the accretion rate using equation 11 which in turn determines the rate at which the envelope mass increases.

3 Calculation of the evolution of envelopes with free surfaces

Once the protoplanetary disc has dispersed, the surface boundary conditions have to be changed from those used above so as to represent a free surface. The energy radiated into space is as before supplied by gravitational settling but at fixed envelope mass which requires the radius to contract. To obtain models belonging to the contracting sequence, we proceed as above but implement the different boundary conditions. Firstly the radius $R_{\rm surf}$ is not fixed to be the radius of the Hill sphere but is determined in the course of the solution. As for stellar models we take the radius to define the photosphere where the optical depth is taken to be $2/3.$ At this location the pressure is given by

[TABLE]

The temperature is given by

[TABLE]

where $T^{4}_{\rm effp}=L_{\rm tot}/(4\pi\sigma R^{2}_{\rm surf})$ with $L_{\rm tot}$ denoting the total luminosity escaping from the surface, and $T^{4}_{*}=L_{*}/(16\pi\sigma R^{2}_{p})$ with $L_{*}$ denoting the luminosity escaping from the central star. Here $T_{\rm effp}$ is the effective temperature associated with the luminosity of the planet and $T_{*}$ corresponds to the equilibrium temperature of a black sphere in the external radiation field of the star. Thus we assume that reprocessing of this radiation leads to a heating of the atmosphere that can be approximated as being spherically symmetric. The interior boundary conditions at the surface of the core are as above.

The calculation of the models proceeds as above except that instead of the radius being fixed the envelope mass is. Thus instead of obtaining a sequence for which the envelope mass is a function of the luminosity, we obtain one for which the radius is a function of luminosity. In this way we get a sequence of models with fixed $M_{\rm p}$ and varying $R_{\rm surf}$, each model being associated with a determined luminosity.

An evolutionary sequence of contracting models can be determined by equating the rate of increase of binding energy with the rate of radiation of luminous energy. To do this the total energy of the envelope for a given model is found using equation 10. Its derivative with respect to $R_{\rm surf}$ is found numerically from considering two models with slightly different $R_{\rm surf}$. Then the radius can be evolved forwards in time having obtained $dR_{\rm surf}/dt$ from

[TABLE]

3.1 End states

The sequence of models described above have decreasing luminosity as they contract. The envelope thus becomes denser and cooler. As the only energy source considered is that due to gravitational settling of the envelope, the end state will be an isothermal envelope heated by stellar radiation with zero internal luminosity. Attainment of this state could of course be retarded if a cooling core provides a heat source. Accordingly, the lifetimes of inflated states might be underestimated by our calculations as we neglect this effect.

As the envelopes become cool and dense we found that in some cases the physical parameters could become out of the range of the applicability of the equation of state of Saumon et al. (1995). Accordingly our end state models were calculated using the equation of state of Becker et al. (2014). Tests showed that very similar results were obtained for models using the different equations of state when the parameters were in the range of validity of both.

4 Background disc model

We utilise the disc model presented in Coleman & Nelson (2014, 2016a) where the standard diffusion equation for a 1D viscous $\alpha$-disc model is solved (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). Temperatures are calculated by balancing black-body cooling against viscous heating and stellar irradiation. The viscous stress parameter $\alpha=2\times 10^{-3}$ throughout most of the disc, but increases to $\alpha=0.01$ in regions where $T\geq 1000$ K to mimic the fact that fully developed magnetohydrodynamic turbulence can develop in regions where the temperature exceeds this value (Desch & Turner, 2015).

The final stages of disc removal occur through a photoevaporative wind. A standard photoevaporation model is used for most of the disc evolution (see Dullemond et al., 2007, for a review), corresponding to a photoevaporative wind being launched from the upper and lower disc surfaces. Direct photoevaporation of the disc is switched on during the final evolution phases when an inner cavity forms in the disc (due to the inner parts viscously accreting onto the star), corresponding to the outer edge of the disc cavity being exposed to the stellar radiation (Alexander & Armitage, 2009). We note that when a giant planet begins to form in a disc, and gap opening criteria are met, then tidal torques from the planet are able to open a gap using the impulse approximation from Lin & Papaloizou (1986).

Figure 1 shows the evolution of a 1 $\times{\rm MMSN}$ (Minimum Mass Solar Nebula, Hayashi, 1981) disc model. Disc surface density profiles are shown in the left-hand panel, temperature profiles are shown in the middle panel, and $H/r$ profiles are shown in the right-hand panel. The times corresponding to each profile are shown in the middle panel, expressed as a percentage of the disc lifetime ($\sim$ 4.6 Myr for a 1 $\times{\rm MMSN}$ disc). The active turbulent region where the local disc temperature exceeds 1000 K can be seen in the inner region of the disc, gradually moving in towards the star as the surface density decreases, reducing the viscous heating rate and the opacity. This region disappears when the temperature no longer exceeds 1000 K anywhere in the disc, as occurs at the evolutionary stage illustrated by the yellow curves in Fig. 1. This happens in all of the disc models when the gas surface density in the inner disc drops below $\sim$ 100 $\rm g/cm^{2}$, which typically occurs 0.5 Myr before complete dispersal of the gas disc.

On longer time-scales the removal of gas by the photoevaporative wind causes the disc to disperse. The loss of mass at large radius results in the inner disc emptying viscously onto the star leaving a gap that separates the remnant inner disc from the outer disc which is clearly seen in the surface density profile at the evolutionary stage illustrated with yellow curves in Fig. 1. Finally the remnant outer disc is removed by the wind (Clarke et al., 2001).

When the disc mass is increased to larger values, as will be considered in this work, the evolution is qualitatively similar to that of the 1 $\times{\rm MMSN}$ disc. The major difference is in the lifetime of the disc. We find that 1 $\times{\rm MMSN}$ discs have lifetimes $\sim 4.6$ Myr, whereas 3 and 5 $\times{\rm MMSN}$ discs last for $\sim$ 7.8 and 9.4 Myr respectively. Apart from existing for longer, the main effect that increasing the disc mass will have on gas accreting planetary cores is to change the physical conditions at their outer boundaries where their gaseous envelopes merge with the gas disc. For more massive discs there will be an increase in the viscous heating rate, leading to slightly higher disc temperatures, and the local gas density and pressure will also be increased. These increases will affect the structure of the gaseous envelope and modify the gas accretion rate. In practice, however, we find that cores of a given mass located at a particular orbital radius evolve very similarly in the different disc models until they reach a critical state, because the local conditions in the models do not vary by a large amount. Once the planets have reached a critical state, they undergo runaway gas accretion and then begin to accrete gas from the disc at the viscous supply rate. The evolution of these planets from this point then depends on the remaining disc mass and lifetime, which will lead to different evolution paths for similar planets embedded in discs with different masses.

Table 1 shows the disc parameters that are used in all of the simulations. The value of $f_{41}$ defines the value of the photoevaporative UV flux that is incident on the disc, as described by Dullemond et al. (2007).

5 In situ accretion of gaseous envelopes

We now present one fiducial model in detail, before discussing the results for cores of different masses that are placed at different locations in the three disc models that we consider.

5.1 Envelope evolution

We solve equations 1–5 using the boundary conditions described in Section 2.2 to obtain the envelope structure models. In this work we only consider single planets that are embedded in protoplanetary discs. We use core masses of 5, 10 and 15 $\,{\rm M}_{\oplus}$, and we place them in discs of mass 1, 3 and 5 $\times{\rm MMSN}$. The planets are placed at different locations in the disc: 0.1, 0.2, 0.5, 1, 2, 5 and 10 $\,{\rm au}$. We embed the planets at different locations in the disc to examine the effects that the different disc locations and conditions have on the accretion of gaseous envelopes. The combination of these parameters results in a total of 63 different simulations. To allow the disc time to settle into a quasi-steady state, we inject the cores into the disc after a time of 0.1 Myr. From this point onwards the cores are able to accrete gas from the disc following equation 11, so that the rate of liberation of accretion energy balances the planet’s luminosity. All gas that is accreted by the planet is removed from the surrounding disc.

When calculating a sequence of hydrostatic models for which the mass of the gaseous envelope increases, one may reach a point where a hydrostatic model can no longer be obtained. Here, the luminosity associated with the liberation of gravitational energy by the settling envelope is insufficient to allow the envelope to be supported, and it is expected that the planet will then enter a period of runaway gas accretion. In this paper we refer to this moment as being when the envelope reaches a critical state. If the planet’s envelope reaches a critical state, we change the gas accretion routine to that used in Coleman & Nelson (2016b), which utilises fits based on models from Movshovitz et al. (2010) to calculate the gas accretion rate, including during the early stages of runaway accretion. In our models, once in this phase of runaway gas accretion, the planet rapidly accretes material from within its feeding zone. Once the material in the feeding zone has been accreted, the planet opens up a gap in the protoplanetary disc as indicated by the torque formulae of Lin & Papaloizou (1986). At this point the planet starts to accrete gas at the viscous supply rate. If the protoplanetary disc disperses before the planet’s envelope reaches a critical state, however, we transition the envelope to a free boundary model and allow it to cool and settle over time as discussed in Section 6.

We will now present the results of the simulations. To begin, we will examine a single case where we show the envelope structure and present its evolution over time. We will then examine the effects that placing the planetary cores in different regions of the discs has on the accretion and evolution of the envelopes, before considering the effects that different initial core masses have on their evolution.

5.2 A fiducial envelope accretion calculation

We begin by describing the time evolution of a gaseous envelope accreting onto a 10 $\,{\rm M}_{\oplus}$ core located at 10 $\,{\rm au}$ in a 1 $\times{\rm MMSN}$ disc. We note that the radius of the core is $1.7\times 10^{9}$ cm. Figure 2 shows the internal structure of this envelope at four moments in time: blue line – 0.1 Myr; red line – 1 Myr; black line – 2.5 Myr; green line – 4.2 Myr. The temporal evolution of the mass of this envelope is shown by the blue line that starts at $10\,{\rm M}_{\oplus}$ in the lower right panel of Fig. 3. The evolution of the luminosity, energy liberated per unit mass, and the mass accretion rate are shown by the black lines in Fig. 4.

We initialise the calculations with an envelope mass that is equal to 1 % of the core mass, i.e. a 0.1 $\,{\rm M}_{\oplus}$ envelope is assumed to be present initially surrounding the 10 $\,{\rm M}_{\oplus}$ core. The model then calculates the planet luminosity that is required to support this envelope according to the assumed equation of state, opacity and boundary conditions. The blue lines in Fig. 2 show the structure of this initial envelope. As expected, the innermost part of the atmosphere is convective. This convection zone extends from the core surface to $\sim 10^{10}$ cm, and contains $\sim 90\%$ of the envelope’s mass. There is also a second convection zone that extends between 7 $\times 10^{10}$ cm to $10^{12}$ cm. These convection zones occur for temperatures $T>1860$ K and $55$ K $\leq T\leq 570$ K, respectively.

As time progresses, the gaseous envelope settles onto the planet’s core as more material is accreted from the surrounding disc. One effect of this is to cause the luminosity of the planet to decrease as the optical depth through the atmosphere increases, as shown in the left panel of Fig. 4. After 1 Myr, shown by the red lines in Fig. 2, the planet has a mass of 13 $\,{\rm M}_{\oplus}$, with the gaseous envelope containing $3\,{\rm M}_{\oplus}$. There are two convection zones, one in the deep interior of the planet extending out to $1.2\times 10^{10}$ cm, and one further out between the radii $3\times 10^{10}$ cm and $5\times 10^{11}$ cm. These two convection zones differ from those in the initial envelope because of the temperatures at which they occur. Compression has heated the inner envelope, and the temperatures in the innermost convective zone now exceed T $>2340$ K. The outer convection zone occurs for $55$ K $\leq T\leq 1280$ K. At this time, approximately 66 % of the envelope’s mass lies within the inner convection zone, with a further 21 % being situated within the outer convective zone.

The black lines in Fig. 2 show the envelope structure after 2.5 Myr. Up until this point, the luminosity of the planet has been decreasing due to the increasing opaqueness of the envelope arising from its early compression (see Fig. 4). After 2.5 Myr, however, the luminosity and gas accretion rates start to increase as the self-gravity of the envelope helps to compress it. Here the gaseous envelope has a mass of 6.2 $\,{\rm M}_{\oplus}$, and hence comprises $\sim 40\%$ of the total planet mass. The inner convective zone again contains $\sim 66\%$ of the envelopes mass, and extends out to 1.7 $\times 10^{10}$ cm. The outer convective zone extends between $3.4\times 10^{10}$ cm and $5.1\times 10^{11}$ cm and contains $\sim 25\%$ of the envelope’s mass. The temperatures in the outer and inner convection zones are $55$ K $\leq T\leq 1400$ K and T $>2160$ K, respectively.

After a further 1.7 Myr, corresponding to a total of 4.2 Myr of evolution, the planet reaches a critical state and is about to enter runaway gas accretion. The green lines in Fig. 2 show the envelope structure at this point. The planet has a total mass of 30.6 $\,{\rm M}_{\oplus}$ with the envelope contributing $67\%$ of this. The inner convective zone now extends out to $\sim 2.8\times 10^{10}$ cm and contains $\sim 73\%$ of the envelope’s mass, whilst the outer convective zone lies between $6\times 10^{10}$ cm and $7.7\times 10^{11}$ cm and contains $\sim 20\%$ of the envelope’s mass. These convective zones occur at similar temperatures to those found for the convective zones of the atmosphere after 2.5 Myr.

We recall that the disc model in which the planet is embedded disperses after 4.6 Myr. As described in Sect. 5.1, when the planet envelope reaches the critical state then the assumption of hydrostatic equilibrium breaks down and another prescription for mass evolution of the planet is required. Our adoption of fits to the mass versus time plots presented in Movshovitz et al. (2010), followed by assuming that the planet grows at the viscous supply rate of gas by the protoplanetary disc when the planet becomes gap forming, do not allow us to calculate the detailed structure of the planet, so the green lines in Fig. 2 represent the final time at which the detailed envelope structure can be computed. The longer term mass growth of the planet, however, is shown in the bottom right panel of Fig. 3, which is discussed in more detail in the next section.

5.3 Effect of location

As shown in Fig. 1, the disc models that we employ have strongly varying temperatures as a function of radius, and this variation in temperature has a significant influence on the accretion histories of the planets that are placed at different locations in the discs, a point that has been explored in previous work (e.g. Piso & Youdin, 2014). For a given time interval over which accretion occurs, gaseous envelopes that accumulate onto planets close to the central star are expected to be less massive than those accreted by planets at large orbital radii. The gas orbiting close to the star is much hotter than that further away, reducing gas accretion rates because of its higher thermal energy content. The colder gas located at large orbital radii is accreted more efficiently because its lower thermal energy allows it to become more deeply bound to the accreting cores. From the point of view of constructing a sequence of hydrostatic models, as undertaken here, increasing the disc temperature increases the surface temperatures of the models through equation 8, and this in turn increases the temperatures and entropies throughout the envelopes, leading to longer cooling times and slower accretion rates.

Figure 3 shows the mass evolution of 5, 10 and 15 $\,{\rm M}_{\oplus}$ cores situated at 0.1 $\,{\rm au}$ (top left-hand panel), 0.5 $\,{\rm au}$ (top right-hand panel), 5 $\,{\rm au}$ (bottom left-hand panel) and 10 $\,{\rm au}$ (bottom right-hand panel) in their respective discs. The blue lines show the evolution in 1 $\times{\rm MMSN}$ discs, with the red and black lines denoting 3 and 5 $\times{\rm MMSN}$ discs respectively, and the black crosses indicate that a planet has reached a critical state and is about to undergo runaway gas accretion.

5.3.1 Accretion at $a_{\rm p}\leq 0.5\,{\rm au}$

The top left-hand panel of Fig. 3 shows the evolution of the masses of the planets orbiting at 0.1 $\,{\rm au}$. The evolution of the luminosity, change in energy liberated per unit mass, and accretion rate for the planet with a $10\,{\rm M}_{\oplus}$ core located in a $1\times{\rm MMSN}$ disc located at $0.1\,{\rm au}$ is shown by the blue line in Fig. 4. These planets are in a hot and dense environment, with temperatures moderately exceeding 1000 K for the majority of the disc lifetime (we note that dust grains remain an important source of opacity in the models at these temperatures).

It is clear from Fig. 3 that none of the planets located at $0.1\,{\rm au}$ accrete large amounts of gas and reach the critical state that precedes runaway gas accretion. As such, our models predict that none of these planets can grow to become gas giants. The most massive envelope that could be accreted was in a 5$\times{\rm MMSN}$ disc, where a 15 $\,{\rm M}_{\oplus}$ core accreted a 7.2 $\,{\rm M}_{\oplus}$ envelope, representing $\sim 33\%$ of the planet’s total mass. The core was able to accrete this significant envelope due to the long lifetime of the disc ($\sim 9.4$ Myr), and the high core mass which enhances the gas accretion rate as discussed in Section 5.3.3. The $10\,{\rm M}_{\oplus}$ cores located at $0.1\,{\rm au}$ were only able to accrete between $1.2\,{\rm M}_{\oplus}$ and $1.6\,{\rm M}_{\oplus}$, with more gas accreting onto cores in the more massive and longer lived disc models. The $5\,{\rm M}_{\oplus}$ cores accreted very small amounts of gas, with the maximum amount being $0.2\,{\rm M}_{\oplus}$ which accreted in the long-lived $5\times{\rm MMSN}$ model.

For a given core mass, we find that the outcomes of the calculations for planets located at $0.1\,{\rm au}$ are broadly similar when varying the disc model, and this is to a significant degree because in all models the disc temperature at this location is $\sim 1000$ K for a large part of the evolution, due to our inclusion of a viscosity switch at $T\geq 1000$ K to mimic the development of fully developed magnetised turbulence in the hot inner disc. As the disc model evolves, the location of the highly viscous region moves inwards, and it eventually moves interior to $0.1\,{\rm au}$ at which point the planet finds itself in a cooler environment. This normally occurs approximately 2 Myr from the end of the disc lifetime in all models, and causes there to be a modest increase in the gas accretion rate as shown by the blue lines in Fig. 4 for the planet with a $10\,{\rm M}_{\oplus}$ core.

Moving the cores further out into cooler disc regions results in enhanced accretion rates, and we can see from the top right panel of Fig. 3 that placing cores at $0.5\,{\rm au}$ leads to the accumulation of more massive envelopes. In particular, the models with $15\,{\rm M}_{\oplus}$ cores all reach the critical state and undergo runaway gas accretion to become gas giants. For a given core mass, we can see that the planet evolution models are rather insensitive to the disc models. The planets with $15\,{\rm M}_{\oplus}$ cores all reach the critical state after $\sim 2.6$ Myr when their total masses are $\sim 35$ - $40\,{\rm M}_{\oplus}$, with the gaseous envelopes contributing approximately $60\%$ of their mass. After reaching the critical state these planets undergo runaway gas accretion and become gas giants by rapidly accreting the gas from their feeding zones before accreting gas from the disc at the viscous rate. These planets are finally able to reach masses $0.8M_{\rm J}\leq m_{\rm p}\leq 2M_{\rm J}$.

Considering the planets with 10 $\,{\rm M}_{\oplus}$ cores, it was found that they are unable to accrete enough gas to reach the critical state. In the $1\times{\rm MMSN}$, with a lifetime of 4.6 Myr, the 10 $\,{\rm M}_{\oplus}$ core was only able to accrete $\sim 2.8\,{\rm M}_{\oplus}$ of gas. The planet in the $5\times{\rm MMSN}$ disc accreted a $\sim 4.3\,{\rm M}_{\oplus}$ gas envelope, comprising $30\%$ of the planet’s total mass. This increase in gas accretion is due to the disc having a longer lifetime of 9.4 Myr. The 5 $\,{\rm M}_{\oplus}$ cores located at $0.5\,{\rm au}$ were still unable to accrete significant gaseous envelopes on account of their low masses. Envelope masses for these cores ranged between 0.3 $\,{\rm M}_{\oplus}$ and 0.4 $\,{\rm M}_{\oplus}$, with the latter envelope being accreted in the more massive, longer lived disc.

To understand the evolution of the planets as a function of orbital radius, it is useful to consider Fig. 4. We recall from equation 11 that the mass accretion rate of the planet is controlled by the luminosity of the hydrostatic model and the quantity $dE/dM_{\rm p}$, which measures the rate at which the total energy of the planet changes as a function of its mass. During the early evolution, the left panel in Fig. 4 shows that the total luminosity drops rapidly as the envelope contracts and the optical depth through the envelope increases, but thereafter the luminosity is approximately constant throughout most of the evolution. We also see that for models that do not undergo runaway gas accretion, $dE/dM_{\rm p}$ is also essentially constant during the evolution. So gas accretion occurs at an almost constant rate throughout most of the evolution. We note, however, that the luminosity required for the models to achieve hydrostatic balance is small at small orbital radii and increases for planets located further from the star, and $dE/dM_{\rm p}$ decreases in magnitude as one moves the models further from the star. This results in the gas accretion rate increasing with increasing stellar distance for a fixed core mass, as observed in the models.

5.3.2 Accretion at $a_{\rm p}\geq 5\,{\rm au}$

The accreting planets located at 0.1 and 0.5 $\,{\rm au}$ discussed above demonstrate how gaseous envelopes settle on to planetary cores of various masses in hot and warm environments close to the central star. Further out in the disc, however, the temperature is significantly lower, and this allows for more rapid accretion of massive gaseous envelopes, resulting in the occurrence of runaway gas accretion after the envelopes have reached the critical state.

The bottom left panel of Fig. 3 shows that at 5 $\,{\rm au}$, planets with 15 $\,{\rm M}_{\oplus}$ cores reach critical states after between 1.2–1.35 Myr, and rapidly accrete the gas from their feeding zones before accreting from the disc at the viscous supply rate, reaching masses in the range $1.5M_{\rm J}\lesssim M_{\rm p}\lesssim 3M_{\rm J}$ at the end of the disc lifetime. Planets with 10 $\,{\rm M}_{\oplus}$ also undergo significantly more rapid gas accretion, such that those placed in the 3 and 5 $\times{\rm MMSN}$ discs undergo runaway gas accretion. A common feature of all of the calculations is the similarity of the quasi-static phase of evolution as a function of the disc model. The disc lifetime is primarily responsible for the variation observed between the disc models. For example, all of the 10 $\,{\rm M}_{\oplus}$ cores placed in 1, 3 and 5 $\times{\rm MMSN}$ discs follow almost nearly identical tracks until the 1 $\times{\rm MMSN}$ disc reaches the end of its lifetime, leaving this 10 $\,{\rm M}_{\oplus}$ core with a sub-critical gaseous envelope containing $\sim 8\,{\rm M}_{\oplus}$. The cores placed in 3 and 5 $\times{\rm MMSN}$ discs, however, reached a critical state after 6.3 and 6.6 Myr respectively, when their masses were 26.1 and 26.6 $\,{\rm M}_{\oplus}$, and grew to become gas giants with masses between 0.5–1.5 $M_{\rm J}$. As with the 5 $\,{\rm M}_{\oplus}$ cores located at 0.1 and 0.5 $\,{\rm au}$, cores of this mass located at 5 $\,{\rm au}$ failed to accrete significant gas envelopes, with the most massive forming in the 5 $\times{\rm MMSN}$ disc with a final envelope mass of 1.3 $\,{\rm M}_{\oplus}$.

Finally, the bottom right panel of Fig. 3 shows the evolution of cores placed at 10 $\,{\rm au}$ in the disc. The planets with 15 $\,{\rm M}_{\oplus}$ cores all reach the critical state after 0.93–0.95 Myr with total masses $\sim 35\,{\rm M}_{\oplus}$. Due to being located further out in the disc, where the planet’s feeding zones are larger, and undergoing runaway gas accretion at earlier times, these giant planets reach masses of between 2–4 $M_{\rm J}$ at the end of the disc lifetime. The planets with 10 $\,{\rm M}_{\oplus}$ cores also reached the critical state once their masses were between 27 and 31 $\,{\rm M}_{\oplus}$. These planets all underwent runaway gas accretion. The planet in the $1\times{\rm MMSN}$ disc entered the runaway state near the end of the disc lifetime, when there was a limited supply of gas available, so this planet did not manage to reach a Jovian mass but instead achieved a final mass of $44\,{\rm M}_{\oplus}$. The planets with 10 $\,{\rm M}_{\oplus}$ cores in the more massive discs finished the simulations with masses between 2–3 $M_{\rm J}$.

As with the other locations in the disc, the 5 $\,{\rm M}_{\oplus}$ cores again failed to accrete enough gas to undergo runaway gas accretion.

To summarise, cores that accrete gaseous envelopes close to the central star are unable reach a critical state due to the high gas temperatures there, but accretion in the cooler disc regions further from the star is much more efficient. The curves shown in the various panels of Fig. 4 explain this behaviour. Over most of the evolution, the luminosity, energy change per unit mass, and the mass accretion rates are more or less constant. Planet models located close to the star in hot regions of the disc have lower luminosities, larger values (in magnitude) of ($dE/dM_{\rm p}$), and hence smaller accretion rates than planet models located further from the star but with the same core mass. Significant variation of these quantities as a function of time is only really observed at the very beginning of the calculations, and during the time when a planet model approaches the critical state. During this latter phase, Fig. 4 indicates that the luminosity increases as a result of the increasing and more strongly bound envelope mass having a higher temperature at the surface of the core which leads to an increased rate of outward heat transport, and ($dE/dM_{\rm p}$) decreases in magnitude sharply, resulting in the dramatic increase in the accretion rates observed in all calculations of this type (e.g. Pollack et al., 1996).

5.3.3 Influence of core mass

The previous sections examined how the accretion and evolution of gaseous envelopes depended on the core’s location in the disc. Figure 5 shows how the accretion of gaseous envelopes is dependant on the initial core mass, and the general trend of higher core masses resulting in faster gas accretion rates can be seen. Each panel shows the initial core mass, linked by a vertical dashed line to a final core and envelope mass. When a core is denoted by a cross, the envelope has reached a critical state, and we plot the mass at this time (not the final mass after runaway gas accretion has occurred). The blue, red and black colours denote the accretion in 1, 3 and $5\times{\rm MMSN}$ discs respectively, and we show the results of cores accreting gaseous envelopes at 0.1, 0.5, 5 and 10 $\,{\rm au}$. The left panel shows the evolution of 5 $\,{\rm M}_{\oplus}$ cores. Only at 10 $\,{\rm au}$ in the most massive disc considered here, is such a core able to accrete a sizeable envelope, comprising $\sim 36\%$ of the total mass. The middle panel shows the evolution of 10 $\,{\rm M}_{\oplus}$ cores, which are able to accrete much more massive gaseous envelopes than their 5 $\,{\rm M}_{\oplus}$ core counterparts. The right panel shows the evolution of 15 $\,{\rm M}_{\oplus}$ cores. The increase in final gas envelope mass with increasing core mass is evident here, with only the planets located close to the star at $0.1\,{\rm au}$ failing to undergo runaway gas accretion. Figure 5 shows that for planets to accrete significant gaseous envelopes and reach a critical state, they need to have masses of $10\,{\rm M}_{\oplus}$ or higher. This is in line with results obtained in some different contexts to ours, such as in prescribed nebula models with planetesimal accretion included (e.g. Pollack et al., 1996).

The trend for increasing gas accretion rate with core mass can be easily understood from inspection of Fig. 6, which shows the evolution of various quantities for planets with different core masses located at $0.5\,{\rm au}$ in the $1\times{\rm MMSN}$ disc. Planets with more massive cores require larger luminosities to support their envelopes due to the larger binding and thermal energies per unit mass, and collectively the planets have values of $dE/dM_{\rm p}$ which are similar in magnitude during the early and middle phases of gas accretion, leading to larger gas accretion rates and shorter evolution times for planets with more massive cores.

6 Long term evolution

In order to have results that are comparable with observations on appropriate time-scales, it is necessary to consider the evolution of protoplanets after disc dispersal. During this phase the envelope contracts in order to release gravitational energy to balance the energy that is radiated away and the radius shrinks while the planet surface is heated by radiation from the central star. We consider models that have not reached a critical state before disc dispersal and thus only objects that may become Super-Earths or Neptunes. For initial models we take a sample of protoplanets that were embedded in their nascent disc during its lifetime just after its dispersal. Their total mass, core mass and orbital location at that stage are listed in table 2. Models A1, A2 and A3 have $5\,{\rm M}_{\oplus}$ cores and are of increasing envelope mass. Models B1, B2 and B3 are of increasing envelope mass, each having $10\,{\rm M}_{\oplus}$ cores.

The calculations are performed as outlined in Section 3. The evolution of the models with initial parameters listed in table 2 is illustrated in Fig. 7. Radii and internally generated luminosities are shown as a function of time. It will be seen that after an initially fairly rapid contraction from the vicinity of the Hill sphere, and reduction in internal luminosity, the evolution rate slows significantly. The evolution rate decreases with both increasing initial core mass and increasing envelope mass, the radii of models A1, A2 and A3 reaching $3\times 10^{9}$ cm (i.e. $\sim 4.6\,{\rm R}_{\oplus}$) after $\sim 2,4$ and 8 Gyr, respectively. To investigate the dependence on the initial envelope mass we ran the same models with their initial envelope masses reduced by a factor of $5.$ The results are illustrated in Fig. 8. It will be seen that the evolution times to attain a radius of $3\times 10^{9}$ cm are reduced to $\sim 0.5$, 3 and 4.5 Gyr, respectively. This type of evolution indicates that the models may not be fully relaxed on Gyr time scales, and as a consequence observed systems may exhibit residual inflation.

6.1 Effect of reducing the opacity

We have found that the main quantity determining the rate of evolution is the form of the opacity. This is because this is determined by the rate at which radiation can diffuse out of the optically thick envelope. To illustrate this we have considered the evolution of models for which the opacity was reduced by a factor of $100$. The models considered were A1, A2 and B3, though in the latter case we also reduced the initial envelope mass by a factor of $10$.

The evolution of these models is illustrated in Fig. 9. Black and blue solid curves correspond to models A1 and A2 with the opacity reduction. The red dashed curve corresponds to model B3 without opacity reduction but with initial envelope mass reduced by a factor of 10. The red solid curve corresponds to the same model but with the opacity reduction. By comparing Figs. 7 and 9 we see that the evolution rate for the models with the opacity reduction is sped up by a factor of $100$. This indicates that the evolution rate $\propto 1/\kappa.$

This may be understood as a consequence of most of the envelope mass being convective and as a consequence almost isentropic. At the boundary of the interior convection zone we have $\nabla_{\rm rad}=\nabla_{\rm ad},$ which is equivalent to

[TABLE]

In the exterior radiative region only the first equality in the above holds. For an opacity law $\kappa\propto T^{\eta},$ with $\eta<4,$ this implies that the radiative region tends to a state for which $P\propto T^{4-\eta}.$ Convection will then arise if

[TABLE]

In our case $\Gamma_{2}\sim 1.4$ so convection ensues if $\eta\gtrsim 1/2.$ This is found to become satisfied with the onset of convection always occurring when $T\sim 1.8\times 10^{3}K$ (see above). Given convection is efficient, the model can be assumed to be isentropic to a good approximation. Equation (15) and the above discussion will then imply that the structure depends only on the product of the luminosity and the opacity scaling. As the rate of evolution is $\propto L_{tot}$ (see equation 14) which is determined as a function of the envelope mass in the course of solution, this must then be $\propto 1/\kappa$ as we have found.

Note too that the model radii at corresponding times, taking into account the opacity scaling, are almost independent of this scaling. Thus we find that model B3 with initial envelope mass reduced by a factor of 10 and the same model with opacity reduced by a factor of 100 differed by < 2% at corresponding times ( $3\times 10^{9}$ yr in former case and $3\times 10^{7}$ yr in the latter). This indicates insensitivity to the surface boundary condition.

We have seen that the opacity scaling is important for determining the rate of evolution of the isolated planets. At late stages the envelope becomes cool and dense. For most of the calculations presented above, as has been commonly done (see e.g. Valencia et al., 2013), we have effectively employed an extrapolation of results obtained for much lower densities to obtain opacities. This is clearly uncertain. However, in spite of this we can determine the end state of the evolution with the opacity scaling then determining the time taken to attain this. Consequently, the proximity of any configuration to its end state at a given time is related to the magnitude of the opacity.

6.2 End states

As we have seen the above models contract at a progressively lower rate with decreasing luminosity. The end state is expected to be one with zero internally generated luminosity for which the envelope will be isothermal with a temperature determined by the external radiation field. We shall assume the envelope is optically thick so that this radiation can be effectively absorbed and provide heating. As indicated in Section 3.1 we have determined the end state configurations for the models listed in table 2.

Results are given in table 3 which gives the mean density in the final state for each of the models. This is seen to range between 0.67 and 1.13 g cm*-3*. We also give the ratio of the radius after 2 Gyr to that of the end state radius. This number ranges between 1.37 and 2.02, giving an indication of the degree of inflation that may still exist after that time. Thus evolving protoplanets may exhibit a range of mean densities and radial inflation which depend on the initial envelope mass and the form of the opacity, the latter determining the proximity to the end state.

7 Discussion and Conclusions

We have presented the results of simulations of gas accretion onto planetary cores that are embedded in the viscously and thermally evolving protoplanetary disc models described in Coleman & Nelson (2016a). The calculations consist of a temporal sequence of hydrostatic planetary envelope models in which the gas accretion rate at each time step is determined by the luminosity of the planet (Papaloizou & Nelson, 2005). We have considered in situ gas accretion onto cores with masses of 5, 10 and 15 $\,{\rm M}_{\oplus}$ located at orbital radii 0.1, 0.2, 0.5, 1, 2, 5, 10 $\,{\rm au}$ within disc models whose masses are equivalent to 1, 3 and $5\times{\rm MMSN}$, giving rise to a total of 63 models. The aim is to examine how much gas can accrete onto planetary cores as a function of their mass and orbital radius within evolving disc models whose lifetimes span the range that is inferred from observations of young stellar clusters (e.g. Haisch et al., 2001). For simplicity we adopt standard interstellar opacities (Bell & Lin, 1994; Bell et al., 1997) and a standard equation of state (Saumon et al., 1995) and leave a more complete exploitation of the effects of varying these for future work. For a subset of the planet models that have not undergone runaway gas accretion at the time of protoplanetary disc dispersal, we evolve the planets forward in time for Gyr time scales to follow their long term cooling and contraction. Our main results can be summarised as follows.

(i) For a fixed core mass and disc model, gas envelope accretion occurs more rapidly further from the central star where the gas temperature is lower. This is in agreement with recent studies of gas accretion onto cores located at large distances from the star which also reached the same conclusion (Piso & Youdin, 2014; Piso et al., 2015). We also find that during the slow gas accretion phase, which occurs prior to runaway gas accretion, the mass growth histories of planets at fixed orbital radii only vary very slightly between the different disc models, with higher mass discs generally favouring slightly more rapid gas accretion. Only once the planets undergo runaway gas accretion and begin to accrete gas from the disc at the viscous supply rate, do the mass evolution tracks differ substantially. These differences arise due to the remaining disc mass and lifetime, where giant planets undergoing runaway gas accretion in more massive discs have more time to increase their masses than those that start to undergo runaway gas accretion at roughly the same time in lower mass discs.

(ii) None of the planets that we considered were able to undergo runaway gas accretion at an orbital radius of $0.1\,{\rm au}$, due to the long accretion times that we obtained in this hot region of the disc where the temperature remains close to 1000 K throughout most of the evolution. At this temperature the opacity of the outer radiative zone of the planetary envelope remains dominated by dust grains, and so the Kelvin-Helmholtz contraction of the envelope takes longer than the disc model lifetimes of between 4.6 Myr and 9.4 Myr. There have been other recent studies of in situ gas accretion that have resulted in planets undergoing runaway gas accretion close to the star that at first sight appear to contradict our results. Batygin et al. (2016) presented models in which $15\,{\rm M}_{\oplus}$ cores located at $0.05\,{\rm au}$ underwent runaway gas accretion within $\sim 1$ Myr, but these models were computed under static nebula conditions quite different to ours. Their disc temperature remained at 1500 K, above the sublimation temperature of grains such that they did not contribute to the opacity, and it is well known from previous studies that opacity reductions increase gas accretion rates onto planetary cores (e.g. Hubickyj et al., 2005; Papaloizou & Nelson, 2005; Movshovitz et al., 2010; Mordasini et al., 2014). A significantly larger local gas disc density was also assumed. We note, however, that Batygin et al. (2016) also considered gas accretion onto 4 and $10\,{\rm M}_{\oplus}$ cores, and in these cases found that only moderate gas accretion occurred, in accord with our results. In situ gas accretion at $0.1\,{\rm au}$ has also been considered by Lee et al. (2014), resulting in more rapid gas accretion onto $10\,{\rm M}_{\oplus}$ cores than we observe, again under nebula conditions very different to ours. In their study a static version of the so-called minimum mass extrasolar nebula (Chiang & Laughlin, 2013) was adopted, which has a much higher density and pressure at $0.1\,{\rm au}$ than our viscously and thermally evolving models, so it would appear that the differences in outcomes and conclusions in their study arise because of the different model assumptions. In particular we remark that the large optical depths of cores in such very massive discs tends to increase the temperature at the disc-protoplanet interface, which in turn acts to drive cores to criticality more easily than those indicated above.

(iii) As we move further from the star, the planets with 10 and $15\,{\rm M}_{\oplus}$ cores start to undergo runaway gas accretion and become gas giant planets. The planets with $15\,{\rm M}_{\oplus}$ cores all undergo runaway gas accretion for orbital radii $\geq 0.5\,{\rm au}$, and once the orbital radius reaches $5\,{\rm au}$ all but one of the $10\,{\rm M}_{\oplus}$ cores experience runaway gas accretion. Taking these results at face value, without considering how reductions of the opacity or other factors could modify the results, we conclude that in situ formation of gas giant planets at $0.1\,{\rm au}$ without migration probably requires the presence of a core that exceeds $15\,{\rm M}_{\oplus}$, but formation at orbital radii $\sim 0.5\,{\rm au}$ can occur for core masses $\sim 15\,{\rm M}_{\oplus}$. Viewed within the context of our models, the existence of hot Jupiters orbiting with short periods can either be explained as a consequence of formation further from the star via gas accretion and further solid accretion onto cores with masses $\sim 10$–$15\,{\rm M}_{\oplus}$, followed by migration to their current locations, or by in situ gas accretion onto cores more massive than $15\,{\rm M}_{\oplus}$. Both of these scenarios are consistent with the typical core/heavy element masses of $\geq 20\,{\rm M}_{\oplus}$ inferred to exist within hot Jupiters by the numerous studies that have examined this point (Guillot et al., 2006; Burrows et al., 2007; Miller & Fortney, 2011; Thorngren et al., 2016).

(iv) None of the $5\,{\rm M}_{\oplus}$ cores in our simulations accreted enough gas to undergo runaway gas accretion at any location in the disc. An obvious conclusion to draw from this is that the systems of short-period Super-Earths and Neptunes that have been observed in abundance (e.g. Fressin et al., 2013) probably accreted moderate gas envelopes without undergoing runaway gas accretion, either in situ or further from the star followed by inward migration, because their solid cores have masses $\leq 10\,{\rm M}_{\oplus}$, or because the cores accreted to larger values late within the disc lifetime. If this conclusion is correct, then there is likely to be an as yet undiscovered population of Super-Earths/mini-Neptunes orbiting at large radii with moderate gas envelopes sitting on cores with masses $\lesssim 5\,{\rm M}_{\oplus}$.

(v) The mean densities of those planets that were evolved for Gyr time scales after the dispersal of the gas disc reach values of $\bar{\rho}\sim 0.1$–0.3 g cm*-3* after contracting for 2 Gyr. The final end states to which these planets evolve on longer time scales have mean densities in the range $\bar{\rho}\sim 0.7$–1 g cm*-3*. This range of values encompasses the range of densities inferred for the Solar System gas and ice giants, and the values inferred for the super-Earths and Neptunes among the exoplanet population (Marcy et al., 2014; Weiss et al., 2016).

7.1 Future Work and Direction

Whilst this paper addresses the issue of accreting gaseous envelopes on to pre-existing cores in situ, there are many areas of parameter space and model improvements that we intend to explore in future work. These include:

(i) Allowing the planets to migrate whilst accreting gaseous envelopes. Planets of the masses described here can have short migration times relative to disc lifetimes, which could result in planets not having time to undergo runaway gas accretion before they migrate into the star, unless they form at very large radii and/or experience some kind of migration stopping mechanism. Another effect of considering migration would be the ever changing disc conditions as the planet migrates from a cold region to a hot one (or vice versa). Rapidly changing conditions could have significant consequences for the evolution of the atmosphere and the subsequent accretion rates.

(ii) Including a stochastic accretion rate for the core’s mass through the introduction of planetesimals/boulders/pebbles into the models. With the accretion of macroscopic solids, there will also be an accretion luminosity which will contribute to the luminosity of the planet. This contribution to the luminosity will affect the gas accretion rate, and subsequently the structure of the envelope. But note that, as indicated in 2.3, this will lead to lower envelope masses than those found here because we assumed planetesimal accretion to be inoperative. Some of the very large heavy element abundances inferred for gas giant exoplanets such as CoRoT 10-b, HD 80606-b and HD 149026-b (Guillot et al., 2006; Burrows et al., 2007; Miller & Fortney, 2011; Thorngren et al., 2016) suggest that heavy elements are likely to be present in both the cores and the envelopes of some planets, such that this model extension is required to account for at least some fraction of the giant exoplanet population.

(iii) Combine the models presented here with ab-initio models of planet formation similar to those presented in Coleman & Nelson (2014, 2016a, 2016b) to examine the types of planets and planetary systems that emerge with more realistic gas accretion models.

(iv) Improve the modelling of the long term evolution of planets after dispersal of the protoplanetary disc through consideration of opacity changes and the heat capacity of the core that may provide a source of long term heating for the planets as they cool and contract.

Acknowledgements

The simulations presented in this paper utilised Queen Mary’s MidPlus computational facilities, supported by QMUL Research-IT and funded by EPSRC grant EP/K000128/1. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915. We acknowledge the referee, Kaitlin Kratter, whose comments helped to improve this paper.