The diverse lives of massive protoplanets in self-gravitating discs

We use the SPH code seren (Hubber et al., 2011a, b) to treat the disc thermodynamics. The code uses an octal tree to compute gravity and find neighbours, multiple particle timesteps for optimisation, and a 2^nd Runge-Kutta integration scheme. We use a time-dependent artificial viscosity (Morris & Monaghan, 1997) with parameters $\alpha_{\min}=0.1$, $\alpha_{\max}=1$ and $\beta=2\alpha$, so as to reduce artificial shear viscosity. The chemical and radiative processes that determine the disc temperature are treated with the diffusion approximation of Stamatellos et al. (2007) and Forgan et al. (2009). The net radiative heating rate for an SPH particle $i$ is

where $u_{i}$ is the specific internal energy of the particle, and $\rho_{i}$, $T_{i}$, its density and temperature, respectively. The positive term on the right hand side represents heating by the various radiation sources (star, protoplanet, and pseudo-background radiation field; see Section 3.4), and ensures that the gas and dust cannot cool radiatively below the pseudo-background temperature $T_{{}_{\rm A}}$. $\sigma_{{}_{\rm SB}}$ is the Stefan-Boltzmann constant, $\bar{{\Sigma}}_{i}$ is the mass-weighted mean column-density, and $\bar{\kappa}_{{}_{\rm R}}(\rho_{i},T_{i})$ and ${\kappa_{{}_{\rm P}}}(\rho_{i},T_{i})$ are suitably adjusted Rosseland- and Planck-mean opacities (see Stamatellos et al., 2007, for details). The method takes into account compressional heating, viscous heating, heating by the background radiation field, and radiative cooling/heating. The method has previously applied to disc studies (Stamatellos & Whitworth, 2008, 2009b; Stamatellos et al., 2011a) and has given similar results with grid-based computational methods (Boley et al., 2006; Cai et al., 2008). The method used has been shown that may overestimate the column density through which the disc gas cools (Wilkins & Clarke, 2012; Young et al., 2012; Lombardi et al., 2015) by a factor of a few; however, considering that the dust opacity in discs is not known accurately enough, this limitation of the method is not critical for the results presented in this paper. The detailed treatment of radiative transfer is important when considering how the gas fragments but it is not expected to be crucial (at least qualitatively) for the issues relating to the mass growth and migration of protoplanets already formed in the disc that are discussed in this paper.

We use two prescriptions for the opacity the disc (see Figure 1):

1. 1.

   The Bell & Lin (1994) opacity which is parameterized as

   $${\kappa}(\rho,T)=\kappa_{{}_{\rm 0}}\ \rho^{l}\ T^{m}\,,$$

   (4)

   where $\kappa_{{}_{\rm 0}}$, $l$, $m$ are constants that depend on the species and the physical processes that contribute to the opacity at different densities and temperatures (note that Planck mean and Rosseland mean opacities are assumed to be the same). For the temperatures in the discs we simulate here (up to $\sim 1,000$ K) the opacity is due to dust grains that at low temperatures are coated with ice, which evaporates at higher temperatures. The dust grains start to evaporate at around 1,000 K.

2. 2.

   The Semenov et al. (2003) model includes similar physical processes but the dust considered includes spherical and aggregate particles of various sizes that consist of ice, organics, troilite, silicates and iron. This composition is thought to be more appropriate for the conditions found in protoplanetary discs.

The Semenov et al. Rosseland mean opacity is similar to the Bell & Lin one at temperatures up to $\sim 100$ K but higher at larger temperatures by up to a factor of 2. At such temperatures the optically thick parts of the disc cool less efficiently when using the Semenov et al. opacity. The Semenov et al. Planck mean opacity is larger by a factor of 5 than the Bell & Lin one, which means that the optically thin parts of the disc cool more efficiently.

The above two sets of opacities are widely used in the literature. However, the composition of dust in discs is rather uncertain, especially since significant dust growth may occur, lowering the opacity. It is unclear whether dust growth is significant at the early disc phase that we study here.

Sink particles are used to represent the central star and the protoplanet (Bate et al., 1995). Sink particles interact with the rest of the computational domain only through their gravity (and their luminosity if their radiative feedback in considered, see below). The sink radius of the central star is $R_{\rm sink,\star}=0.2$ AU, and the sink radius of the protoplanets is $R_{\rm sink,p}=0.1$ AU. This value is chosen to be smaller than the Hill radius of the protoplanet, i.e. the region around it where the its gravity dominates over the gravity of the star:

For example, when a Jovian planet is at 50 AU away from the central star then its Hill radius is $R_{\rm H}\sim 3.5$ AU. The Hill radius may increase as the protoplanet accretes material from the disc or it may decrease as the protoplanet moves closer to the central star.

Gas particles accrete onto a sink when they are within the sink radius and bound to the sink (Hubber et al., 2011b). The accretion rate onto the sink is computed by calculating the increase of the sink mass over a time interval that is comparable to the dynamical time of each particle at the edge of of the sink.

The radiation feedback from the protoplanet is taken into account (in a subset of the simulations) using the method described in Stamatellos (2015). We use the method of Stamatellos et al. (2011b) & Stamatellos et al. (2012) that invoke a pseudo-ambient radiation field with temperature $T_{{}_{\rm A}}^{\rm planet}({\bf r})$ due to radiation from the protoplanet (see also Stamatellos et al., 2007; Stamatellos & Whitworth, 2009a). This pseudo-ambient temperature sets the minimum temperature that the gas can attain if it cools radiatively. The contribution from the protoplanet is set to

