© ESO, 2017

1. Introduction

A relative newcomer to the toolbox of exoplanet discovery and characterization techniques, direct imaging occupies an important niche with a potential for growth in the near future. Indeed, direct detections are particularly sensitive to planets at large separations (≈10–100 AU) from their host star, thus probing the outer architecture of planetary systems and complementing radial-velocity, transit, and microlensing observations. This can inform models of the migration and build-up of planets also in the inner regions (≲10 AU), with orbital and compositional consequences for instance for hot Jupiters or close-in super-Earths. Moreover, the detection of photons originating in the (non-irradiated) atmosphere of these objects affords precious information about their composition, chemistry (e.g., Lavie et al. 2017), and various microphysical processes such as dust growth and sedimentation, cloud formation, or possibly even lightning (e.g., Buenzli et al. 2015a,b; Bonnefoy et al. 2016; Bailey et al. 2014). Spectroscopy can also serve to measure the temporal variability and, at high resolution, the rotation period of these young (≲100 Myr) objects (e.g., Snellen et al. 2014; Zhou et al. 2016; Allers et al. 2016; Schwarz et al. 2016), which are connected to their atmospheric dynamics. A yet unexplored facet of exoplanet formation, the spin may also reflect the angular-momentum accretion history.

Several past and ongoing surveys have revealed a scarce but interesting population of young gas giants on wide orbits (e.g., Chauvin et al. 2010; Tamura 2016; Rameau et al. 2013b; Galicher et al. 2016). These surveys, which have various observing strategies and target different stellar masses¹, and other searches as well (e.g., Durkan et al. 2016; Bryan et al. 2016), consistently find that low-mass companions are rare. In a recent meta-analysis, Bowler (2016) derived a companion fraction of less than roughly one percent around all stars for masses from 5 to 13 M_♃ and semi-major axes of 10 to 100 or 1000 AU (see also Lannier et al. 2016). Ongoing surveys, principally the SPHERE GTO program, GPIES, LEECH, and Project 1640 (Beuzit et al. 2008; Macintosh et al. 2014; Skemer et al. 2014; Hinkley et al. 2011), as well as the upcoming SCExAO system (Jovanovic et al. 2015), should be able to increase the sample size and statistical significance of the results.

Direct imaging discoveries include the well-studied β Pic b (Lagrange et al. 2010) and HR 8799 bcde objects (Marois et al. 2008, 2010), but also a handful of exciting recent additions: HD 95086 b and 51 Eri b with small lower mass limits of a few Jupiter masses, and HIP 65426 b (Rameau et al. 2013a,c; De Rosa et al. 2016; Macintosh et al. 2015; Samland et al. 2017; Chauvin et al. 2017). As advances both in observational and data-reduction techniques make it possible to probe closer in to the host stars, one could expect an increasing number of discoveries N if the radial-velocity result that dN/ dP ≈ −0.7 < 0 (Cumming et al. 2008), derived for masses M = 0.3–10 M_♃ and periods P = 2–2000 days, holds also at separations relevant for direct imaging.

Direct observations measure an object’s flux but the fundamental quantity relevant for understanding its formation is its mass. Therefore, interpreting these observations intrinsically relies on models. In principle, one can derive mass and radius from the spectroscopy or photometry by comparing to theoretical spectra. In practice, however, there are several objects whose spectral appearance currently cannot be explained, with non-equilibrium chemistry (linked to vertical mixing) and patchy or variably thick cloud decks of ill-known composition (leading to a dusty photosphere) thought to play a role. Also, constraining the surface gravity of these objects has proven difficult since atmospheric properties are relatively insensitive to it. A more robust approach consists of using model atmospheres to derive a bolometric luminosity or, as in Morzinski et al. (2015) for β Pic b, deriving it most empirically by combining photometry of a sufficiently large part of its spectrum. However, the issue is that to derive a mass from the bolometric luminosity, the post-formation luminosity needs to be known; theoretical studies of the last decade have made clear that this is still an open issue, and in fact the single most important one in the field.

Predicting the mass-post-formation luminosity relationship of gas giants requires modeling the energetics of gas accretion. In particular, as first pointed out by Marley et al. (2007), one of the key aspects is the nature of the planetary gas accretion shock sitting on the surface of the planet during giant planet formation. The details are reviewed in Sect. 2.1, and the extreme outcomes are “cold starts” and “hot starts” (“start” referring to the beginning of the cooling, i.e., the end of the accretion). The issue is that the difference between the two, which in the first several Myr can be as large as a factor ~1000 in the luminosity (Marley et al. 2007) or up to 8 mag in the K band (Fortney et al. 2008; Spiegel & Burrows 2012), persists for the initial Kelvin-Helmholtz (cooling) timescale; for cold starts, this is often larger than the age of the observed systems. This is an important issue because, as Marleau & Cumming (2014) explore in detail, a given brightness can be either explained by a low mass with a hot start, a high mass with a cold start, or an intermediate combination. Without further information on the mass, this degeneracy between the initial luminosity and the mass derived from a measurement of the current luminosity cannot be lifted. Depending on the system, this leads to an uncertainty of a few to several M_♃. Therefore, any prediction of the post-formation luminosity of planets would represent an important step.

In only the last few years, cutting-edge observations have caught a small number of putative low-mass companions in their mass-assembly phase (around LkCa 15, HD 142527, HD 100546, and HD 169142; Kraus & Ireland 2012; Sallum et al. 2015; Close et al. 2014; Quanz et al. 2013; Reggiani et al. 2014). These fascinating observations of planet formation as it happens provide an additional motivation to study the luminosity of planets already during the formation phase, whether or not their radiation can leave relatively unimpeded their natal circumstellar and -planetary disk(s). Indeed, a part of the accretion luminosity could be directly visible in H_α (e.g., Sallum et al. 2015), or in the infrared (van Boekel et al. 2017) and the planet’s radiative feedback is expected to change the local thermal and density structure of the disk (Montesinos et al. 2015; Klahr & Kley 2006) and thus its chemistry (Cleeves et al. 2015). In the near future, this should be accessible to the unprecedented resolution and exquisite sensitivity of the Atacama Large Millimeter Array (ALMA; Wolf 2008; Cleeves et al. 2015). Future instruments like METIS at the ELT (Brandl et al. 2014) may be able to detect the vast population of low-mass planets during formation as found by comparing the detection capability of METIS with the luminosity of forming low-mass planets predicted here (van Boekel et al. 2017).

Previous theoretical studies have looked at the evolution of a planet’s luminosity, from its beginnings as a rocky core to its cooling on Gyr scales, for different masses (Marley et al. 2007; Mordasini et al. 2012c; Bodenheimer et al. 2013; Mordasini 2013). However, this was only for often artificial simplifications such as a prescribed gas accretion rate and a fixed semi-major axis, and for a limited set of conditions (e.g., a given stellar mass, nebula temperature, and planetesimals surface density), while Mordasini (2013) showed that the post-formation luminosity depends critically on these assumptions, in particular on the planetesimal surface density. Therefore, to make realistic predictions and enable a statistically meaningful comparison to observations, one needs to sample the whole range of input parameter values and to analyze the outcome statistically. This is precisely what population synthesis makes possible and the subject of this work. The resulting fundamental statistical predictions like the planetary luminosity distribution in time may be compared in future with the results of direct imaging surveys.

The contents of this paper are as follows. In Sect. 2 we shortly present our global planet formation and evolution model. We also list the parameters that were used to generate the synthetic populations. Section 3 describes the three populations studied in this work, the cold-nominal, the hot, and the cold-classical population. In Sect. 4 the simulation of the formation and evolution of a 5 M_♃ planet taken from the cold-nominal synthetic population is discussed. Many of the findings for this individual planet are useful to understand the statistical population-wide findings in Sect. 5, which is the main part of this paper. In this part we discuss three fundamental statistical properties, which is first the planetary mass-luminosity relation both during formation and evolution (Sect. 5.1), the mass-entropy diagram at the moment when the protoplanetary disk disappears (Sect. 5.2), and finally the luminosity distribution as a function of time (Sect. 5.3). We furthermore revisit the mass-radius relation that was extensively discussed in Paper II in Sect. 5.4, now including the effect of deuterium burning. In Sect. 6 we summarize our findings and give the conclusions. Appendices A and B contain fits to the post-formation properties of the synthetic planets which are of interest as initial condition for evolution models. In Appendix C we study the energy deposition into a protoplanet by planetesimal impacts which is important for the post-formation entropy of the planets.

2. Formation and evolution model

Our global planet formation and evolution model was described in detail in several earlier publications (in particular Alibert et al. 2005; Mordasini et al. 2012c, hereafter Paper I; Mordasini et al. 2012b, hereafter Paper II). Therefore we only give here a short summary. Based on the core accretion paradigm, our planet formation model combines three standard elements:

• 1.

  A classical giant planet formation model very similar to the onedescribed in Pollack et al. (1996)and Bodenheimer et al. (2000). Itcalculates the growth of the solid core by the accretion ofplanetesimals based on the Safronov equation as well as the gasaccretion rate and structure of the gaseous envelope by solving the1D internal structure equations. The envelope structure iscalculated both in the attached (or nebular) phase for planets withmasses less than 10–100 M_⊕ and in the detached (or transition) phase. The internal structure calculations in particular also yield the luminosity. In the attached phase the planet’s gaseous envelope smoothly transitions into the background nebula. The planet’s outer radius in this phase is on the order of the planet’s Hill sphere or Bondi radius, whichever is smaller. At the beginning of the detached phase, the planet’s outer radius detaches from the disk as the disk can no more deliver enough gas given the contraction of the envelope to keep the envelope and disk in contact. This contraction becomes increasingly rapid as the core grows (runaway accretion). The outer radius then rapidly decreases to a value that is much smaller than the Hill sphere (Bodenheimer et al. 2000; Paper I), and the planet now has an actual surface. In the attached phase, the Kelvin-Helmholtz contraction of the envelope sets the gas accretion rate, while in the detached phase, the accretion rate is controlled by the protoplanetary disk rather than the planet itself (disk-limited gas accretion phase, see Paper I). After a short transition phase with Bondi-limited accretion, the gas accretion rate in this phase is set equal to 0.9 the non-equilibrium gas accretion rate in the local protoplanetary disk (Paper I).

• 2.

  A standard 1+1D α model for the protoplanetary gas disk (Lynden-Bell & Pringle 1974; Papaloizou & Terquem 1999) including stellar irradiation (Fouchet et al. 2012) and internal and external photoevaporation (Paper II).

• 3.

  Orbital migration due to angular momentum exchange with the protoplanetary gas disk (Goldreich & Tremaine 1980; Lin & Papaloizou 1986), including both non-isothermal type I and type II migration (Dittkrist et al. 2014).

After the protoplanetary disk has disappeared, the planet evolution model which was described in detail in Papers I and II calculates the long-term thermodynamical evolution (cooling and contraction) up to an age of 10 Gyr. It is self-consistently linked to the formation model by directly using the internal structure of the planets at the end of the formation phase. This direct link yields in particular the post-formation entropy s_pf in the deep convective zone of a giant planet. In the evolutionary model simple gray outer boundary conditions are employed which yield however cooling tracks that are generally in good agreement with non-gray atmospheric models (Bodenheimer et al. 2000; Paper I). In the convective regions standard zero-entropy gradient convection is assumed. Deuterium burning is included (Mollière & Mordasini 2012) while atmospheric evaporation is neglected as we are interested in massive planets which are mostly unaffected by this process (Jin et al. 2014). The most important settings and parameters of the formation and evolution model can be found in Table 1.

2.1. Thermodynamics of the accretion process

The thermodynamics of the accretion process of gas and planetesimals set the post-formation entropy and therefore luminosity, radius, and magnitudes of a planet (Bodenheimer et al. 2000; Fortney et al. 2005b; Marley et al. 2007; Fortney et al. 2008; Spiegel & Burrows 2012; Mordasini et al. 2012b; Mordasini et al. 2012c; Mordasini 2013; Bodenheimer et al. 2013; Owen & Menou 2016; Berardo et al. 2017; Szulágyi & Mordasini 2017; Marleau et al. 2017).

Regarding the accretion of gas, during the disk-limited gas accretion phase (M_planet> 10–100 M_⊕) where giant planets gain most of their mass, gas falls, in the 1D spherically symmetric picture, in near-free fall from the Hill sphere onto the planet’s surface where it shocks. The recent 3D global radiation-hydrodynamic simulations of Szulágyi & Mordasini (2017) show that (polar) accretion shocks on the circumplanetary disk and protoplanet are crucial for a planet’s thermodynamic evolution also in a more realistic accretion geometry.

In the current absence of detailed radiation-hydrodynamic simulations of the planetary gas accretion shock with a realistic EOS (see Marleau et al. 2017, for 1D results with an ideal EOS), we consider the limiting cases of either completely cold accretion, where the entire gas accretion shock luminosity is radiated away in a supercritical shock, or of completely hot accretion where no radiative losses occur (subcritical shock), as described in Paper I. As for stars (Hartmann et al. 1997; Baraffe et al. 2009), such cold (hot) accretion lead to the accretion of low (high) entropy gas which in turn leads to giant planets with low (high) post-formation luminosities and small (large) radii. It is interesting to note here that for stars, radiation-hydrodynamic simulations show that the accretion shock onto the first Larson core (M ~ 200 M_♃, R ~ 10 AU, Ṁ ~ 10^-5M_⊙/yr) is supercritical (Commerçon et al. 2011; Vaytet et al. 2012), while the accretion shock of the second Larson core is subcritical and almost all the energy from the infalling material is absorbed into it (Tomida et al. 2013; Vaytet et al. 2013). The second Larson core initially shares some properties with an accreting giant planet (M ~ 1 M_♃, R ~ 6 R_♃, compare with Fig. 20), but with Ṁ ~ 10^-1M_⊙/yr its accretion rate is about seven orders of magnitude higher. The recent study of the planetary gas accretion shock in the 1D spherically symmetric approximation by Marleau et al. (2017) has shown that for an ideal EOS, the accretion shock is supercritical, but that it may still lead to a significant advection of thermal energy into the planet due to energy recycling by the infalling gas. Therefore, both limiting cases of completely hot and cold accretion are currently of interest.

Regarding the accretion of planetesimals, we assume in all except one population (see below) that planets continue to normally accrete planetesimals in the detached phase as found from the Safronov equation, and that the planetesimals always sink to the solid core (sinking approximation, Pollack et al. 1996). As demonstrated in Mordasini (2013) and Bodenheimer et al. (2013), a high planetesimal accretion luminosity associated with the formation of a massive core leads to a high post-formation entropy (and luminosity) even for a completely cold accretion of the gas due to the self-amplifying mechanism found in Mordasini (2013). For high core masses (~100 M_⊕), the post-formation luminosities even approach those in classical hot start simulations. The very low post-formation luminosities of Marley et al. (2007) were in contrast obtained under the assumption that the accretion of planetesimals artificially stops when the planet detaches, resulting in low-mass cores.

As an additional source of luminosity, deuterium burning in massive planets (Baraffe et al. 2008) is also included as described in Mollière & Mordasini (2012). The (small) radiogenic core luminosity is included as well (Paper I), but the thermal cooling of the core is currently neglected. During the formation phase, the luminosity resulting from the accretion of planetesimals is anyway usually much bigger than the core’s cooling. But it leads to an underestimation of the luminosity of low-mass, core-dominated planets with masses ≲30 M_⊕ during the evolutionary phase at constant mass (Lopez & Fortney 2014). We therefore restrict this work mostly to giant planets where the luminosity is strongly dominated by the cooling of the H/He envelope also during evolution.

In this study we work in a strictly 1D spherically symmetric approximation, meaning in particular that we do not consider the presence of a circumplanetary disk. Multi-dimensional hydrodynamic simulations have found more complex accretion geometries involving preferential accretion at high latitudes (e.g., Ayliffe & Bate 2012; Tanigawa et al. 2012; Gressel et al. 2013; Szulágyi et al. 2014; Szulágyi & Mordasini 2017). While the total accretional luminosity should remain similar in such more complex geometries for a fixed planetary mass, radius, and gas accretion rate, the characteristic surface over which it is radiated may change, meaning that the SED would be different than when assuming that the accretional luminosity originates homogeneously from the planet’s entire surface. This could have important observational consequences (e.g., Zhou et al. 2014; Eisner 2015; Zhu 2015; Szulágyi & Mordasini 2017) and must therefore be critically kept in mind.

2.2. Simplifications and limitations of the internal structure model

In this work we have modeled the interiors of the planets and their evolution under the standard assumptions of fully convective adiabatic interiors without net rearrangements of matter. In the solar system, the evolution to the present-day luminosity of 2 of the 4 giant planets (Jupiter and Neptune) can be well reproduced with such simple models (Fortney et al. 2011). For the other two, namely Saturn which is brighter relative to these models by about 60% (Fortney et al. 2011) and Uranus which is fainter by at least one order of magnitude (Guillot & Gautier 2014), additional physical effects must play a role that we do not consider here. These effects could be: a demixing of different chemical species followed by gravitational settling causing a net rearrangement of matter like a helium rain, which could explain Saturn (Stevenson & Salpeter 1977), compositional gradients (e.g., Vazan et al. 2016) that can lead to semiconvection and non-adiabatic interiors, which could explain both Saturn (Leconte & Chabrier 2013) and Uranus (Nettelmann et al. 2013), and finally core erosion where a part of the planet’s luminosity is used to dredge up core material (e.g., Guillot et al. 2004).

2.2.1. Simplifications for the EOS and statistical imprints

Furthermore, we have assumed that the envelope can be describe by the H/He equation of state of Saumon et al. (1995), i.e., we have neglected the enrichment by heavy elements. For low-mass planets, this enrichment mostly by water (for a formation outside of the iceline) can be very high (Fortney et al. 2013). For giant planets that are at the focus of this paper because of their better detectability, the planet metallicity Z should usually be rather low at least during the evolution phase with Z ~ 0.1 (Thorngren et al. 2016).

These processes, if occurring frequently in a statistical sense, imprint into the statistical quantities like the luminosity distribution or mass-luminosity relation that we study here. Our theoretical results can then serve as a baseline comparison sample the difference to which could be used as a diagnostics of the aforementioned processes. A challenge could be that in contrast to the solar system, the data on extrasolar planets is often incomplete. For example, to understand whether a luminosity differs from the prediction of simple models, a sufficiently exact knowledge of the planetary age or mass is necessary.

2.2.2. Uncertainties related to the core-mass effect

Another simplification in the model is the sinking approximation that is used for all solids that are accreted into the planets.

During the formation of a giant planet, the accretion of planetesimals can go through a (second) maximum at the moment when the planet starts to accrete gas in a runaway fashion. At this moment, envelope growth leads via the increasing Hill sphere to an extension of the feeding zone of planetesimals (Zhou & Lin 2007; Shiraishi & Ida 2008; Mordasini 2013). This process can significantly increase the total mass of solids that is accreted into the planet (Helled & Lunine 2014). It also leads, at least in the sinking approximation used here, to a significant energy input into the planet during the runaway and early disk-limited gas accretion phase. This has, as mentioned, important effects on the post-formation entropy via the core-mass effect: as found independently by Mordasini (2013) and Bodenheimer et al. (2013), the stronger pressure support in a planet with higher planetesimal accretion luminosity L_pla leads during the early disk-limited gas accretion phase to a larger radius of the planet and therefore to a lower gas accretion shock luminosity L_shock. A lower L_shock means for cold gas accretion that gas of higher entropy is incorporated into the planet which in turn leads to a slower decrease of the radius, such that the mechanism is self-amplifying (see also Berardo et al. 2017, for the importance of the entropy value at detachment).

This core-mass effect means that high luminosities are found for large core masses even for fully cold gas accretion where the entire gas accretion shock luminosity is radiated away. The high luminosities are almost comparable to those found for hot gas accretion (see Sect. 5.2). This could mean that the expected outcome for core accretion are quite bright planets, and not the very low luminosities found in Marley et al. (2007) where only low core masses ≲17 M_⊕ were considered. In the nominal population synthesis discussed below, we find core masses that are usually significantly higher (see Fig. 14). These high heavy element contents agree at least in an approximate way with the observationally inferred values (Thorngren et al. 2016) as discussed in detail in Sect. 5.1.8. So from this point of view, the core-mass effect could indeed be at work in reality.

However, there is an important caveat to this: high luminosities because of the core-mass effect require not only high heavy element contents nowadays, but also that the solids are accreted rapidly during the late attached and early detached phase, and that these solids sink quickly deep into the potential well to provide a strong accretional heating of the planet’s interior (Mordasini 2013).

This has not yet been studied with giant planet formation models tracking the thermodynamical and compositional evolution of the interior (planetesimal accretion, dissolution, sinking, mixing), representing an important open question. First works have recently started to study the compositional aspect (Venturini et al. 2016; Lozovsky et al. 2017), but they do not yet address the consequences for the luminosity. Because of this, it is currently also not clear if, for the core-mass effect to work efficiently, a high total heavy element content today (in the envelope and/or core) is sufficient as implicitly assumed in the argument above. This is because it is currently unclear whether a centrally concentrated distribution of heavy elements (in a solid core or a strongly enriched inner region) at the beginning of detachment could shift to a more homogeneous distribution over time as heavy elements are soluble in H/He under conditions typical for giant planet interiors (e.g., Soubiran & Militzer 2015), or alternatively, that heavy elements sink over long timescale to the center. The results of Vazan et al. (2016) indicate that no full mixing occurs if the initial compositional gradient is strong enough. If the heavy element mass fraction as a function of enclosed mass Z(M) in a planet today (see the Juno spacecraft data, Wahl et al. 2017) therefore reflects at least partially the structure during build-up, this would open an interesting avenue to constrain the ratio of the solid accretion to the gas accretion rate as a function of planet mass, Ṁ_Z/Ṁ_XY(M). This would obviously be of high interest in the context of the core-mass effect.

To investigate the efficiency of the core-mass effect if the solids do not sink to the core, we present in Appendix C a preliminary study based on analytical considerations using polytropic models and numerical results. We find that a homogeneous mixing of the solids into the envelope instead of sinking still provides a significant energy source, at least for the cases we studied. The reduction factors of the heating relative to the sinking case (2–3 for homogeneous mixing, 3–8 for no sinking) are not very large. This could imply that the highest luminosities caused by the core-mass effect are reduced, but that the luminosities still do not become very low as found for no impact heating at all as in Marley et al. (2007). This is the observationally most relevant question, as hot and “hotter” starts converge rapidly, whereas cold starts differ for a long time. However, it is clear that future work must address this important question more thoroughly. In that sense, the high luminosities because of the core-mass effect found here must also be confirmed by future work.

2.3. Nomenclature for the luminosities

In the following sections, several different sources of luminosity are addressed, for which we use the following nomenclature: L_HHe, the usual luminosity resulting from the cooling and contraction of the gaseous H/He envelope (material that is already part of the planet); L_radio, the radiogenic luminosity in the solid core; L_D, the deuterium burning luminosity; L_pla, the luminosity due to the accretion of planetesimals during the formation phase; and L_shock, the luminosity due to the gas accretion shock during the formation phase. With internal luminosity L_int = L_HHe + L_radio + L_D + L_pla we denote the luminosity that is generated in the planet’s interior. The total luminosity L which we address most frequently additionally includes the accretion shock luminosity that is generated at the planet’s surface, so that for cold accretion L = L_int + L_shock. It should be noted that for hot accretion L_shock = 0, as no shock luminosity escapes in this case. The accretional energy is rather deposited into the planet and radiated later on. In the evolutionary phase, i.e., after the dissipation of the protoplanetary disk, when L_shock = 0 and L_pla = 0, L = L_int.

3. Synthetic planetary populations

For the statistical study of the planetary luminosities, and in particular the global impact of cold and hot gas accretion and of massive cores, we conduct planetary population syntheses. For them, four fundamental initial conditions are varied in a Monte Carlo way: the disk metallicity, the initial disk gas mass, the external photoevaporation rate which controls together with viscous dissipation the disk lifetime, and the initial starting position of the planetary embryo. The Shakura & Sunyaev (1973) disk viscosity parameter α is fixed to 7 × 10^-3. The population synthesis method and the probability distributions from which the four Monte Carlo variables are drawn were described in Mordasini et al. (2009a). We consider the following three synthetic planetary populations:

• 1.

  The cold-nominal population CD753. It assumes completelycold gas accretion, i.e., supercritical shocks radiating allaccretional luminosity, and includes the accretion ofplanetesimals in the detached phase, assuming thatplanetesimals always sink rapidly to the core. This population isidentical to the nominal one presented in Paper IIwith the difference that deuterium burning is now included.

• 2.

  The hot population CD752. It is identical to the cold-nominal population with the only difference that completely hot gas accretion is assumed (subcritical shocks, all accretional energy is added to the planets’ interior). In Paper I it was found that under this assumption, giant planet formation by core accretion leads to high post-formation entropies that are traditionally often associated only with a formation via the direct collapse (gravitational instability) mechanism (e.g., Galvagni et al. 2012). This shows that for the post-formation entropy the structure of the accretion shock is as important as the fundamental formation mechanism. In our model the treatment of the accretion shock energetics does not influence the gas accretion rate, therefore some giant planets in the cold and hot populations can be cross-matched (see Fig. 11). Note however that the total number of planets differs in the two populations. The accretion rate of planetesimal is in contrast slightly influenced by the thermodynamics: for hot accretion, the planet’s capture radius for planetesimals remains larger during the detached phase, allowing the planets to accrete more planetesimals. The impact is however very small, with differences in the final core masses of 5% or less.

• 3.

  The cold-classical population CD777. As for the cold-nominal population, it assumes that the gas accretion is completely cold. But in contrast to the cold-nominal population, it is here assumed that planetesimal accretion stops artificially once a giant planet enters the disk-limited gas accretion (detached) phase. This causes the maximal core masses to be lower than in the cold-nominal population. Therefore the effect that massive cores make hot planets (Mordasini 2013) is of reduced importance in this population. Furthermore, orbital migration is not included in this population so that planets form in situ. Both these assumptions are the same as in the classical work of Marley et al. (2007) where planets with a very low entropy were found to form. This population was already studied in Mordasini et al. (2014), where a further description can be found.

The populations CD752 and CD753 are publicly available on DACE, the Data Analysis Centre for Exoplanets of the NCCR PlanetS reachable at https://dace.unige.ch. Additional populations will be added in future. DACE yields both interactive snapshots of the entire population at a given moment in time like the population-wide a−M or M−R diagrams and formation tracks of individual planet (e.g., M(t), R(t), etc.) for all synthetic planets.

An important idealization of the syntheses presented here is that they assume that only one embryo per disk forms (see Alibert et al. 2013, for the impact of the concurrent formation of many embryos). While this should not usually affect the thermodynamics of the accretion process itself, it can change the formation tracks in particular in systems where several giant planets form. Future work will study the impact of multiplicity on the planetary luminosities and address in particular the predictions for the number and luminosities of giant planet scattered to large orbital distances.

4. Formation and evolution of a 5 M_♃ planet

We now turn to the results obtained with the framework introduced in the previous two sections. But before we address the statistical population-wide results, we study the formation and evolution of an individual planet taken from the cold-nominal population. This detailed analysis is useful to understand the statistical results presented later in Sect. 5.

4.1. Temporal evolution of the mass and radius

The left panel of Fig. 1 shows the temporal evolution of the core, envelope, and total mass of a giant planet starting to form at 4.6 AU. As typical for the core accretion mechanism, in a first step, the core is formed. Thanks to migration, no lengthy “phase II” occurs as no isolation of the core happens (Alibert et al. 2004). Instead, due the large feeding zone and the high surface density in this disk (about 12 times the MMSN) the core continues to grow and reaches a mass of 45 M_⊕ at the crossover point which occurs at 1.96 Myr. At the crossover point, core and envelope mass are equal. In contrast to the simulations of Pollack et al. (1996), because of the large core mass, the crossover point occurs slightly after the moment when the planet enters the disk-limited gas accretion phase at 1.94 Myr. The gas accretion rate reaches a maximum of about 0.0012 M_⊕/yr at this point, about an order of magnitude smaller than the often assumed 0.01 M_⊕/yr in the disk limited phase (e.g., Hubickyj et al. 2005). The total mass of the planet is then 66 M_⊕. Afterwards, the gas accretion rate decreases roughly as 1 /t² as the disk dissipates because of accretion onto the star and photoevaporation. These rather low gas accretion rates have important implications as discussed below. As in earlier simulations, we do not reduce the planetary gas accretion rate due to gap formation because of the eccentric instability found by Kley & Dirksen (2006), but assume that the planetary accretion rate is simply 0.9 times the local gas accretion rate in the protoplanetary disk (Paper I). The disk disappears at about 8.6 Myr so that the accretion of gas and solids stops. This relatively long disk lifetime (Haisch et al. 2001; but see also Pfalzner et al. 2014) reflects the correlation of long disk lifetimes and a high probability of the formation of a massive giant planet (Mordasini et al. 2012a). The final total mass of the planet is about 1578 M_⊕ (4.97 M_♃), while the final core mass is 49.4 M_⊕.

The temporal evolution of the planetary radius is also shown in the left panel. During the attached phase, the radius is approximately the smaller of the Bondi and of one third of the Hill sphere radius (Lissauer et al. 2009) that is given for a planet of mass M and semimajor axis a around a star mass of mass M_⋆ as (1)It therefore initially increases as the planet grows in mass. But after about 1.8 Myr, it decreases as the reduction of the Hill sphere by the planet’s inward migration becomes dominant. At 1.94 Myr the planet detaches from the disk as it enters the disk limited gas accretion phase, which is visible as a change in slope in the curve. The radius then decreases rapidly (Bodenheimer et al. 2000) to reach a value of about 1.4 R_♃ at the end of the formation phase, which agrees very well with the value given by Spiegel & Burrows (2012) for a planet of this mass and entropy. This radius is intermediate between the radii predicted by Marley et al. (2007) which are about 1.25 R_♃ for a classical (very) cold start, and about 2.1 R_♃ for a hot start.

4.2. Temporal evolution of the luminosity

The right panel of Fig. 1 shows the luminosity as a function of time during both the formation and subsequent evolution phase at constant mass. The internal, shock and planetesimal accretion luminosity are shown separately. For this planet both the radiogenic and deuterium burning luminosity are of negligible importance. At the beginning, the planet’s luminosity is dominated by the accretion of the planetesimals L_pla. Since we assume that planetesimals sink to the core (2)where G is the gravitational constant, M_core the core mass, Ṁ_core the protoplanet’s planetesimal accretion rate, and R_core the radius of the solid core. As discussed in Appendix C this expression is strictly speaking only applicable when the envelope mass is negligible compared to the core mass. The planetesimal accretion luminosity reaches a maximal value of about 10^-4L_⊙, about a factor 10 higher than in Pollack et al. (1996) which is due to the higher planetesimal surface density and therefore higher Ṁ_core, and the higher M_core itself. After about 2 Myr, L_pla decreases again because of two reasons: first, it is a consequence of the decrease of the planetesimal capture radius in the detached phase (Paper I) which enters quadratically in Ṁ_core. Second, at about this moment, the planet passes into slower type II migration which also reduces the core accretion rate because the planet now migrates slower into parts of the disk with a high planetesimal surface density.

The high energy input by planetesimal accretion during the runaway and early disk-limited gas accretion phase has via the core-mass effect (see discussion in Sect. 2.2) the important consequence that despite the cold accretion of gas, the planet has a post-formation entropy and therefore luminosity that is significantly higher than in the classical (very) cold start simulation of Marley et al. (2007). In this latter paper, only low core masses of 17–19 M_⊕ were considered. In the specific case, the post-formation luminosity is approximately log (L/L_⊙) = −4.5, whereas Marley et al. (2007) had for 5 M_♃ about log (L/L_⊙) = −5.6 or about a factor 13 less, an observationally relevant difference.

Figure 1 also shows that during the disk-limited gas accretion phase where the planets gains almost all of its mass, the gas accretion shock luminosity dominates over the internal luminosity by almost one order of magnitude (factor 8.5). The accretion shock luminosity of a planet with mass M and radius R is calculated as (3)where Ṁ_XY is the planet’s gas accretion rate. For cold accretion, this dominance of L_shock is typical for giant planets in their main growth phase as we will see in the statistical analysis below (Sect. 5.1.5).

The maximum gas accretion rate in this simulation is, as mentioned, only 0.0012 M_⊕/yr at the beginning of the disk-limited phase at 1.94 Myr, followed by even lower Ṁ_XY during the subsequent 7 Myr during which the planet grows from 66 M_⊕ to its final mass. This means that the time-averaged gas accretion rate is only about 2 × 10^-4M_⊕/yr. Such a gas accretion rate is clearly lower than the value of 10^-2M_⊕/yr assumed in Marley et al. (2007). In the syntheses, it is found that such lower disk-limited gas accretion rates (a few 10^-4 to 10^-3M_⊕/yr) are actually the typical case (see Fig. 7). This has important implications for the observability of the planets in this phase: instead of short spike of very high luminosity as in Marley et al. (2007), the planets are typically less bright in our simulations, but retain this level during several million years. For a 5 M_♃ planet for example, Marley et al. (2007) find a luminosity of up to log (L/L_⊙) = −1.5 using an arbitrary gas accretion rate of 10^-2M_⊕/yr, but the high luminosity phase (log (L/L_⊙) > −4) only lasts about 3 × 10⁵ years. In our simulation where the gas accretion rate is given self-consistently by the disk model, the planet only reaches a peak luminosity of about log (L/L_⊙) = −3.15, but retains a log (L/L_⊙) > −4 during more than 5 Myr. The specific form of L_shock as a function of time (Fig. 1, right panel) is given by the interplay of the gas accretion rate and the planet mass as L_shock ∝ MṀ_XY. Initially L_shock is lower because the planet’s mass is still low. As its mass grows, also L_shock increases. But at the same time Ṁ_XY decreases in time as the disk gradually dissipates which becomes the dominant effect after some time, such that L_shock again diminishes after the maximum at about 3.1 Myr is reached. Finally, when the gas disk disappears, which marks the end of the formation phase, L_shock vanishes on a relatively short timescale, and the evolutionary phase starts. Also this temporal behavior of the luminosity during the formation phase is typical, as found again in the syntheses.

As an additional effect not seen in this simulation, for planets more massive than about 13 M_♃ (Baraffe et al. 2008; Spiegel et al. 2011), deuterium burning delays the luminosity decrease (for hot accretion) respectively re-increases it (for cold accretion) after the planets have stopped accreting during a timescale of 10⁷ to a few 10⁸ years, depending on the mass of the planet (Mollière & Mordasini 2012; Bodenheimer et al. 2013).

4.3. Post-formation entropy s_pf

A practical measure to quantify the thermodynamical state of a giant planet is the specific entropy in its deep convective zone (e.g., Marleau & Cumming 2014). In this simulation we find a post-formation entropy s_pf of about 9.1 k_B/baryon². In the classical (very) cold starts studied by Marley et al. (2007), a value of 8.5 k_B/baryon was found, which leads as seen above to a difference in the luminosity of more than one order of magnitude. The difference is due to the mentioned core-mass effect. But the value of 9.1 k_B/baryon is also clearly lower than the 10.2 k_B/baryon found in Mordasini (2013) for a planet with M = 5M_♃ and M_core = 49M_⊕, i.e., for a very similar core mass as in the simulation here. However, in Mordasini (2013), the gas accretion rate in the disk limited phase was (arbitrarily) 10^-2M_⊕/yr, while here it is more than one order of magnitude less. As shown by Spiegel & Burrows (2012), lower gas accretion rates lead to lower s_pf because the planet has more time to cool while it is still accreting. This means that the gas accretion rate not only directly affects L_shock as seen before, but indirectly, and even more importantly, also the post-formation luminosity. Spiegel & Burrows (2012) show that a reduction of Ṁ_XY by a factor 30 leads to decrease of s_pf by about 1 k_B/baryon. In our simulations, we also see this effect: if the simulation of Mordasini (2013) for M = 5M_♃ and M_core = 49M_⊕ is repeated with a disk limited gas accretion rate of 10^-3 instead of 10^-2M_⊕/yr, then a s_pf = 8.9k_B/baryon results. This is relatively close to the value found in the simulation here. An exactly identical result cannot be expected, because in the in situ simulations in Mordasini (2013), the overall temporal evolution of both Ṁ_core and Ṁ_XY differ significantly from the case simulated here. We thus see that the post-formation entropy can only be predicted reliably if all aspects of a planet’s formation tracks (temporal evolution of the gas and solid accretion rate, semimajor axis, opacity, etc.) are taken into account.

This dependency of s_pf and thus the post-formation luminosity on the individual formation track of a planet suggests that we should expect a significant spread in post-formation properties for giant planets, even for a fixed radiative efficiency of the gas accretion shock. With the population syntheses presented below, we confirm this, and can partially quantify it. In Fig. 12 we for example present the population wide mass-s_pf relation, finding that even for completely cold gas accretion, there is a wide spread of cold, warm, and even relatively hot starts. Since the radiative efficiency of the accretion shock will likely also vary depending on the planet’s properties (for example via its Mach number, Marleau et al. 2017), an even larger spread should exist in reality, and additional factors like the opacity will also play a role. This shows that statistical constraints will be among the most useful (and necessary ones) to better understand giant planet formation from direct imaging observations.

5. Statistical results

We now turn to main subject of this paper which is the statistical study of the luminosity of forming and young giant planets. We discuss the mass-luminosity relation, the statistics of the post-formation properties, the luminosity distribution, and the mass-radius relationship.

5.1. The planetary mass-luminosity relation

Figures 2 and 4 show one of the most important results of this study. It is the mass-luminosity relation of giant planets during the formation phase for cold and hot gas accretion, respectively.

While the mass-luminosity relation during the evolutionary phase at constant mass after formation has been studied in many works like the classical models of Burrows et al. (1997) or Baraffe et al. (2003), the mass-luminosity relation during formation has in contrast not been studied in a systematic way. Recently, several discoveries of (candidate) low-mass companions (low-mass stars and protoplanets) embedded in the circumstellar disk of young stars were announced: around LkCa 15 (Kraus & Ireland 2012), HD 142527 (Biller et al. 2012), HD 100546 (Quanz et al. 2013), and HD 169142 (Reggiani et al. 2014; Biller et al. 2014). While the exact nature of some of these sources needs to be further investigated, some of these objects are probably indeed observed in the act of formation as indicated by the Hα emission of LkCa 15 b Sallum et al. (2015; see also HD 142527 b, Close et al. 2014). This warrants an extension of the M−L relation to the formation phase.

5.1.1. Formation phase: cold-nominal population

Figure 2 shows the total luminosity L = L_int + L_shock as a function of mass for all giant synthetic planets in the cold-nominal population at 1, 3, 5 and 10 Myr after the start of the simulation, which can approximately be associated with the stellar age. We note the problem of a potential time delay between stellar age and planetary age, as encountered when comparing observations with purely evolutionary models not considering formation (Fortney et al. 2005a; Kraus & Ireland 2012), does not exist here, because the time needed to form the planets is self-consistently included. In the plot, planets for which the accretion shock luminosity is dominant over the internal luminosity (L_shock>L_int) are marked with a black point in the middle. The separate contributions of L_int and L_shock are discussed in Sect. 5.1.5. The colors indicate the importance of the deuterium burning luminosity. We now discuss each panel in turn.

1 Myr. In the first panel at 1 Myr, the luminosity of almost all giant planets is clearly dominated by the accretion shock luminosity as visible from the black points. Almost all giant planets are in the detached, disk-limited gas accretion phase. No protoplanetary disk has yet dissipated. At this relatively early time, the accretion rate in the disks is still high, therefore (see Sect. 5.1.5) also the accretion rates of the planets is high (up to 5 × 10^-3M_⊕/yr), and thus also L_shock. High accretion rates in the disks are even more likely because to form a giant planet on a short timescale of 1 Myr, a high disk mass is required. Within the constant α viscosity disk model, such massive disks also have high accretion rates (Lynden-Bell & Pringle 1974).

Because of the dominance of L_shock, and because R ≪ R_H (except for a short time directly after detachment) the luminosity can here be approximated simply as (4)For cold accretion, the radius of giant planets at 1 Myr is almost independent of mass between about 1 and 10 M_♃ (see Fig. 8 in Paper II and Fig. 19 below) and approximately equal to 1.8 R_♃. This facilitates estimating L. Figure 3 compares the luminosity estimated with Eq. (4)assuming for all planets a fixed radius of 1.8 R_♃ (but the accretion rate from the simulation) with the luminosity found in the simulation. One sees that the agreement is relatively good for masses between about 3 to 17 M_♃ with a typical difference of 20% or less. For less massive planets, the actual luminosity is larger than L_est, because the contribution from L_int is not negligible (see Sect. 5.1.5). For more massive planets, the actual luminosity is in contrast smaller then L_est. This is due to the increase of the radius R for these massive planets that burn deuterium in their interior as can be seen in Fig. 19. This reduces L_shock ∝ 1 /R.

When calculating L_est in Fig. 3, we have taken the Ṁ_XY of the planet from the simulation. In reality, the gas accretion rate of a planet is usually unknown. However, the accretion rate of the host star is more often known (e.g., Rigliaco et al. 2012). If the planetary and stellar accretion rates are of the same order of magnitude as suggested by hydrodynamic simulations (Lubow & D’Angelo 2006), then the independency of R on M for cold accretion can be used to estimate, via L_est, the planet’s mass provided that L_shock can be observationally obtained. This approach (see also Kraus & Ireland 2012; Close et al. 2014) then yields a constraint on the planet’s mass besides other indicators such as the gap morphology (e.g., Kanagawa et al. 2015) or the local disk temperature (Montesinos et al. 2015) and scale height (Klahr & Kley 2006). For hot accretion it is more complex to estimate the planet’s mass in this way, as the radius becomes then also a function of mass (Fig. 20).

The dominance of L_shock in this phase, the independency of R on M, and our assumption that the planet’s gas accretion rate is proportional to the disk accretion rate, i.e., a priori independent of the planet mass (no reduction due to gap formation) mean that L ≈ L_shock should roughly increase linearly with planet mass for strongly accreting planets. The slope of the envelope of points in Fig. 2 indeed increases approximately as M^1.2. The slightly larger exponent is due to a weak positive correlation of the gas accretion rate and planet mass. It is an indirect consequence of the fact that more massive planets form in more massive disks that have a higher gas accretion rate. One furthermore notes a spread of about one order of magnitude in L at a given mass.

Regarding deuterium burning, L_D/L_int is equal to 0.05, 0.1, 0.5, and 1 at 12.6, 13.2, 15.8, and 16.7 M_♃. There is thus no single mass where deuterium burning starts but a smooth transition to higher L_D contributions (Spiegel et al. 2011). Planets with masses above 16.7 M_♃ have L_D/L_int> 1, meaning that part of the deuterium burning luminosity is spent on re-inflating their radius (Mollière & Mordasini 2012).

Near the left border of the panel, one notes six planets of low mass (between about 0.15 and 0.7 M_♃) that are quite luminous with log (L/L_⊙) ≳ −2.5. A similar case is also seen in the panel at 3 Myr. The absence of black dots shows that their high luminosity is not due to the gas accretion shock. Instead, these are planets which are just in the phase where the envelope contracts rapidly (but still quasi-hydrostatically, Paper I) immediately after the detached phase starts. The high planetesimal accretion rate in this phase (Paper I,D’Angelo & Podolak 2015) plus the fast contraction powers their luminosity. Since this phase typically lasts for only a few 10⁴ yr, only a low number of planets are caught in this stage at any given moment in time. But it is an interesting phase where even low-mass planets (typical masses of 50–200 M_⊕) can reach high L for a short period.