where ${\bf r}$ is the position on the disc midplane, and $L_{p}$ and ${\bf r}_{p}$ are the luminosity and position of the protoplanet, respectively. The above approximation may overestimate the effect of the feedback in optically thick regions of the disc, that are well-shielded from the protoplanet (see discussion in Mercer & Stamatellos, 2017).

The luminosity of the protoplanet is set to

where $M_{p}$ is the mass of the protoplanet, $\dot{M}_{p}$ is the accretion rate on to it, and $R_{\rm acc}$ the assumed accretion radius. $f=0.75$ is the fraction of the accretion energy that is radiated away at the photosphere of the protoplanet, rather than being expended driving jets (Machida et al., 2006; Offner et al., 2009), or deposited within the protoplanet (e.g. Baraffe et al., 2017).

The exact amount of the energy radiated away from the protoplanet depends on the detailed properties of the accretion shock around the protoplanet (Zhu, 2015; Marleau et al., 2017; Szulágyi & Mordasini, 2017; Mordasini et al., 2017). If this protoplanet has formed by gravitational instabilities in the disc then the accretion happens onto the second hydrostatic core (e.g. Larson, 1969; Stamatellos & Whitworth, 2009b). The radius of the second core is uncertain: $\sim~1-20~{\rm R}_{\sun}$ (Masunaga & Inutsuka, 2000; Tomida et al., 2013; Vaytet et al., 2013; Bate, 2014; Tsukamoto et al., 2015). Here we assume $R_{\rm acc}=3~{\rm R}_{\sun}$, but the results of this work are not sensitive on this assumption. The accretion luminosity of the protoplanet is significant due to the relatively high accretion rate onto it during the initial stages of its evolution and can dominate over the stellar luminosity in the disc region around the protoplanet (Owen, 2014; Stamatellos, 2015; Montesinos et al., 2015), especially for the Jupiter-like planet-seeds in high-mass discs that we study here.

If the protoplanet has formed by disc fragmentation, on top of the significant accretion luminosity, it will also radiate energy due to its contraction. This energy release is on the order of 0.1 L_☉ (Inutsuka et al., 2010), which is similar to the one expected from accretion. Therefore, the protoplanet’s luminosity is deemed to play an important role in its evolution (Nayakshin & Cha, 2013; Stamatellos, 2015; Benítez-Llambay et al., 2015). We do not take this radiation explicitly into account in the simulations presented in this paper but we have performed runs with smaller $R_{\rm acc}$ resulting in higher luminosity for the protoplanet. These runs give qualitatively similar results to the ones we discuss in this paper.

The radiation feedback from the star is taken into account (in a subset of the simulations) by invoking an additional pseudo-ambient temperature,

where $R$ is the distance from the central star on the disc midplane, and $T_{0}$ the temperature at $R=1$ AU from the star. The total pseudo-ambient temperature (due to the central star and the protoplanet) is then set to

where $T_{{}_{\rm A}}^{\rm planet}({\bf r})$ is the contribution from the protoplanet. Note that the radiative feedback from the star is fixed in time for simplicity, whereas the radiative feedback from the protoplanet is variable and depends on the accretion rate of gas onto it as it moves within the disc.

The protoplanet grows in mass initially very fast as it opens up a gap but its mass increase slows down once the gap is opened up (see Figure 4). The accretion rate onto the protoplanet (see Figure 4, top) is relatively high ($\sim~10^{-4}-10^{-2}\ {\rm M_{J}\ yr^{-1}}$) during the gap opening phase but then it drops down; nevertheless, accretion continues through streams within the gap (Lubow & D’Angelo, 2006). Similarly high accretion rates are also seen in previous studies (D’Angelo & Lubow, 2008; Ayliffe & Bate, 2009; Zhu et al., 2012; Gressel et al., 2013). The resulting luminosity of the protoplanet (see Figure 4, bottom) could rival the luminosity of the young star making its detection easier in terms of the sensitivity required. However, this phase of high accretion lasts only for a relatively short time, making such a detection difficult. These luminosity estimates are only up to a few times higher than the luminosities estimated from hot-start models of planet formation (Marley et al., 2007; Mordasini, 2013; Mordasini et al., 2017). The value of the luminosity depends on the accretion radius, which in our model is set assuming that gas accretion happens onto the second core formed after the temperature at centre of the protoplanet-precursor clump reaches $\sim 2000$ K and molecular hydrogen dissociates (e.g. see Stamatellos & Whitworth, 2009b). The accretion rate onto the protoplanet decreases after the gap is opened up but thereafter it shows many spikes, some with periodicity similar to the orbital period of the planet at each specific time, indicating periodic interactions with spiral structures in the disc that drive gas accretion onto the protoplanet. Such spikes are absent when the radiative feedback of the planet is taken into account (Run 2) as strong spiral features are suppressed. There is a delay of $\sim 2-4$ kyr between the gap opening and the gap edges becoming gravitationally unstable, driving accretion onto the protoplanet. The delay is longer for protoplanets that are closer to the central star, i.e. in a warmer region of the disc, as more gas needs to accumulate for the gap edges to become unstable. The accretion rate onto the protoplanet is generally lower when the protoplanet’s radiative feedback is included in the simulation as seen in previous studies (Nayakshin & Cha, 2013; Stamatellos, 2015; Gárate et al., 2017).