3 Myr. In the second panel at 3 Myr, the maximal gas accretion rate has fallen to about 1 × 10^-3M_⊕/yr, and, in contrast to 1 Myr, a significant number of protoplanetary disks have already dissipated. This has the following interesting consequence: there are now two classes of planets. First those that are still accreting in the detached disk-limited regime, and second those that have already stopped accreting because the disk has dissipated. They are already in the evolutionary phase at constant mass. The two groups can be identified in the figure by the presence/absence of black dots. The planets that are still accreting (the “accreting sequence”) have higher luminosities because of L_shock, while the non-accreting planets (the “evolving sequence”) are already on standard cooling tracks, and have lower L. Between the two sequences, there is a region which is less populated. This is due to the fact that because of photoevaporation, the final dissipation of the protoplanetary disk, and with it the decrease of the planetary Ṁ_XY and thus L_shock occurs on short timescale (“two-time-scale” behavior, Clarke et al. 2001). This means that at a specific moment in time not many planets are in the final phase where L_shock falls to zero.

The two sequences are characterized by different M−L scalings: as at 1 Myr, the accreting planets have roughly speaking L ≈ L_shock ∝ M, while the non-accreting planets follow roughly speaking L = L_int ∝ M², as expected analytically (Burrows & Liebert 1993; Marleau & Cumming 2014). To illustrate these scalings, we have added in the panel two gray lines which encompass the cloud of points: for planets without significant D-burning, the envelope of points is approximately bounded by an upper limit given by those planets of the accreting sequence that accrete at the highest rate in the most massive disks. A fit to the upper limit (by eye) gives (5)This is shown by the gray dotted line in the figure. This is about a factor 5 smaller than the highest luminosities at 1 Myr because of the decrease of the maximal Ṁ_XY by the same factor, a consequence of disk evolution. The lower limit (the least luminous planet in the evolving sequence) is given by the dashed line which roughly follows (6)The lines show that at a given mass, the spread of associated total luminosities is very large, about two orders of magnitude. This is a consequence of the fact that the planets not only have different entropies depending on their formation history (see Sect. 5.2) which causes a spread in L_int (for the evolving sequence), but also that planets of a given mass can have different Ṁ_XY, which leads to a spread in L_shock. This obviously means that if only the sum of L_int and L_shock can be measured, it is very difficult to derive the mass from the observed luminosity.

Concerning D-burning, at 3 Myr the boundary for a L_D/L_int equal to 0.05, 0.1, 0.5, and 1 is at 11.4, 11.9, 13.4, and 14.2 M_♃, about 1–2 M_♃ lower than at 1 Myr. This decrease is due to the cooling of the planets which has the well-known but counterintuitive effect (Kippenhahn & Weigert 1994; Mollière & Mordasini 2012) that their central temperature increases (the interiors are not yet strongly degenerate). Planets more massive than about 20 M_♃ are now in their D-burning main sequence, defined as the phase where L_int = L_D (Mollière & Mordasini 2012). We also see that for D-burning planets, a weaker decrease of L occurs when accretion stops. There is narrower depleted region between the “accreting” and the evolving sequence than for M ≲ 10 M_♃. The reason for this is twofold: first, L_shock approaches a constant value as a function of mass for D-burning planets at a given gas accretion rate as they have an increasing radius with mass (see Fig. 20), and second, L_int is larger because of L_D.

5 Myr. The global structure at 5 Myr is quantitatively similar to the one at 3 Myr, with the difference that the accreting sequence is increasingly less populated, while the evolving sequence becomes more populated, as more and more protoplanetary disks dissipate. We also see that the luminosities have further decreased. For the accreting planets because of a decrease of the gas accretion rates (now less than 5 × 10^-4M_⊕/yr), while for the evolving sequence it is because of cooling, at least for those planets that have already stopped accreting a certain time ago. Additionally, at late times, planets accrete at a lower gas accretion rate, which allows them to (partially) cool already while they accrete (Spiegel & Burrows 2012). The mass limits for a contribution of L_D/L_int of 0.05, 0.1, 0.5, and 1 has again decreased slightly and is now at 10.8, 11.4, 12.7, and 13.5 M_♃. As indicated by the blue and green color of the most massive planets (~40 M_♃), these planets which start intense D-burning already during formation (Mollière & Mordasini 2012) have already burnt most of their deuterium because of their high central temperatures and the extreme dependency of the D-burning rate approximately ∝T^11.8 (Stahler 1988). Thus, they have already started to simply cool. The only planets for which the luminosity has not systematically decreased between 3 and 5 Myr are the those in the lower arc-like structure at about 20 M_♃: they are still in the D-burning main sequence, where the thermostatic nature of D-burning (see, e.g., Mollière & Mordasini 2012) keeps their luminosity approximately constant for some time.

10 Myr. At 10 Myr, almost all protoplanetary disk have disappeared, so that almost all planets have entered the non-accreting evolving sequence. Only a low number of planets in long-lived disks are still accreting, giving them luminosities up to 2 orders of magnitude higher than their non-accreting counterparts. The minimum mass for D-burning has again very slightly decreased, while the planets more massive than about 20 M_♃ have now already almost fully consumed their D-reservoir, so that they simply cool.

5.1.2. Comparison with classical hot and cold start models: with sinking, core accretion makes hot planets

In the panels at 3 at 5 Myr, the gray solid line shows the mass-luminosity relation in the classical hot start models of Burrows et al. (1997) at these times. The Baraffe et al. (2003) would give very similar results. These purely evolutionary models start at t = 0 with a fully formed planet of arbitrarily high specific entropy of the interior adiabat (≳10–12 k_B/baryon). The black dotted line shows the mass-luminosity relation predicted by the formation and evolution models of Marley et al. (2007), also at 3 and 5 Myr. They assumed as in our simulation here cold gas accretion and found low specific entropies of ≲8–9 k_B/baryon. As it is well known, the Marley et al. (2007) cold start models predict luminosities that are much lower then the hot start ones, by up to 2 orders of magnitudes at 10 M_♃.

Comparing our non-accreting “finished” synthetic planets with these lines shows that they are, at first sight surprisingly, almost as luminous as in the hot start simulations. We thus find that despite the completely cold accretion of gas, the luminosity of the synthetic planets in the evolving sequence is very similar to the classical hot start models of the same age. This already holds as early as at 3 Myr. The planets are thus much more luminous than originally predicted by the classical core accretion models of Marley et al. (2007). As will be further illustrated below in Sect. 5.1.8, this is a consequence of the core-mass effect (Mordasini 2013; Bodenheimer et al. 2013): especially the more massive planets (several M_♃) in our simulations have core masses clearly larger (≳50 M_⊕, see Mordasini et al. 2014) than assumed in Marley et al. (2007) who only had core masses of 17–19 M_⊕. This has a large effect on the luminosity because of the self-amplifying mechanism described in Mordasini (2013). This has the following important consequence: if the core-mass effect is efficiently acting in reality (see Sect. 2.2), then core accretion actually predicts hot start planets even if the gas accretion shock is completely cold. This has the important implication that based on a high luminosity alone, core accretion cannot be excluded as the formation mechanism.

5.1.3. Formation phase: hot gas accretion

Figure 4 shows the total luminosity L as a function of mass for the hot population, analogous to Fig. 2 for the cold-nominal population. We recall that L = L_int and L_shock = 0 in this population, either because there is no shock, or because the shock is radiatively inefficient. Rather, it is assumed that newly accreted gas is directly added to the planet’s gravothermal energy reservoir without radiative losses, assuming a radially constant luminosity as described in Paper I. It is clear that the actual accretional process could be much more complicated, involve strong luminosity gradients in the accreting envelope, and depend in particular on the relation between the planet’s initial internal entropy and the entropy of the material added by the accretion shock, as recently shown by Berardo et al. (2017). The results of this study, combined with those of works predicting the post-shock entropy (for shocks on the circumplanetary disk and then on the planet itself; see respectively Szulágyi & Mordasini 2017; Marleau et al. 2017) should allow us to obtain more realistic theoretical models for accreting young planets in the future. These can then in turn be used for future statistical studies as conducted here.

The panels show that the general trend of decreasing luminosities in time and of higher luminosities for more massive planets found for the cold-nominal population also holds for this population, as expected. However, one can also identify three important differences that are best visible in the panel at 3 Myr:

• First, there is no split into two distinct groups of planets in an“accreting” and “evolving” sequence as in the cold-nominalpopulation. This is readily understood: in the cold-nominalpopulation, at the end of the disk lifetime, the total luminosity(L = L_int + L_shock) decreases rapidly and significantly, as L_shock disappears at this moment such that only L_int is left, as visible in Fig. 1 here and Fig. 10 of Paper I. In the hot population, no such sudden drop of L occurs, first because there is no L_shock during accretion, and second because the L_int in the hot population is higher. The L_int in the cold-nominal population is also high when compared to the cold-classical population considered below because of the core-mass effect, but still less than in the hot population (see also Fig. 13). Rather, the envelope of points now covers without much substructure a triangular region in the M−L space with strongly accreting planets in the upper part, and weakly or non-accreting planets close to the lower boundary of the triangle.

• Second, one sees that when comparing the highest total luminosities at a given mass (planets that are accreting strongly), planets in the cold-nominal have a higher total luminosity (L = L_int + L_shock) than the brightest planets in the hot population of the same mass, which only have L_int. This inversion (cold-accretion planets being brighter than hot-accretion planets during formation) relative to the situation during the evolution is also visible in Fig. 10 of Paper I. This can be readily be understood from energy conservation arguments (see also Lissauer et al. 2009): during formation, the hot-accretion planets are heated by the kinetic energy influx originating from the accreting gas; this leads to larger radius and therefore to a higher gravothermal energy reservoir that can be radiated later during the evolutionary phase. Consequently, they are less luminous during formation, as the material gets accreted at a larger radius (see Sect. 5.4). For cold accretion, it is the opposite: these planets are not heated by gas accretion and accrete the gas therefore onto a smaller radius, resulting in a high L_shock that is radiated away. This makes them brighter during formation but dimmer during evolution because their gravothermal energy reservoir is smaller. For hot accretion, the same small radius is only reached at some later moment during evolution. This also means that the effective temperature of accreting planets is higher for cold accretion than for hot accretion at a given planetary mass. Note finally that during formation, if the L_int and L_shock components of the total luminosity L can be measured separately via a tracer of L_shock such as Hα emission, we should see that the L_int of cold-accretion objects tends to be lower than the L_int of hot-accretion objects.

On the other hand, when comparing the planets that have already stopped accreting and have the lowest luminosity at a given mass (at the bottom of the triangular region in the hot population and in the “evolving” sequence in the cold-nominal population), we see the opposite, namely that the planets of the hot population have a higher L_int than those of the cold-nominal population. The typical difference is however not very large, as the core-mass effect heats the planet in the cold-nominal population almost as strongly as hot gas accretion, as discussed in the previous section. The difference is best visible when comparing the envelope of points with the gray lines at 5 Myr. The difference can be much larger for planets in the cold-nominal population as will be shown below.

To illustrate these findings, we have added two empirical M−L scalings in the panel at 3 Myr for the cold-nominal population that enclose the envelope of points and form the aforementioned triangular region. The upper dotted line is given by (7)While the lower boundary given by the dashed line approximately follows: (8)We see that the upper envelope thus corresponds to a lower luminosity than in the cold-nominal population, while it is the opposite for the lower envelope. The differences are, however, not very large. But we also note that the simple power-law fits lead to a certain underestimation in the difference of L_int between the cold-nominal and the hot population for massive planets with masses between about 5–7 M_♃ and the lower limit for D-burning. Here, the actual difference in L_int at 3 Myr can actually be half an order of magnitude.

Third, there is a difference in the extent and timing of deuterium burning and also the shape of the M−L relation in the D-burning regime. Regarding the extent and timing when D-burning occurs, we see that at 1 Myr in the cold-nominal population, a significant number of planets is already intensively burning deuterium as shown by the red points, while no significant D-burning occurs in the hot population. At 3 Myr, strong D-burning is ongoing in all planets more massive than about 13 M_♃ in the cold-nominal population, whereas in the hot population only a handful more massive planets already burn significant deuterium. At even later times, strong D-burning exists in both populations, but in the hot population it is restricted to higher masses relative to the cold-nominal population. In summary we thus see that in the hot population, deuterium burning sets in later, as can also be seen in the simulations of Mollière & Mordasini (2012). This is not related to a different mass of the planets as a function of time, as the accretion of the planets is unaffected from the hot/cold assumption in our model, as mentioned above. It is rather related to their different thermodynamic state and in particular central temperature, because of the following mechanism: for cold accretion, the planets have a smaller radius at a given mass. Crudely estimating the (central) pressure P and density ρ of a planet with mass M and radius R (e.g., Cox & Giuli 1968) (9)where G is the gravitational constant, and because the planets are at this time not yet strongly degenerate so that the ideal gas law applies (Mollière & Mordasini 2012), one estimates a central temperature as a function of planetary radius (ℛ is the gas constant and μ is the mean molecular weight), (10)The smaller radius thus means that the central temperature of the cold accretion objects is higher in the first phase of D-burning when the thermostatic D-burning main sequence (where L ≈ L_D, see Mollière & Mordasini 2012) is not yet reached. As the deuterium nuclear energy generation rate is approximately ε ∝ ρT^11.8 (Stahler 1988) this means that in the cold start objects D-burning sets in earlier during formation, and uses up some of the deuterium which they then lack later on, during evolution. For hot accretion, sufficiently high central temperatures are, in contrast, only reached later on when the planetary radii have decreased, and at this moment also more deuterium is still available, so that deuterium burning occurs later and then also more intensively. The latter point is manifested by the fact that the radii of the planets during the D-burning main sequence is somewhat larger for hot accretion than for cold accretion despite the thermostatic effect of D-burning. This is found in the tracks of the planets in the syntheses and is visible in Fig. 7 of Mollière & Mordasini (2012).

Finally, regarding the shape of M−L of D-burning planets we see that for the cold-nominal population, the M−L of such planets takes an arc-like form bending upwards with a higher dL/ dM, whereas for the hot population, the M−L relation is more an approximately straight extension of the M−L relation of lower mass planets not burning deuterium with a lower dL/ dM. The reason for this is the following: in the cold-nominal population, the mass interval where D-burning occurs must connect non-deuterium burning planets with a lower L than in the hot population (compare the right end of the gray line in the panel at 10 Myr) with D-burning planets that have a similar, higher L in the two populations. This is because deuterium burning tends to remove the hot-cold difference because of its thermostatic nature, although it does not do so completely, as discussed above and in Mollière & Mordasini (2012).

5.1.4. Formation phase: cold-classical population

Figure 5 shows the mass-luminosity relation for the third population, the cold-classical population. As no accretion of planetesimals occurs in this population after detachment (as in Marley et al. 2007), giant planets can have much lower core masses than in the cold-nominal population (see Fig. 14 and Mordasini et al. 2014, for the core masses). As the gas accretion is fully cold, this has, via the core-mass effect, the consequence that some planets are born with very low luminosities. The cold-classical population thus departs from the cold-nominal population in the opposite direction than the hot population does. This has the following consequences:

• (1)

  First, generally speaking, this is the population with the lowestluminosities L_int after formation. This is best seen when comparing the envelope of non-accreting planets (no black dots) with the gray solid line in the panel at 10 Myr. During the intense accretion phase, the total luminosity L_int + L_shock is in contrast even (slightly) higher than in the cold-nominal population (and thus also in the hot population). The reason is the same as discussed in the previous section (total energy conservation).

• (2)

  In contrast to the two other populations, there is a group of planets that have (essentially) stopped accreting and do not exhibit a positive correlation between planet mass and luminosity. In the other two populations, this positive correlation holds both during formation and evolution, despite a significant scatter around the mean M−L relation. In the cold-classical population there is, in contrast, a significant group of planets with masses between about 1 to 7 M_♃ that have an L_int that is concentrated around log (L/L_⊙) ~ −6 independent of mass (see the panel at 5 and 10 Myr). This is very reminiscent of the Marley et al. (2007) result who also found such a nearly mass-independent post-formation luminosity, which is not surprising as the cold-classical population is built on the same key assumptions. The M−L relation of this work is shown in the 5 and 10 Myr panels with the black dotted line. We see that some of the synthetic planets have luminosities that are even lower than in Marley et al. (2007). Most planets have, however, not such extremely low luminosities, they rather fall into the warm start regime (see Fig. 6 below). Taken at face value, this population thus predicts that most giant planets have a warm start, while a smaller but significant group also has very low luminosities (cold start).

• (3)

  The low L_int and high L_int + L_shock, as well as the very large spread in L_int alone means that this is the population with the largest variation in the (total) luminosity at a given mass. The spread in L_int is caused by the spread in the core mass of the planets. This is shown in Sect. 5.2.3. A large spread could also result if the core-mass effect is not as efficient as assumed here (leading to lower L_int, see Sect. 2.2), but if the efficiency of the accretion shock in radiating away the accretional energy is not 100%. It will be interesting to repeat the statistical analysis once predictive models for the shock’s efficiency (e.g., Marleau et al. 2017) are coupled to planet formation and evolution models.

• (4)

  These points have the two following secondary consequences: during the main formation phase (at 3 and 5 Myr) there is even a bigger separation between the accreting and evolving sequence compared to the cold-nominal population, as there is even a bigger decrease of the total luminosity at the moment when L_shock vanishes.

• (5)

  There is an even steeper L−M relation in the deuterium burning mass domain at 10 Myr compared to the cold-nominal population, as in this mass domain an even lower L_int (compared to the cold-nominal population) for the non-burning planets must be connected with the high L ≈ L_D of the planets in the D-burning main sequence which is approximately equally high as in the cold nominal population because of the thermostatic effect.

To quantify the difference in the luminosities we show in Fig. 6 the luminosity distribution of giant planets at 10 Myr. The plot shows the L_int of giant planets with masses between 2 and 10 M_♃ for the cold-nominal, hot, and cold-classical population. The plot also indicates roughly the luminosity domains that may be associated with cold, warm, and hot starts. It is clear that there is, in principle, not a single luminosity that corresponds to the cold/warm/hot transitions as the transition at least for warm and hot starts depends on the mass (e.g., Spiegel & Burrows 2012). Therefore we have taken approximative limiting L_int values for a planet of about 5 M_♃. Given that the mass distribution is similar in the three populations (see Fig. 18), the plot nevertheless makes clear that there is an obvious shift in the luminosity distribution between the three populations that comes from the different thermodynamics during accretion. One sees the shift from the cold-classical population that contains a significant number of cold cases, to the intermediate cold-nominal population, and finally to the hot population with a number of planets that have very high luminosities.

5.1.5. Accretion shock and internal luminosities

The separate detection of Hα emission (a typical tracker of accretion) from LkCa 15 b (Sallum et al. 2015) indicates that it is possible to spectroscopically separate L_int from L_shock. We therefore plot the separate contributions to the total luminosity in Fig. 7, namely L_int (originating from the cooling and contraction of matter already in the planet, planetesimal accretion, and potentially deuterium burning) in the top left panel, and L_shock originating from the accretion of the gas in the top right. The cold-nominal population is shown at an age of 3 Myr. In the bottom left we also show an estimate of L_Hα, while the bottom right panel shows the ratio L_shock/L_int.

The H-α luminosity is calculated from L_shock via the scaling relation (Rigliaco et al. 2012) for low-mass T Tauri stars (typical M_∗ ~ 0.2 M_⊙) (11)where a = 1.49 and b = 2.99. We stress that it is currently unknown whether these scaling relations also hold for objects that have masses that are ~2 orders of magnitudes lower (see also Zhu 2015).

The colors in the top left panel showing L_int indicate that for accreting giant planets well in the detached phase (M ≳ 1 M_♃) there is a correlation between L_int and the accretion rate Ṁ_XY. For planets below the D-burning limit, this can be analyzed by writing (Hartmann et al. 1997; Paper I) (12)where (13)with E_tot = E_grav + E_int the total energy that we determine numerically from the internal structures. We see that for ξ> 1, the first term is positive and will introduce, if dominant over the two other terms, a L_int ∝ MṀ_XY. For a n = 1 polytrope that can been used to model Jupiter nowadays (e.g., Hubbard 1975), ξ = 1.5, while in the simulations, during the disk-limited gas accretion phase, we rather see ξ ≈ 1.2. Thus this term is indeed positive. But if L_int is plotted as a function of M and Ṁ_XY, in the simulations, roughly speaking, rather a is found for strongly accreting giant planets, indicating that also the other terms are important. This scaling is clearly different than the L ∝ M² for non-accreting giant planets (Burrows & Liebert 1993) and should be investigated further in dedicated work. In the panels one furthermore sees the clear imprint of D-burning, strongly increasing L_int for sufficiently massive planets with M ≳ 14 M_♃.

The top right panel shows L_shock that can be well approximated with Eq. (4)for planets between about 1 to 14 M_♃ assuming a constant radius of 1.5 R_♃ (see Fig. 19). This is illustrated by the gray lines that show Eq. (4)for a fixed R = 1.5R_♃ and Ṁ_XY = 10^-5, 10^-4, 10^-3, and 10^-2M_⊕/yr. For lower masses, L_shock is smaller than this expression because the radii are larger. The plot also shows that for the most massive planets, L_shock does not further increase. This is a consequence of their increased radii because of deuterium burning (Fig. 2).

The bottom right panel shows L_shock/L_int. This is an interesting quantity in the context of the detectability of accreting planets. One would maybe expect that from L_shock ∝ M, and L_int ∝ M² for non-accreting planets, the ration should scale as 1 /M. Instead, it increases with mass, in the 1 to 10 M_♃ domain approximately as M^0.5. This is the slope of the gray solid line in this panel. This slope can, however, again roughly speaking, be understood with the previously mentioned results because of (14)For the second proportionality we have used that in our model, the gas accretion rate in the disk-limited phase is simply a fixed fraction of the local gas accretion rate in the protoplanetary disk, which is independent of planet mass. If the mass of the growing giant planet influences its gas accretion rate in this phase, a different scaling would result. This is interesting as it may provide observational hints about the process of disk-limited accretion: if the planet remains on a circular orbit in the disk-limited phase, the gas accretion rate is a decreasing function of mass as the gap becomes increasingly wide and deep (Lubow et al. 1999). If the planet-disk system becomes instead eccentric because of the eccentric instability (Kley & Dirksen 2006), then the gas accretion rate is high also for massive planets.

At even higher masses, once deuterium burning sets in, the ratio decreases as L_int increases because of L_D, whereas L_shock decreases because of the larger radii.

These results indicate that the accretion shock luminosity of accreting giant planets with masses between about 1 to 14 M_♃ could be between a factor of a few to more than an order of magnitude higher than their internal luminosity, which would suggest that searching for accretion trackers like H-α or Paschen-β line emissions is a promising way to detect forming, detached giant planets.

However, three caveats must be added: (1) concerning the energetics of the accretion, our model is a simple one-zone model (although it does take the entire radial structure via the numerical calculation of ξ into account) with a radially constant luminosity. While numerous comparisons with other works that solve the full set of equations (like Burrows et al. 1997; Bodenheimer et al. 2000; Marley et al. 2007; Bodenheimer et al. 2013, to name just a few) showed excellent agreement (Paper I), it cannot be excluded that for some scenarios this approximation leads to different results than a more detailed treatment of the energetics of the accretion process (as in Berardo et al. 2017). This could, in turn, influence the predicted L_int and L_shock. (2) Our calculations are for a spherically symmetric accretion flow and neglect the influence of a circumplanetary disk. (3) The shock luminosity of embedded object may not escape unaltered the surrounding protoplanetary disk.

The bottom left panel shows L_Hα, estimated from L_shock with the aforementioned scaling relation. As expected from the relation, it follows a similar morphology as L_shock itself. The comparison with LkCa 15 b is found in the next section. The L_Hα of HD 142527 b (Close et al. 2014), which is a low-mass M dwarf (Lacour et al. 2016), is also shown for comparison.

5.1.6. Comparison with observed forming companions: HD 100546 b and LkCa 15 b

It is interesting to see how the luminosities found theoretically for the accreting or just formed protoplanets compare to the recent observations of the several (candidate) companions that are presumably still forming in their parent disk. In this simple explorative analysis in the context of a statistical study, we will only use the inferred luminosities and ages of the observed objects to search for synthetic counterparts. We thus do not use other available constraints like the metallicity, accretion rate, and mass of the host star (which is for some objects significantly different than the 1 M_⊙ adopted here), the properties of the protoplanetary disk (mass, gaps, SEDs, etc.), and the specific orbital, photometric, and morphological properties of the embedded companion here. Ignoring these constraints means that the results must be understood as only illustrative in nature. At the same time this long list of observational data demonstrates all the available constraints for future dedicated object-by-object comparisons, or in other words the high constraining power of observations of planet formation as it happens for theoretical models. This is especially the case as all these quantities are also available from global planet formation models as presented here, meaning that comprehensive comparisons covering many different aspects are possible.

HD 100546 b. This object was discovered by Quanz et al. (2013), recovered by Currie et al. (2014), and confirmed by Quanz et al. (2015; but see also Rameau et al. 2017). The Herbig Ae/Be host star is a young, actively accreting intermediate-mass B9Vne star (M = 2.4 ± 0.1 M_⊙, L ≈ 32 L_⊙, T_eff ≈ 10 500 K). The star has a complex transition disk, and its age is estimated to be 3–10 Myr (Brittain et al. 2014), with a poorly constrained stellar accretion rate of 5 × 10^-8M_⊙/yr (corresponding to 0.017 M_⊕/yr) according to Wright et al. (2015), while other estimates give a ~ten times lower value (Grady et al. 2005). Here we assume a fixed age of 5 Myr, but this could be extended to cover the full allowed age range in an in-depth analysis.

According to Quanz et al. (2015), the co-moving companion detected in L′ and M′ is embedded in the disk and orbits the star at a deprojected distance of about 53 AU. The observed emission consists of a point source component surrounded by spatially resolved emission. By fitting a single-temperature blackbody to the point source, an effective temperature of about T = 932 ± 200 K, a radius of 6.9 ± 2.8 R_♃ and a luminosity of log (L/L_⊙) between −3.72 and −3.54 is inferred. This range is indicated in Figs. 2 and 4 at 5 Myr. Here we associate this luminosity with L_shock + L_int. The large radius of the emitting zone potentially rather indicates emission from a circumplanetary disk, or an early post-detachment state at a lower mass than shown in Fig. 2 and 4 (see below). For an accreting circumplanetary disk in Keplerian rotation, half of L_shock would be emitted from the disk (e.g., Frank et al. 2002), and the other half from the viscous boundary layer or magnetospheric accretion. This would change the SED (Eisner 2015; Zhu 2015), but the total bolometric luminosity would remain the same.

For cold accretion, synthetic planets within the observed luminosity range have the following properties (ignoring for the moment all other constraints like T_eff): a mass range between about 0.2 and 9.7 M_♃ for the planets still actively accreting where typical values are 1–4 M_♃, and between 9.7 and 13.9 M_♃ for non-accreting planets (which seems however unlikely in the current embedded case). The planet can thus either have a high mass without (or very little) accretion, or a low (or even very low) mass down to <1 M_♃ if it is undergoing strong gas accretion. For comparison, at 5 Myr, Quanz et al. (2015) find a best-fit mass a of non-accreting, hot start planet of about 10 M_♃. The gas accretion rates of the synthetic planets with observed luminosities are between 10^-5–4 × 10^-4M_⊕/yr, with the lowest value for the most massive, still accreting planets (for higher Ṁ_XY they would have a too high L). For the lowest mass planets, the accretion of planetesimals also contributes significantly, meaning that potentially very luminous impacts may occur frequently on this planet. The radii are 1.9 R_♃ for the smallest masses (which are still in the process of significant contraction), decreasing to about 1.4 R_♃ towards higher masses ≲10 M_♃. For even more massive, non-accreting planets, the radius increases again to 1.5–1.6 R_♃ because of D-burning. The effective temperatures are between 1600–2000 K. Regarding the radius, the values of the synthetic analogs are clearly smaller, and consequently larger for T_eff than the ones estimated by Quanz et al. (2015), indicating a more extended region of emission. However, it has to be kept in mind that the effective temperature of the HD 100546 b is derived from fitting a blackbody spectrum to 4 photometric points, one of which is only an upper limit. The spectral shapes of cool objects may deviate from a blackbody, due to the presence of absorbing molecules like water or methane. This could affect the estimated total luminosity and effective temperature. On the other hand, a significant dust envelope could weaken the effects of molecules because dust opacities are gray absorbers in comparison to the molecules. Clearly, observations at a higher spectral resolution would be very helpful.

The ratio L_shock/L_int of the subset of synthetic planets with the same luminosity as HD 100546 b is about 0.1 at the lowest mass, then increases to about 15 at 3 M_♃, and then decreases again towards more massive but still accreting planets, qualitatively similar to Fig. 7. The accretion luminosity is thus clearly dominant for intermediate mass planets, meaning there is enough luminosity that could be liberated over an extended region around the planet, instead of directly on the planet’s surface. This could be a circumplanetary disk or an otherwise heated zone where the hard accretion shock luminosity gets absorbed and then re-emitted at longer wavelength. L_D/L_int is about 0.15 for the most massive, but still accreting planets, and completely negligible for the more typical lower masses. For the non-accreting planets, L_D/L_int ranges between 10^-3 and 0.8. The planets have entropies between 9 and 11 k_B/baryon, and core masses between 20 and 400 M_⊕, with typical values of around 70 M_⊕. One important constraint can clearly not be met by the synthetic planets: the semimajor axes are just 0.2–16 AU, much less than the 53 AU in the observation. Here, mechanisms currently not included in the formation model must be invoked, which could be scattering (Alibert et al. 2013) or a faster core growth for example via pebbles (Ormel & Klahr 2010; Bitsch et al. 2015).

Due to the core-mass effect, even cold gas accretion leads to quite warm planets in our model, therefore considering hot instead of cold gas accretion (Fig. 4) does not strongly change the properties of the synthetic planets with the same luminosities as HD 100546 b. Still, the following numbers change somewhat: the masses are now between 0.2 and 8.5 M_♃, i.e., the highest mass is smaller because of the still higher entropy. The accreting planets have masses between about 0.2 and 6.5 M_♃, mainly uniformly distributed between 2 and 5 M_♃. The radii of the accreting planets are larger than for cold accretion, usually between 1.6 and 2 R_♃, but still much smaller than derived from the observation. Because of the lower masses, deuterium burning is of negligible importance, while the entropies are between 10 and 11 k_B/baryon.

Up to this point we have mainly concentrated on synthetic planets that have the same luminosity as observed, assuming that the observed luminosity originates at least partially from the accretion shock. We have seen that such objects have radii smaller than observed, making it necessary to invoke elements not present in the simple 1D framework, like the presence of a circumplanetary disk or a heated gas blob around the planet (Klahr & Kley 2006), which is certainly possible.

An alternative hypothesis within the 1D picture that takes, besides the luminosity, also the effective temperature into account is that the object is in an earlier evolutionary state. As mentioned above (Sect. 5.1.1), protoplanets can also be very bright for a short period during the fast contraction phase just after detachment. Indeed, in the cold-nominal population one finds for example a synthetic planet that has at 4.96 Myr a log (L/L_⊙) = −3.54, a radius of about 6.5 R_♃ and an effective temperature of 915 K, all in agreement with the observed values. This planet has just detached 30 000 year earlier, and is still in the rapid contraction phase which is together with planetesimal impacts the main source of its luminosity. The gas accretion shock luminosity ∝M/R becomes important only later because of the still extended state of the planet and its low mass. The planet has a mass of just 78 M_⊕ and accretes gas at a rate of about 5 × 10^-4M_⊕/yr which would be about a third of the lower estimates for the stellar accretion rate. While with 5 AU the semi-major axis of the synthetic planet is much less than the one of the actual object such that this is not a real analog of HD 100546 b, it indicates that such an earlier stage can explain at least some of the observations. The critical point is here that this phase only lasts a few 10⁴ years, while the phase were the accretion shock yields the luminosity lasts 10⁵–10⁶ years. We would thus observe HD 100546 b in a special moment, which is not very likely (but neither impossible). Because of its low mass the planet would likely not have yet opened a clear gap in the disk, which is in agreement with the absence of a gap in the polarized light images of Avenhaus et al. (2014), and it would also not be a significant source of hard radiation like Hα. Similar objects can also be found in the hot population. These considerations also point towards the fascinating possibility to directly observe the fundamental phases of giant planet formation predicted by the theoretical models (attached, detached, during the rapid contraction, etc.). This will be further addressed in dedicated work.

In summary we see that because of accretion, a wide range of masses is compatible with the observed luminosity. In particular planets with low, potentially sub-Jovian masses could generate the observed luminosity, either due to accretional luminosity, or due to rapid contraction of the gaseous envelope.

LkCa 15 b. LkCa 15 is a young (2–5 Myr) solar analog (L ≈ 0.74 L_⊙, M = 1.01 ± 0.03 M_⊙) in Taurus-Auriga (e.g., Isella et al. 2014; Piétu et al. 2007) that hosts a gas-rich circumstellar disk (Andrews & Williams 2005; Andrews et al. 2011). Despite a large cavity in the dust extending to about 45 AU, the star is accreting at rate of about 1.3 × 10^-9M_⊙/yr as inferred from the stellar UV excess (Hartmann et al. 1998).

Kraus & Ireland (2012) have reported the discovery of an unresolved source in K′ and a surrounding extended emission in L′ that they interpret as a forming protoplanet surrounded by heated circumplanetary material like a circumplanetary disk. If this candidate protoplanet is coplanar with the disk and on a low-eccentricity orbit, a semimajor axis of 16 AU results. They also estimated a total bolometric luminosity of ~2 × 10^-3L_⊙, and that the unresolved emission could be reproduced naively with the photospheric emission of a non-accreting giant planet of 6 to 10 M_♃, but stress that the luminosity could rather result from accretion. Sallum et al. (2015) detect two sources, one of which is detected both in Ks and L′ and in Hα. They explain the infrared emission by a circumplanetary disk, and the Hα by emission from a hot accretion shock with L_Hα ≈ 6 × 10^-5L_⊙ that they convert via T Tauri scaling relations (Rigliaco et al. 2012) into an accretional luminosity of ~4 × 10^-4L_⊙.

The bottom left panel of Fig. 7 compares the inferred L_Hα of LkCa 15 b with the synthetic population at 3 Myr (a comparison at 2 Myr yields similar results). It suggests that planets with masses between 1 and several tens of Jupiter masses (i.e., a wide range) can have such accretion luminosities if the luminosity is due to the planetary gas accretion shock and the Rigliaco et al. (2012) relation extends to these masses, with low masses corresponding to high gas accretion rates and vice versa. However, if it is additionally requested that the planetary gas accretion rate (color coded in Fig. 7) is as least as high as the observed stellar accretion rate as suggested by hydrodynamic simulation (Lubow & D’Angelo 2006), then only a narrow mass interval of giant planets with masses between 1 to 2 M_♃ concurrently fulfills these two constraints. These synthetic planets have a ratio L_shock/L_int of 3–10, radii of about 1.6–2 R_♃, and escape velocities between 30 to 70 km s^-1. Their effective temperature to radiate L_int is about 1100 to 1800 K, while the effective temperature to radiate both L_int + L_shock is between 1900 and 2400 K if the accretion shock covers the entire surface of the planet (strictly spherical accretion).

In the context of magnetospheric accretion onto CTTS one rather observes covering factors f of the accretion columns on the stellar surface that are much less than unity, f ~ 0.1–30% (e.g., Calvet & Gullbring 1998; Bouvier et al. 2007; Ingleby et al. 2013). If this picture also applies to planets, for f = 1%, these planets would have shock-heated spots on the photosphere with temperatures of about 6000–8000 K, in good agreement with the estimates of Zhu (2015). We can finally use Eq. (22) of his work (which assumes that the temperatures in a planet’s magnetosphere – if there is one – is 8000 K as for CTTS) to very roughly estimate L_Hα. One finds L_Hα/L_⊙ ~ 1–2 × 10^-5 × (R_T/R_♃)² where R_T is the magnetospheric cavity truncation radius, which depends on the unknown magnetic field strength of the planet. For R_T ~ R_♃ this is at least to order of magnitude consistent with the observationally determined L_Hα ≈ 6 × 10^-5L_⊙ that stood at the beginning of this analysis. We conclude that giant planets with rather low masses of 1–2 M_♃ could be responsible for the observed L_Hα, and that observational determinations of L_int, as well as of the magnetic field strength of young planets would be important to further understand the process of planetary gas accretion. We also add the caveat that our synthetic analogs have semimajor axes between 1–5 AU, clearly less than LkCa 15 b, such that they cannot meet also this observational constraint.

5.1.7. Evolution phase: a spread in the M – L relation at young ages

As evident from the last section, it is difficult to constrain the mass from a measured total luminosity (L_int + L_shock) alone during formation because of the possible contribution of accretion. A comprehensive view with separate (spectroscopic ) determination of the different contributions (planetary internal, planetary accretion shock, and circumplanetary disk luminosity) combined with other indicators (like accretion rate or disk structure) should – if possible – be used to get a clearer picture.

But also during the evolution phase, the diversity of planetary formation histories can still influence a planet’s thermodynamical state, meaning that there is no one-to-one mapping from one luminosity to one mass. This has the implication that even if a young planet’s luminosity could be measured perfectly without any error bar, one could still only infer its mass from it within a certain interval. This is well known from the cold vs. hot start dichotomy. But even if all planets have a quasi-hot start as it is the case in the cold-nominal population, there is an intrinsic spread in the M−L relation at early ages. We show this in Fig. 8 displaying the M−L relation of the cold-nominal population at four moments in time during the evolutionary phase at constant mass. We discuss in the following sections the main reason for the scatter, which is the scatter in the core mass, and exemplify the implications of the non-unique M−L relation by applying it to β Pic b and 51 Eri b. The even much bigger spread introduced by core masses that vary much more is discussed in Sect. 5.1.11 which described the M−L relation of the cold-classical and hot population.

At 20 Myr, the spread in the M−L relation leads to a L that can vary by 50% or more at a given mass. The panel for 50 Myr shows that at this time, the cold-nominal population already follows well the general hot start M−L relation indicated by the gray line, but still it is not unique (see Fig. 9). This panel, together with the one at 0.5 Gyr, also illustrates the effects of deuterium burning: planets undergoing D-burning show up as a clear bump in the M−L relation that in time moves to lower masses and becomes less prominent. The green symbols at 5 Gyr demonstrate that low-level D-burning (5% of the total L) continues to this age, but only with a very minor imprint on the M−L. One can also observe a “rejuvenating” effect of D-burning. In the 50 Myr panel, the planets below the D-burning limit and those where the D-burning is already over (blue symbols) both follow the identical L ∝ M² scaling. But the planets that previously underwent D-burning have an offset to higher L, as if they would be younger objects.

One furthermore notes the thickening and flattening of the M−L relation at log (L/L_⊙) ≈ −6.5 that is visible at 50 Myr, 0.5, and 5 Gyr (where it coincides with the D-burning planets). This is in contrast to the aforementioned spread not an echo of the formation process, but is caused by a change of the internal structure of the planets (a detached radiative zone). We discuss this feature, which is additionally distance (irradiation) dependent, and its impact on the planetary luminosity distribution in Sect. 5.3.

5.1.8. Impact of M_core for the cold-nominal population

The left panel of Fig. 9 again shows the internal luminosity (excluding L_shock for the handful planets still accreting at 10 Myr) as a function of mass at 10, 20, and 50 Myr for the cold-nominal population. The plot shows three things:

• (1)

  First, as the plot is zoomed in into a smaller mass andluminosity range than Fig. 8, it allows to betterquantify the intrinsic luminosity spread at a given mass and age.We see that at an age of 20 Myr which is well into theevolutionary phase, there is still a variation in luminosity by afactor of about 1.5 to nearly 2 at a fixed mass. At anage of 50 Myr, the spread is as expected reduced asthe hotter planets cool faster leading to convergence, but even thenthe intrinsic spread is still between 20 to 50%. This should becritically kept in mind when deriving masses from measuredluminosities even at ages of ~50 Myr.

• (2)

  Second, it shows via the color code the planets’ (high) core masses. These high core masses are the reason for the high luminosities even for cold gas accretion, via the aforementioned core-mass effect (Mordasini 2013; Bodenheimer et al. 2013). The core masses in our simulations are, in particular, usually clearly higher than in Marley et al. (2007). While these authors had core masses of about 18 M_⊕, all synthetic planets with M> 0.3 M_♃ in this population have core masses higher than that except for three planets. For 1 M_♃ planets, e.g., the typical heavy element masses are around 40–60 M_⊕. Because of the high importance of the core mass for the post-formation luminosity, it is important to see how the heavy element masses in the synthetic planets compare to observational constraints. In the solar system, the highest total heavy element content of Jupiter allowed by internal structure models is about 42 M_⊕ for adiabatic models and the SCvH EOS (Guillot & Gautier 2014), while for semi-convective models up to 63 M_⊕ are possible (Leconte & Chabrier 2012). Concerning exoplanets, thanks to recent analyses (Miller & Fortney 2011; Thorngren et al. 2016) of the mass-radius relation of warm Jupiters without significant bloating, it is now possible to infer observationally estimates of the heavy element mass contained in planets also outside the solar system. While these analyses are (as those for Jupiter) affected by many uncertainties, e.g., regarding the equation of state and cannot yield (in contrast to the solar system) information about the distribution of the heavy elements within the planet (in the core or mixed throughout the envelope) they still significantly increase the number of planets with heavy element estimates and allow to make statistical inferences. Accepting these caveats, and the fact that these analyses apply to planets inside of 1 AU, while in the context of direct imaging one is mainly interested in planets further out, then the aforementioned analyses also indicate that rather high heavy element contents seem to be quite common in giant exoplanets. Specifically, for a sample of 47 planets with masses between about 0.1 and 10 M_♃, Thorngren et al. (2016) derive a mean heavy element mass M_Z as a function of total mass M of M_Z/M_⊕ ≈ (58 ± 7)(M/M_♃)^{0.61 ± 0.08}. A least-squares fit to the synthetic planets in the cold nominal population with 0.1 ≤ M/M_♃ ≤ 10 yields, for comparison, M_Z/M_⊕ ≈ 45(M/M_♃)^0.24 if only planets with a semimajor axis of less than 1 AU are included as in the sample of Thorngren et al. (2016), and M_Z/M_⊕ ≈ 63(M/M_♃)^0.26 if all orbital distances are included. This shows that the amount of planetesimals accreted in our simulations seems to be for approximately Jovian mass planets in rough agreement with or, for high planetary masses, even lower than the amount indicated by currently available observational constraints. This in turn could indicate that the high luminosities we find in the cold-nominal population (almost comparable to those found for hot accretion) are the outcome for core accretion, and not the very low luminosities in Marley et al. (2007) or in the low-core mass population where the accretion of planetesimals is artificially shut off. However, there is an important caveat to this which is the currently poorly known efficiency of the core-mass effect (see Sect. 2.2).

• (3)

  Besides the generally high core masses and associated high luminosities, the plot also directly shows the main reason for the intrinsic spread of L at fixed total mass and age. As is visible from colors of the symbols at a given total mass and age, the planets with a higher core mass usually also have a higher luminosity. The correlation between core mass and luminosity found by Mordasini (2013) and Bodenheimer et al. (2013) for idealized conditions (like a fixed semi-major axis, non-evolving disk, externally prescribed gas accretion rate and final planet mass) is thus also recovered in the population syntheses where these quantities are obtained self-consistently from the disk and migration model, covering also a much larger parameter space. The plot, however, also shows that there are planets that have a rather low core mass, but still a relatively high luminosity, and vice versa. This illustrates that while the spread in core masses is the main cause for the spread in luminosity, it is not the only one. Also the rate of gas accretion or the moment when most of the accretion happened relative to the disk lifetime influence the luminosity. This means that while there is a general trend of increasing post-formation L with total and core mass, the exact post-formation properties of a planet can only be understood when the complete formation history is considered. In Mordasini (2013) it was found that the post-formation luminosity L_pf increases with core mass at fixed total mass roughly like . It is interesting to explore whether a similar scaling also exists in the syntheses by studying L_pf(M_core) in a narrow bin of total mass (say of 1 M_♃ width). It is found that there is indeed such a scaling, maybe like , but the spread around it is large (see Fig. 14).

5.1.9. A correlation of stellar [Fe/H] and planetary L?