The way that gas is accreted onto the protoplanet depends on the dynamical state of the protostellar disc. In Figure 6 we plot the difference between the angular momentum of the gas entering the protoplanet’s Hill sphere from inside the protoplanet’s orbit ($L_{\rm in}$), and the angular momentum of the gas entering the Hill sphere from outside the protoplanet’s orbit ($L_{\rm out}$). To calculate these we use the mass of gas within distance $1<r<1.2~R_{\rm H}$¹¹1 The same result is obtained for a different outer limit, e.g. if we consider the gas within distance $1<r<1.5~R_{\rm H}$ from the protoplanet. from the protoplanet that moves towards the protoplanet from inside its orbit ($m_{\rm in}$), and the gas that moves towards the protoplanet from outside its orbit ($m_{\rm out}$) (Figure 6, bottom). During the initial gap opening phase, in both cases, gas is accreted onto the protoplanet mostly from the inner disc. After the initial gap opening phase ($\sim 5$ kyr), we see that for the run including the protoplanet’s radiative feedback (Run 2, red), in which the disc is gravitationally stable because of the effect of the protoplanet’s feedback, most of of the accreted gas comes from inside the protoplanet’s orbit. The total accreted gas in this case has lower angular momentum than the protoplanet, so that the protoplanet’s angular momentum decreases. On the other hand, in the case without the feedback (Run 1, black), the disc is gravitationally unstable and a high fraction of the accreted gas (up to almost 50%) comes from outside the protoplanet’s orbit. This gas has higher angular momentum than the protoplanet.

This can also be seen in Figure 6 where we plot the change of angular momentum (top) and mass (bottom) of the protoplanet ($\Delta L_{p}$, $\Delta M_{p}$), the inner disc (i.e. the disc inside the protoplanet’s orbit; $\Delta L_{\rm inner\ disc}$, $\Delta M_{\rm inner\ disc}$), and the outer disc ($\Delta L_{\rm outer\ disc}$, $\Delta M_{\rm outer\ disc}$). In Run 1 the protoplanet increases in mass while the outer disc mass decreases, whereas in Run 2, in which the protoplanet’s feedback is included, the outer disc mass remains almost steady (after 5 kyr) whereas the inner disc mass decreases as the mass of the protoplanet increases.

Apart from the run with the protoplanet’s radiative feedback, there are only relatively subtle differences in the accretion onto the protoplanets in the other runs. As expected, for a higher disc viscosity the gap opens up at a later time (Crida et al., 2006), but thereafter the accretion rate into the protoplanet is higher, resulting in a higher final mass.

The gap opens up faster in the run with the Semenov et al. (2003) opacities, which are generally larger than the Bell & Lin (1994) opacities, and the accretion rate thereafter is larger by a factor of $\sim 2$ than in the other runs. This is because the Semenov et al. (2003) opacities correspond to more efficient cooling (see Nayakshin, 2017b), and therefore: (i) the disc is cooler, enabling the quick opening of the gap, and (ii) the thermal pressure of the circumplanetary disc is lower and gas accretion is not opposed as much as in the other runs. The dependence of the accretion rate and the final mass of the object on the assumed opacity seems to be different from what is seen in Nayakshin (2017b). However, what matters is not the opacity itself but the cooling of the gas around the protoplanet. Nayakshin (2017b) employs an opacity that regulates gas cooling both in optically thin and optically thick regimes (i.e. the Planck-mean and Rosseland-mean opacities are the same), whereas the Semenov et al. (2003) opacity account of the differences between these regime. The higher Planck-mean opacity of Semenov et al. (2003) results in more efficient cooling in the optically thin regions (see Equation 3; in the optically thin limit the denominator on the right hand side is dominated by the $\kappa_{P}^{-1}$ term). Therefore as in Nayakshin (2017b) our results mean that the gas accretion onto the protoplanet with more efficient cooling is higher. Similar results are also found by Ormel et al. (2015a, b) when studying accretion onto massive solid planet cores.

We see that in all cases, apart from the one in which the protoplanet radiative feedback has been included, the protoplanet’s final mass (i.e. the mass at the end of the hydrodynamic simulation; $t=20$ kyr) is in the brown dwarf regime ($>16.3~{\rm M_{J}}$). The protoplanet’s mass in the case in which its radiative feedback is taken into account is at the borderline between the planetary and brown dwarf regime. However, it is expected that the protoplanet will continue to accrete gas and eventually become a brown dwarf (unless the disc dissipates by other processes, like e.g. photo-evaporation, on a faster timescale). Nevertheless, when the protoplanet’s radiative feedback is taken into account the mass growth of the protoplanet has been suppressed and, at the end of the simulation, its mass is lower by a factor of $\sim 2$ than in the other runs.