In view of these results, it is interesting whether a link can be made between the stellar [Fe/H] which is an observable quantity, and the planetary luminosity which can likely also be observed. The reason is that the observed mass-radius relation of extrasolar hot and warm Jupiters indicates that there is a positive correlation between stellar [Fe/H] and planetary heavy element content (Guillot et al. 2006; Burrows et al. 2007), even if this correlation appears to be less clear in a bigger more recent data set (Thorngren et al. 2016). This could lead, via the core-mass effect, to a stellar [Fe/H] − planetary L correlation. Is there such a correlation in the synthetic population?

Therefore, in the right panel of Fig. 9 the luminosity at 20 Myr of synthetic planets more massive than 3 M_♃ is shown as a function of host star’s [Fe/H] around which they formed. The stellar [Fe/H] is one of the Monte Carlo initial conditions that is varied in the population syntheses, where it is assumed that the stellar [Fe/H] and disk [M/H] (dust to gas ratio) are proportional to each other (see Mordasini et al. 2009a).

In the plot, several features can be seen:

• (1)

  First, there is an approximately horizontal over-density ofplanets with a luminosity of about 10^-3L_⊙. These are planets with masses between about 15 and 25 M_♃ which are in the process of burning deuterium or which have burnt it recently, and that are now cooling off from it. The tracks of Mollière & Mordasini (2012) show that at this age (both for cold and hot accretion), the luminosity is approximately independent of mass in the ~15–25 M_♃ range. The reason is the following: within this mass range, more massive objects reach higher luminosities because of burning deuterium (about 10^-2L_⊙), but burn most of it earlier (at ≲10⁷ years), so that at 20 Myr they have again cooled down by about one order of magnitude in L. Lower mass objects in contrast reach less high luminosities because of D-burning (about 10^-3L_⊙), but they do it later on, so that at 20 Myr, they are still strongly burning D or have not yet cooled down so much. The combination of the intensity of D-burning and the moment when it happens thus leads to an approximately mass-independent L. The luminosity value of the horizontal bar (and the mass of the objects in it) decreases in time which is visible in the distribution of luminosities (Sect. 5.3).

• (2)

  Second, we see that there is a paucity of very luminous planets (corresponding to masses ≳10 M_♃) at low metallicities. In general, in core accretion, for more common giant planets with lower masses, the final mass of a planet is not correlated with the stellar [Fe/H] (Mordasini et al. 2012a). Rather, a higher [Fe/H] leads to a higher fraction of disks that can form a giant, in agreement with the observed metallicity effect (Santos et al. 2001; Fischer & Valenti 2005). This effect is also clearly visible in the panel by the high density of points at high [Fe/H]. In other words, [Fe/H] controls whether a giant planet can form, but usually not its mass, which is given by the disk’s gas mass (Mordasini et al. 2012a). Only for the highest masses, it is different: to grow to a very high mass, a critical core must form before disk evolutions has had time to significantly deplete the disk mass. Such an early start is impossible at low metallicity where growing a critical core take longer.

• (3)

  Third, the colors of the points show the positive correlation of disk/stellar [Fe/H] and the heavy element content that was mentioned for observations earlier (Guillot et al. 2006; Burrows et al. 2007). As discussed already earlier (Mordasini et al. 2009b), this is also recovered in the synthetic population. The plot however also shows that the correlation only holds in a statistical sense, as there are also planets at high [Fe/H] without a particularly high heavy element content.

• (4)

  However, we see that apart of the absence of the highest L at low [Fe/H], there is no correlation of [Fe/H] and L. The reason is that the most important quantity determining the luminosity is the total mass, which is, as mentioned, not directly dependent on [Fe/H] for most planets. In this panel, and in contrast to the left panel, the total mass of the planets is not directly visible. But it is more fundamental for L than the effect of M_core at given total mass. Only if the population is divided in fine bins in total mass of about 0.5 M_♃ in width, a positive correlation between L and [Fe/H] becomes visible, but even then, the correlation is weak with a lot of scatter. This is a consequence of the other three Monte Carlo variables which are the disk gas mass, disk photoevaporation rate (corresponding to different disk lifetimes), and the starting position of the embryo. They introduce so much variation in the formation tracks that there is only a weak imprint.

So, in summary we see that the answer to the question posed at the beginning of the section, whether there is a correlation of [Fe/H] and L, is no, except for an absence of very luminous (and massive) planets at low [Fe/H]. Regarding them, one should kept in mind that such planets are in the syntheses actually a rare outcome: only about 0.6% of the embryos reach a mass higher than 13 M_♃. For comparison Marcy & Butler (2000) have estimated a frequency of companions more massive than 13 M_♃ within 3 AU of ≲0.5%. In the plots shown here they only appear numerous because we are dealing with a very large synthetic population of about 50 000 planets in total.

5.1.10. Comparison with β Pic b and 51 Eri b

In the top left panel of Fig. 8, the gray lines and shaded regions indicate synthetic planets that have a luminosity compatible with measurements of two important directly imaged extrasolar planets at 20 Myr, namely β Pictoris b (Lagrange et al. 2009, 2010) and 51 Eridani b (Macintosh et al. 2015). Both stars are members of the β Pic moving group with an age estimate of 21 ± 4 Myr (Binks & Jeffries 2014), while Macintosh et al. (2015) adopt 20 ± 6 Myr for 51 Eri. It is interesting to study the properties of the synthetic analogs of the two planets, but it should be kept in mind that we only selected them based on a compatible luminosity. Additional constraints such as the planetary semimajor axis, the stellar mass (here 1 M_⊙, in reality about 1.75 M_⊙ for both stars) or the stellar [Fe/H] were not considered as constraints. This can, however, be included in specific population syntheses, as demonstrated for β Pic (but for older observational constraints) in Bonnefoy et al. (2013), Mordasini et al. (2015).

βPic b: from recent spectroscopic and photometric observations, Bonnefoy et al. (2014) derive a revised luminosity of log (L/L_⊙) = −3.90 ± 0.07, a T_eff of 1650 ± 150 K and a log g ≤ 4.7 dex. Synthetic planets compatible with the luminosity have the following properties: a mass of 10.4–13.3 M_♃, which is compatible with, but as expected a somewhat wider range than earlier estimates of about 10.7–12.3 M_♃ from hot start evolutionary models (Bonnefoy et al. 2014). As visible from the figure, even if the luminosity could be determined observationally without any error bar, an uncertainty of about 2 M_♃ would remain due to the aforementioned intrinsic spread in the M−L relation at this early age. The synthetic planets have a T_eff of 1600–1700 K, and a log g of 4.18–4.26 dex which is also compatible with observations. Their radius is 1.28–1.38 R_♃, in agreement with the observationally inferred values of 1.2–1.7 R_♃ (Bonnefoy et al. 2014, 2013; Currie et al. 2013). The planets have a current specific entropy of the interior adiabat of 9.35–9.55 k_B/baryon³, which is a quite narrow range, and a central temperature of 3.7–4.2 × 10⁵ K. Mixing length theory predicts that there are about 21 layers of eddies transporting the heat in the convective zone of the planet, assuming a mixing length parameter of α_ml = 1. As can be seen from the colors of the β Pic b analogs in Fig. 8, this planet is at the very interesting transition zone between planets that do not burn deuterium and deuterium-burning planets. At the current epoch, the relative contribution of D-burning to the planet’s luminosity is L_D/L_int of 0.2 to 0.4 for the lower masses ≤12.5 M_♃, and up to 0.8 if the mass is higher (also depending on the core mass). Up to now, the planet has however only consumed about 3 to 10% of its initial deuterium, assumed to initially have a standard interstellar medium abundance of D:H of 2 × 10^-5 (Spiegel et al. 2011). During the main sequence lifetime of β Pic of about 2 Gyr, the planet will typically burn about 10% of its deuterium if its mass is on the lower side of the allowed range, and up to 90% if its mass is rather 13 M_♃. In view of atmospheric spectra, it is interesting to consider the heavy element content in the planets. It is found that the synthetic analogs have typically accreted about 30 to 100 M_⊕ of solids, but a handful have up to 350 M_⊕ of metals. This gives them a wide range of heavy element mass fractions between 0.005 and 0.1, with a typical value of about 0.015. This corresponds to an enrichment relative to solar of about 0.3 to 7, with a typical value of around 1.

51 Eridani b: according to Macintosh et al. (2015), this planet has a luminosity of log (L/L_⊙) of −5.4 to −5.8, which corresponds to synthetic planets with masses between 1.7 and 3.6 M_♃. This is in agreement with, but covers again (due to the intrinsic scatter), a somewhat wider range than the hot-start masses reported by Macintosh et al. (2015), 1.8–2.9 M_♃. The plot also shows that even if the luminosity were perfectly known, the intrinsic spread in the M−L relation would lead to a uncertainty of almost 1 M_♃ in the mass. The synthetic planets have an effective temperature of 570–740 K (630 K at log (L/L_⊙) = −5.6), similar to the values derived from the observations (600–750 K), and most radii are between 1.12 to 1.25 R_♃, but a handful have lower R, down to about 0.85 R_♃. This is related to the heavy element content of the planets. Most synthetic planets contain between 30–100 M_⊕ of metals, but a handful have accreted up to 500 M_⊕ of solids. Checking the metallicity of the initial conditions that lead to these strongly enriched synthetic planets shows that they only form in clearly supersolar metallicity disks, in contrast to the approximately solar metallicity of 51 Eri A (Macintosh et al. 2015). Thus these are not self-consistent analogs, meaning that for an in-depth analysis all observational constrains should be taken into account. The typical bulk heavy element mass fraction then ranges from 0.04 to 0.2, corresponding to about 3 to 16 times solar. If 51 Eri b’s atmosphere reflects this bulk composition, it is thus predicted to have supersolar abundances (which new results indicate, see Samland et al. 2017). As discussed in Mordasini et al. (2014), Mordasini et al. (2016), a higher metal content relative to the star is a typical prediction of the core accretion paradigm for a rather low-mass giant planet like 51 Eri b but not necessarily for a massive planet like β Pic b. The specific entropy in the planet’s convective zone is about 8.5–8.7 k_B/baryon. From mixing length theory it is estimated that about 16–18 layers of eddies swirl in the planet’s convective zone. Deuterium burning is of course completely negligible, with log (L_D/L_int) ≈ −11.

5.1.11. Evolutionary phase: cold vs. hot accretion

Figure 10 shows the mass-luminosity diagram of the cold-classical population. Compared to the cold-nominal and hot population, the cold-classical population contains a characteristic “bulge” of planets with masses between about 1 to 7 M_♃ that have very low luminosities and “hang” at early times below the M−L relation of the hot and cold-nominal population. Its temporal evolution can be understood considering Marley et al. (2007): because of the long Kelvin-Helmholtz timescale of these very cold planets, this “bulge” remains almost static in time, until the planets with warm and hot starts have themselves cooled down to these low luminosity values. At this moment the evolution converges and the memory of the initial conditions is lost. Put differently, planets roughly stay at their initial entropy until they have reached an age equal to their initial cooling time, at which point they join the hot-track sequence (see Eq. (13) of Marleau & Cumming 2014). The timescale on which this happens increases with planetary mass.

Figure 11 compares the luminosities of the giant planets in the three populations during evolution. In the panel at 10 Myr, it is possible to directly cross-match some accreting planets in the hot and cold-nominal population. This is not possible for the planets in the cold-nominal population as no migration is included in this population, leading to different growth tracks. This cross-matching shows the effect mentioned above that for strongly accreting planets, the cold accretion planet are brighter, at least if the total luminosity L_int + L_shock is considered. The plot also shows the offset in the timing of deuterium burning.

In the panel at 20 Myr, the luminosities compatible with those of 51 Eri b and β Pic b and the corresponding mass intervals are again shown, analogous to Fig. 8 where they were compared with the cold-nominal population only. We see that the situation is quite different for the two planets.

For 51 Eri b the allowed mass interval is much bigger than when compared with the cold-nominal population only where a mass between 1.7 and 3.6 M_♃ was found. Compatible planets from the cold-classical population have in contrast masses between about 1.7 and almost 10 M_♃. Macintosh et al. (2015) estimated for comparison cold-start masses of 2 to 12 M_♃. The plot furthermore shows that the planets with core masses below 17 M_⊕ (black dots) cannot (except for two synthetic planets) reproduce the observed luminosity. This means that the lower limit on the heavy-element content is anti-correlated with the (true) total mass.

The compatible mass interval for β Pictoris b is in contrast only marginally increased relative to comparison with the cold-nominal population only, namely by less than 1 M_♃ to higher masses. This is a consequence of D-burning that tends to erase the hot vs. cold accretion imprint (Sect. 5.1.3).

At 50 Myr, the plot shows that the hot and cold-nominal populations have converged to virtually identical luminosities, while for the cold-classical this takes clearly longer. At 500 Myr, some of the very cold planets just below the deuterium burning limit still have somewhat smaller luminosities. This panels shows that interestingly, some D-burning planets in this population are more luminous than in the hot or cold-nominal population. This can be understood as a mild form of a late D-flash (Salpeter 1992; Marleau & Cumming 2014).

We finally also comment that for low-mass planets, there is no hot vs. cold accretion difference in the same sense as for giant planets. The more massive a giant planet, the higher the fraction of its final mass it accretes after detachment (at ~100 M_⊕) through an increasingly strong, potentially entropy reducing shock. Detachment occurs when the planetary gas accretion rate as controlled by the planet’s Kelvin-Helmholtz contraction rate becomes higher than the rate at which the protoplanetary disk can deliver gas to the planet. As low-mass planets have very long KH-timescales and thus low gas accretion rates (e.g., Ikoma et al. 2000), they only detach from the nebula when the nebula itself has already almost completely dissipated. This means that only a tiny amount of gas is accreted after detachment through a potentially entropy lowering shock. In this sense, these planets have a hot start.

However, the post-formation entropies of the low-mass planets (M ≲ 1 M_♃) can still vary significantly depending on the formation history because of other mechanisms. As an example, a planet that underwent strong solid accretion not long before the disk disappeared will have a higher entropy at the start of the evolutionary phase than a planet that accreted strongly at the beginning of the disk lifetime, and then already cooled for several million years. The entropy and luminosity of young low-mass planets and their observational consequences will be quantified in future work (Linder et al., in prep.).

5.2. Planetary properties at the moment of disk dissipation

In the previous sections, we have analyzed the luminosity of the populations at given uniform ages. An age that is of special interest but different for each planet (or protoplanetary disk) is the moment when the protoplanetary disk disappears. The properties of the planets at this moment are interesting as they set the initial conditions for the evolutionary phase, and are thus important for planet evolution models. Additionally, the properties of the planets at this moment have been previously predicted for specific initial conditions by planet formation models like Marley et al. (2007), Mordasini (2013), Bodenheimer et al. (2013).

5.2.1. Post-formation entropies s_pf

In particular the planetary mass–specific entropy relation at the end of the formation phase is of interest, as the specific entropy of the planet’s interior adiabat is the ideal quantity to describe the gravothermal heat content of a planet, in the sense that the mass and the specific entropy of the planet fully describe the interior structure of the planet, at least for a given composition, and under the assumption of a fully convective planet (e.g., Marleau & Cumming 2014). In the original work of Marley et al. (2007) it was shown that if the specific post-formation entropy s_pf is plotted as a function mass for cold and hot start planets (their Fig. 2), then the curves take the form of a tuning fork, with an upper prong of high entropies that increase with mass for the hot start planets, and a lower prong of low entropies that decrease with mass for the cold start planets. Spiegel & Burrows (2012) then showed that intermediate radiative efficiencies of the accretion shock lead to additional prongs of intermediate entropies (warm starts). They also found that lower gas accretion rates lead to lower s_pf. Mordasini (2013) then showed that such intermediate warm states also result from a fully radiatively efficient shock (cold accretion) but a core mass that is higher than in Marley et al. (2007). This mechanism was also found by Bodenheimer et al. (2013).

These results indicate that given the diversity of protoplanetary disks, one expects a significant spread of planetary post-formation entropies (and thus luminosities and radii), and not only two extreme values like in Marley et al. (2007). To quantify this, we plot the specific entropy of synthetic planets more massive than 0.3 M_♃ at the moment when their protoplanetary disk disappears in Fig. 12.

In the population synthesis calculations (see Mordasini et al. 2009a), the distribution of external photoevaporation rates is chosen in such a way that the synthetic disks have a lifetime distribution that is approximately in agreement with observations (e.g., Haisch et al. 2001; but see also Pfalzner et al. 2014). Consequently, the mean lifetime of all synthetic disks is 3.05 Myr. The mean lifetime of disks which form a giant planet (M> 0.3 M_♃) is, however, longer, 5.1 Myr. This is expected, since giant planet formation by core accretion is facilitated in long-lived disks, as there is more time to build the critical core and to accrete gas (Mordasini et al. 2012a). Thus in Fig. 12 the age of the stars will on average be 5.1 Myr, but it can be as short as 1.6 Myr for the most short-lived disks producing a giant planet, and 16 Myr for the longest one.

Considering first the hot (red stars) and the cold-classical mass population (black points), we see the same general picture as found in earlier work. In the hot population, the entropy increases with mass, demonstrating again (Paper I) that core accretion with a radiatively inefficient shock leads to high entropies comparable to classical hot start simulations (but see Paper I, and Berardo et al. 2017, for a possible issue when simulating hot accretion with the ∂L/∂r = 0 simplification employed here). This is traditionally rather associated with gravitational instability (Galvagni et al. 2012), but we show here that this is not necessarily the case. For the cold-classical population, and for core masses less than 17 M_⊕ which mimics the simulations of Marley et al. (2007), s_pf in contrast decreases with mass. Thus, the fundamental shape of the classical “tuning fork” is recovered. However, and this is a new aspect, due to the different formation histories, there is a large spread of about 1 to 1.5 k_B/baryon in s_pf at a given mass in the hot population. In the cold-classical population, the s_pf even have a ~2 k_B/baryon scatter at about 5 M_♃. One also sees that there is a significant overlap in the covered entropy range between the populations, especially between the hot and the cold-nominal populations, but also between the cold-nominal and the cold-classical populations.

If finally also the cold-nominal population (blue empty circles) is taken into account, the complete covering of a wide entropy range is even more apparent. This population is characterized by a mean entropy that is approximately constant⁴ at about 9.2 k_B/baryon for masses ≲5 M_♃, followed by an increase of s_pf at higher masses. Also in this population there is an intrinsic scatter in s_pf at fixed mass of about 1 k_B/baryon at 1 M_♃ and 1.5 k_B/baryon at 10 M_♃. Such a spread in entropy is associated with important luminosity differences.

The mass range in Fig. 12 extends to higher masses compared to earlier studies and now includes also D-burning planets. One sees that D-burning leads in the cold population at about 12–13 M_♃ to a sharp and well-defined upturn in the entropy, as expected (Mollière & Mordasini 2012). Also in the hot population the entropy starts to increase stronger with mass at this point, but the transition is smoother. For the most massive planets (≳20 M_♃) the entropy is again smaller, as these planets have already burnt most their deuterium when the disk disappears (Sect. 5.1.1).

It is interesting to compare the hot and cold start entropies in Fig. 12 with the “tuning fork” in Marley et al. (2007). One finds that the black points are at even lower entropies, by 0.5–0.7 k_B/baryon depending on the mass. Also the red points are lower than the upper line in Marley et al. (2007), but these entropies were arbitrarily set in Marley et al. (2007) and not obtained from a formation calculation as here, so that this does not pose an issue. The reason for the difference between the lower prong in Marley et al. (2007) and the black points is not completely obvious. It is likely a consequence of several different settings in the two models. First, in Marley et al. (2007), a grain opacity reduced to 2% of the ISM opacity was used, while here it is reduced to 0.3%. As demonstrated by Spiegel & Burrows (2012) and Mordasini (2013), a lower opacity leads to lower s_pf because the planets can cool more efficiency already during formation. Second, on a related note, the atmospheric boundary conditions used here differ from those used in Marley et al. (2007). Marleau & Cumming (2014), who also use an Eddington atmosphere, found that the entropy in their models needed to be increased by 0.14 k_B/baryon⁵ to match a given luminosity in Marley et al. (2007). Finally, in the synthesis the typical gas accretion rates are (as also seen in the example of Sect. 4) up to two orders smaller than in Marley et al. (2007) where the Ṁ_XY is (arbitrarily) constant at a high 10^-2M_⊕/yr. This is quite different from the behavior here where the disk evolution leads to a gradual reduction of Ṁ_XY down to very low values. As found by Spiegel & Burrows (2012) and Mordasini (2013), lower gas accretion rates during formation lead to lower s_pf.

In Appendix A we present the post-formation entropy s_pf as a function of H/He envelope mass M_env for the cold and hot population, covering a very wide range in envelope masses. Such a relation is of interest especially as initial conditions for evolutionary models. We therefore provide a least-squares fit to the numerical results in Appendix A.

5.2.2. Post-formation luminosities L_pf

A more observational quantification of the thermodynamic state of new-born giant planets is the post-formation luminosity. In analogy to Fig. 12, Fig. 13 shows the luminosity of giant planets as a function of mass at the moment when the disk disappears. Again, the cold-nominal, hot, and cold-classical populations are shown. The general shape of the M−L_pf is a consequence of the M−s_pf relation studied before and the way the luminosity depends on mass and entropy, L(M,S), as studied in detail in Marleau & Cumming (2014). In the plot, one can identify four different interesting M−L_pf relations.

• (1)

  First, the hottest planets. The planets with the highest post-formation luminosity are found in the hot population. Theupper envelope of points approximately follows the relationL_pf/L_⊙ = 7 × 10^-5(M/M_♃)^1.4 shown by the short-dashed curve in the figure with a somewhat lower values at intermediate masses.

• (2)

  Second, the cold-nominal planets. A least-squares fit to the cold-nominal planets (blue circles) gives L_pf/L_⊙ = 1.2 × 10^-5(M/M_♃)^1.3. This is the dashed-dotted curve in the figure. The (slightly) lower exponent compared to the first case reflects the fact that planets with masses higher than about 3 M_♃ and up to the D-burning limit have still somewhat lower luminosities because of cold gas accretion which an increasing difference in luminosities with increasing mass, even though the impact of cold vs. hot gas accretion is strongly reduced because of the core-mass effect.