The protoplanet initially interacts strongly with the disc resulting in fast inward migration, similar to the Type I migration for low-mass planets in low-mass discs. The migration timescale is only $10^{4}$ yr and it is similar for all 5 runs (see Figures 7, 8 and Tables 7, 7). However, once a gap opens up the migration stalls; this happens at different times for each run. When the planet orbit plane is not the same with the disc midplane (Run 3) it is more difficult for the protoplanet to open up a gap. In the same way, the more viscous is the disc (Run 4) the more difficult is to open up a gap (e.g Crida et al., 2006; Malik et al., 2015). In both cases there is delay in opening up a gap and the protoplanet migrates further inwards (but only by $\sim 4$ AU compared with the other runs). In the four runs without radiative feedback from the protoplanet, the protoplanet switches to a slow outward migration. At the end of the simulations (20 kyr) the protoplanet’s semi-major axis has moved to 40, 38, and 36 au, for Runs 1, 2, and 3, respectively, from an initial semi-major axis of 50 au. In Run 5 with the higher opacity (faster cooling) (Semenov et al., 2003) the gap opens up quicker and the protoplanet starts migrating outwards reaching a semi-major axis of 53 AU at the end of the simulation.

The case that it is distinctly different is the one in which the radiative feedback from the planet is taken into account (Run 2). In this case the opening up of the gap is more difficult as the protoplanet heats its circumplanetary disc but also the host protostellar disc (Marley et al., 2007; Nayakshin & Cha, 2013; Benítez-Llambay et al., 2015; Montesinos et al., 2015; Stamatellos, 2015; Szulágyi & Mordasini, 2017). Therefore the gap opens at a later time ($\sim 1$ kyr later) and the protoplanet migrates closer to central star. Its inward migration continues after the gap is opened up but at a much longer timescale ($\sim 10^{5}$ yr; typical of Type II migration). At the end of the simulations (20 kyr) the protoplanet is at 18 AU from the central star.

Lin & Papaloizou (2012) and Cloutier & Lin (2013) have shown that when a planet resides in a gap that has gravitationally unstable edges, gas is brought into the planet’s co-orbital region by the spiral arms, resulting in a positive torque that pushes the planet outwards. This is indeed what we see here in the runs without the radiative feedback from the protoplanet. In these runs, a significant amount of the gas accreted onto the protoplanet is high angular momentum material coming from outside the protoplanet’s orbit (see Figs. 6-6). On the other hand, when the radiative feedback from the protoplanet is taken into account the temperature at the edges of the gap is raised, they are not gravitationally unstable anymore, and the protoplanet accretes a higher fraction of lower angular momentum gas from the disc inside its orbit.

Contrary to previous studies of (Baruteau et al., 2011; Michael et al., 2011; Malik et al., 2015), we allow the protoplanet to accrete gas and grow in mass. By growing in mass, it is more capable to open up a gap. In previous studies mass is accumulated in the circumplanetary disc and therefore contributed to strongly couple the protoplanet with the gaseous disc in which it is formed. This behaviour is also seen in the simulations of Nayakshin & Cha (2013). However, the crucial difference with previous studies is that we use a more detailed prescription for the radiative transfer in the disc. Our prescription allows the circumstellar and circumplanetary material to moderate the efficiency of its cooling as it becomes more or less dense, whereas the studies of Baruteau et al. (2011) and Malik et al. (2015) use the $\beta$-cooling approximation, in which the cooling time is proportional to the local orbital period of the gas, i.e. ${t_{\rm cool}}=\beta\ \Omega_{K}^{-1}(R)$. In Section 5 we discuss this issue in detail.

The estimated migration timescales ($\tau_{\rm mig}={\alpha_{p}}/|\dot{\alpha_{p}}|$) are shown in Table 7, and the migration rates ($v_{\rm mig}=\dot{\alpha_{p}}$) inTable 7. These values are calculated at 5 different times during the simulations. The migration timescales for all 5 runs are plotted against the protoplanet mass in Figure 8. On the same graph we plot the estimated migration timescales for Type I and Type II migration as calculated by Ward (1997) and Tanaka et al. (2002) (see also Bate et al., 2003).

The Type I migration rate estimated by Tanaka et al. (2002) is

where $\alpha_{\Sigma}$ is the exponent of the disc surface density profile ($\Sigma\propto R^{-\alpha_{\Sigma}}$), $R_{p}$ the orbital radius of the planet, $\Sigma_{p}$ the surface density of the disc at the position of the planet, $H_{p}$ the disc thickness at the position of the planet, and $\Omega_{p}$ the planet’s angular frequency. The rate of migration of the Type II case is set by the viscous evolution of the disc (e.g. Bate et al., 2003):

where $\alpha$ is the disc viscosity parameter. We note though that the actual Type II migration timescale may vary from the above approximation (Crida & Morbidelli, 2007; Duffell et al., 2014; Dürmann & Kley, 2015, 2017). For example, Dürmann & Kley (2015, 2017) find that due to gas crossing through the gap Type II migration may be faster or slower. The total rate of migration, combining the Type I and Type II regimes, is therefore

$M_{t}$ is the transition mass between Type I and Type II migration (Ward, 1997) and is set to

The migration timescale is therefore

We calculate the migration timescales from the above equations assuming a planet at $R_{p}=5.2$ AU in a disc with $\alpha_{\Sigma}=0.5$, $\Sigma_{p}=75\ {\rm g\ cm^{-2}}$, and $H_{p}/R_{p}=0.05$. In Figure 8, we plot 3 cases that correspond to different disc viscosity ($\alpha=0.1,0.01,0.002$). We note that these calculations are for planets in low-mass discs at specific regimes: Tanaka et al. (2002) study 3D discs including both Lindbland and corotational resonances assuming an isothermal disc, whereas Ward (1997) study a 2D disc ignoring corotational resonances. These calculations are not applicable for the cases we examine in this paper. Therefore, these analytic solutions are plotted in the graph in Figure 8 merely for reference, so as to put our results in context with the classic picture of Type I and Type II migration.

The migration timescales in Figure 8 are plotted for 5 different times in the simulation: 0.7, 1.5, 2.5, 5 and 18 kyr. As the mass of the protoplanet increases with time, the same graph may be used to track the evolution of the migration timescale with time. In all runs the protoplanet initially migrates inwards on timescales similar to the Type I migration timescale of a $1-{\rm M_{J}}$ planet, i.e. $\sim(1-2)\times 10^{4}$ yr. As the protoplanet moves towards the host star and grows in mass, the migration timescale becomes shorter (down to $\sim 3\times 10^{3}$ yr), but once the gap opens up the migration slows down and eventually the protoplanet starts moving outwards in all runs apart from the case when its radiative feedback is taken into account. In this case, inward migration continues albeit at a much longer timescale, similar to the Type II migration case ($\sim 2\times 10^{5}$ yr).

The eccentricity of the protoplanet grows as the protoplanet interacts with the gravitationally unstable disc (see Figure 9). The growth pattern is rather stochastic with the eccentricity growing up to $\sim 0.15$. Even if the eccentricity after this initial growth period dampens with time due to secular interactions with the disc, these initial interactions can provide the seed for subsequent eccentricity growth (Goldreich & Sari, 2003; Duffell & Dong, 2015). Therefore, disc-planet interactions, while the disc is still relatively massive, may help explain the observed high eccentricities of giant planets in systems that are thought to contain only one planet (Wright et al., 2011; Dunhill & Stamatellos, 2018). In the case when the radiative feedback of the protoplanet is taken into account (Run 2), the eccentricity initially grows during the gap opening phase but then it is dampened quickly, as there is no strong stochastic driving, so that the protoplanet ends up in a nearly circular orbit. Therefore, interactions within a gravitationally unstable disc is necessary for eccentricity growth to occur.

The protoplanet in the run in which its orbit is inclined (Run 4, blue lines) shows a mass growth that resembles closely the standard run (see Figure 4). This is because the protoplanet’s orbit inclination gets smaller quickly as the protoplanet crosses through the disc midplane. Within $\sim 2$ kyr ($\sim 5$ orbits) the protoplanet’s orbit has been aligned with the disc midplane (inclination has become zero, see Figure 10). Therefore, the long-term evolution of a protoplanet is not significantly different if the protoplanet initially forms on an inclined orbit in a relatively massive protostellar disc.

The surface density, the temperature profile and the Toomre parameter Q for the region around the protoplanet for all five runs are shown in Figure 11 (on the protoplanet’s corotational frame). The snapshots correspond to $t=8$ kyr, i.e. after the gap has been opened up in the disc. The gap size is a few Hill radii and rather similar for all runs (note though that the size of the Hill radius is different in each run, as the protoplanet is on different orbits in different runs). The increase of surface density towards the protoplanet within the Hill radius corresponds to the circumplanetary disc (see section below).

The surface density at the boundaries of the Hill radius is asymmetric (Figure 11a). In Run 2 (which includes the protoplanet’s radiative feedback) the surface density at the inner boundary of the protoplanet’s Hill radius (e.g. at $-1.5~{\rm R_{H}}$) is higher than the surface density at outer boundary of the Hill radius (e.g. at $+1.5~{\rm R_{H}}$) whereas in all other runs the outer boundary surface density is higher. This is consistent with the point made in the previous subsection, i.e. that a large fraction of gas is accreted onto the protoplanet from outside its orbit when the gap edges are gravitationally unstable. The asymmetry is more pronounced for Run 5 (higher opacity run) which also exhibits shorter outward migration timescales (see Table 7 and Figure 8). The disc temperature (Figure 11b) at the run with the protoplanet feedback is higher than in the other runs, which results in stabilizing the disc edges as the Q parameter is larger than $\sim 3$ (Figure 11c). The temperature of the disc around the protoplanet is lower for the higher opacity run (Run 5) which helps opening up the gap fast and results in higher gas accretion onto the protoplanet.

The properties of the circumplanetary discs in all 5 runs are shown in Figure 12 and the time evolution of the Hill radius and the mass within it for the protoplanet in each run are shown in Figures 14 and 14. The typical circumplanetary disc mass is about $0.1~{\rm M_{J}}$, which means that the disc is resolved by $\sim 1,000$ SPH particles, i.e. approximately 7 smoothing lengths (assuming a nearly 2D disc). The sink radius of the protoplanet is set to 0.1 AU so it is always smaller than the Hill radii of the protoplanets at all times, in all runs (see Figure 14). The circumplanetary discs are therefore just resolved and the presented properties near the planet are probably resolution dependent. We expect that the surface density and the temperature near the planet are underestimated. In Run 2 (with protoplanet feedback) the circumplanetary disc is smaller (with mass of $0.01~{\rm M_{J}}$; see Figure 14) as the protoplanet is closer to the parent star; therefore the circumplanetary disc is just resolved in this case (only by $\sim 2$ smoothing lengths).

Despite the above drawback, by comparing the derived properties we see that in Run 2, with the protoplanet radiative feedback, the circumplanetary disc is hotter than in the other runs, as expected. Radiative feedback from the protoplanet together with the gas thermodynamics (Gressel et al., 2013) are important in determining the properties of circumplanetary discs. Further studies with higher resolution are needed for more secure results. Additionally, the presence of magnetic field may also play a significant role (Gressel et al., 2013; Fujii et al., 2014) but we do not examine this case here.