• (3)

  Third, cold-classical population planets with an intermediate luminosity, shown by the gray symbols. They can be associated with an intermediate warm start, and the dotted line in the plot shows L_pf/L_⊙ = 4 × 10^-6(M/M_♃)^0.5. In this population there is the largest spread in the post-formation luminosity (a factor 200 at 5–10 M_♃), with L_pf correlating strongly with the core mass. One furthermore sees that once deuterium burning kicks in at sufficiently high masses, it reduces the spread of the post-formation thermodynamic conditions, with the variation in L_pf for planets more massive than about 15 M_♃ reduced to about a factor ≤5.

• (4)

  Finally, the coldest planets. The coldest planets in the cold-classical population with M_core ≤ 17M_⊕ have a very low luminosity around L_pf = 4 × 10^-7L_⊙ (as shown by the long-dashed line in the figure). The luminosity is approximately independent of mass, with a slight trend of a decreasing luminosity with increasing mass. For comparison, Marley et al. (2007) had also found a nearly mass-independent post-formation luminosity, at a typical value of about 2–3 × 10^-6L_⊙.

These groups represent physically motivated initial values for the post-formation luminosities of giant planets obtained from a much higher ensemble of possible formation tracks compared to earlier works. Cooling curves and magnitudes for these different initial conditions employing non-gray atmospheric boundary conditions are in preparation (Marleau et al., in prep.).

5.2.3. Cold-classical population: impact of M_core on L_pf

The influence of the core mass for the cold-nominal population was discussed above in Sect. 5.1.8. However, if planetesimal accretion is assumed to continue freely after detachment, as is the case in the cold-nominal population, the core masses are effectively high enough to result in relatively high post-formation luminosities, which in part increase with increasing mass, as shown by the blue dots in Fig. 13. The population where the impact of the core mass is even much larger is the cold-classical population. In this population, giant planets can have a very wide range of core masses from less than 10 M_⊕ to more than 100 M_⊕. This leads to a very wide spread of post-formation entropies and luminosity, as shown by the black and gray points in Figs. 12 and 13. To directly illustrate the effect of the core mass on L_pf in this population, we again show in Fig. 14 the luminosity at the moment when the protoplanetary disk disappears L_pf as a function of mass, but now color-coding the planets’ core mass. The plot shows first the very clear and strong correlation between M_core and L_pf at a given total mass. It leads to a very large spread in L_pf at one mass, with the spread increasing with total mass for 0.1≲ M/M_♃ ≲ 7. At the mass showing the biggest variation (about 7 M_♃) L_pf covers almost three orders of magnitude. It is a consequence of the spread in M_core that varies at this mass by more than a factor 20.

At even higher masses than about 7 M_♃, the spread in L_pf becomes smaller again as all these very massive planets also have massive cores, and for even higher masses deuterium burning further reduces the spread in L_pf, bringing it down to just about a factor 2 at M ≈ 17 M_♃. For the most massive planets in this population (M ≈ 30 M_♃) L_pf varies again by about a factor 10 because of different timings of the earlier D-burning.

The plot furthermore illustrates the correlation between core mass M_core and total mass M that was extensively studied in Mordasini et al. (2014). There it was shown that for planets forming in situ and accreting all planetesimals in the feeding zone, one finds analytically that M_core ∝ M^1/3, and furthermore that M_core ∝ Σ_P where Σ_P is the surface density of planetesimals. In the population synthesis, Σ_P is given as the product of the (initial) disk gas mass and the dust-to-gas ratio, which are both Monte Carlo variables (Mordasini et al. 2009a). Their distributions cover about one (for the dust-to-gas ratio) or two orders of magnitude (for the initial disk gas mass), such that the large spread in M_core and thus L_pf is a direct consequence of the varying initial conditions in protoplanetary disks. We however stress again (Sect. 2.2) that in reality the core-mass effect might be not as efficient as assumed in the present model, in which accreted planetesimals instantaneously sink to the central core (Mordasini 2013). This would then naturally also reduce its impact on L_pf. This should be quantified with future work treating self-consistently the thermodynamics and compositional evolution during formation (e.g., Venturini et al. 2016).

5.3. The planetary luminosity distribution

In the previous sections, we have mostly studied scatter plots involving the luminosity which give an insight into the diversity of it and its dependencies on the mass. However, also the simple distribution of L is of interest, e.g., for statistical comparisons of the theoretically predicted luminosity distribution with observational data as soon as there will be a higher number of directly imaged planets, or to estimate the necessary sensitivity to enter a certain discovery parameter space. Compared to the mass-luminosity relation, the advantage here is that the mass is not necessary for a comparison. The interest in the luminosity distribution is analogous to the situation of radial velocity surveys in the past which allowed to compare the observed planetary mass distribution (Howard et al. 2010; Mayor et al. 2011) with theoretical predictions from population syntheses (Mordasini et al. 2009b; Benz et al. 2014).

5.3.1. Luminosity distribution as a function of time for the cold-nominal population

In Fig. 15 we show the distribution of the luminosities in the cold-nominal population at six moments in time. The histogram extends down to a luminosity of about 0.1 present time Jovian luminosities. As expected, the general trend is that the distribution shifts to lower L as time goes on because of cooling. At a fixed moment in time, the L distribution mainly reflects the mass distribution, and how the luminosity depends on it. The distribution of log (M) in the synthetic populations is characterized by a relatively flat distribution in the giant planet regime (see Fig. 18 and Mordasini et al. 2009b), corresponding to a distribution scaling approximately with M^-1 in linear units, in good agreement with observations of extrasolar giant planets inside of 6 AU (e.g., Marcy et al. 2005). Below 10–40 M_⊕, there is a sharp upturn of the mass function, reflecting the transition from giant planets that underwent gas runaway accretion to the much more numerous lower mass planets that did not do so. Such a transition is also seen in the observational data (Howard et al. 2010; Mayor et al. 2011). However, there are also a number of features that are specific to the luminosity distribution, especially at young ages. We now discuss each age in turn.

2 Myr. In the middle of the formation phase, one sees an approximately flat part (in log) of high luminosities, and strong upturn starting at log (L/L_⊙) ≈ −6. This upturn divides lower mass subcritical planets (masses of less than a few 10 M_⊕) that are in the attached phase from forming giant planets that are detached. The luminosity of these lower mass planets is mainly powered by the accretion of planetesimals L_pla, whereas the luminosity of the forming giant planets originates predominately from the gas accretion, as visible from the green dashed line showing L_shock. To order of magnitude, one can estimate the luminosity resulting from the accretion of a core of mass M_core, radius R_core with a rate Ṁ_core as (15)where we have assumed rocky cores following the mass-radius relation of Valencia et al. (2006) and a constant accretion timescale of t_form ≡ M_core/Ṁ_core for the right hand side. In reality, the accretion rate depends on the surface density of planetesimals, the distance from the star, etc., and varies in time (see Sect. 4). The equation nevertheless illustrates that for the accretion of solid cores of 1–10 M_⊕ on a timescale of a few Myr, one expects luminosities between log (L/L_⊙) ≈ −8 and −6. At even lower luminosities, there is also an important contribution by low-mass planets (M ≲ 5 M_⊕) that are not significantly accreting planetesimals. Their luminosity is mainly given by the cooling and contraction of the low-mass gaseous envelope.

Both high precision RV surveys (e.g., Mayor et al. 2011) and the Kepler satellite (e.g., Fressin et al. 2013) have shown that Neptunian and super-Earth planets are very abundant around solar like stars, a result recovered in the syntheses. These planets have luminosities ≲10^-6L_⊙ during formation. If it would be possible to detect and disentangle them form other disk features despite the fact that these planets are still embedded in the disk (attached phase) and cold, one would expect to find them in most protoplanetary disks. Their indirect impact on the disk, e.g., heating it up locally, which increases the vertical scale height (Klahr & Kley 2006), and which changes the disk’s local molecular abundances (e.g., in HCN) could possibly be detectable with ALMA (Cleeves et al. 2015).

5 Myr. At this age, the general shape is similar to the one at 2 Myr, except for a general reduction of the luminosities that are related to accretion. The peak of the luminosity generated from gas accretion has fallen by about 1 order of magnitude to now around log (L_shock/L_⊙) ≈ −3, while the upturn related to the subcritical planets now rather starts at log (L/L_⊙) ≲ −7.

10, 50 Myr. In these panels, almost all planets have entered the evolutionary phase, meaning that the underlying mass distribution is constant. At 10 Myr the green line indicates that in some long-lived disk, giant planet formation is still ongoing, but this affects only a few planets. Thus the luminosity distribution mainly reflects the mass distribution via the mass-luminosity relation for objects not undergoing accretion. We thus see an upturn at about log (L/L_⊙) ≈ −8.5 and −9.5 at 10 and 50 Myr, respectively, corresponding to luminosities of planets with masses of a few 10 M_⊕ shortly after formation.

At higher luminosities (e.g., log (L/L_⊙) ≈ −8 to −3.5 at 10 Myr), the luminosity distribution is to first order flat⁶, modulated by a local maximum at about log (L/L_⊙) ≈ −6.5. This range contains the “normal” lower mass giant planets that do not burn deuterium, and is therefore of key importance from an observational point of view. The flat distribution and the local maximum can be understood in the following way:

• (a)

  Regarding the flat distribution. The planetary mass function inthe giant planet regime from about 0.3 to 10 M_♃ scales as mentioned approximately as 1 /M followed by a faster decrease at even higher masses, both in the observations (Marcy et al. 2005; Cumming et al. 2008) and the synthetic population studied here (Fig. 18, see also Mordasini et al. 2009b). If the luminosity is a power law in M, like in particular L ∝ M² (Burrows & Liebert 1993) during evolution, then the probability distribution of L will also be proportional to 1 /L, or uniform in log (L) from the fundamental transformation law of probabilities. This explains why in the giant planet regime, the luminosity distributions are to first order flat over almost four orders magnitude. As the L ∝ M² relation is independent of time, this flat distribution can be seen in several panels of Figs. 15 and 17 below. This is an important results result of this study, and has important implications for example for the expected yield of direct imaging searches, especially once they start to probe closer-in giant planets where radial velocity surveys have found the aforementioned mass distribution.

• (b)

  Regarding the maximum at about log (L/L_⊙) ≈ −6.5. This feature is seen at several ages in both Figs. 15 and 17 and is not related to the underlying mass distribution. It rather related to a feature in the luminosity evolution as a function of time of giant planets which is in turn controlled by the microphysics (opacity) as we will see now. It can be understood from considering Fig. 16, which shows the luminosity as a function of time for giant planets of 1 to 10 M_♃ during the evolutionary phase. The plot, which is a zoomed in version of Fig. 10 of Paper I, shows the cooling curves of the model used in this paper as solid lines. In the figure one sees that at a luminosity of about log (L/L_⊙) ≈ −6.5, there is a phase where the luminosity decreases less rapidly as a function of time than otherwise. This feature is very similar to the one in Figs. 3 and 8 of Marleau & Cumming (2014) who use very similar atmospheric boundary conditions as in this work (gray atmosphere with Freedman et al. 2008, opacities). Since the value of L where this flattening of the L(t) curves occurs is independent of mass (but it occurs at different moments in time for different masses), it has the population-wide consequence that at a give time (like 50 Myr), there is a flattening of the mass-luminosity relation (a lower dL/ dM). This is, e.g., visible in Fig. 8 where there is a flattening of the M−L relation at log (L/L_⊙) ≈ −6.5 at 50 Myr, 0.5, and 5 Gyr. This flattening means that compared to otherwise, a larger mass interval maps approximately into the same luminosity interval, which means for the luminosity histogram that a larger number of planets falls into the same luminosity bin of log (L/L_⊙) ≈ −6.5 (given the underlying flat mass distribution). From this consideration we can also see that at 50 and 500 Myr, the maximum in N(L) is formed by planets of masses around 1 and 5 M_♃, respectively.

The underlying physical reason for the change in slope of L(t) is a change in the atmospheric structure and thus the efficiency of cooling, caused by the existence of a detached radiative zone (Fortney et al. 2011). A related, but much deeper detached radiative zone was originally proposed to exist in Jupiter to this day (Guillot et al. 1994a,b). The later addition of the opacity provided by alkali metals, that are important opacity sources where other elements (hydrogen, helium, water, methane and ammonia) have opacity windows, has since much reduced the possible p−T domain where radiative transport is possible (Guillot et al. 2004). Figure 16 indeed shows that for a Jupiter-mass planet, the detached radiative zone already disappears at around 20 Myr, when the opacity provided by the alkali metals is included as it is the case in the Freedman et al. (2008) opacities we use. More generally, it is found that when the effective temperature of a planet has fallen to about 400 K (and the temperature at the radiative-convective boundary is ~1500 K), the detached radiative zone in the planet disappears, and the radiative convective boundary jumps from about 10 bar to 0.3 bar, much closer to the photosphere. This is visible in the temporal evolution of the p−T diagram of Jupiter shown in Fig. 5 of Paper I or for other masses in Fig. 1 of Marleau & Cumming (2014). As discussed by Marleau & Cumming (2014), the disappearance of the detached convective zone is due to the change of the behavior of the Rosseland opacity that changes from decreasing along the p−T profile to increasing with it once the effective temperature has fallen below about 400 K. This is in turn due to the appearance of CH₄ which is crucial for the opacity in this regime. Interestingly, we thus see an imprint of the microphysics into the planetary luminosity distribution.

However, Fig. 16 also shows that for the atmospheric models of Burrows et al. (1997) and Baraffe et al. (2003), there is no comparable bump in the L(t) at around log (L/L_⊙) ≈ −6.5. These models use both non-gray atmospheres and different underlying line lists for the opacities. The absence of the bump is also seen in the more modern hybrid models of Saumon & Marley (2008). These models include the transition from cloudy L dwarfs to cloudless T dwarfs at a T_eff of 1400–1200 K⁷. If the evolution of the population would be calculated with these boundary conditions, there would not be a maximum in the luminosity distribution at log (L/L_⊙) ≈ −6.5. The recent models of Fortney et al. (2011) on the other hand contain a small bump in L(t) at a comparable, but somewhat lower luminosity of log (L/L_⊙) ≈ −7. Interestingly enough, these non-gray models also find a detached radiative zone that disappears at T_eff ≲ 400 K, or about 30 Myr for Jupiter. This feature thus depends on the specific abundances and opacity tables. We can thus speculate that the diversity of atmospheric compositions, which is identical (solar composition) in all synthetic planets in an artificial way, could lead to a smearing out, or even disappearance of the peak. We therefore stress that the key result for the luminosity distribution of giant planets is that it is, except for the “D-peak” discussed below, roughly flat in log (L) (and not the potential maximum at log (L/L_⊙) ≈ −6.5).

It should be noted that the maximum at log (L/L_⊙) ≈ −6.5 would similarly not be present if one were to consider only closer-in planets in the synthetic population where irradiation prevents the temperature in the atmosphere to fall below about 400 K, i.e., where the equilibrium temperature that is given by the stellar irradiation is larger than this value. The temporal evolution of the atmospheric structure proceeds differently in this case, in particular no detached radiative zone forms (see Fig. 8 in Mordasini et al. 2016). Rather, there is just one exterior radiative zone that deepens in time. Therefore, there is also no change of the slope of L(t) as seen for the less irradiated planets when the detached radiative zone disappears. This in turn leads to very smooth L(t) curves similar to the ones seen in the Burrows et al. (1997) and Baraffe et al. (2003) models. This is exemplified by the luminosity as a function of time after formation of the planet shown in Fig. 1, which has a final semimajor axis of 0.24 AU, and thus an equilibrium temperature of about 570 K. This distance dependence explains why the M−L relation (e.g., in Fig. 8) not only flattens at log (L/L_⊙) ≈ −6.5, but also widens.

Coming back to the luminosity distribution shown in Fig. 15, at luminosities larger than log (L/L_⊙) ≈ −5 at 50 Myr, there is a decrease in the luminosity distribution with increasing L, reflecting the further decrease of the number of planets with increasing mass. Finally at even higher luminosities (log (L/L_⊙) ≳ −3), we see another feature that is not visible in the mass distribution. As discussed in Sect. 5.1.9, for massive companions, D-burning causes the luminosity to be similar for a relatively wide mass range, as more massive planets burn D earlier and at a higher L, while less massive ones burn it later and at a lower luminosity, leading to overlapping L(t) tracks (see the tracks in Mollière & Mordasini 2012). This manifests in the L distribution with a “D-peak”, i.e., a narrow local maximum in the luminosity. At 10 Myr, the peak is at log (L/L_⊙) ≈ −2.5 formed by planets with masses of about 18 to 35 M_♃ (Fig. 2), while at 50 Myr, it is at log (L/L_⊙) ≈ −3.8, containing planets with masses between about 12 and 25 M_♃ (Fig. 8). There is a paucity of luminosities just below these values, as the usual planetary cooling is delayed by the D burning.

0.5, 5 Gyr. In this panels the luminosity of the non-giant planets has fallen below the lowest plotted value, so that the increase of the luminosity distribution at the transition from solid to gas-dominated planets is not visible anymore. The flat part in the distribution formed by giant planets extends from about log (L/L_⊙) ≈ −10 to −6 at 0.5 Gyr. The “D-peak” is still well visible at 0.5 Gyr at log (L/L_⊙) ≈ −5.5 consisting of planets with masses between about 11 and 19 M_♃, whereas at 5 Gyr, the contribution of D-burning has decreased so much that the “D-peak”, which would at this moment be at log (L/L_⊙) ≈ −7 (see Fig. 8), is no more clearly visible. The luminosity distribution therefore simply reflects the giant planet mass distribution and the upper end of the planetary mass function.

5.3.2. Luminosity distribution: comparison of the cold-nominal, hot, and cold-classical populations

After the detailed study of the distribution of L in the cold-nominal population as a function of time, we show in Fig. 17 a comparison of the luminosity distributions of the cold-nominal, hot, and cold-classical populations at four moments in time during the evolutionary phase. It is important to note that the cold-nominal and hot populations have nearly identical mass distributions (as mentioned above there is no impact of hot/cold accretion for the mass accretion rate included in the formation model). Also the cold-classical population has a very similar mass function in the giant planet range, except for a somewhat lower number of very massive planets with M ≳ 10 M_♃. This is shown in Fig. 18 which compares the mass distributions in the giant planet domain after accretion has essentially stopped. These similar mass functions mean that the differences in the luminosity distributions between the three populations are only due to difference in the thermodynamic state of the planets.

As expected from Fig. 11, there is only a small difference between the luminosity distributions of the cold-nominal and hot population, namely a somewhat higher number of planets with high log (L/L_⊙) ≈ −4 in the hot population at 10 Myr. Otherwise, one recognizes the same features as seen in the distribution of the cold-nominal population (Fig. 15). The “D-peak” is, e.g., particularly well visible in both cases. The cold-classical population has, in contrast, a very different luminosity distribution: its distribution contains a strong local maximum at about log (L/L_⊙) ≈ −6 at 10 Myr as all giant planets with masses between about 1 and 10 M_♃ and core masses of less than 10–20 M_⊕ in this population approximately have a mass-independent luminosity after formation (Marley et al. 2007, Fig. 11). Luminosities exceeding this value correspond to planets with higher core masses (Fig. 14). This strong local maximum remains visible in the distribution up to an age of at least 100 Myr, additionally amplified by the bump in L(t) at log (L/L_⊙) ≈ −6.5 discussed in the previous section. This very different luminosity distribution with a peak can clearly serve as a statistical diagnostic of the thermodynamics of the formation process once a higher number of directly imaged exoplanets is known.

5.4. Mass–radius relation with deuterium burning for cold and hot accretion

The planetary mass-radius relation for cold accretion was extensively studied in Paper II, but these calculation neglected deuterium burning. After a (giant) planet has finished its rapid contraction occurring after it detaches from the disk, its radius only changes by a factor of a few (about 3–5) during the remainder of its lifetime, while the luminosity still changes by many orders of magnitude. This weak dependency of the radius is a consequence of the partially degenerate interiors of the giant planets that are characterized by an EOS that is only weakly temperature dependent. The M−R relation nevertheless reflects the thermodynamic state of a planet and depends on the cold/hot accretion mode and on the occurrence of D-burning.

5.4.1. Cold gas accretion (cold-nominal population)

In Figure 19 we therefore show the impact of deuterium burning on the temporal evolution of the M−R relation of the cold-nominal population. The plot shows the mass-radius relation during formation at 3 Myr and at three times during the evolutionary phase, at 20, 50, and 500 Myr. Green, yellow, and red points correspond to planets with L_D/L_int of at least 0.05, 0.1, and 0.5, as in Fig. 2. Small black dots additionally show planets with L_D/L_int> 1.

The plot illustrates the population-wide impact of a main finding of Mollière & Mordasini (2012), namely that for cold-accretion, D-burning does not merely delay the contraction of the radius as it is the case in classical hot start simulations (e.g., Burrows et al. 1997), but leads to a re-inflation of the radius. This is illustrated by the example of one M−R track of a deuterium burning planet shown by the brown line in the top left panel. After the rapid decrease of the radius occurring after the detachment from the protoplanetary starting here at about 0.3 M_♃, the planet⁸ then grows from 1 to 17 M_♃ while the radius decreases from about 3.5 to 1.7 R_♃, a typical sign of cold gas accretion. At this moment, the deuterium burning becomes so strong that it starts to drive an expansion of the outer radius. In this phase, the luminosity radiated at the surface is lower than the internally produced deuterium luminosity. This excess is used to increase the planets radius and thus gravothermal energy, which is then again radiated at later times. The planet’s radius reaches a maximum of almost 3 R_♃, after which is evolves vertically down in the M−R plane. We thus see that for cold accretion, D burning has a clear and quite severe impact on the radii during formation, with an increasing radius with mass at 3 Myr for masses higher than about 15 M_♃.

The top left panel of Fig. 19 furthermore shows that for cold accretion, the typical radii of fully collapsed giant planets not undergoing D-burning (masses between about 1 and 13 M_♃) in the flat part of the M−R relation cover at 3 Myr a quite narrow range of about 1.4 to 1.6 R_♃. At 1 Myr, the situation is qualitatively similar, but the typical radius is now rather 1.8 R_♃, and the spread around this value is larger. At 8 Myr, the typical radius is 1.3–1.4 R_♃. These nearly mass-independent radii are interesting because the first facilitate estimating the accretion shock luminosity (Sect. 5.1.1) and second because this is different from hot accretion, as we will see below. This could be interesting from an observational point of view if future spectroscopic observations of accreting planets can determine the surface gravity observationally.