We note that we do not find temperatures in the circumplanetary discs higher than a few hundred Kelvin in contrast to the simulations of Szulágyi (2017), in which they find that temperatures close to the protoplanet rise up to a few thousand Kelvin. However, their simulations are able to resolve the region down to $10^{-3}\times{\rm R_{Hill}}$, which can indeed become very hot. We note though that the Szulágyi (2017) simulations do not include the effect of molecular hydrogen dissociation at $\sim 2,000$ K (nor the ionisation of hydrogen and the first and second ionisation of helium), so they may overestimate the temperatures in the inner region. When we compare the temperatures farther out in the circumplanetary disc (e.g. at $0.5\times{\rm R_{Hill}}$) we see that the temperatures that we find here are lower than the ones in Szulágyi (2017).

The presence of a protoplanet in the disc regulates the accretion rate onto the star (see Figure 15). Interestingly, the accretion rate onto the protoplanet during the initial stages of its evolution that we model here, is comparable to the accretion rate onto the star itself. During the gap opening phase it can be even ten times higher. This is consistent with previous studies of planets embedded in lower-mass discs (Nelson et al., 2000; Lubow & D’Angelo, 2006; Owen, 2014). During this phase the accretion rate onto the star decreases but as the protoplanet moves closer to the star and once the gap opens up, the accretion rate onto the star increases by a factor of a few, as the protoplanet acts to drive accretion onto the star; the closer the protoplanet orbit to the star and the higher its mass, the higher the accretion rate onto the star becomes. At this stage, and for a few thousand years, the accretion rate onto the star may become higher than the accretion rate onto the protoplanet. In Run 5 (with the Semenov opacities) the protoplanet stays sufficiently away from the star so that the accretion rate onto the star is not affected.

In Run 2 (with protoplanet feedback) the protoplanet migrates closer to the star than in the other runs, and the accretion rate onto the star increases significantly until a cavity forms around the star; thereafter, the accretion rate drops considerably as the presence of the protoplanet starves the star from gas. At the same time the protoplanet experiences similar gas starvation as it resides within the same cavity. This stage can be thought as similar to the transition disc phase (see review by Espaillat et al., 2014). However, the accretion rate onto the protoplanet is still higher than the accretion onto the star by a factor of 2. In this case the star-protoplanet system behaves as a low mass-ratio binary system with a circumbinary disc, in which the secondary component (i.e. the protoplanet) accretes more than the primary component (i.e. the star) (Artymowicz & Lubow, 1994, 1996). In such systems, secondaries increase in mass faster than primaries, resulting in an almost equal-mass binaries (e.g. Satsuka et al., 2017). However, in the case of a planetary-mass companion, its mass cannot become comparable to that of the star; even if an unlikely high accretion rate of $\sim 10^{-4}{\rm M_{J}\ yr^{-1}}$ is maintained for $\sim 10^{5}$  yr, then the protoplanet’s mass will increase only by $10~{\rm M_{J}}$.

The actual accretion rate onto the protoplanet is important as it regulates the strength of its radiative feedback, which it turn determines whether the edges of the gap opened by the protoplanet are gravitationally unstable (so that the protoplanet migrates outwards), or gravitationally stable (so that protoplanet migrates inwards). The luminosity of the protoplanet, $L_{p}$, is proportional to the accretion rate, ${\dot{M}_{p}}$, onto the protoplanet (see Equation 7), therefore the temperature due to the presence the protoplanet (see Equation 6) scales as $T_{{}_{\rm A}}^{\rm planet}(r)\propto{\dot{M}_{p}}^{1/4}$. The Toomre parameter, $Q$ that determines whether the gap edges are gravitationally unstable is $Q(r)\equiv{\kappa c_{s}(r)}/{\pi G\Sigma(r)}$, where $\kappa$ is the epicyclic frequency, and $c_{s}$ the sound speed (we assume that the distance, $r$, is measured from the protoplanet). Hence, $Q\propto\left[T_{{}_{\rm A}}^{\rm planet}(r)\right]^{1/2}$ or equivalently

Therefore, there is only a weak dependence of the Toomre parameter on the accretion rate onto the protoplanet. If we assume that the actual accretion rate is 10 times lower than the one we estimate in the simulations we present here, then the Q value is lower only by a factor of $\sim 1.3$ (assuming that the other parameters remain the same). Then the Q value in the outer edge of the gap in the run with the protoplanet radiative feedback (see Figure 11, red line) will drop from 3 to 2.3, i.e. the gap edge will still be gravitationally stable ($Q>1.5$). We conclude the value of the accretion rate onto the protoplanet is not critical, at least qualitatively, regarding the effect of radiative feedback on the migration of the protoplanet.

Gas accretion is significant during the initial stages of the protoplanet’s evolution, during the gap-opening phase. In most of the simulations presented here, the protoplanet grows to the planet-brown dwarf mass-limit ($\sim$13 M_J) within a few thousand years. Once the gap is opened up, accretion occurs mainly from inside the protoplanet’s orbit. However, a fraction of it is also accreted from outside the protoplanet’s orbit; when the gap edges are gravitationally unstable (low Toomre-Q, i.e $Q\stackrel{{\scriptstyle<}}{{{}_{\sim}}}2$), then an almost equal amount of gas is accreted from outside the protoplanet’s orbit. The accretion rate onto the protoplanet is higher for higher viscosity, and for more efficient cooling (note that cooling is modulated by the opacity). The role of radiation feedback from the protoplanet is critical as it considerably decreases the accretion rate of gas onto it and therefore its final mass. The effect of radiation feedback from the central star is rather minimal (Appendix B). We also found a dependence of the protoplanet’s final mass on its initial position within the disc: a planet that is relatively closer to the central star ($\stackrel{{\scriptstyle<}}{{{}_{\sim}}}20$ AU) quickly opens up a cavity around its parent star and resides near the outer edge of this cavity. In this case, the star-protoplanet system behaves as a low mass-ratio binary system attended by a circumbinary disc. In such a system gas flows to the star and the protoplanet only through streams that pass through the Lagrangian points of the system (Artymowicz & Lubow, 1994, 1996), and therefore accretion is slow. These are the only cases (with or without radiative feedback from the protoplanet) that the protoplanet’s mass remains within the planetary regime at the end of the simulations. As the mass loss of the disc due to accretion onto the star and the protoplanet is lower in this case, the disc is expected to live longer, allowing other processes, e.g. photo-evaporation (Alexander et al., 2006; Owen et al., 2010) or disc winds (Suzuki & Inutsuka, 2009; Suzuki et al., 2010; Suzuki & Inutsuka, 2014; Fromang et al., 2013; Lesur et al., 2013; Bai & Stone, 2013; Gressel et al., 2015; Bai, 2016), to disperse the disc fast enough for the protoplanet to survive as a planet, not a brown dwarf. Another way to suppress mass growth is to eject the protoplanet from the disc (in multiple-planet systems) so that it ends up as a free-floating planet (Li et al., 2015, 2016; Mercer & Stamatellos, 2017).

Disc-planet interactions result in an initial phase of fast inward migration that lasts for a few thousand years, apart from the runs in which the protoplanet is very close to its parent star ($\stackrel{{\scriptstyle<}}{{{}_{\sim}}}10$ AU); in these runs the protoplanet’s orbit changes only slightly as it quickly opens up a cavity around the star. The migration timescale for this initial phase is $\sim 10^{4}$ yr, i.e. similar to the Type I migration timescale established for planet migration in low-mass discs. This is consistent with previous studies (Michael et al., 2011; Baruteau et al., 2011; Malik et al., 2015). However, we find that the protoplanet is always able to open up a gap in the disc, and thereafter the migration pattern diverges from that in previous studies.

We find that if the gap edges are gravitationally unstable (low Toomre-Q, i.e $Q\stackrel{{\scriptstyle<}}{{{}_{\sim}}}2$), then the protoplanet starts migrating outwards on a timescale $\sim 10^{5}$ yr, as a high fraction of the gas it accretes (almost about 50%) is higher angular momentum gas from outside its orbit. If the gap edges are gravitationally stable (i.e. in the runs with protoplanet radiative feedback, or when the protoplanet is closer to the star, in a hotter disc region) then the migration continues to be inwards but on a longer, Type II migration timescale, i.e. $\sim 10^{5}$ yr, with the planet accreting mainly lower-angular momentum gas from inside its orbit.

Contrary to previous studies we have allowed the protoplanet to grown in mass, and, critically, we have used a more sophisticated treatment of the radiative transfer rather than the commonly used $\beta$-cooling approximation (Baruteau et al., 2011; Malik et al., 2015). The method we use allows the gas to modulate its cooling/heating depending on its properties (density, temperature); this enables gap opening and facilitates enhanced mass growth of the circumplanetary disc and of the planet.

In the present calculations this “accretion from outside" that is responsible for initiating outward migration is driven by the gravitational instability in the outer disc. However, it is possible this may be promoted by any mechanism that makes large amplitude fluctuations in the outer disc, such as the magneto-rotational instability (MRI, but see Gressel et al., 2012) or other instabilities.

The eccentricity of a protoplanet grows up to 0.05 during the initial gap opening phase. Thereafter, two patterns are seen: protoplanets that reside within a gap with gravitationally unstable gap edges show further eccentricity growth up to $0.1-0.2$, whereas the eccentricity of protoplanets within a gap with stable edges dampens, making their orbits nearly circular. The eccentricity growth is not monotonic but shows a stochastic pattern, characteristic of both periodic and random interactions with condensations present near the gap edges. In the long term runs ($10^{5}$ yr) we see further eccentricity growth up to $\sim 0.3$ (see Appendix C). Almost in all cases that we examined here we see a minimum eccentricity growth of a few 0.01 which may provide the seed for further eccentricity growth (Goldreich & Sari, 2003; Duffell & Dong, 2015). Therefore, early protoplanet interactions with a massive, but gravitationally stable disc may explain in some cases the high eccentricities of giant planets in single-planet systems (Wright et al., 2011; Dunhill & Stamatellos, 2018). On the other hand, we see that protoplanets ending up in the inner disc region ($\stackrel{{\scriptstyle<}}{{{}_{\sim}}}20$ AU), tend to have circular orbits.

We find that in most cases the accretion rate onto the protoplanet is higher or at comparable to the accretion rate onto the star. If the protoplanet is close enough to the central star, it creates a cavity, resides within it (at its outer edge), and seemingly starves off the central star from gas, accreting most of the incoming gas. This stage shows similarities with the transition disc phase (Espaillat et al., 2014).