In the panel at 20 Myr, gas accretion has stopped. The specific shape of the M−R reflects now the following effects (Mollière & Mordasini 2012; Bodenheimer et al. 2013): because of their higher central temperature, the most massive planets with M ≳ 20 M_♃ have already fully burned their deuterium at this time, and are back at cooling and contracting, while lower mass planets only now undergo significant D-burning, with the highest D-burning luminosities occurring at a mass of about 17 M_♃. The radius distribution now bends upwards at mass of about 12.5 M_♃, less than the often used 13 M_♃ D-burning boundary, as the presence of massive cores shifts the limit downwards. This means that a “wave” of planets undergoing D-burning and expanding moves to lower masses in time. Lower mass planets burn the deuterium later, on a longer timescale, and with a lower L_D. Once the deuterium in the planet is consumed the planet re-contracts.

This leads to two local maxima of radius values in the population at 20 and 50 Myr, at 14.5 M_♃ and for the most massive planets at about 45 M_♃. Finally, at 500 Myr, only two subtle imprints of deuterium burning remain in the M−R relation: first, there is a barely noticeable bump of about 0.02 R_♃ at a mass of about 13 M_♃ formed by planets that still undergo mild burning. Second, the most massive planets (M ≳ 35 M_♃) still bear the imprint of their strong deuterium burning during formation. This causes a flat M−R relation diverging from the usual pattern of a decreasing radius with mass for M ≳ 5 M_♃. At even later times, also these most massive planets follow this trend that stems from the higher compressibility of more degenerate objects (Chabrier et al. 2009).

5.4.2. Hot gas accretion

Figure 20 shows the mass-radius relation for hot gas accretion during the formation phase at 1 and 3 Myr when most planets are accreting, and during the early evolutionary phase at constant mass at 20 Myr. Due to the deposition of the shock luminosity into the planet, the radius is, for planets that have completely undergone the fast contraction after the detachment, an increasing function of the mass, in contrast to the cold accretion case. Such an increase of R with M is also seen in classical hot start models which do not calculate the formation of the planet but assume a very hot state of the planet as an arbitrary initial condition.

In the plot we have also included three empirical mean M−R relations. For planets with 1 <M/M_♃< 10, the mean radius is approximately R ≈ 2.3 R_♃ × (M/M_♃)^0.22 at 1 Myr which means that these planets have quite large radii. Even larger radii are seen for planets in the fast contraction phase (M< 1 M_♃), but this phase is short, a few 10⁴ years. At 3 Myr, we find a mean radius of about R ≈ 1.6 R_♃(M/M_♃)^0.16. Interestingly, this is quite similar to the radii in the hot start simulations with arbitrary initial conditions of Marley et al. (2007) at 1 Myr (their Fig. 4), but for already fully formed planets. At 20 Myr, the radii are about 1.3 R_♃, already weakly dependent on mass as it is the case for mature planets (Paper II), with a difference in R between 1 and 10 M_♃ of only about 0.2 R_♃. This is again similar to the classical hot start models. The flat M−R in this phase is a consequence of a competition between the higher entropy in the more massive planets increasing R and their higher compressibility due to a weaker ionic contribution in the EOS, decreasing R. During formation, there is a spread of about 1 R_♃ around the mean value.

During formation, no very clear imprint of D-burning is seen in contrast to the cold accretion case. This is not surprising, because for hot accretion, D burning only delays the contraction, but does not lead to a re-increase of the planetary radius after the collapse (Mollière & Mordasini 2012). Therefore, one notes that for planets undergoing D-burning, there is an absence of small planetary radii, but in general, the M−R is relatively smooth for M ≳ 1 M_♃, at least during the formation phase. During evolution, at 20 Myr, the M−R is, in contrast, similar to the cold case with the characteristic upturn of the radii at about 10–12 M_♃. Compared to the cold case, the radii of these planets are slightly bigger, by about 0.1 to 0.2 R_♃. We thus see that for planets that are sufficiently massive, D-burning tends to partially erase the effect between cold and hot accretion (Mollière & Mordasini 2012) as in both cases the planets approach the gravothermal state obtained after the complete deuterium reservoir has been burned. The way this state is approached, and at which moment it happens, is however affected by the hot/cold accretion mode.

6. Summary and conclusions

In this paper, we studied for the first time the statistics of planetary luminosities during formation and during evolution in the framework of the canonical core accretion scenario for giant planet formation (Pollack et al. 1996; Bodenheimer et al. 1980) using the tool of population synthesis. Since the effects of the accretion shock are thought to be important, but are currently not definitely predicted by theory (see Szulágyi & Mordasini 2017; Marleau et al. 2017) nor constrained by observations, we considered both a cold-nominal and a hot population, in which the kinetic energy of the incoming gas is respectively fully absorbed by the planet or fully radiated away. We also calculated a cold-classical case where, as in the pioneering work of Marley et al. (2007), it was additionally assumed that planetesimal accretion stops artificially once a giant planet enters the disk-limited gas accretion (detached) phase and that the planets form in situ.

We discussed three fundamental statistical properties of the formed planets: the planetary mass-luminosity relation during both formation and evolution (Sect. 5.1), the mass-entropy diagram at the moment when the protoplanetary disk disappears (Sect. 5.2), and finally the luminosity distribution as a function of time (Sect. 5.3). We furthermore revisited the mass-radius relation that was extensively discussed in Paper II in Sect. 5.4, now including the effect of deuterium burning as implemented in Mollière & Mordasini (2012). Our main findings are as follows:

• 1.

  The planetary mass-luminosity relation during the formationphase is shown for the cold-nominal, hot, and cold-classicalpopulations in Figs. 2, 4,and 5, respectively. Important features include:

  • (a)

    At a stellar age of 1 Myr, the accretionluminosity of giant cold-nominal planets almost alwaysdominates over the internal luminosity. If the planetary shockluminosity L_shock can be measured observationally from accretion trackers, and the gas accretion rate Ṁ_XY can be constrained from disk and stellar observations like the stellar accretion rate, the planet’s mass can be estimated by (Eq. (4)) (16)with a radius R ≈ 1.8 R_♃ (at 1 Myr; 1.5 R_♃ at 3 Myr; and 1.3 R_♃ at 8 Myr) approximately independently of mass. For known Ṁ_XY, this is accurate to about 20 % for M = 3–17 M_♃ (Sect. 5.1.1).

  • (b)

    At ≈3–5 Myr, the L−M diagram is characterized by two groups: an upper accreting sequence of protoplanets (with L ∝ M), whose parent disk has not yet dissipated, and an evolving sequence (with L ∝ M²), which is already on standard cooling tracks. Over time, all luminosities drop due to decreasing accretion rates and normal cooling, respectively, and planets move onto the evolving sequence (Sect. 5.1.1).

  • (c)

    Comparing the highest total luminosities at a given mass (strongest accretors), planets in the cold-nominal population have during formation a higher total luminosity than the brightest planets in the hot population of the same mass, which only have the interior luminosity. This situation is inverted relative to the subsequent evolutionary phase at constant mass. This also means that the effective temperature of accreting planets is higher for cold accretion than for hot accretion at a given planetary mas. This inversion comes from energy conservation (Sect. 5.1.3).

  • (d)

    In the cold-classical population, we recover for a significant group of planets the Marley et al. (2007) result of a low luminosity (cold start) independent of planet mass (Sect. 5.1.4). However, most planets exhibit a warm start because of the core-mass effect, from the disk mass and metallicity distribution. The cold-classical population is the population with the greatest variation in the (total) luminosity at a given mass. A large spread could also result if the core-mass effect is not as efficient as assumed here (leading to lower L_int) but if the efficiency of the accretion shock in radiating away the accretional energy is less than 100%.

  • (e)

    For giant planets in their main accretion phase, the shock luminosity L_shock is higher than the internal luminosity L_int by a factor 2 to more than one order of magnitude. The ratio L_shock/L_int increases with mass up to the deuterium burning limit (Sect. 5.1.5).

• 2.

  We performed simple comparisons to observations ofembedded (and potentially accreting) planets in the formationphase and older purely cooling planets in the evolutionary phaseand found the following:

  • (a)

    Comparing the cold-nominal population to HD 100546 b, the observed luminosity could be matched by a mass of 0.2–9.7 M_♃, depending on the (unknown) relative importance of the accretional and internal luminosity. For the lowest mass planets, the accretion of planetesimals contributes significantly to the luminosity, meaning that potentially very luminous impacts may occur frequently on this planet. Results are similar for hot starts due to the core-mass effect (Sect. 5.1.6).

  • (b)

    The same exercise for LkCa 15 b yielded the result that masses of 1–2 M_♃ are consistent with the inferred Hα luminosity and the (estimated) minimum gas accretion rate (Sect. 5.1.6).

  • (c)

    We found that it is difficult to constrain the mass from a measured total luminosity (L_int + L_shock) alone during formation because of the possible contribution of accretion. A comprehensive view with spectroscopic determination of the different contributions (planetary internal, planetary accretion shock, and circumplanetary disk luminosity) combined with other indicators (accretion rate, disk structure, etc.) should be used to get a clearer picture (Sect. 5.1.6).

  • (d)

    Looking at β Pic b, we find it is with 10–13 M_♃ at the interesting transition near the deuterium-burning limit. During the main-sequence lifetime of β Pic of about 2 Gyr, β Pic b will typically burn about 10% of its deuterium reservoir if its mass is on the lower side of the allowed range and up to 90% if its mass is rather 13 M_♃. The synthetic analogs have a wide range of heavy element mass fractions of 0.005–0.1, with a typical value of about 0.015. This corresponds to an enrichment relative to solar of about 0.3 to 7, with a typical value of around 1 (Sect. 5.1.10).

  • (e)

    For 51 Eri b, the luminosity measurement implies M = 1.7–3.6 M_♃, a somewhat wider range than the hot-start masses reported by Macintosh et al. (2015) (Sect. 5.1.10). This wider range is a consequence of the intrinsic scatter in the M−L relation at young ages. It stems from different formation histories. Planets in the cold-classical population have in contrast masses between about 1.7 and almost 10 M_♃ (Sect. 5.1.11).

• 3.

  One of the key findings of this work is theimportance of the core-mass effect (Mor-dasini 2013; Bodenheimeret al. 2013): Due to it and to thelarge mean amount of heavy elements expected in giant plan-ets (Thorngren et al. 2016),the luminosities in the cold-nominal population are almostcomparable to those found for hot accretion. Hot or at leastwarm planets could thus indeed be the expected outcomefor core accretion (but see also Sect. 2.2 andAppendix C), and not the very low luminosi-ties in Marley et al. (2007)where the accretion of planetesimals is artificially shut off(Sect. 5.1.7).

  Nevertheless, even at an age of 20 or 50 Myr, well into the evolutionary phase, there is still an intrinsic scatter in luminosity in the cold-nominal population by a factor of respectively ≈1.5–2 or 20–50% at a given total mass, mainly due to the different core masses (Fig. 9). The actual “spread” in the evolution could in reality be even wider than found in our model. Relevant factors could be different opacities, envelope metal enrichments, non-solar stellar masses, and other effects that are not included in this current generation of models (such as sub-unity shock efficiencies, complex infall geometries, or planetary magnetic fields). This should be critically kept in mind when deriving masses from measured luminosities. It means that even if the luminosity could be determined observationally with a vanishing error bar, there is no one-to-one translation into a single planet mass (Sect. 5.1.8).

• 4.

  In Fig. 12, we presented for all three populations a key diagnostic outcome of a formation model, the entropy “tuning fork” diagram of exoplanets. This shows the specific entropy in the convective zone as a function of planet mass at the end of the formation phase when the protoplanetary disk disappears. Fig. 13 shows the corresponding luminosities.

  A new aspect is that due to the different formation histories, there is a large spread of about 1 to 1.5 k_B/baryon at a given mass in the hot and cold-nominal population. In the cold-classical population, the post-formation entropies even have a ≈2 k_B/baryon scatter at about 5 M_♃ which results primarily from the spread in core masses as demonstrated by Fig. 14. This means that the entire phase space from cold through warm to hot initial states is populated even if only the limiting cases of completely hot and cold gas accretion are considered. Appendices A and B contain fits to the post-formation properties of the synthetic planets which are of interest as initial condition for evolution models.

• 5.

  Next, we examined an important statistical result of this work, the planetary luminosity function (i.e., the M–L distribution marginalized over mass) during formation and evolution. The main results are (Figs. 15 and 17):

  • (a)

    At a given time, the L distribution mainly reflects the mass distribution. The distribution of log (M) in the synthetic populations is relatively flat in the giant planet regime (see Fig. 18 and Mordasini et al. 2009b), corresponding to a distribution scaling approximately with M^-1 in linear units, in good agreement with observations.

  • (b)

    Up to ca. 5 Myr, one sees an approximately flat part (in log L) of high luminosities, and strong upturn starting at log (L/L_⊙) ≈ −6. This upturn divides lower mass subcritical planets that are in the attached phase from forming giant planets that are detached. The luminosity of these lower mass planets (M ≲ 10–40 M_⊕) is mainly powered by the accretion of planetesimals (Sect. 5.3.1).

  • (c)

    In the pure cooling phase (after ~10 Myr), the luminosity distribution for non-deuterium-burning giant planets is, to first order, uniform in the logarithm of the luminosity over four orders of magnitude. This flatness comes from the 1 /M mass function from about 0.3 to 10 M_♃ combined with the fixed-time L ∝ M² scaling (Burrows & Liebert 1993). This is an important results result of this study and has important implications for example for the expected yield of direct imaging searches, especially once they start to probe closer-in giant planets, which are more likely to be formed by core accretion (Sect. 5.3.1).

  • (d)

    At even higher masses, there is a clearly visible, ≈0.5 dex-wide peak in the luminosity distribution that is separated from the remaining distribution. The peak is located at log (L/L_⊙) ≈ −4 and −5.5 at 50 and 500 Myr, respectively (Sect. 5.3.1). It is due to deuterium burning, which slows the cooling or even reverses it and thus leads to a flattening of L(t).

  These points pertained to the cold-nominal and hot population, with no significant difference between the two (Sect. 5.3.2).

  The cold-classical population has a very different luminosity distribution, with a strong local maximum near log (L/L_⊙) ≈ −6 at 10 Myr instead of a flat distribution, as all giant planets with masses between about 1 and 10 M_♃ and core masses of less than 10–20 M_⊕ in this population approximately have a mass-independent luminosity after formation, as seen in Marley et al. (2007). This strong local maximum remains visible in the distribution up to an age of at least 100 Myr. This very different luminosity distribution with a peak can clearly serve as a statistical diagnostic of the thermodynamics of the formation process once a higher number of directly imaged exoplanets is known (Fig. 17).

• 6.

  Finally, we studied the impact of deuterium burning on the mass-radius relationship (Sect. 5.4). We saw, as in Mollière & Mordasini (2012) and Bodenheimer et al. (2013), that for cold accretion, deuterium burning does not merely delay the contraction of the radius as for classical hot starts but also brings about a re-inflation of the planet. This leads to an increasing radius with mass at 3 Myr for masses higher than about 15 M_♃. The M–R relation was additionally seen to (i) be nearly flat from 1 to 10 M_♃ up to ≈50 Myr, and (ii) show a local peak (higher radii) near 15 M_♃.

  In the hot population, the radii increase in contrast during formation with mass and can reach 4 to 5 R_♃. D burning (which sets in later) here only delays the contraction, but does not lead to a re-increase of the planetary radius; therefore, no very clear imprint of D-burning is seen in the M−R relation during formation in contrast to the cold-accretion case.

One important point to take away from this work is that core accretion cannot be excluded as the formation mechanism based on an observed high luminosity alone, as seen in Sects. 5.1.1 and 5.1.8. However, one should bear in mind that the self-amplifying “core-mass effect” which is responsible for this as described in Mordasini (2013; see also Bodenheimer et al. 2013) requires that the planetesimals be accreted rapidly during the late attached and early detached phase and that the solids sink quickly deep into the potential well (Mordasini 2013). This has not yet been studied with giant planet formation models tracking both the thermodynamical and compositional evolution of the interior (Sect. 2.2). A preliminary analysis (Appendix C) indicates that relative to sinking, the heating by impacting planetesimals could be reduced by factors 2–3 for homogeneous mixing into the envelope, and 3–8 for no sinking.

Furthermore, our results indicate that during formation a rough estimation of the planetary mass may be possible if the planetary gas accretion rate and its accretion shock luminosity can be determined, at least for cold gas accretion where the radius is nearly independent of mass. If it is unknown whether the planet still accretes gas, the total luminosity (accretional and internal) spread at a given mass may be as large as two orders of magnitude, therefore inhibiting the mass estimation in this way. Due to the core-mass effect even planets which underwent cold gas accretion can have large post-formation entropies and luminosities. This means that alternative formation scenarios such as gravitational instabilities do not need to be invoked from a (high) luminosity point of view alone. Once the number of self-luminous exoplanets with known ages and luminosities increases thanks to future survey with small inner working angles, the observed luminosity distribution may be compared with our theoretical predictions. This comparison will eventually allow to develop a better understanding of the thermodynamics of giant planet formation, and the planet formation process overall.

Looking ahead, it will be interesting to repeat the statistical analysis as the theoretical description of the underlying physical processes in the planet formation and evolution models are improved. Of primary importance will be the coupling of predictive models (e.g., Marleau et al. 2017) for the accretion shock’s efficiency in radiating away L_shock to structure calculations. This will eliminate the current need to arbitrarily assume a fully cold or hot gas accretion. Another point is to include the radial non-constancy of the luminosity in accreting objects (Stahler 1988; Berardo et al. 2017), or the replacement of the gray atmospheric boundary conditions with more realistic non-gray models (Marleau et al., in prep.). This will lead to a direct prediction of the magnitudes instead of the total luminosity only. Also, considering multiple embryos per disk to take their dynamical interactions into account (Alibert et al. 2013) will enable us to also look at the luminosity–separation diagram, of crucial importance for direct detections. Finally, it will be interesting to perform dedicated simulations for specific recently discovered systems, tuning in particular the stellar mass and metallicity, as was done for β Pic b (Bonnefoy et al. 2013).

We finish by commenting on two recently observed extrasolar protoplanet candidates. They are potentially among the first examples where we may directly observe planet formation as it happens. Comparing HD 100546 b and LkCA 15 b, it seems possible that the two protoplanets represent two different evolutionary stages of giant planet formation. HD 100546 b could still be in the earlier, more extended, and cooler phase during the transition from the attached to the detached phase at a lower mass without significant hard accretion-driven radiation, and without a (deep) gap in the protoplanetary disk. LkCa 15 b would in contrast correspond to the later, fully detached phase characterized by a higher mass, higher temperatures, gap formation, and hard accretion shock radiation. To firmly establish this, more observations seem necessary, also to disentangle the impact of the circumplanetary disk. But directly observing the different theoretically predicted main phases of (giant) planet formation would represent an important step forward to a better understanding of this process. We note that the model predicts (see Sect. 5.3) that protoplanetary disk should contain a much higher number of still fully embedded, even lower mass protoplanets that are in the even earlier attached phase. Past radial velocity and transit surveys have shown that such low-mass planets are frequent. They are characterized during formation by low temperatures (a few 100 K) and low luminosities ≲10^-6L_⊙, making them difficult to detect and disentangle from disk features. As shown by van Boekel et al. (2017), comparisons of the synthetic populations presented here with the expected detection limits of the Mid-Infrared ELT Imager and Spectrograph (METIS; Brandl et al. 2014) indicate that at least the more massive planets in this stage could potentially nevertheless be detected with this instrument. This opens up the fascinating possibility to directly observe the “when” and “where” of this crucial phase of planet formation, which would yield observational constraints of uttermost importance for planet formation theory.

---

Acknowledgments

3M thank Y. Alibert, W. Benz, B. Biller, M. Bonavita, M. Bonnefoy, E. Büenzli, K.-M. Dittkrist, R. Helled, H. Klahr, A.-M. Lagrange, M. Meyer, S. Quanz, and R. van Boekel for useful discussions. We also thank the referee J. Fortney for a helpful report. C.M. and G.-D.M. acknowledge the support from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. Parts of this work have been carried out within the frame of the National Center for Competence in Research PlanetS supported by the SNSF.

References

Appendix A: Empirical fit to the post-formation entropy s_pf

As an initial condition for evolutionary models and cooling tracks we provide a least-squares fit to the relation of the planetary envelope mass M_env to the post-formation specific entropy in the inner convective zone at the moment when the disk disappears, s_pf. For giant planets, the envelope mass is to first order equal to the total mass, or can be estimated by using the observationally inferred relation of core and total mass of Thorngren et al. (2016). Theoretically predicted M_core−M_env relations for low-mass planets can for example be found in Mordasini et al. (2014) or Lee & Chiang (2015). Considering s_pf as a function of M_env is found to reduce the scatter compared to considering it directly as a function of the total mass.

The M_env−s_pf relation is shown in Fig. A.1. For the fit, the entropy at the end of the formation phase s_pf in units of k_B/baryon can written as a polynomial of degree 9 in terms of χ = log ₁₀(M_env/M_⊕) as (A.1)The coefficients a_i for the cold-nominal and hot populations are given in Table A.1. It is important to note that these fits can only be used in the domain −2 ≤ χ ≤ 4.25 and diverge outside. The rms of residuals around the fit is in both cases 0.46 k_B/baryon, reflecting a total spread of typically about 1 k_B/baryon around the mean value for a given mass due to the different formation histories. Some low-mass planets (M_env ≤ 1 M_⊕) have clearly larger entropies than predicted by the fit. These are planets at large distances which were in the process of significant planetesimal accretion at the moment when the disk disappears. They have not yet accreted all planetesimals in the feeding zone as it is otherwise usually the case for planets at smaller separations. The associated high L_pla leads to the high s_pf.