The luminosity of a protoplanet can be quite high during the initial stages of its evolution, as the accretion rate on it is relatively high (typically $\sim 10^{-3}{\rm M_{J}\ yr^{-1}}$, but could be up to 10 times higher during the gap opening phase). Assuming a relatively large radius for the protoplanet of $R_{\rm acc}=3~{\rm R}_{\sun}$ (which is expected for a protoplanet that has formed by disc fragmentation), we find a typical accretion luminosity of a few $0.1~L_{\sun}$ (see also Inutsuka et al., 2010). Combined with the fact that the protoplanet starves off its central star from gas, it may make its detection easier.

The high accretion phase lasts only for a few $10^{3}-10^{4}$ yr and therefore only a small fraction of transition discs (assuming that they all are due to planetary companions) are expected to be seen in this high-luminosity phase. If transition discs live for 1 Myr, then, from timescale arguments, we suggest that only a small percentage, $0.01\%$ to $0.1\%$, of transition discs may have observable, high-luminosity protoplanets (or brown dwarfs).

Protoplanets are attended by circumplanetary discs that are fed with gas from the disc. We find that these discs are hotter than their surroundings and therefore easier to observe. Their temperatures are a few hundred Kelvin and they are hotter when radiative feedback from the protoplanet is taken into account. If the protoplanet migrates outwards its circumplanetary disc becomes more massive ($\sim 1~{\rm M_{J}}$) and larger $\sim 10$ AU.

A protoplanet could be detected indirectly through its circumplanetary disc (Szulágyi et al., 2017). The circumplanetary disc of a young luminous protoplanet is hotter than in the case of non-luminous protoplanet (see Figure 12), therefore it would emit significant long-wavelength radiation even if it is not particularly dense. This long-wavelength radiation would escape mostly unattenuated from the system, whereas short-wavelength radiation that is emitted from the protoplanet may be significantly attenuated, especially if the protoplanet orbits a young embedded protostar (e.g. in a Class 0 object).

Our study shows that the final fate of a Jupiter-like planet-seed formed early-on during the life of a protostellar disc depends both on the initial properties of the seed (e.g. its initial orbit), and the physical processes that are at play (e.g. the role of the protoplanet’s accretion luminosity). The protoplanet may end up as a massive giant planet on a circular orbit close to its parent star or as a low- or high- mass brown dwarf on an eccentric wide orbit.

The parameter space has not been fully explored. It is expected that the initial disc mass will play an important role as it will determine to a wide extent whether the disc is close to being gravitationally unstable or not. Additionally, the initial mass of the protoplanet is important; lower-mass planet-seeds (i.e. a few tenths of the mass of Jupiter) will grow in mass slower and may have a better chance to survive as planets (Nayakshin, 2017a; Malik et al., 2015).

We have assumed that a protoplanet has fully formed (i.e. it is a bound object) at the beginning of the simulation. However, if the planet has formed by disc fragmentation then the formation timescale of the clump and how this evolves to a bound object should be taken into account (Boley et al., 2010; Zhu et al., 2012). Previous studies have shown that a large fraction of these pre-protoplanet clumps may be disrupted, with only a few of them surviving (Boley et al., 2010; Tsukamoto et al., 2015; Hall et al., 2017).

The protoplanets is our simulations are represented by sink particles that accrete gas as this approaches within $0.1$ AU from them. The use of sinks in hydrodynamic simulations is needed to avoid small timesteps but it could lead to enhanced gas accretion that may depend on how well the circumplanetary disc is resolved. The study of how the accretion rate onto the protoplanet varies with resolution and with respect to the hydrodynamic method used (e.g. SPH vs grid-based) is important but outside the scope of this paper. However, we note that the accretion rates onto the protoplanets that we obtain in this work are similarly high (within a factor of a few) as the accretion rates obtained from previous studies (D’Angelo & Lubow, 2008; Ayliffe & Bate, 2009; Zhu et al., 2012; Gressel et al., 2013), which use different numerical methods and disc-planet initial setups.

The physical processes involved are also important. We have shown that the disc thermodynamics, i.e. how fast the disc heats and cools, affects the accretion rate onto the protoplanet (see also Nayakshin, 2017b), its ability to open up a gap in the disc, and its final mass and orbital radius. The $\beta$-cooling approximation (e.g. Baruteau et al., 2011), in which the cooling rate is proportional to the local orbital period, does not properly capture the thermodynamics of the circumstellar and circumplanetary disc, and may underestimate the mass flow onto the protoplanet and its ability to open up a gap. We also found that the protoplanet accretion luminosity is significant in shaping the properties of the disc in the vicinity of the planet and may stabilize the edges of the gap in which the planet resides, altering its migration pattern.

Therefore, there are many uncertainties regarding the final fate of a young protoplanet that needs to be taken into account in population synthesis models (e.g. Mordasini et al., 2009; Forgan & Rice, 2013; Nayakshin, 2017a) in order for numerical results on planet formation at an early phase during the protostellar disc evolution to be compared with the observed properties of exoplanets. The observational evidence that only 1-10% of star host gas giant planets on wide orbits (Brandt et al., 2014; Galicher et al., 2016; Vigan et al., 2017), corresponds to the final outcome of the planet formation process which does not necessarily reflect the properties of newly formed planet-seeds in young protostellar discs.