The two fits overlap as expected for low-mass planets with M_env ≲ 30 M_⊕. For such low-mass planets no distinction of cold or hot starts exists because they do not accrete gas through a potentially entropy reducing shock, as discussed in Sect. 5.1.11. In the giant planet mass regime, the fit for the cold population predicts a slight decrease of s_pf, before it increases again due to deuterium burning. In the hot population, no such decrease occurs (see Fig. A.1). It should be noted that the s_pf for low envelope masses is likely affected by the lack of including in our model the core’s thermal contribution to energy budget.

Appendix B: Empirical fits to L, s, and [D/H] at 10 Myr

For some low-mass planets with a envelope mass of less than ~1 M_⊕, the post-formation entropy is as mentioned still heavily increased as the planet was undergoing an intense accretion of planetesimals near the end of the disk lifetime. These planets shows up in Fig. A.1 with a s_pf ≳ 10k_B/baryon. We therefore also give empirical fits to the planetary M_env−L, M_env−s, and M_env− [D/H] relations at 10 Myr when accretion has stopped for almost all planets (only a small fraction of synthetic disks has a lifetime exceeding 10 Myr, namely 390 out of 50 968), such that no exceptionally high entropies exist any more. A time of 10 Myr has often been used as a starting time for evaporation in evolutionary calculations (e.g., Lopez & Fortney 2013; Chen & Rogers 2016).

Figure B.1 shows the total luminosity L₁₀ (which in the evolutionary phase is equal to the internal luminosity), the specific entropy s₁₀ at the bottom of the gaseous envelope, and for massive planets, the fraction of remaining deuterium f_RD10 for the cold-nominal population at 10 Myr. The black lines are non-linear least-squares fits. The three quantities are plotted as a function of envelope mass rather than total mass, because this is found to reduce the scatter around the fit compared to using the total mass.

The logarithm of the luminosity can be fitted as a function of χ = log ₁₀(M_env/M_⊕) with a polynomial of degree 5, (B.1)Two sets of coefficients are needed depending on the envelope mass. For lower envelope masses in the range −3.0 ≤ χ ≤ 3.55, the coefficients in the second column of Table B.1 are used, while for 3.55 ≤ χ ≤ 4.1 those in the third column apply.

The entropy at 10 Myr s₁₀ as a function of envelope mass has a similar shape as the luminosity, but the upturn at higher masses is more prominent. Below this upturn at about M_env ≈ 1000 M_⊕, the entropy roughly lies on a straight line, indicating a logarithmic dependency. For these planets, a simple fit is given as s₁₀ ≈ 7.5 + 0.21ln(M_env/M_⊕) k_B/baryon.

The full fit to s₁₀ in units of k_B/baryon as function of χ = log ₁₀(M_env/M_⊕) for −3 ≤ χ ≤ 4.1 can be written as (B.2)where for χ ≤ 3.55 the coefficients in the second columns of Table B.2 must be used, and for higher masses the ones in the third column.

In the panel showing the entropy, one notes that for a given envelope mass between about 0.5 and 20 M_⊕, there is an upper and a lower group of points, and a depletion of points between them. It is found that this is related to the semimajor axis of the planet, with the upper group of points being formed by close-in planets (a ≤ 0.1–1 AU) while the lower entropy group belongs to planets at larger semimajor axes. The fact that the fit runs through the depleted region reflects that such a dependency cannot be captured by the simple one parameter fit presented here.

Finally, we also give a fit of the remaining fraction of deuterium f_RD10 as a function of χ = log ₁₀(M_env/M_⊕), since for deuterium burning planets, also this quantity must be specified for a complete description of the initial conditions for the evolutionary phase. The initial deuterium number fraction [D/H] is 2 × 10^-5 (Mollière & Mordasini 2012). As one would expect, the right panel of Figure B.1 shows that below a certain χ, no deuterium has been burned, while above a certain limit, no deuterium remains. This suggest a functional form of the fit of (B.3)The fit parameters a and b are again determined by the least-squares method which yields a = 0.101561 and b = 36.616373.

Appendix C: Considerations on the efficiency of the core-mass effect

In this section we present an exploratory study of the efficiency of the core-mass effect in the case that planetesimals do not sink to the core. It is clear the future work will have to investigate this in a self-consistent fashion (see Venturini et al. 2016; Lozovsky et al. 2017, for works that do this for the composition of protoplanets without addressing the luminosity and the solid accretion during gas runaway, though).

First, we assess whether much less potential energy is liberated as L_pla,mix by planetesimals that mix throughout the envelope instead of sinking, i.e., whether we are greatly overestimating the heating by impacting planetesimals L_pla that in the sinking approximation we use in the syntheses is given approximately by (C.1)In the equation G is the gravitational constant, M_core, R_core the core mass and radius, and Ṁ_Z the accretion rate of solids. This expression is strictly speaking only valid if the envelope mass is negligible. Otherwise, the numerically obtained value of the difference of the gravitational potential between the surface of the core and infinity must be used to obtain the luminosity. This is done in the numerical calculations, but for the phase we are mostly interested in (total mass M ≲ 100 M_⊕) where core and envelope mass are comparable, Eq. (C.1)only differs by less than 30% from the numerically found value for the case we studied. Therefore we use it here because of its simple analytical form.

If in reality the planetesimals do not reach the core but get dissolved into the envelope less potential energy gets liberated. No direct core hits are expected for, e.g., ≲100 km-sized rocky planetesimals plunging radially into planetary envelopes of ≳3 M_⊕ (Mordasini et al. 2006, 2015).

We finally also consider that solids do not sink at all, but stay at the place where the planetesimals were destroyed in the envelope. This should give a lower limit on the impact heating, L_pla,nosink.

Appendix C.1: The case of a n = 1 polytrope

Before we study more realistic (numerical) models, we calculate the difference in energy deposition in a n = 1 polytrope when matter is added at the outer surface and at the center. While an n = 1 polytrope can clearly not capture all aspects of Jupiter like the presence of a core, it nevertheless represent a useful analytical representation of Jupiter’s basic structure nowadays (Hubbard 1974; de Pater & Lissauer 2010). The density as a function of radius r in a n = 1 polytrope is given as (C.2)where ρ_c is the central density and R the total radius. From this, the gravitational potential Φ is found by integrating the Poisson equation in spherical symmetry (C.3)Setting the zero point of Φ at infinity, and requesting that Φ and its spatial derivative are continuous at R, we find a potential (C.4)So if mass is added at the surface of the planet at a rate Ṁ, an accretion luminosity of GMṀ/R is generated, as expected. If the mass is instead added at the center, we have to consider lim_{r → 0}Φ(r). As lim_{x → 0}sin(x) /x = 1, one finds lim_{r → 0}Φ(r) = −2GM/R. This means that if material is sinking to the center, twice as much luminosity is generated. This result that the two expressions only differ by a factor of a few is repeated also in the numerical results below. It reflects that for the acceleration of the sinking material, the enclosed (and thus accelerating) mass diminishes, but also 1 /r² becomes larger.

Appendix C.2: Uniform mixing: L_pla,mix

For the full mixing case, we calculate the change in potential energy of the planet due to increasing its metallicity uniformly in every layer, yielding L_pla,mix. We address the assumption of uniformness afterwards by considering several timescales describing the system.

Let ρ_XY(r) and ρ_Z(r) be the time-dependent local densities of hydrogen and helium and of metals, respectively, with total density ρ = ρ_XY + ρ_Z, and let the total mass of hydrogen and helium and of metals in the envelope only be M_XY and M_Z,env, respectively. The local and global metallicities are defined as For concurrent accretion of gas and solids, the changes in the metallicity are therefore given by where the total envelope mass is M_env = M_XY + M_Z,env. It should be noted that since Eq. (C.7) is written at fixed radius (and not in the Lagrange frame at fixed mass coordinate), homologous contraction without global accretion (Ṁ_XY = Ṁ_Z = 0) can be sufficient to have . When Ṁ_XY = 0, Eq. (C.8) reduces to (C.9)Similarly, ignoring contraction and if the hydrogen–helium component of the mixture is not redistributed, the local change in metallicity reduces to (C.10)In the case that the metals are mixing uniformly throughout the planet, it must hold that Z(r) = Z at all radii.

Thus in particular (C.11)implying, when Ṁ_XY = 0, with Eqs. (C.9) and (C.10) that (C.12)This result will be needed shortly.

We now compute the luminosity generated by adding Ṁ_Z homogeneously into the envelope. To simplify, let us assume that the potential Φ(r) inside the planet is fixed and that the core does not contract. For total mass M = M_env + M_core and radius R = R_core + R_env, the luminosity is given by These quantities are obtained numerically in the formation model at each timestep, and will be used below to compare L_pla,sink, L_pla,mix, and L_pla,nosink.

Appendix C.3: No sinking: L_pla,nosink

Besides the case of homogeneous mixing, we can also consider the limiting case that the solids do not sink at all from the point where they are deposited into the envelope by the impact. From the calculation of the trajectories of impacting planetesimals (Mordasini et al. 2006) we know the (enclosed) mass and radius in the envelope of maximum energy deposition, M_maxE and R_maxE. For the large impactors we consider, this is also the place of maximum mass deposition (e.g., Mordasini 2014). When the planet is still attached or early after detachment, M_maxE is similar to the total mass, but R_maxE is much smaller than the total radius, as the planet has a very voluminous thin outer atmosphere through which the planetesimal penetrate. In late stages when the planet has fully contracted, both the mass and radius of maximum energy deposition are similar to the total mass and radius.

Figure C.1 shows the luminosity caused by radial planetesimal impacts at 1.95 Myr shortly after detachment for the planet discussed in Sect. 4. The properties of the planet at this moment are shown in Table C.1.

The plot shows the two contribution to L_pla separately. The first results directly from the impact of the planetesimal itself. The second contribution is due to the subsequent settling of the debris to the core after the planetesimal was destroyed. In the sinking approximation, the sum of the two is assumed to heat the planet. The shape of the direct contribution shows that the energy and mass is deposited essentially in one narrow location at 2.2 R_♃ because of the terminal explosion that is typical for large aerodynamically disrupted bodies as it was also the case for Shoemaker-Levy 9 (Mac Low & Zahnle 1994). This location corresponds to the locus of maximum energy deposition M_maxE and R_maxE and also the location where the planetesimal’s mass is deposited. The plot and Table C.1 show that in terms of radius, the planetesimals penetrate deep into the planet, in particular into the convective zone. But in terms of enclosed mass, they are destroyed relatively close to the surface “below” about 2 M_⊕ of envelope. We also see that the settling contribution is about 6 times as large as the direct contribution (cf. Fig. C.2).

With M_maxE and R_maxE we obtain the planetesimal luminosity in the no sinking case as (C.19)In Fig. C.1 this corresponds to the violet curve, i.e., the direct contribution. The reduction of the luminosity in the no sinking case relative to the sinking case is (C.20)When the planet has not yet accreted a significant envelope and planetesimals penetrate to the core, this expression becomes unity, as can be seen in Fig. C.2 below. At crossover, M_maxE/M_c ~ 2, and the syntheses show that at that point R_core/R_maxE lies between 0.08 and 0.24, with a typical value of about 0.15. This corresponds to a planetesimal heating that is reduced relative to the sinking case by a factor ~3 with a spread between about 2 to 6. When M_env ≫ M_core, M_maxE ≈ M and R_maxE ≈ R. Using the result from the syntheses that for typical giant planet during formation with masses between 100 and 1000 M_⊕, R_core/R_maxE ~ 0.2, we see that L_pla,nosink/L_pla,sink could become larger than unity. This is however an artefact of approximating the potential at the surface of the core as GM_core/R_core which breaks down exactly if M_env ≫ M_core. In reality, the ratio is given by the ratio of the actual gravitational potentials, Φ_maxE/ Φ(R_core) ≈ Φ(R)/Φ(R_core). This expression is always smaller than unity, as otherwise, parts of the planet would be accelerated outward, which is of course not the case.

Appendix C.4: Comparison of sinking, uniform mixing, and no sinking

Figure C.2 show L, L_shock, L_pla,sink, L_pla,mix and L_pla,nosink during the early formation phase of the planet discussed in Sect. 4. This includes the phase when the planet goes into gas runaway and then disk-limited gas accretion. The planet detaches from the disk at about 1.94 Myr. This is the time when the core-mass effect could be acting. During the time interval shown, the planet’s core grows from about 2 to 48 M_⊕, almost the final core mass, while the total mass increases from 2 to 335 M_⊕.

In the left panel we see that after the initial build-up of the core, L_pla,sink, L_pla,mix and L_pla,nosink follow a similar temporal pattern. We also see that before significant gas accretion starts at around 1.9 Myr, L_pla,sink provides the dominant luminosity source of the planet. The assumption of sinking has therefore important consequences for the efficiency of gas accretion and the moment when runaway gas accretion occurs, as already demonstrated by Pollack et al. (1996). After detachment, the gas accretion shock luminosity L_shock grows quickly.

The plot shows that for homogeneous mixing and no sinking, the heating of the planet by impacting planetesimals is reduced relative to sinking, as expected. The reduction is stronger in the no sinking case relative to mixing, which is expected as well, as in the mixing case, some material still sinks deep into the potential well. In the right panel, we have divided the different luminosities by L_pla,sink so that we can see how strongly the energy input is reduced. One sees that relative to sinking, the heating is reduced by a factor 2 to 3 for homogeneous mixing, and by a factor 3 to 7 for no sinking. We thus find that also a homogeneous mixing of the solids into the envelope instead of full sinking provides a significant energy source, at least for this case.

To give an impression how general these results are, we show in Fig. C.3 the ratio L_pla,nosink/L_pla,sink (Eq. (C.20)) as a function of planet mass for the CD777 synthetic population at 2 Myr. Only planets where M_env ≤ 2M_core are included such that the approximation of the core potential as GM_core/R_core still approximately holds. The color code gives the envelope mass fraction, M_env/M. We see that in the no sinking case, the heating by planetesimals is about 2.5 to 8 times less strong than for sinking when the planets are near runaway, i.e., when M_env/M ≳ 0.5. This is comparable to the result of Fig. C.2.

An interesting point to address in future work is to check whether such a non-central energy input could shut down convection. For this it is however necessary to consider that impacts are not a spherically symmetric phenomenon, as indicated by the timescale estimates below.

Appendix C.5: Mixing timescale of metals in the envelope

We have showed in the previous section that the uniform mixing of the metals provided by thermally and mechanically destroyed planetesimals represents a source of luminosity comparable to the case when the planetesimals hit the core. The question is now whether this mixing can happen fast enough for it to be relevant for the cooling, i.e., whether the mixing timescale for the whole envelope τ_{mix, env} is much smaller than the Kelvin-Helmholtz timescale τ_KH. It should be noted that contrary to what one might intuitively imagine, the mixing timescale throughout the envelope does not need to be shorter than the time between impacts at a given location (as it nevertheless is, as Fig. C.4a will show); in the case it is not, a “pile-up” of metals in one layer still only needs to mix with the rest of the envelope faster than the envelope cools in order to contribute significantly to the luminosity.

The time needed for the composition to become uniform throughout the planet can be calculated from the mixing time within one convective eddy, equal to its lifetime, and the picture of diffusion across the eddies. Thus, as implemented in the Bern planet formation code by Mollière & Mordasini (2012), which implicitly defines the average mixing time in one eddy τ_{mix, loc}; the integral is computed over the convective region and weighs by mass, v_conv is the convective velocity in every layer, the mixing-length parameter is α_MLT = 1, and the total number of eddies in the radial direction is given by (C.23)Equation (C.21) represents an approximate upper bound since in general R_maxE could be closer to the middle of the planet, reducing the number of eddies through which the metals must diffuse.

Figure C.4a shows the relevant timescales as a function of time for the 5 M_♃ example studied in the previous sections. The global mixing timescale varies from τ_{mix, env} ~ 1 yr at 1.8 Myr, when the planetesimals do not reach the core anymore, to τ_{mix, env} ~ 100 yr at 2.3 yr, when the planet has grown to 1 M_♃. Over this period, there are always between 15 and 20 convective eddies. The reciprocal of the mass-weigthed inverse convective velocity is 0.003–0.01 km s^-1 decreasing with time⁹, the mass-weighted pressure scale height is of order 0.2 R_♃, and the mass-weigthed average mixing time for a single eddy is τ_{mix, loc} ~ few days at 1.8 Myr, going up to τ_{mix, loc} ≈ 2.5 months at 2.3 Myr and 6 months at the end of formation. Thus the global mixing time according to Eq. (C.21) is between τ_{mix, env} ~ 1 yr at 1.8 Myr and 100 yr at 2.3 Myr. However, the Kelvin-Helmholtz timescale¹⁰τ_KH = E_tot/L_int, with E_tot and L_int the total energy and luminosity, goes over this period from around τ_KH ~ 0.1 to 1 Myr, for a ratio τ_KH/τ_{mix, env} ~ 10⁵–10⁴. Thus, it is a robust result that the mixing should happen much faster than the planet cools.

We note that we ignore the possibility that a sufficiently strong metallicity gradient could inhibit convection on the largest scales (Leconte & Chabrier 2012; Leconte & Chabrier 2013; Vazan et al. 2016), preventing a homogenization of the composition. See also Lozovsky et al. (2017), Berardo & Cumming (2017).

Appendix C.6: Impact timescale of planetesimals

We now estimate whether the impacts could shut off or at least impair convection by generating a spherically symmetric hot layer of material at R_maxE, which is close to the outer radius but within the convective zone (see Fig. C.4a). Heating by impacts will hinder convection if there is not enough time between impacts for the effect of a planetesimal to become washed out before the next one arrives. This means that the convection would “see” a spherically-symmetric blanket of planetesimals falling down. We thus need to determine the global impact timescale, the size of the region affected by an impact, and from this then the local impact timescale.

The global average time between impacts is τ_{impact, glob} = m_pmal/Ṁ_Z for a monodispersive planetesimal mass distribution of m_pmal. Inside of the iceline, where our example planet remains during its formation, the planetesimals are rocky, with a material density ρ_• = 3.2 g cm^-3. With a diameter d_pmal = 100 km, the mass of each planetesimal is m_pmal = 1.7 × 10²¹ g, implying (C.24)The results of Eq. (C.24) are shown in Fig. C.4b. We note in passing that these impacts are thus quite frequent and might be observable by dedicated monitoring campaigns. The caveats are however that at this stage the protoplanet might be hidden by the circumstellar (and possibility also the circumplanetary) disc(s) but mostly that the phase in which these impacts are expected is very short compared to the typical age of young stars.

One can estimate the size of the region affected by a single impact in at least two ways. One measure ℓ_ρ is given by the volume at which the mean density of the dissolved metals drops to the background density at the impact location ρ_amb: (C.25)Another estimate of the size of the explosion region is given by the Sedov blast wave solution for a point explosion, (C.26)where E_impact is the energy deposited by the impact of one planetesimal (which is very localized; see Fig. A.1b of Mordasini 2014, even though for icy planetesimals) and e_{int, amb}, γ_{eff, amb}, and P_amb are respectively the ambient internal energy density, effective adiabatic exponent γ_eff = P/ (ρe_int) + 1 ≈ 1.1–1.4 at the relevant conditions, and pressure at the blast location. The explosion or disruption energy is (C.27)with v_esc(R_maxE) the escape velocity at R_maxE. For completeness, the density at R_maxE varies from ρ_amb = 10^-3 at the onset of detachment to 10^-2 g cm^-3 at 2.3 Myr, while P and T range typically from 200 to 10⁴ bar and 5000 to 10⁴ K, respectively.

The energy deposition is typically around 5 × 10³³–10³⁴ erg per planetesimal and is at any time a factor 10⁹–10¹⁰ smaller than the planet’s binding energy.

As Fig. C.4b shows, the density-based and the energy-based lengthscales (ℓ_Sedov and ℓ_ρ) agree quite well and yield a typical ℓ ≈ 300–1000 km.

The eddy size at R_maxE is also shown and is roughly one order of magnitude larger than ℓ. Also, the mass contained in the background gas in the impact volume m_impact = 4πℓ³/ 3 × ρ_amb is a factor of 2 to 10 larger than the planetesimal mass, with m_impact/m_pmal = γ_amb(γ_{eff, amb}−1)ℳ² where γ is the heat capacity ratio and ℳ the Mach number of the impact¹¹. Therefore, by all counts the effects of the impact are very localized in the sense that ℓ ≪ H_P ≪ R_maxE. This also justifies a posteriori our use of the point explosion expression in Eq. (C.26) instead of the expression for a line charge.

We can now estimate the time between impacts into the same region of size ~ℓ³ at a height R_maxE. Since there are –10⁶ (from 1.8 to 2.3 Myr) such regions in a shell at height R_maxE, the timescale is approximately (C.28)This is also shown in Fig. C.4a and is roughly τ_{impact, loc} ~ 1000 yr at crossover to 0.1 Myr at 2.3 Myr (within a factor of ten since ℓ_ρ and ℓ_Sedov differ by at most a factor of four, at 2.3 Myr). The local impact timescale is thus 10³–10⁶ times longer than the eddy mixing time. Also, the local cooling time τ_cool = E_impact/L_int = 4πℓ³/ 3 × e_int is 1–10 h. Therefore, the impacts are not spherically symmetric as far as the convection is concerned, in the sense that each impact can be locally “forgotten” before the next one takes place. As a consequence, one would not expect the energy deposition of planetesimals to be a thermal impediment to convection.

Appendix C.7: Summary

In summary we find with this preliminary analysis that relative to the sinking case, the planetesimal impact heating during runaway could be reduced by a factor 2 to 3 for homogeneous mixing, and a factor 3 to 8 for no sinking. How this translates into post-formation luminosities needs to be assessed with future work. It seems that the reduction by less than an order of magnitude would still lead to luminosities higher than in the classical models of Marley et al. (2007) where no planetesimal impact heating at all occurs during this phase. A more definitive answer will however require a fully self-consistent treatment of the problem which is not trivial (non-spherical symmetry of the impacts and thus mass deposition, modification of the EOS, local variation of the opacity, transport processes under the influence of compositional gradients, etc.). However, such studies are important given the significant implications for the formation and evolution itself (e.g., Venturini et al. 2016; Vazan et al. 2016; Lozovsky et al. 2017) as well as for the detectability via direct imaging.

All Tables

All Figures