© E. F. Linder et al. 2019

1 Introduction

During the last few years, the Kepler satellite has detected numerous exoplanets of which many are in the sub-Neptunian or super-Earth mass range (e.g. Batalha et al. 2012; Burke et al. 2014; Fressin et al. 2013; Petigura et al. 2013). Different from the solar system, especially in close-in orbits, sub-Neptunian planets seem to be very abundant in the solar neighbourhood (Howard et al. 2012). Various studies on sub-Neptunians and super-Earths have been conducted. For example, Bodenheimer & Lissauer (2014) studied the origin and evolution of low-density planets in the mass range from 1–10 M_⊕, which are within 0.5 au from their star. In particular, they wanted to find out if these planets formed in situ or further out and then moved inwards. Another analysis conducted by Chen & Rogers (2016) was dedicated to computing mass-radius-composition-age relations for low-mass planets and how these depend on the evolution history of the planets. Finally, Jin & Mordasini (2018) studied whether through photoevaporation a certain planetary composition is revealed.

In the literature, evolutionary calculations for gas giants are abundant (e.g. Burrows et al. 1997; Baraffe et al. 2003; Podolak et al. 1995; Nettelmann et al. 2013; Fortney & Nettelmann 2010; Fortney et al. 2011). However, studies of the thermodynamic evolution of low-mass planets have been scarce so far (e.g. Howe & Burrows 2015; Nettelmann et al. 2013; Fortney et al. 2011; Beichman et al. 2010). In this paper we want to extend calculations of the thermodynamic evolution and cooling curves to lower mass planets in a small parameter study. An important initial condition for planetary cooling calculations is the post formation entropy (e.g. Marley et al. 2007; Marleau et al. 2017; Mordasini et al. 2017). For low-mass planets, the core-envelope mass ratios are also important. We study non-irradiated planets in a mass range of 5–636 M_⊕ (2 M_♃) and provide magnitudes corresponding to filter bands from various instruments and systems. While non-irradiated planets are simpler in the sense that the three-dimensional redistribution of insolation energy through the planetary atmosphere does not have to be modelled, they still present a unique set of challenges when trying to model their atmospheric structures and spectra. For one, we expect the atmospheres of the objects studied here to be heavily enriched in metals, because the bulk enrichment of a planet appears to be a function of its mass (Miller & Fortney 2011; Thorngren et al. 2016), and so may the atmospheric enrichment (see Fig. 3 in Mordasini et al. 2016). The degree to which the the bulk enrichment of a planet is visible in its atmosphere is not straightforward to assess and hence remains an open challenge (see the discussion in Sect. 2.4.4 in Mordasini et al. 2016). In addition, the question of when and how clouds form in self-luminous planets is far from being understood. The challenge lies in understanding and quantifying the multitude of microphysical processes that lead to cloud formation and evolution (Rossow 1978). Because of this, most cloud models currently in use are heavily simplified or parametrized (Tsuji et al. 1996; Ackerman & Marley 2001; Allard et al. 2001, 2003; Zsom et al. 2012; Mollière et al. 2017), and remove certain cloud species as a function of temperature in an ad hoc fashion in order to mimic the settling of these species below the planet’s photosphere (Morley et al. 2012; Mollière et al. 2017). Even with the use of a sophisticated micro-physical model, the cloud properties depend strongly on the unknown vertical mixing and the detailed atmospheric chemistry and hence remain a priori under-determined without further observational constraints (Gao et al. 2018; Ohno & Okuzumi 2018). Moreover, the recovery of the optical properties of such a vast variety of potential condensates is still an ongoing endeavour (Kitzmann & Heng 2018). Finally, since some self-luminous sub-stellar objects cannot be reasonably fitted with current cloud models, an altogether different process has been suggested to affect the atmospheres and spectra of planets (Tremblin et al. 2017).

In terms of imaging observations of exoplanets, the soon-to-be-launched James Webb Space Telescope (JWST) as well as the next generation of ground-based optical and near-infrared telescopes with 30–40 m primary mirrors will probe currently uncharted parameter space in terms of exoplanet mass and luminosity. The Near Infrared Camera (NIRCam) and the Mid Infrared Instrument (MIRI) on the JWST will provide unprecedented sensitivity to cool and/or low-mass objects at near- and mid-infrared wavelengths, respectively, and instruments like the Mid-infrared E-ELT Imager and Spectrograph (METIS; Brandl et al. 2016), to be installed at European Southern Observatory’s (ESO’s) 39 m European Extremely Large Telescope (E-ELT), with its unparalleled spatial resolution and superior sensitivity compared to current ground-based instruments, will be able to search for low-luminosity objects in the solar vicinity (e.g. Crossfield 2013; Quanz et al. 2015). Therefore, evolutionary models extending to smaller masses (ice giants, super-Earths) are needed to inform these future observations and interpret possible detections.

The structure of this paper is as follows. In Sect. 2, the model and the improvements made in the planetary evolution code for calculating the evolution of (low-mass) planets are described. Following in Sect. 3 is the description of the atmospheres used in this work. Section 4 shows various example and benchmark calculations, where the results obtained in this work are compared to measurements and earlier theoretical evolution calculations. The initial conditions for the final calculations are presented in Sect. 5. After this, the results and discussion are presented in Sect. 6. Section 7 summarizes the findings and major conclusions.

2 Evolutionary and internal model

The evolutionary calculations presented here were obtained with the evolutionary model described in Jin et al. (2014), which is itself based on the model of planetary evolution of Mordasini et al. (2012a,b). This model describes the planets as consisting of three distinct homogeneous layers, namely a H/He envelope (using the equation of state (EoS) of Saumon et al. 1995), an ice layer (for planets which have accreted outside of the iceline), and a rocky core, which itself consists of silicates and iron. To address the cooling and contraction of very low-mass planets, we have extended the model in regard to two aspects.

In our previous simulations, the source of luminosity of the planets were the cooling and contraction of the H/He envelope, the radiogenic luminosity due to long-lived radionuclides, and the luminosity generated from the contraction of the solid core when the external pressure exerted by the envelope changes. The contributions resulting from the non-zero temperature of the core were, in contrast, not considered. Therefore, the contraction of the core due to a change of its mean density, because of a decrease in its mean temperature, as well as the decrease of the core’s internal energy, were neglected. While the contributionof the core to the total luminosity is negligible for H/He dominated giant planets (e.g. Baraffe et al. 2008), neglecting itfor core-dominated low-mass planets leads to inaccurate cooling sequences, as demonstrated by Baraffe et al. (2008) and Lopez & Fortney (2014). We have therefore added a first order temperature correction of the mean core density to take into account the temperature dependency of the core radius, which is described in more detail in Appendix A. We also take into account additional terms in the energy equation that is used to calculate the temporal evolution and thus the luminosity of the whole planet. This addition is described below.

As described in Mordasini et al. (2012a), the calculation of evolutionary sequences in our model is based on the fundamental relation between the change of the total energy of the planet, and its luminosity because of energy conservation, d E_tot∕dt = −L (for other energy based approaches, see Leconte & Chabrier 2013; Piso & Youdin 2014). In previous versions of our planet evolution model (except for Linder & Mordasini 2016), the thermal energy of the core was, however, neglected. As the second modification of the code, we include it here, considering both the isothermal and adiabatic case.

The total energy E_tot of the planet consists of four terms, $$E_{\text{tot}}=E_{\text{grav},\text{e}}+E_{\text{int},\text{e}}+E_{\text{grav},\text{c}}+E_{\text{int},\text{c}},$$(1)

which are the gravitational potential energy of the gaseous envelope E_grav,e, its thermal (internal) energy E_int,e, the potential energy of the core E_grav,c, and finally its thermal energy E_int,c. Additional sources of luminosity we include are the radiogenic decay (calculated as in Mordasini et al. 2012b, assuming a chondritic abundance of radionucleides), which is important for low-mass planets, as well as deuterium burning (Mollière & Mordasini 2012; Bodenheimer et al. 2013), which is in contrast negligible for the planets we study here.

The potential energy of the gaseous envelope is found as $$E_{\text{grav},\text{e}}=-{\displaystyle {\int }_{M_{\text{core}}}^M\frac{Gm}{r}}\text{\hspace{0.17em}d}m,$$(2)

where M_core is the core mass, M the total planet mass, r the distance from the planet’s centre, and G is the gravitational constant. The internal energy is $$E_{\text{int},\text{e}}={\displaystyle {\int }_{M_{\text{core}}}^{M}u}\text{\hspace{0.17em}d}m,$$(3)

where u is the specific internal energy of the H/He gas, which is directly given by the Saumon, Chabrier, and van Horn (SCvH) EoS (Saumon et al. 1995).

For the core’s gravitational energy, we assume for simplicity a (mean) density that is constant within the core, as density contrasts in the cores are smaller than in the gaseous envelope, even if this is strictly speaking not self-consistent with the internal structure model of the core. We then have for the potential energy of the core $$E_{\text{grav},\text{c}}=-\frac{3G{M}_{\text{core}}^2}{5R_{\text{core}}},$$(4)

where R_c is the core’s radius that is found as described in the previous section. To relate the pressure to the density in the core, we use the modified polytropic EoS of Seager et al. (2007) for iron, rock (perovskite: MgSiO₃), and ice.

Finally, for the core’s internal energy we consider two cases, reflecting the uncertainty in the heat transport mechanism in the core (Baraffe et al. 2008). First, as in Lopez & Fortney (2014), we consider the isothermal case, where the internal energy is given as $$E_{\text{int},\text{c},\text{iso}}=c_{\text{v}}{M}_{\text{core}}{T}_{\text{ceb}},$$(5)

where c_v is the specific heat capacity that is set to 10⁷ erg g⁻¹ K⁻¹ for rocky material (Guillot et al. 1995). As noted by Baraffe et al. (2008), this value is compatible with the predictions by the Analytic Equations of State (ANEOS) in the relevant pressure and temperature range. For (water) ice we assume c_v = 6 × 10⁷ erg g⁻¹ K⁻¹. For cores consisting of both rocky material and ice, we use the mass weighted average.

We also consider an adiabatic thermal structure of the core to estimate the core’s thermal energy content. Under the (rough) approximation that in the core the density ρ, heat capacity at constant pressure c_p, and thermal expansion coefficient α are constant with radius, one can find the temperature as a function of radius r from the adiabatic temperature gradient $$\frac{\text{d}T}{\text{d}P}=\frac{T\alpha }{\rho c_{\text{p}}}$$(6)

and the pressure as a function of radius $$P(r)=P_{\text{ceb}}+\frac{3G{M}_{\text{core}}^2}{8\pi R_{\text{core}}^6}(R_{\text{core}}^2-r^2),$$(7)

where P_ceb is the pressure at the core-envelope boundary, that is, at R_core. Integrating the first equation and replacing the pressure using the second one yields the temperature structure as $$T(r)=T_{\text{ceb}}\exp \left[\frac{\alpha G{M}_{\text{core}}}{2c_{\text{p}}{R}_{\text{core}}^3}(R_{\text{core}}^2-r^2)\right].$$(8)

This temperature structure can be integrated to find E_int,c,adia via $$E_{\text{int},\text{c},\text{adia}}={\displaystyle {\int }_0^{R_{\text{core}}}4}\pi r^2\rho c_{\text{v}}T(r)\text{d}r.$$(9)

This integral is evaluated to $$\begin{array}{l}{E}_{\text{int},\text{c},\text{adia}}=4\pi c_{\text{v}}{T}_{\text{ceb}}\rho [\sqrt{\frac{\pi c_{\text{p}}^3{R}_{\text{core}}^9}{2G^3{M}_{\text{core}}^3{\alpha }^3}}\exp \left(\frac{G{M}_{\text{core}}\alpha }{2c_{\text{p}}{R}_{\text{core}}}\right)\hfill \\ \times \text{\hspace{0.17em}erf}\left(\sqrt{\frac{G{M}_{\text{core}}\alpha }{2c_{\text{p}}{R}_{\text{core}}}}\right)-\frac{c_{\text{p}}{R}_{\text{core}}^4}{G{M}_{\text{core}}\alpha }\sqrt{\frac{\pi c_{\text{p}}^3{R}_{\text{core}}^9}{2G^3{M}_{\text{core}}^3{\alpha }^3}}].\hfill \end{array}$$(10)

In the simulations presented below, this is used as the nominal expression for the core thermal energy, since the core’s energy transport is assumed to be convective. However, as already noted by Baraffe et al. (2008) and Linder & Mordasini (2016), we find that using E_int,c,adia instead of E_int,c,iso has a much smaller impact on the results compared to neglecting the temperature dependence of the core material or its energy contribution in general. The influence on the luminosity of an isothermal versus an adiabatic core is biggest for the 5, 10, and 20 M_⊕ planets. The simulations with an isothermal core are up to 21% smaller in luminosity, whereas for simulations where the core is excluded the luminosity is up to 86% smaller than in the adiabatic case. For the 50 M_⊕ planet, the luminosity with an isothermal core is 11% smaller compared to the luminosity with an adiabatic core, and the luminosity without any core contribution is 57% smaller than in the adiabatic case. For higher masses the differences become smaller, but the general trend that simulations without any core contribution included have a much smaller luminosity compared to those with an isothermal or adiabatic core stays the same.

In both the isothermal and adiabatic cases, we assume that the temperature is continuous across the core-envelope boundary. As discussed by Lopez & Fortney (2014), this should be the case for planets with an envelope sufficiently massive that the surface of the rocky core is partially or completely molten (T_ceb ≳ 2000 K), allowing an efficient heat transfer, as also discussed by Ginzburg et al. (2016) in the context of icy cores. For the smallest planetary masses studied in this work for ages younger than 1 Gyr, T_ceb is always above 2000 K. For the heavier masses this is true for even later times. These times are much later than the time when the planets could be detected, such that this does not pose a problem. In this paper we deal with planets at large orbital distances, therefore a mechanism that is in principle also included in the evolution model, namely atmospheric escape (Jin & Mordasini 2018), can be neglected.

3 Atmospheric models

The atmospheric models provide boundary conditions for the cooling of the planet and also determine the spectral appearance of the planet. In this work we use the approach of Chabrier & Baraffe (1997) to couple externally calculated atmospheres to the interior calculation. Briefly, for a given logg and T_eff reached in the course of evolution, this simply means looking up (interpolating) the pressure and temperature in the convective very deep layers of the atmospheric model, and using this (P, T) pair as the starting point for the inward structure integration (see Mordasini et al. 2012a for the structure calculation). Further details will be given in Marleau et al., in prep. For petitCODE and HELIOS we use a pressure of 50 bar as the connecting level but we verified that taking another pressure does not change the cooling curves and that the error in the radius from neglecting the layers above50 bar is negligible (at roughly the percent level, smaller than the effects of model uncertainties such as cloud types or metallicities). For the AMES-Cond model, we took the structures available on F. Allard’s website ¹ and extracted the layer at τ_std ≈ 100, where all models are convective (this is almost but not quite the case at τ_std ≈ 30), where the index “std” stands for “standard” and refers to2.15 μm (F. Allard, 2014, priv. comm.).

Additionally, we also calculate evolutionary tracks using so-called Eddington atmospheres (Eddington atmospheres as always simply use T = T_eff at P = P_phot = (2 g)∕(3κ); Mordasini et al. 2012a). Since Eddington atmospheres as relevant for this work have been described in Mordasini et al. (2012a), we only summarize the main features of AMES-Cond and briefly describe the more recent models petitCODE and HELIOS.

When the planet leaves the atmospheric grid, we extrapolate linearly from the last two grid points for the AMES-Cond grid. For the petitCODE and HELIOS grid, we extrapolate using splines when the planet leaves the atmospheric grid. We have verified that this extrapolation is reasonable and similar to a linear extrapolation.

3.1 AMES-Cond grid

The AMES-Cond grid (Allard et al. 2001) consists of cloud-free atmosphere calculations that have been obtained with the PHOENIX code (Hauschildt & Baron 1999). The models are calculated assuming radiative-convective and chemical equilibrium. While clouds are not included in these calculations, the sequestration of elements into condensates is treated. The AMES-Cond models treat condensation in strictly local chemical equilibrium, which means that the condensated particles do not rain out and are therefore still available for chemical reactions (Allard et al. 2001). The models contain opacities important for young, high temperature brown dwarfs (e.g. oxides such as TiO, VO, as well as hydrides such as FeH and MgH). Molecules that are important at intermediate to low temperatures, such as H₂ O and CH₄, are also included, where the CH₄ line list with 47 415 lines is likely to be very incomplete when compared to modern CH₄ line lists with ~10¹⁰ lines (Yurchenko & Tennyson 2014). Convection in the AMES-Cond models is treated using mixing length theory. AMES-Cond models have been widely used in the literature for the calculation of planet evolutionary tracks, but also for the spectral fitting of brown dwarf and dwarf star atmospheres (see e.g. Baraffe et al. 2003).

3.2 petitCODE grid

We calculated another grid using the petitCODE, a 1D model that self-consistently calculates atmospheric structures and spectra of exoplanets. The main assumptions of the code are radiative-convective and chemical equilibrium. petitCODE treats condensation of solids in chemical equilibrium. This means that gas phase chemistry will be affected by the depletion of elements into solids. Rain-out is neglected, however. The choice to neglect feldspar condensation in the equilibrium chemistry calculations effectively mimics the rain-out of silicon atoms into silicates (which are included). Hence the alkalis, which would otherwise condense into feldspars, will stay in the atmosphere until Na2S and KCl condense, which seems to be confirmed by observations (Line et al. 2017). Indeed petitCODE spectra and structures agree with calculations including rain-out (Baudino et al. 2017). The equilibrium condensate mass fractions are used as an input to the Ackerman & Marley cloud model, as described in Mollière et al. (2017). This cloud model includes rain-out of the cloud particles, but this is not coupled back to the chemical equilibrium calculations.

Convection is modelled applying adiabatic adjustment: petitCODE solves the radiative temperature structure of the atmosphere from top to bottom. If, during that process, the temperature gradient from one layer to the next is found to be steeper thanthe adiabatic temperature gradient, the temperature of the bottom layer is corrected to follow the adiabatic temperaturegradient instead. This process is repeated until the bottom of the atmosphere is reached, and the full atmospheric temperaturestructure has been found. Further details can be found in Mollière et al. (2015) on how the adiabatic adjustment is implemented. The radiative transfer implementation treats absorption, emission, and scattering. Clouds can be added self-consistently, making use of the Ackerman & Marley (2001) model, or a model that simply parametrizes the cloud particle size and maximum cloud mass density (Mollière et al. 2017). In the calculations presented here, we used the Ackerman & Marley (2001) model. The gas absorbtion of the following species is considered: CH₄, HCN (ExoMol, see Tennyson & Yurchenko 2012), H₂O, CO, CO₂, OH (HITEMP, see Rothman et al. 2010), H₂, H₂ S, C₂ H₂, NH₃, PH₃ (HITRAN, see Rothman et al. 2013), and Na, K (VALD3, see Piskunov et al. 1995). Ultraviolet electronic transitions are included for H₂ and CO (Kurucz 1993). The code also includes collision induced absorption (CIA) of H₂ –H₂ and H₂–He (Borysow & Frommhold 1989; Borysow et al. 1989; Richard et al. 2012). Lastly, Rayleigh scattering is included arising from H₂, He, CO₂, CO, CH₄, and H₂ O. The cross sections are taken from Dalgarno & Williams (1962; H₂), Chan & Dalgarno (1965; He), Sneep & Ubachs (2005; CO₂, CO, CH₄), and Harvey et al. (1998; H₂O). petitCODE is described in Mollière et al. (2015, 2017). Recently, petitCODE was benchmarked against the ATMO and Exo-REM codes (Baudino et al. 2017).

The grid presented here is an extension of the grid of self-luminous atmospheres calculated for spectral fitting of 51 Eri b (Samland et al. 2017). Identical to Samland et al. (2017), we assume the clouds to consist of Na₂ S and KCl. The following parameter values were considered, spanning the grid in a rectangular fashion: T_eff = [150, 1000], ΔT_eff = 50; log₁₀(g) (cgs) = [1.5, 4.0], Δ log₁₀(g) = 0.5; [Fe/H] = [−0.4, 1.4], Δ[Fe/H] = 0.2; f_sed = [0.5, 3.0], Δf_sed = 0.5, cloud-free. As described in Ackerman & Marley (2001), f_sed is the cloud settling parameter. In contrast to the original Ackerman & Marley (2001) model, the Mollière et al. (2017) implementation assumes the cloud mixing length to be always equal to the pressure scale height. This effectively lowers the f_sed value when compared to Ackerman & Marley (2001), as described in Mollière et al. (2017) and Samland et al. (2017).

3.3 Absence of water clouds in atmosphere models

The major limitation of the petitCODE cloudy grid for the temperature range to which it is applied in this work is the absence of water clouds. Water clouds are expected to form at temperatures from 300–400 K (see e.g. Morley et al. 2012, 2014). This is especially relevant for planets with masses smaller than 20 M_⊕, which start their evolutions at temperature below 300 K. Water clouds can heavily impact the spectra by absorbing flux in the 4.5 μm region at higher temperatures (300 K) and heavily absorbing across the full spectral range for even cooler planets (Morley et al. 2014). Thus, the models are not very realistic at low temperatures, until very low temperatures are reached where the water cloud has disappeared below the photosphere and cloudless models become relevant again (until methane and ammonia condense). This has a potentially large effect on the surface boundary conditions as well as the predicted magnitudes.

3.4 HELIOS grid

A third atmosphere grid is produced with the HELIOS code (Malik et al. 2017), which is an open-source 1D radiative transfer code developed specifically for exoplanetary atmospheres. As in the petitCODE, HELIOS self-consistently calculates the temperature structures and resulting emission spectra in radiative-convective equilibrium through radiative iteration and convective adjustment. The chemistry model FASTCHEM is used (Stock et al. 2018), which includes 550gas-phase reactants to calculate the atmospheric abundances in chemical equilibrium. In the current version of FASTCHEM, condensation is not taken into account for the gas phase abundances. This may lead to an overestimation of the gas absorption in the cool temperature regime, for example of water vapour for T ≲ 300 K. The following opacities are included from these line lists – EXOMOL: H₂ O (Barber et al. 2006), CH₄ (Yurchenko & Tennyson 2014), NH₃ (Yurchenko et al. 2011), HCN (Harris et al. 2006), H₂S (Azzam et al. 2016); HITEMP (Rothman et al. 2010): CO₂, CO; HITRAN (Rothman et al. 2013): C₂ H₂. Also, added are the resonance lines for Na, K as described in Heng et al. (2015) and Heng (2016) and CIA H₂ –H₂ H₂ –He absorption (Richard et al. 2012). Also, included is isotropic Rayleigh scattering of H₂ (Sneep & Ubachs 2005). With HELIOS we calculated a grid of cloud-free, self-luminous atmospheres with the following parameters: [Fe/H] = 0, 0.6; T_eff = [100, 1200], ΔT_eff = 50; log₁₀(g) (cgs) = [1.6, 4.0], Δ log₁₀(g) = 0.1.

4 Examples of evolutionary calculations

In this section examples of cooling curves of simple models for a Neptune-, Uranus-, Jupiter-, and Saturn-like planet, for a planet like GJ 436 b, and for a close-in, core-dominated, sub-Neptunian planet are presented. The results are compared with other thermal evolution models to validate our evolutionary model.

The thermal evolution of the planets is modelled with a three layer interior of the planet, namely an iron/silicate core, potentially an ice layer, and a (pure) H/He envelope, as was described in more detail in Sect. 2. This is similar to but simpler than the approach in Fortney & Nettelmann (2010), who include water mixed into the H/He layer, and H/He mixed into the water layer above the rock core. The assumption is that a mixed envelope composition is favoured from a planet formation point of view, because planetesimals might get dissolved in the envelope of the accreting planet (e.g. Podolak et al. 1988; Mordasini et al. 2006). In contrast to the calculations further down, a grey atmosphere is assumed here as boundary condition for the inward integration, as introduced in Sect. 3. The metallicity [M/H] enters the opacity calculation, using the Freedman et al. (2014) Rosseland mean opacity. Our aim is not to present detailed models for the evolution of the giant planets of the solar system for which numerous observational constraints exist. Rather, we want to understand how our simplified model that is used for exoplanets compares to existing more detailed simulations. For exoplanets we usually only have little observational constraint, for example a rough age and magnitudes in some bands. This makes a simplified approach appropriate, as more complexity would in any case remain unconstrained.

The evolution of the planet is calculated in the following way. First, with static interior structure calculations, the envelope-to-core mass ratio is determined for the planet today, given its measured mass, luminosity, and radius. Then, with the derived core-to-envelope mass ratio and a starting luminosity that is several orders of magnitudes higher than the one measured today, the planets’ temporal evolution is calculated. Since the Kelvin-Helmholtz timescale is short in the beginning, the exact starting luminosity respectively starting entropy is no longer important at the present time. The evolution is then calculated by taking into account the contraction of the envelope, as well as its cooling, and the contraction, cooling and radioactive energy production in the core, as explained in Sect. 2.

4.1 Neptune and Uranus

In Fig. 1 the bottom right panel shows the luminosity as a function of time for our simplified Neptune model. The planet has a total mass of 1 M_"???" = 17.147 M_⊕. The atmospheric opacity corresponds to a [M/H] = 1.8 (Guillot & Gautier 2014). An ice mass fraction of 50% in the core is assumed for all four giant planets in the solar system. A similar value is expected for a condensation of water ice in the solar nebula (Min et al. 2011; Lodders 2003) and is also motivated from planetary formation models (Guillot & Gautier 2014). Much higher ice mass fractions have sometimes been used in interior structure models, which is difficult to understand from a formation point of view (Guillot & Gautier 2014). For this fixed ice mass fraction, the core and envelope mass is determined with static interior structure calculations so that the planet has for its observed intrinsic luminosity L_"???" = 9.85 × 10⁻³ L_♃ (Guillot & Gautier 2007) a radius of 1 R_"???" = 3.87 R_⊕. We find a composition with a H/Heenvelope mass of 1.80 M_⊕. The central part, built of an ice layer wrapped around an iron/silicate core, thus has a mass of 15.347 M_⊕. An overview of these values is given in Table 1. This composition can be compared with the models of Podolak et al. (1995). Their Neptune model 1 (2) has a H/He envelope of 2.2 (0.9) M_⊕, and a total heavy element mass of 14.9 (16.2) M_⊕. These values bracket the ones found in our model. The bulk mass of heavy elements in Neptune found in Nettelmann et al. (2013) is 14− 14.5 M_⊕, which is also similar to our value.

We then followed the thermal evolution of this planet, starting with a high luminosity of 50 L_♃. The total luminosity is split in the contributions coming from the core cooling, core contraction, envelope cooling and contraction, as well as radiogenic heating. The total luminosity at the age of the solar system agrees well (difference of 3%) with the observed value, which is shown as a black dot. Our cooling curve for Neptune overlaps especially at later times with the one presented in Fortney et al. (2011). In Fig. 2, the change of the radius in time is shown, where the measured radius of today of the respective planet is given as a coloured dot. Neptune’s measured radius of today is by construction well reproduced (within 0.5%) by our simulations. The change in intrinsic luminosity over time of Neptune (bottom right panel in Fig. 1) follows a t^−4∕3 slope, which is analytically expected for the cooling of a planet where the dominant energy source is the thermal cooling of the core (Ginzburg et al. 2016).

We also simulated the cooling of an Uranus-like object. It is well known that Uranus is much fainter than expected from fully convective cooling models (for example Podolak et al. 1991, 1995; Fortney et al. 2011; Nettelmann et al. 2013). We therefore expect that it is not possible to find a correct cooling age with our model. This is indeed the case: for the simulation of a simple Uranus model, we assume that the planet has a total mass of 1 M_"???" = 14.536 M_⊕ with an opacity corresponding to [M/H] = 1.8 (Guillot & Gautier 2014). In order to reproduce the radius of Uranus today, R_"???" = 3.98 R_⊕ with the current upper limit of its intrinsic luminosity of L_"???" = 1.01 × 10⁻³ L_♃ (one order of magnitude less than Neptune, Guillot & Gautier 2014), we need an interior consisting of a H/He envelope of 2.470 M_⊕, and a solid part of 12.066 M_⊕, which is split into 50% ice and 50% silicate/iron (see Table 1). These values are similar to the Uranus model of Podolak et al. (1995) with 1.5 M_⊕ for the envelope and 13 M_⊕ for the ice and silicate/iron part, as well as to the Uranus model computed in Nettelmann et al. (2013), who find a bulk composition of heavy elements of 12.5 M_⊕.

With this interior composition of Uranus, we simulated the cooling of the planet, starting from a high initial luminosity (50 L_♃ as for Neptune, but the precise value is not important as long as it is high). Our evolution calculation of Uranus agrees at later times very well with the one in Fortney et al. (2011). We find that at the current age of the solar system the simulated planet has a luminosity that is eight times too high relative to the upper limit of the observations, which can also be seen in the bottom left panel of Fig. 1. We therefore recover the result (e.g. Fortney et al. 2011; Nettelmann et al. 2013) that standard fully convective cooling models fail to explain Uranus’s very low luminosity. This is also mirrored in the fact that the simulated radius is 5% too big compared to the measured radius of today (see Fig. 2). So even though we matched a static model of the planet with the present measured luminosity and radius of the planet, these may not be reproduced by the modelled evolution of the planet at 4.5 Gyr.

With the improved gravity field data from the Voyager fly-by of Uranus and Neptune and modified rotation period and shape of the planets, Nettelmann et al. (2013) compute adiabatic three-layer structures for the two planets and find that Uranus and Neptune might differ in their atmospheric enrichment within an observationally significant amount. This could be due to a stable stratification in the interior of Uranus and originate from a giant impact (Nettelmann et al. 2013). We conclude that Uranus and Neptune might have very different internal structures and/or thermodynamic states despite their similar masses and radii.

Fig. 1 Temporal evolution of the luminosity of the simplified solar system gas giants. The different contributions to the total luminosity are indicated with colours, the luminosity of today is shown as a black dot.

4.2 Jupiter and Saturn

Jupiter is simulated with a total mass of M_♃ = 317.83 M_⊕. We assumedagain an ice fraction of 50% in the core as for Uranus and Neptune and an opacity corresponding to [M/H] = 0.5 (Guillot & Gautier 2014). As before, we matched today’s given radius to the observed present-day luminosity by varying the core-envelope mass ratio. With static interior calculations, we obtained the radius of the planet today (R_♃ = 10.97 R_⊕) for a model with a central part containing iron, silicate, and ice of 27.50 M_⊕ and a H/He envelope mass of 290.33 M_⊕ (Table 1). Fortney & Nettelmann (2010) conclude that the current Jupiter models show a range in the core mass of 0− 18 M_⊕ and a heavy element mass in the envelope of 0−37 M_⊕. The more recent analysis of Wahl et al. (2017) based on Juno data finds core masses between 6 and 24 M_⊕, and total heavy element masses (core and metals mixed in the envelope) of 24–46 M_⊕, depending on assumptions concerning the core’s state and the EoS. Thus, keeping in mind that the core mass in our model rather represents the bulk heavy element in the planet, our results lie in a similar interval for the heavy element content of Jupiter.

We thenfollowed the evolution of this planet. The evolution of the luminosity can be seen in the top left panel of Fig. 1, and the evolution of the radius in Fig. 2. Our modelled Jupiter is slightly too bright compared to the measured luminosity (difference of 13%), as was already found for example by Fortney et al. (2011). For the Neptune model, the biggest contribution to the total luminosity came from the core cooling, as is expected for core-dominated planets (Baraffe et al. 2008). However, for our Jupiter model, the biggest contribution comes from the envelope contraction. Our modelled Jupiter is cooling too slowly compared to the real planet, reaching the observed luminosity at 4.91 Gyr. Therefore, also the modelled radius of the planet is slightly too big (1.0014 × R_♃) at the present age, which can be seen in Fig. 2.

We also modelled the evolution of Saturn. The planet was simulated with a mass of  1M_♄ = 95.159 M_⊕. As for the other solar system planets, we assumed an ice mass fraction of 50% in the core, but an opacity corresponding to [M/H] = 1 (Guillot & Gautier 2014). We find a solution for the static interior model with an envelope mass of 72.159 M_⊕ and a central part of 23 M_⊕ (Table 1). Following the evolution of the planet over time, we expect a cooling time for our homogeneous adiabatic models of 2–3 Gyr, as in Püstow et al. (2016). Our modelled Saturn reaches the luminosity of today at 3.1 Gyr, while at 4.6 Gyr the simulated Saturn is 48% too dim (note the logarithmic scale). The difference between today’s simulated luminosity and the measured one for Saturn is much larger than for Jupiter, by a factor of 3.7. Also the radius of the simulated planet is thus smaller than the measured radius. The modelled Saturn is therefore cooling too fast. Hence, we reproduce the common result in the literature that adiabatic homogeneous models underpredict Saturn’s current luminosity, or in other words that Saturn exhibits a strong excess luminosity (Stevenson & Salpeter 1977). This is commonly attributed to H/He demixing (e.g. Stevenson & Salpeter 1977; Fortney & Hubbard 2003; Püstow et al. 2016), but other explanations exist as well (Leconte & Chabrier 2013). We note that for a Saturnian-mass planet, the demixing sets in at an age of 1–2 Gyr (Püstow et al. 2016), so for objects younger than that it should not pose a big problem.

Fig. 2 Temporal evolution of the radius of the simplified solar system gas giants. The radius of today is given in the respective colour of the planet as a cross. On this scale, the symbols for the radii of Uranus and Neptune partially overlap.

4.3 GJ 436 b: comparison with Baraffe et al. (2008)

The transiting planet GJ 436 b (Gillon et al. 2007) orbits a 0.44 M_⊙ M-star at 0.028 au. The mass of the planet is M_GJ436b = 22.6 ± 1.9 M_⊕ and its radius is determined to be 0.386 ± 0.016 R_♃ (Gillon et al. 2007; Deming et al. 2007). The age of the system is unconstrained by observations. Baraffe et al. (2008) assume a system age of 1–5 Gyr and do not take irradiation on the planet into account because of the low luminosity of the parent star. They find a good match with the observed radius within the uncertainty of the age of the system for a 22.6 M_⊕ planet with a 21 M_⊕ water core surrounded by a 1.6 M_⊕ H/He envelope. We want to test whether we can match the observations of this planet with the same composition which is what we find.

Due to its core-dominated nature, GJ 436 b is a good example to study the contributions of the core and envelope to the planet’s temporal evolution. The contributions to the gravothermal energy release for the water core look qualitatively similar if the structure of the planet is calculated with the EoS SESAME or ANEOS (Baraffe et al. 2005, Baraffe 2015, priv. comm.). Figure 3 shows the luminosity from the core cooling and contraction as a sum and separate relative to the total luminosity L over time, as well as the contribution from the envelope relative to the total L. Shown are the luminosities predicted by our model assuming a water core together with the results from Baraffe (2015, priv. comm.).

In constrast to Fig. 10 in Baraffe et al. (2008), which used SESAME, the black and blue lines in Fig. 3 were obtained using the ANEOS EoS on which we base our expressions for the heat capacities and the thermal expansion coefficient as well. Looking in more detail at the calculation employing ANEOS, which is shown in Fig. 3 (blue and black lines), the biggest contribution to the gravothermal energy release comes from the core with a fraction of ~ 0.85 from the total luminosity averaged over time. The total contribution from the core is split into a bigger part coming from the thermal cooling of the core with a fraction of ~ 0.55 between 10 and 100 Myr and then rising to reach ~0.8 around 1 Gyr. The smaller part of the total core contribution comes from the contraction of the core with a fraction of ~ 0.3 from 10 to 100 Myr and then sinking to ~0.1 at 1 Gyr. The envelope contribution to the total luminosity of the planet is the smallest with a fraction of ~ 0.15.

The relative contributions do not strictly sum up to 1, as is expected. This is due to the fact that the total luminosity is calculated using the entropy that is given by the EoS, L_tot = −∫ T(dS∕dt) dm. If the applied EOS were thermodynamically coherent, one would have −TdS∕dt = − PdV∕dt −dU∕dt. However, this is not exactly the case with ANEOS (from Baraffe 2015, priv. comm.). Therefore, the sum of the volume work and the internal energy terms, which are shown separately in Fig. 3, do not sum up to the total luminosity obtained from the entropy change.

The gravothermal energy contributions from this work are similar (shown in Fig. 3 as well, olive and orange lines), namely a fraction of ~ 0.9 for the biggest contribution to the total energy from the core which is split into a contribution of ~0.73 from the cooling of the core and ~0.17 from the contraction of the core. The smallest contribution to the total luminosity originates again from the envelope (fraction of ~ 0.1). Different to the calculation by Baraffe 2015, the energy contributions sum up to 1.

In this calculation, we aim to reproduce the model of Baraffe et al. (2008) and use the same model parameters. Accordingly, the core is composed of pure water and thus there is no radiogenic contribution to the luminosity, as the radioactive material content of the planet is assumed to be proportional to the rock mass fraction of the planet. Simulating a planet with the same core and envelope mass fraction as above, but containing 0.8 M_⊕ of iron and rock in the core (corresponding to an ice mass fraction in the core of 50%), so containing 0.53 M_⊕ rock, the radiogenic luminosity contributions can be estimated. For example, it accounts for 0.01% (0.5%) at 1 (100) Myr of the total luminosity.

As we discussed above, the overall agreement between the two simulations regarding the luminosity contributions in the planet is good, as expected. The agreement especially concerns the fact that the core makes the dominant energy contribution to the total luminosity budget of the planet. Therefore, neglecting the core’s contribution in the evolution calculation of such a planet has asignificant impact.

Fig. 3 Relative contributions to the total luminosity L over time for the water core simulation of GJ 436 b. The blue and black lines are from Baraffe (2015; priv. comm.), the red and brown lines are from this work and represent the following: solid blue/ brown line: energy release from the core, which is the sum of core cooling and core contraction; dashed blue/ brown line: core cooling; dash-dotted blue/ brown line: core contraction; dashed black/ red line: envelope cooling and contraction. A discussion can be found in Sect. 4.3.

Fig. 4 Luminosity as a function of time for a 5 M⊕ planet with a 1% H/He envelope at a distance of 0.1 au from a 1 M⊙ star. The core is assumed to be isothermal. The different contributions are indicated in the plot. The solid black line shows the total luminosity as found by Lopez & Fortney (2014).

4.4 A 5 M_⊙ planet with a 1% H/He envelope: comparison with Lopez & Fortney (2014)

As an example of the evolution of a strongly core-dominated planet, and for comparison with another independent evolutionary model, we present in Fig. 4 the cooling curve of a planet with the same properties as the one simulated in Fig. 3 of Lopez & Fortney (2014). It is a 5 M_⊕ planet with a rocky core and a 1% H/He envelope. The planet is located at 0.1 au from a solar-like star. Such close-in sub-Neptunian planets have been found in high numbers by the Kepler satellite (e.g. Fressin et al. 2013; Petigura et al. 2013). The envelope opacity corresponds to a 50-times-solar heavy element enrichment ([M/H] = 1.7). In contrast to our simpler model, Lopez & Fortney (2014) directly use the ANEOS (Thompson 1990) and SESAME (Lyon et al. 1992)equations of state to model the interior of the solid core, and they employ fully non-grey atmospheric models.

Different from all other simulations in the present paper, we assume for this comparison that (i) the core is isothermal, (ii) that the heat capacity of the core is 7.5 × 10⁶ erg K⁻¹ g⁻¹, (iii) that the core does not shrink due to the reduction of its temperature, which means that the release of internal energy (d u∕dt), but not the associated release of gravitational potential energy (pdV∕dt) of the core, contributes to its luminosity. All these three settings mimic the ones made by Lopez & Fortney (2014). The adiabatic model for the core (Eq. (10)) leads in this case to total luminosities that are about 66% higher compared to the isothermal core. This comes from the fact that the thermal energy of the core is more than two times higher for the adiabatic core model than for the isothermal core model. For planets with a less extreme core-envelope mass ratio, the impact is correspondingly smaller, as discussed in Sect. 2. For giant planets, the difference is completely negligible.

A comparison of the two simulations in Fig. 4 shows that the overall agreement is very good, with a tendency of our model to predict a somewhat dimmer planet at earlier times (≤100 Myr). After 300 Myr, the luminosity found in our model is slightly higher than in Lopez & Fortney (2014), with a difference ofabout 9% at 10 Gyr. A feature that is shared in both models is the dominance of the radiogenic heating mainly caused by 40 K decay (Mordasini et al. 2012b) over the cooling of the core at intermediately late times (0.7–5 Gyr). Further comparison shows that our model predicts a somewhat higher contribution of the core at later times, whereas it is the opposite at earlier times. The radiogenic luminosity is virtually identical in the two models.

5 Initial conditions

For the simulation of the thermodynamic evolution of the planets and the calculation of their magnitudes we need to specify the post formation properties of these planets. As exemplified by the large impact of cold versus hot starts for giant planets (e.g. Marley et al. 2007; Mordasini et al. 2017), the post formation properties can have a very strong influence on the predicted detectability and observational characteristics. To find such initial conditions for the simulations, the output of planetary population syntheses based on the core accretion paradigm (see e.g. Mordasini et al. 2017) was studied. Population synthesis is the attempt to model planet formation globally. Therefore, a model containing many sub-models coming from specialized studies on one aspect of planet formation (for example gas and dust disc dynamics, type I & II migration, ice-line behaviour, planetary accretion) is constructed. By simplifying and putting these specialized models together, a global planet formation model can be built. The initial conditions, such as the total disc (gas) mass, the dust-to-gas-ratio, and the lifetime of the disc are sampled in a Monte Carlo way from probability distributions that are derived as closely as possible from observations. In the particular population synthesis used for this work, ten planetary embryo were inserted into a disc around a solar-like star. More model settings used in the population synthesis are described in Mordasini et al. (2012b, 2017). Formally, the synthesis was conducted under the assumption of cold gas accretion. However, for low-mass planets studied here, there is no hot versus cold accretion difference in the same way as for giant planets. For giant planets, the dominant fraction of their mass is accreted after detachment from the protoplanetary nebula (at ~ 100 M_⊕) through a potentially entropy-reducing shock. Low-mass planets only detach from the nebula when the nebula itself has already almost completely dissipated. Because of this, only a very small amount of gas is accreted after detachment through a potentially entropy-reducing shock. Different formation histories of the individual planets (e.g. moment when gas and solids are accreted, surrounding disc conditions), however, still induce a diversity in post formation properties. Additional physical mechanisms such as solid accretion at late times, giant impacts, enriched envelope composition, or semi-convection might change the outcome of our planet formation model, but are currently not considered.

From the planet population synthesis, the planet’s core and H/He envelope mass fractions as well as the luminosity at the end of formationwas estimated. The output from the population synthesis was studied at an age of 3 Myr, the typical life-time of the synthetic discs. We do not include planets that are still undergoing strong planetesimal accretion. For the case of a hot protoplanet after collisional afterglows, the reader is referred to for example Schaefer & Fegley (2009) and Miller-Ricci et al. (2009).

The population synthesis results for envelope mass as a function of total planet masses from 5–30 M_⊕ are shown in the left panel of Fig. 5 as blue circles. The data were fit by eye, which is given by Eq. (11): $$\frac{M_{\text{enve}}}{M_{\oplus }}=3{\left(\frac{M_{\text{tot}}}{10M_{\oplus }}\right)}^{1.8},\text{for}\${\text{M}}_{\text{tot}}\$\text{textless}\$ \text{30\hspace{0.17em}}{\text{M}}_{\oplus }\$,$$(11)

and shown as an orange-red line in Fig. 5, together with the three lowest mass planets simulated in this workindicated as orange dots. We note that this fit relation represents mainly planets with a high M_enve for a given total mass (see Fig. 5), and that studies of the composition of close-in, low-mass sub-Neptunes indicate rather lowerenvelope mass fractions (e.g. Wolfgang & Lopez 2015). However, the envelope mass of many of these planets was potentially reduced by atmospheric escape (e.g. Owen & Wu 2013; Jin & Mordasini 2018), and the planets we consider here are at larger semi-major axes.

For heavier planets, the gas accretion rate changes from being limited by the planet’s Kelvin-Helmholtz contraction to the disc-limited regime (e.g. D’Angelo et al. 2011). Thorngren et al. (2016) derive from the mass-radius relation of observed planets a fit for the total heavy element content as a function of the total planetary masses, which is given in Eq. (12): $$\frac{M_{\text{z}}}{M_{\oplus }}=(57.9\pm 7.03){\left(\frac{M_{\text{tot}}}{318M_{\oplus }}\right)}^{0.61\pm 0.08},\text{for}\${\text{M}}_{\text{tot}}\$\text{textgreater}\$ \text{30\hspace{0.17em}}{\text{M}}_{\oplus }.\$$$(12)

For planets that are more massive than 30 M_⊕, their fit was used. In the left panel of Fig. 5 for reference the envelope mass of a 30 M_⊕ is given (green square) as computed with their expression. A satisfyingly smooth transition between the two fits is found.

From the population syntheses, also the post formation luminosity (L_pf) and thus entropy (S_pf) was obtained. The L_pf are shown as blue circles in the right panel of Fig. 5. For planets between 5 and 30 M_⊕, the population syntheses output was fit by eye, which is given in Eq. (13): $$\frac{L}{L_{\text{mboxmohinichar88}}}=42{\left(\frac{M_{\text{tot}}}{10M_{\oplus }}\right)}^{2.1},\text{for}\${\text{M}}_{\text{tot}}\$\text{textless}\$ \text{40\hspace{0.17em}}{\text{M}}_{\oplus },\$$$(13)

and shown as an orange-red line in the right panel of Fig. 5.

For giant planets, L_pf increases more slowly with M_tot than for M_tot ≤ 100 M_⊕ because giant planets go through an entropy-reducing shock, therefore the slope becomes shallower. In Mordasini et al. (2017) a fit for the luminosity of giant planets (M ≥ 100 M_⊕) depending on the total planetary mass was provided. In this work, the fit for the cold-nominal planets is applied and given also in Eq. (14): $$\frac{L}{L_{\text{mboxmohinichar88}}}=1.378\times {10}^4{\left(\frac{M_{\text{tot}}}{318M_{\oplus }}\right)}^{1.3},\text{for}\${\text{M}}_{\text{tot}}\ge \text{\hspace{0.17em}50\hspace{0.17em}}{\text{M}}_{\oplus }.\$$$(14)

This fit is shown as a dark-green line in the right panel of Fig. 5 and used for the heavier planets in our dataset (orange dots in the same figure).

Table 2 finally gives an overview of the properties of the eight planets the evolution of which we simulate in this study, spanning a wide mass range from 5 M_⊕ to 2 M_♃. The post formation entropy S_pf is calculated at the bottom of the convective zone with a grey atmosphere after the bulk composition and luminosity of the planet is given.

Fig. 5 Results from population synthesis calculations (blue circles) of combined planet formation and evolution that are used for the initial conditions. Orange dots give our final dataset. Left panel: envelope mass as function of total mass with the fit from this work (orange-red, Eq. (11)). For more massive planets the envelope mass as calculated from the relationship found by Thorngren et al. (2016; Eq. (12)) is shown for the 30 M⊕ planet with a green square. Right panel: post formation luminosity as a function of total planetary mass, together with the fit from Mordasini et al. (2017; dark-green, Eq. (14)) and from this work (orange-red, Eq. (13)).

6 Results and discussion

In this section, we show cooling curves resulting from our evolution code using the initial conditions given in Table 2 for different planetary atmosphere models and metallicites. Following this, the impact of different post formation entropies on the planet’s evolution is studied. Finally, we calculate magnitudes for various typical filters that can be found for instance in the Nasmyth Adaptive Optics System Near-Infrared Imager and Spectrograph (NACO) on the Very Large Telescope (VLT), in the VLT Imagerand Spectrometer for mid Infrared (VISIR), in the Polarimetric High-contrast Exoplanet Research (SPHERE), or on JWST. The influences on JWST magnitudes from different atmospheric models or post formation entropies are estimated. Also, we carry out a comparison with the JWST sensitivity limits for the simulated planets during their evolution.

For all planets, the ice mass fraction in the core was set to 0.5 as direct imaging is sensitive to planets at large orbital distances. Because of this, stellar irradiation was neglected. It is important to note that in the simulations of this work time zero is when the gas disc disappears. The stellar age could thus be up to ~ 10 Myr higher, depending on the specific disc lifetime.

6.1 Cloud-free solar metallicity models

At first, the evolution with cloud-free solar metallicity model atmospheres was studied. Figure 6 shows the evolution for eight planetary masses from 5 M_⊕ to 636 M_⊕ (2 M_♃) and four different atmosphere types, from 0.1 Myr to 10 Gyr for the heavier planets. The luminosity of Jupiter and Saturn is given for reference as squares in the colour corresponding (roughly, for Saturn) to their masses. Planets with M_tot < 20M_⊕ are not plotted with an AMES-Cond atmosphere, as these masses evolve outside of the atmospheric grid (see Fig. 7). For the AMES-Cond, petitCODE, and HELIOS atmospheres, the tracks are shown as long as the planet is evolving on the atmospheric grid. The dots represent the luminosities, temperatures, and radii from Baraffe et al. (2003).

The differences between Baraffe et al. (2003) and our cooling curves at early times in the luminosity and temperature panel simply reflect the choice of different initial luminosities, and at later times, the two different cooling calculations agree well. For example, at 50 Myr the difference for the 636 M_⊕ planet between the luminosity curve with an AMES-Cond atmosphere versus a HELIOS (petitCODE) atmosphere is 14% (29%). The mean difference is around 25%. The mean was obtained by averaging over the maximal procentual differences in bolometric luminosity for the different atmospheric models at 1, 10, and 100 Myr for the four heaviest planets. These differences are smaller than the error bars today for measured luminosities (e.g. 48% for Eri b in Macintosh et al. 2015). With future more precise measurements, it could be possible to distinguish between different atmospheric models. We conclude that the choice of the atmosphere has, most of the time, only a limited impact on the evolution of the bolometric luminosity, as also found for example in Burrows & Liebert (1993), Burrows et al. (2001), Baraffe et al. (2002), Saumon & Marley (2008), and Dupuy et al. (2015). A posteriori, it is therefore justified to use a grey atmosphere for the comparison calculations in Sect. 4.

In the limit of core-dominated planets with very low M_enve, where the luminosity is dominated by the core cooling, and in the approximation of a constant radius, the analytical model of Ginzburg et al. (2016) predicts L ∝t^−4∕3 for a constant mass. However, the M_enve even of the 5 M_⊕ planet in this work is sufficiently massive so that the radius change cannot be neglected, and numerically a decrease rather like t⁻¹ is found for a fixed mass, even for the low-mass planets. This is similar to what was found in Burrows & Liebert (1993), and is valid for grey and non-grey atmospheres. This can be understood from the fact that at the radiative-convective boundary (RCB), the bottleneck for the transport of the luminosity occurs, independent of the atmospheric model, at the high optical depth of the RCB (Arras & Bildsten 2006; Lee et al. 2018). Radiative transport occurs there by diffusion, meaning that only the Rosseland mean opacity matters, and not the specific wavelength-dependent opacity of the different atmospheric models.

In the top right panel in Fig. 6, the decline of the specific entropy in the inner convective zone with time is shown. The entropy is a good measure of the total gravothermal energy of a planet because it contains both the inner energy content and the gravitational energy (volume work). The entropy contained in our simulated planets ranges from 7.6 to 9.5 k_B baryon⁻¹ at young ages and from 7.0 to 8.7 k_B baryon⁻¹ at 10 Myr. We have fixed the post formation luminosity. This means that the post formation entropy varies for a given mass depending on the atmospheric model (Marleau & Cumming 2014).

We also studied the effective temperature. The effective temperature varies slowly until about one Kelvin–Helmholtz timescale (τ_KH) has passed. The initial τ_KH = E∕L ranges from 8.96 to 9.62 Myr for the 318 M_⊕ to 13.50 to 14.52 Myr for the 5 M_⊕ planet, depending on the atmosphere model. For example, the effective temperature of the 100 M_⊕ planet decreases from 565 to 603 K at 0.11 Myr to 107 to 116 K at 1 Gyr depending on the atmospheric model. The difference between the temperatures from this work and from Baraffe et al. (2003) correspond to the differences in luminosity that were discussed before.

The bottom right panel in Fig. 6 shows the temporal change of the radius. It corresponds to τ = 2∕3 for the Eddington grey atmosphere and to the coupling pressure (50 bar) for the other atmosphere grids. The radii from Baraffe et al. (2003) are constantly bigger than the radii from our simulation of the planets with an AMES-Cond atmospheric grid. When simulating the planets with a luminosity as was chosen in Baraffe et al. (2003), and assuming a core mass of 0.5 M_⊕ to mimic the settings in Baraffe et al. (2003), we find a much better agreement in the radii. For example, at 1 Myr the radius from this work would then be 0.4% bigger (1.9% smaller) than the Baraffe et al. (2008) radius for the 1 (2) M_♃ planet. We see that initially there is a non-monotonic relationship between mass and radius, with the biggest radius occurring for the 20 M_⊕ planet with the AMES-Cond atmosphere. This is a consequence of the initial conditions (L_pf and M_enve versus M_core), which predicts for the 20 M_⊕ planet a rather high envelope mass fraction relative to the total mass of the planet. At later times, the radii computed with different atmosphere models converge, the radius also increases monotonically with mass, even though there is a spread of up to 0.9 R_♃ (for the 20 M_⊕) at early times.

It is important to note here that the AMES-Cond grid was not designed for such low-mass planets. To give an overview, Fig. 7 illustrates the cooling of the planets in the log g–T_eff plane. The top panel shows the cooling in the AMES-Cond grid. Thick lines indicate that the evolution path is inside the AMES-Cond grid, which starts at log g = 2.5 and an effective temperature of 100 K. As can be seen, the 636 M_⊕ planet is evolving on the AMES-Cond grid, whereas all the others fall partially off the grid due to either too small surface gravities at early times or too low temperatures at later times. This means that all AMES-Cond results for planets below 50 M_⊕ must be taken with caution. For the petitCODE and the HELIOS grid, the coverage in log g–T_eff is in contrast good, initially, when T_eff is larger than 150 K (petitCODE) or 100 K (HELIOS). These limits are shown in the T_eff panel. Again, once the planets are below these temperatures, caution must be used when employing the cooling tracks.

Fig. 6 Temporal evolution of fundamental properties of the eight simulated planets with initial conditions given in Table 2. The change in intrinsic luminosity, entropy, temperature, and radius is shown. The colour code for the masses is given in the bottom left panel. The luminosities of Jupiter and Saturn are indicated as squares in the top left panel. The dots in the luminosity, effective temperature, and radius panel show the results from Baraffe et al. (2003). Solid lines show the evolution with a cloud-free atmosphere with solar metallicity from the petitCODE grid, dashed lines show it with an AMES-Cond atmosphere, dash-dotted lines show it with a solar metallicity, cloud-free atmosphere from the HELIOS grid, and dotted lines finally show the evolution calculated with a grey atmosphere.Since the three smallest masses are evolving outside of the AMES-Cond grid (see Fig. 7), these masses are not plotted for an AMES-Cond atmosphere. For the AMES-Cond, petitCODE and HELIOS atmospheres, the tracks are shown as long as the planet is evolving on the atmospheric grid. For better visibility, three resp. two evolution curves are not shown in the entropy resp. radius plot. The horizontal lines at 100 and 150 K in the Teff plot show the lower limit of the HELIOS and petitCODE grid, respectively.

Fig. 7 Evolution of the planets in the logg–Teff space, together with the coverage by the atmosphere grids given as rectangles. The colour code for the masses is given in the middle panel. Thick lines indicate cooling curves that are in the atmosphere grid and thin lines those that are outside of the grid. Top panel: evolution in the AMES-Cond grid. Middle panel: evolution in the petitCODE grid for the cloudy (fsed = 0.5) and solar metallicity atmosphere is shown. This is representative also for other types of petitCODE atmospheres. Bottom panel: evolution in the HELIOS grid for a cloud-free atmosphere with a [M/H] = 0.6 and is representative also for the solar -metallicity evolution paths.

6.2 Non-solar or cloudy atmospheric models

In contrast with the former section, we now study the influence of non-solar metallicities and clouds on the evolution of the planets. In Fig. 8 in the top left panel, the evolution of the intrinsic luminosity from 0.1 Myr to 10 Gyr for planets from 636 M_⊕ (2 M_♃) down to 5 M_⊕ is shown; again the lower mass planets are only shown as long as they are evolving not too far from the boundary of the atmospheric grid. The colour code is given in the top left panel. These evolution calculations were done with the petitCODE and HELIOS grids described in Sect. 3. As a baseline model the evolution of the planets with a cloud-free solar metallicity atmosphere in the petitCODE grid is shown with a solid line; these are the same lines as in Fig. 6. We note that a f_sed = 0.5 is probably not realistic for the coldest planets studied here (T_int ≲ 300−400 K), because the cloud species implemented in our simulations (Na₂S and KCl) form deep in the atmosphere in these cases. A low f_sed potentially mixes the clouds to locations too high up in the atmosphere. However, the evolution with a f_sed = 3.0 was also calculated in the petitCODE grid. Since the difference to the cloud-free evolution is not visible on the scale chosen here (the cloud-free luminosity is less than 1% fainter at 5 Myr than the luminosity calculated with an atmosphere with f_sed = 3.0), it is not shown. This is not surprising, as the high f_sed corresponds to a strong cloud settling or weak cloudiness. As mentioned before, at such low temperatures a more realistic choice for the cloud description would have been to also include the effect of water clouds. Finally, the evolution in the cloud-free atmosphere with a metallicity of [M/H] = 0.6 in the HELIOS grid is given as a dash-dotted line.

In the top left panel of Fig. 8, the change of the bolometric luminosity over time is shown. For example, at 16 Myr the cloud-free HELIOS atmosphere with [M/H] = 0.6 predicts a 1.11 times brighter luminosity than the clear solar metallicity petitCODE atmosphere, whereas the cloudy petitCODE atmosphere predicts a 0.15 times fainter luminosity. Analysing the figure inmore detail, the evolution of the bolometric luminosity of the planets shows average deviations (calculated as explained inSect. 6.1) of about 6% for metallicities of 0.6 instead of 0.0. This is smaller than the deviations seen in Sect. 6.2 where we compared different clear atmospheric models with solar composition, and thus smaller than the error bars in current luminosity measurements. We thus see a non-negligible, but also not very large effect of the atmospheric model on the bolometric luminosity of the planet. We remind the reader, however, that Fig. 8 only shows a maximum enrichment of [M/H] = 0.6. For the maximum enrichment that we simulated ([M/H] = 1.2), the deviation relative to the solar case is around 12% on average (for the clear case). For f_sed = 0.5, the differences to the clear solar metallicity case are similar as those between clear solar metallicity and the [M/H] = 0.6 case. We caution thatwe have not investigated systematically the consequences of all cloud parameters for the bolometric luminosity. Again the result from the literature (e.g. Burrows & Liebert 1993; Burrows et al. 2001; Baraffe et al. 2002; Saumon & Marley 2008; Dupuy et al. 2015) is reproduced, as in the former section, that the choice of the atmospheric model has only a limited impact on the evolution of the bolometric luminosity, at least within the models studied here.

The change of the entropy over time is given in the top right panel in Fig. 8. For clarity the evolution of the 318, 159, and 50 M_⊕ mass planet is not shown. It is again important to keep in mind that the simulations start always with the same luminosity, which corresponds to a different entropy at the beginning depending on the atmospheric model. The biggest postformation entropy spread occurs for the three lightest planets with a difference of up to 0.4 k_B baryon⁻¹. The discrepancy between the atmospheric models can be attributed to a large extent to differences in the employed opacity sources. Since petitCODE and HELIOS use different line lists for some of the absorbing species, such as H₂ O, NH₃, H₂ S, and the alkali metals, and use different scatterers, their calculations may result in somewhat deviating atmospheric temperatureprofiles. This directly impacts the coupling to the deep convective adiabat and the corresponding entropy value. With increasing metallicity this discrepancy is exacerbated as the relative amount of absorbing gases increases.

The bottom left panel shows the evolution of the temperature of the planets and mirrors the luminosity evolution in the top left panel. As an example, the 100 M_⊕ planet has a temperature of 585–605 K depending on the atmosphere model at 0.1 Myr and cools down to 107–125 K at 1 Gyr, respectively.

Finally, the change of the radius over time can be seen in the bottom right panel. The spread in radius at early times can be up to 0.2 R_♃ for the 10 M_⊕ mass planet, up to 0.3 R_♃ for the 20 M_⊕ mass planet, and up to 0.1 R_♃ for the 100 M_⊕ planet. However, they all approximately converge at later times, qualitatively similar to cloud-free, solar metallicity atmospheres.

6.3 Varying the post formation luminosity

Different post formation luminosity or entropy (S_pf or L_pf) can be used to represent different planet formation scenarios (Marley et al. 2007; Mordasini et al. 2017). The post formation Kelvin–Helmholtz timescale is given as $${\tau }_{\text{KH},\text{pf}}\approx \frac{G{M}_{\text{tot}}^2}{R_{\text{pf}}{L}_{\text{pf}}},$$(15)

where G is the gravitational constant, M_tot is the total planetary mass, and R_pf and L_pf are the post formation radius and luminosity of the planet. At late times, much longer than τ_KH,pf, the influence of the initial condition has disappeared, therefore the choice of the post formation entropy (S_pf) is considered to have only a minor influence on the late planetary evolution. However, since τ_KH,pf ∝ 1∕L_pf, a low post formation luminosity (a very cold start) can influence the L(t) over extended and observationally relevant times of up to ~1 Gyr for giant planets (Marley et al. 2007). In contrast, a planet with a “hotter start” will converge on the same evolutionary track as a “hot start”, again because τ_KH,pf ∝ 1∕L_pf. Not surprisingly, very low S_pf (implying low L_pf) can have a strong impact on the predicted magnitudes of giant planets (Fortney et al. 2008; Spiegel & Burrows 2012), with magnitudes that are 1–6 mag fainter than standard hot start models. The very low post formation luminosities found in the original core accretion model of Marley et al. (2007) were originally thought to be a diagnostic distinguishing core accretion from other formation pathways (disc instability, turbulent fragmentation of molecular clouds). However, later core accretion models showed that core accretion can also lead to warm and even hot starts depending on the mass of the planet’s core (Mordasini 2013; Marleau et al. 2017; Berardo et al. 2017). This means that the brightness is likely not a property distinguishing strictly the formation modes, and recent population syntheses based on the core accretion paradigm (Mordasini et al. 2017) instead find warm starts. Therefore, because we want to make predictions about the detectability of young planets, it is important to quantify the impact of the S_pf and also to consider for how long the impact on the evolution of the planets remains. To study the impact of different S_pf, the cooling of planets with a L_pf ten times brighter (hot scenario) or fainter (cold scenario) relative to what is noted in Table 2 (nominal) was simulated. Such a spread in luminosity is suggested by population synthesis calculation of planet formation and can be seen in Fig. 5. The resulting spread in S_pf reaches from 8.4 to 10.5 k_B baryon⁻¹ (difference of 2.1 k_B baryon⁻¹) for the 318 M_⊕ planet, and from 6.9 to 8.5 k_B baryon⁻¹ for the 5 M_⊕ planet.

Figure 9 shows the cooling curves for four different planetary masses (5, 20, 100, and 318 M_⊕) and various L_pf. The evolution is calculated with a cloud-free atmosphere with solar metallicity in the petitCODE grid. The solid lines represent the nominal scenario and are therefore the same as in Figs. 8 and 6. The cooling curves with a ten times higher L_pf are shown as dashed lines, those with a ten times lower L_pf are given as dash-dotted lines.

The left panel in Fig. 9 shows the evolution of the intrinsic luminosity. It can be seen that the cooling curves for different L_pf converge as expected (e.g. Marley et al. 2007; Marleau & Cumming 2014). The planets of the hot scenario have a Kelvin–Helmholtz-time (τ_KH,pf) of ~1 Myr, the ones with nominal L_pf as shown in Table 2 have a τ_KH,pf of ~10 Myr and finally, the planets in the cold scenario have a τ_KH,pf of ~100 Myr. This means that the post formation state (e.g. the initial condition) influences the planets’ properties during about τ_KH,pf ~1, 10, and 100 Myr after formation for cold, nominal, and hot scenario. These phases can be seen in Fig. 9 as those parts of the lines that do not yet follow a L ∝ t⁻¹ behaviour, but which are more horizontal. As an example, at 1 Myr, the hot scenario model for the 100 M_⊕ planet is 1.6 times brighter than the nominal evolution, and the cold scenario model is five times fainter than the nominal model.

In Marley et al. (2007), L_pf is independent of the planetary mass because a larger mass fraction for the more massive planets went through an entropy-reducing shock, and therefore ${\tau }_{\text{KH,pf}}\propto M_{\text{tot}}^2$. Here we find that τ_KH,pf is approximately independentof M_tot (at least for M_tot < 100 M_⊕), as circa $L_{\text{pf}}\propto M_{\text{tot}}^2$, and the planetary radius, which enters linearly into τ_KH,pf, only changes by factors of 2–3 between hot and cold scenarios.

We note that because of the non-uniqueness of the M− L_pf relation, there is a mass-luminosity dependency during ~τ_KH,pf, as already widely discussed in the literature in the context of giant planets (e.g. Mordasini 2013). For example, if the age of a planet could be estimated to be 1 Myr and the planet has a brightness of 2 × 10⁻⁷L_⊙, then this could correspond to either a 20 M_⊕ mass planet with a high post formation luminosity, or to a 100 M_⊕ mass planet with a low post formation luminosity. With a further mass measurement through for example the astrometric or radial-velocity method, it is possible to disentangle this mass-luminosity ambiguity. In principle, if the planetary temperatures can be measured precisely enough, it would also be possible to reduce the degeneracy. This should be possible for differentiating between the hot and cold scenario for planets that are heavier than 2 M_♃ (Samland et al. 2017). Fortunately, at 10 Myr (a more likely observable age), none of the hot versus cold scenario lines representing 5, 20, 100, and 318 M_⊕ from Fig. 9 overlap any more, so that at least a very rough mass estimate should be possible.

The right panel in Fig. 9 shows the evolution of the temperature over time. At 1 Myr the temperature for the 318 M_⊕ planet in the hot, nominal, and cold scenario is 1014, 781, and 536 K, respectively. For comparison, Spiegel & Burrows (2012) have at 1 Myr a T_eff of ~850 and ~550 K for their hot and cold case for a 318 M_⊕ planet, which is similar to our result.

If planet masses and luminosities could be determined observationally at early times (< τ_KH,pf), important constraints about hot versus cold scenarios could be made for planet formation theory, regarding for example the heating by giant impacts (e.g. Anic et al. 2007), the efficiency of heat transport in the interior (Nettelmann et al. 2013), or the magnitude and timing of heating by planetesimal accretion Mordasini (2013).

Fig. 8 Temporal evolution of fundamental properties of the simulated planets with different metallicities and cloudy atmospheres. The colours represent different masses, where the colour code is given in the bottom left panel. The two squares in the top left panel indicate the luminosity of Jupiter and Saturn. Solid lines correspond to a cloud-free atmosphere with solar metallicity, dashed lines correspond to an atmosphere with clouds (fsed = 0.5) and solar metallicity (as in Fig. 6), dotted lines correspond to a cloud-free atmosphere with a metallicity [M/H] = 0.6, all of them using the petitCODE model. Dash-dotted lines indicate a cloud-free atmosphere with a metallicity [M/H] = 0.6 from the HELIOS grid. For better visibility, three resp. two planets are not shown in the entropy resp. radius plot. The horizontal lines at 100 and 150 K in the Teff plot show the lower limit of the HELIOS and petitCODE grid, respectively. The tracks are shown as long as the planet is evolving on the atmospheric grid.

Fig. 9 Influence of varied post formation luminosity (Lpf) on the intrinsic luminosity and temperature for four different masses. The colour code is given in the left panel, the line style is given in the right panel. All cooling curves are calculated with a cloud-free atmosphere with solar metallicity in the petitCODE grid. Solid lines show the evolution with the nominal Lpf of Table 2. Dashed lines show the evolution with a ten times bigger Lpf (hot scenario). Dash-dotted lines show the evolution with a ten times fainter Lpf (cold scenario).

6.4 JWST magnitudes and fluxes

Magnitudes for the most relevant JWST filters for exoplanet imaging, namely, JWST/MIRI (F560W, F770W, F1000W, F1280W, F1500W, F1800W, F2100W, F2550W) and JWST/ NIRCam (F115W, F150W, F200W, F277W, F356W, F444W) were calculated. The wide (W) filters were chosen as they are considered the standard imaging filters and also because we quantify the detection of planets in the background limited regime and not necessarily in the contrast limit, where certain coronagraph-filter combinations may be the preferred choice. The zero points of the magnitudes were obtained with a Vega spectrum ². To check our magnitude calculation, the spectra from F. Allard’s website were downloaded, convolved with the filter profiles, and the resulting magnitudes were compared with the magnitudes given on F. Allard’s website ³ as well. It was found that the magnitudes are the same. The magnitudes of the modelled planets were obtained by convolving their spectra with the filter transmission profiles. For this we interpolated their spectra to the desired log g–T_eff combinationgiven the grid of atmospheres. The filter profiles used in this work are available together with the magnitudes for different planetary masses and ages, metallicities, and post formation luminosities for clear or cloudy atmospheres. To illuminate the effects of clouds, a f_sed = 1.0 was chosen. We are considering the most important cloud species that occur at intermediate temperatures (Na₂S and KCl, for T ≳ 400 K). Consequently, we show the evolution tracks with clouds only down to 200 K (instead of 150 K as is the case for the clear petitCODE atmosphere magnitudes). An example of the tables that can be found at the CDS is given in the Table A.2 for a fixed planetary mass as a function of time, and in Table A.3 for a fixed time for the eight masses used in this work; all the tables can be found at the CDS and at http://www. space.unibe.ch/research/research_groups/planets_ in_time/numerical_data/index_eng.html.

As an illustration, in Figs. 10 and 11, the absolute magnitudes in the F356W filter (centred at 3.45 μm) and the F444W filter (centred at 4.44 μm) are shown for the 5, 20,100, and 318 M_⊕ mass planets and for different atmosphere grids with different parameters, such as metallicity and f_sed, as well as for the variation in post formation luminosity (L_pf) introduced in Sect. 6.2. These two filters were chosen for reasons discussed below.

The clear solar metallicity line in the petitCODE grid is sometimes hard to see due to the other lines. For the case of a ten times lower L_pf, the 5 M_⊕ mass planet is too cool to still be on the atmosphere grid and hence these magnitudes were not calculated. Since the 5 and 20 M_⊕ planets are evolving outside of the AMES-Cond grid, the magnitudes corresponding to this atmosphere are not shown for these two low-mass planets.

In principle, the AMES-Cond, the HELIOS, and the petitCODE magnitudes for a clear and solar metallicity atmosphere should be the same. However, because of different input line lists, this is not the case. We quantified this for the 159 M_⊕ planet. In all the JWST/MIRI filters we calculated, the magnitudes calculated with these atmospheres for the evolution of the 159 M_⊕ planet differ at most by 0.5–2 magnitudes. In the JWST/NIRCam filters we considered, the maximum difference ranges from 0.3 magnitudes up to 6 magnitudes in the F227W and F115W filter at a few 100 Myr. In the F356W filter, which is shown in Fig. 10, and chosen because the impact of the line lists is very strong, the heavier planets are brightest when calculated with an AMES-Cond atmosphere. For example, the 100 M_⊕ planet is 2.9 mag brighter at 100 Myr with an AMES-Cond atmosphere than with a clear solar petitCODE atmosphere. A strong methane absorption feature is located at this wavelength, which makes the corresponding flux emission very sensitive to the employed methane line list. The petitCODE and HELIOS models use the recent EXOMOL line list for CH₄ (Yurchenko & Tennyson 2014) turning the atmosphere more opaque than assumed in the older AMES-Cond models before. Hence there is more flux passing through in the AMES-Cond spectra, which leads to a much brighter planet at shorter wavelengths and is most prominent in some filters at shorter wavelengths. For example in the F356W filter the different methane line lists can have a higher impact than a higher post formation luminosity. A higher L_pf influences the magnitudes only early on, by about 1 mag for the 20 M_⊕ for example. In contrast, the same planet with a lower L_pf is up to 3.2 magnitudes fainter in the F356W band early on. For higher planetary masses, when the atmosphere becomes warmer, the effect of the methane becomes smaller as the relative abundance of methane decreases and hence its spectral effect compared to water diminishes, even in the methane bands. The F444W (Fig. 11) filter is, in contrast to the F345W, not sensitive to the methane abundance. Hence, all atmospheres lead to similar magnitudes. The largest impact is now due to different L_pf.

Figure 12 shows the temporal evolution of the blackbody and non-grey spectrum for four planetary masses. The colour code gives the objects’ age. The spectra are calculated for a cloud-free, solar metallicity petitCODE atmosphere by interpolating to the required temperature and surface gravity and are thus given as long as the planet is evolving on the atmospheric grid. We over-plot the sensitivity limits for the JWST/NIRCam instrument (grey dots) and for the JWST/MIRI instrument (black dots)⁴. These are background limits and do not take the final contrast performance of the instruments into account, which will only be known after commissioning of the high-contrast imaging modes and which will also depend on the selected targets and observing strategy. The prominent features that can be seen, for example at 1–2 μm and the peak at 4.7 μm, originate from the water and methane opacities. Especially the cut-off at 1.6 μm is typical for methane. At longer wavelengths, the calculated spectral emission resembles more closely the theoretical blackbody emission. However, for certain temperatures and surface gravities, the spectral flux at shorter wavelengths of 1–2 μm can be up to orders of magnitude higher than the theoretical blackbody flux. As an example, the 5 M_⊕ object shows a blackbody flux at 1 Myr and at 1.6 μm that is 13 orders of magnitudes lower than the one from the calculated spectrum. This can be understood because at the corresponding low temperature of 182 K there are very few absorbers in the atmosphere; water is condensed to quite high pressures in the atmosphere (0.5 bar), and we do not model water clouds here. Therefore, the depth that is probed in the atmosphere is strongly influenced by collision induced absorption (CIA) in addition to the water and methane opacities. As can be seen in Fig. 12, when the atmosphere of the planet becomes warmer or the planet is more massive, this mechanism is diminished.

Figure 12 also shows the unprecedented sensitivity of JWST at thermal infrared wavelengths. For young nearby stars, such as, membersof the β Pictoris moving group with an estimated average age and distance of 23 Myr and 15 pc, respectively (Mamajek 2016), planets with masses below that of Neptune seem to be within reach at separations from the star where background limited performance is achieved. This is truly uncharted territory in comparison to what has been achievable up to now with exoplanet imaging in terms of mass limits (see e.g. Bowler 2016, for a recent review).

Fig. 10 Absolute magnitudes in JWST filter band F356W (JWST/NIRCam) for four different planetary masses and various atmosphere parameters and grids as well as a variation on post formation luminosity with a clear and solar metallicity petitCODE grid. The magnitudes are shown as long as they are in the atmosphere grid, and for the cloudy models (fsed = 1.0) as long as they are above 200 K.

Fig. 11 Absolute magnitudes in JWST filter band F444W (JWST/NIRCam) for four different planetary masses and various atmosphere parameters and grids as well as a variation on post formation luminosity with a clear and solar metallicity petitCODE grid. The magnitudes are shown as long as they are in the atmosphere grid, and for the cloudy models (fsed = 1.0) as long as they are above 200 K.

6.5 Predictions for space and ground-based observations

While JWST will remain unchallenged in terms of sensitivity in the thermal infrared for many years to come, the currently operational extreme adaptive-optics (AO), high-contrast imaging near-infrared (NIR) instruments, such as VLT/SPHERE or Gemini Planet Imager (GPI; Beuzit et al. 2008; Macintosh et al. 2008) achieve better detection limits at small separations close to the diffraction limit, that is, in the contrast-limited regime. To put our models in the context of the exoplanet imaging surveys presently conducted, we show in Fig. 13 our model predictions for the SPHERE/IRDIS H filter. The magnitude was calculated as described in Sect. 6.4.

For the solar metallicity and clear atmosphere, simulations with a ten times higher and lower post formation luminosity relative to the nominal scenario were calculated, as introduced in Sect. 6.3. This was done for four planetary masses, namely for the 5, 20, 100, and 318 M_⊕ mass planets. Since the 5 M_⊕ is too cool to still be on the atmosphere grid, “fainter” magnitudes are not calculated for this mass. The range in post formation luminosity can lead up to a magnitude difference in the SPHERE/IRDIS H band of about 10, 5, and 4 mag for the 20, 100, and 318 M_⊕ (1 M_♃) mass planets at 1 Myr. At 30 Myr, the different cooling paths for the 100 and 318 M_⊕ are no longer distinguishable, which is in agreement with what was found by Spiegel & Burrows (2012) for masses of 1–2 M_♃.

Also note the mass-magnitude degeneracy in Fig. 13, similar to the mass-luminosity degeneracy already noted in Sect. 6.3 and in the literature (e.g. Spiegel & Burrows 2012). For example, at 3 Myr a 20 mag object could correspond to a 159 M_⊕ planet with nominal post formation luminosity or to a 318 M_⊕ planet with a cold scenario post formation luminosity. From 10 Myr on, however, the lines representing the different masses and post formation luminosities no longer overlap, so that a rough mass estimate should be possible at this more likely observable age.

For comparison, the thin lines in the background show the evolution in a clear petitCODE atmosphere (the same lines as in the first column). It is interesting to note that, depending on the metallicity, the clouds seem to have a dimming (for high metallicities) or brightening (for low metallicites) effect.

For further reference, in Fig. 14, isochrones for the VISIR SiC magnitude are shown for solar -metallicity HELIOS (solid lines) and petitCODE (dashed lines) grid atmosphere for times starting at 1 Myr. Magnitudes are calculated as long as the planet is evolving within the atmospheric grid. For the HELIOS grid, all the masses considered here evolve on the atmospheric grid from 1 to 100 Myr. On the other hand, for the petitCODE grid, all the masses evolve on the grid from 3 to 10 Myr. This difference comes from the different temperature coverage of the grids: the petitCODE grid goes from 150 to 1000 K, whereas the HELIOS grid goes from 100 up to 1200 K (see Fig. 7). At 1 Gyr, masses from 50 M_⊕ on are still in the HELIOS grid. The tables are available for log(age/yr) = 6–10 in steps of 0.1 dex while the planet is evolving on the atmospheric grid.

Our evolutionary models together with the applied filter profiles are available at the CDS for a variety of space and ground-based filter systems: NACO (J, H, K_s, L_p, M_p), Cousin (R, I), Wide-field Infrared Survey Explorer (WISE1, WISE2, WISE3, WISE4), VISIR (B87, SiC), and SPHERE (Y, J, H, K_s, Y23, J23, H23, K12). The magnitudes were calculated as described in Sect. 6.4. The available atmospheric models are shown in Fig. 13: petitCODE, clear, [Fe/H] = −0.4, 0.0, 0.4, 0.8, 1.2; petitCODE, f_sed = 1.0, [Fe/H] = −0.4, 0.0, 0.4, 0.8, 1.2; and HELIOS, clear, [Fe/H] = 0.0. This gives a total of 11 atmospheric models and 35 filters, each for eight planetary masses, plus four (three) masses in the clear solar enrichment petitCODE models with a higher (lower) post formation luminosity. There are three masses for the lower post formation luminosity cases because the fainter 5 M_⊕ planet is too cold to be on the atmospheric grid. An example of the tables is given in Tables A.2 and A.3.

Fig. 12 Spectra for cloud-free solar metallicity atmospheres from the petitCODE grid together with the theoretical blackbody for four planetary masses. The age is given in color, the x- and y-axis are the same for all the figures. The temperature of the blackbodies corresponds to the temperature of the planet at the age given in color. The grey dots show the background sensitivity limits for the JWST/NIRCam filters included in this work, the black dots those for JWST/MIRI. There are 2, 3, 5, and 7 spectra shown for the 5, 20, 100, and 318 M⊕ planet.

7 Summary and conclusion

In this study, we first presented (Sect. 2) the extensions made to our evolution model that was originally designed for gas giants but now we applied it to core-dominated, low-mass planets. In Sect. 4, the updated model was then validated against the solar system gas and ice giants. The results we find are in agreement with the literature (e.g. Fortney et al. 2011; Nettelmann et al. 2013). Comparing our simulations with two independent cooling calculations (Baraffe et al. 2008; Lopez & Fortney 2014), we find a satisfactory agreement as well.

We then turned to the main subject of the paper, which is the extension of classical cooling models like Baraffe et al. (2003) or Burrows et al. (1997) to lower mass planets. For this, we computed initial conditions from formation models (Sect. 5) and applied three different sets of atmosphere models (Sect. 3).

The models used in this work include simple Eddington grey atmospheres, the AMES-Cond atmosphere, and recent petitCODE (Mollière et al. 2017) and HELIOS (Malik et al. 2017) atmospheres. The clouds in the petitCODE models are Na₂ S and KCl. Cloud species that are important at lower temperatures such as water are not (yet) included in the cloudy atmospheric models. The surface boundary conditions and therefore the predicted magnitudes have the potential to be significantly affected by water clouds that could form at low temperatures. When comparing the condensation curves of water with the atmospheric p–T structures, one finds that planets with masses below 20 M_⊕ could have water clouds from the start of their evolution in the outer atmosphere. For the more massive planets, water clouds could appear at about 30 (100) Myr for a 50 (100) M_⊕ planet.

In this first publication we have not considered different C/O ratios in the atmosphere. The C/O ratio gives constraints on a planet’s formation path (e.g. Öberg et al. 2011; Madhusudhan et al. 2014; Mordasini et al. 2016; Lavie et al. 2016), and to study the impact of a varied C/O on magnitudes will be an important next step.

The petitCODE, HELIOS, and the AMES-Cond grids assume chemical equilibrium. Although we expect disequilibrium processes like turbulent mixing and photo-chemistry to have a non-negligible influence on the atmospheric composition (Moses et al. 2016), we focus in this study on the treatment of clouds. We postpone the consideration of disequilibrium chemistry to future work.

The initial conditions together with the atmospheric grids and p–T structures were then used to calculate the evolution of a set of eight planetary masses (Sects. 6.1–6.3). These range from 5 M_⊕ to 2 M_♃. The atmospheric models as well as the post formation luminosities were varied. Following this, magnitudes for 35 filters were calculated for clear and cloudy atmospheres for various metallicities (Sects. 6.4 and 6.5). The magnitudes were calculated as described in Sect. 6.4⁵. We summarize the main findings in the following.

• When simulating hot and cold formation scenarios (Sect. 6.3), we found that the spread in the initial (i.e. post formation) luminosity, as suggested by formation models, has a greater influence on a planet’s bolometric luminosity than the atmospheric model, as already noted by Spiegel & Burrows (2012).

• In the SPHERE/IRDIS H filter for example, the difference between a hot versus cold start can be up to 10 mag for the 20 M_⊕ planet (Sect. 6.5). We also note that there is a mass-magnitude indeterminacy with initial conditions in certain filter bands at young ages, similar to what was found by Spiegel & Burrows (2012). However, after 10 Myr a rough mass estimate should be possible (Fig. 13).

• The atmospheres have a large impact on magnitudes in specific filter bands. For example, in JWST/NIRCam F356W the magnitude calculated with an AMES-Cond atmosphere is up to 2.9 mag brighter at 100 Myr for the 100 M_⊕ planet than with a clear solar petitCODE atmosphere (Sect. 6.4). This can be explained by the differences in methane line lists. We assume that the newer line list ExoMol yields more accurate results because it should be much more complete, especially at higher temperatures.

• Comparing the sensitivity limits of JWST with the emergent flux from the planet, we find that a 20, 100, and 318 M_⊕ (1 M_♃) mass planet should be detectable with JWST/MIRI in the background-limited regime until 10, 100, and 1000 Myr after formation, respectively (see Sect. 6.4 and Fig. 12).

• Filters at wavelengths around 4.7 μm seem to be favourable for exoplanet detection, as there is a prominent window in the water and methane opacities (see Fig. 12) that enhances the emergent flux relative to a blackbody of the same temperature.

While it seems unlikely that (sub-)Jupiter mass planets are within the reach of current instruments like SPHERE or GPI for a large sample of targets, there are a few special cases where our models are applicable. A good example is the nearest pre-main sequence star AP Col (Riedel et al. 2011), where ground-based, high-contrast imaging can actually probe for young Jupiter analogs (Quanz et al. 2012).

More importantly, however, our models will be of relevance when planning future observations with JWST and/or the next generation of 30–40 m ground-based telescopes and their exoplanet imaging instruments. While JWST will provide unprecedented sensitivity at mid-infrared wavelength, future ground-based facilities will remain unchallenged in terms of spatial resolution for many years to come and provide complementary discovery space in comparison to JWST. Most notably METIS (Brandl et al. 2016) and the Planetary Camera and Spectrograph (PCS; formally called Exoplanet Imaging Camera and Spectrograph (EPICS); Kasper & Beuzit 2010) for the ESO E-ELT will be equipped with high-contrast imaging cameras and search for exoplanets around the nearest and nearby young stars. This opens up the exciting perspective of detecting young and forming (e.g. van Boekel et al. 2017) low-mass planets and old gas giants instead of young giant planets only (e.g. Bowler 2016) and will put important new constraints on theoretical models of planet formation and evolution. For these future observations, this paper provides a theoretical framework for interpretation.

Fig. 13 Overview of all evolutionary models that were calculated for this work. As an example, the SPHERE/IRDIS H band magnitude is shown. Some of the lower mass planets might not be visible on the axis range chosen here. The petitCODE grid clear magnitudes for the nominal post formation luminosity that are shown as thick solid lines in the first column are repeated in the other panels as thin lines for comparison between clear versus cloudy cases. In the clear solar metallicity panel (first panel, second row), we also show the magnitudes corresponding to an evolution with a ten times higherpost formation luminosity (brighter, dotted) and with a ten times lower post formation luminosity (fainter, dashed) than the nominal case. This simulations were introduced in Sect. 6.3. The magnitudes are only shown as long as the planets evolve in the atmospheric grid, and for the cloudy models as long as they are above 200 K.

Fig. 14 Isochrones in the VISIR SiC filter for the HELIOS (solid) and petitCODE (dashed) solar metallicity grid. The different lines represent times in years, where the dashed isochronescorrespond to the same times as the adjacent solid isochrones. Masses from 5 to 636 M⊕ are covered. The tables are available for log(age yr−1) = 6–10 in steps of 0.1 dex while the planets are evolving within the grid of the atmospheric models.

Acknowledgements

We thank Eric Lopez for providing the numerical data used in Fig. 4, and Isabelle Baraffe for informative correspondence and the data shown in Fig. 3. We also thank them for interesting discussions. E.F.L thanks Marcus Wernberger Jonsson for a thorough Python introduction. E.F.L., 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. G.-D.M. also acknowledges the support of the DFG priority programme SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KU 2849/7-1). We thank the referee for helpful feedback that greatly improved the manuscript.

Appendix A First-order correction for the temperature

In our previous works, the radius of a solid core (using the astrophysical, not geophysical nomenclature) is determined by numerically solving the internal structure equation of mass conservation and hydrostatic equilibrium assuming a differentiated interior consisting of iron, silicate, and ice. As equation of state, the modified polytropic equation of state (EoS) of Seager et al. (2007) is used. This simple equation of state yields the density as a function of pressure for a wide pressure range including the degenerate limit, but neglects the change of the density with temperature. At a given moment in time, this has only little effect on the resulting radii (Valencia et al. 2007; Grasset et al. 2009), especially for the more massive super-Earth/Neptunian and Saturnian planets considered here (Seager et al. 2007), but for the long-term cooling we need to consider it by adding a first order temperature correction of the mean density.

From the law of thermal expansion, we estimate the variation of the mean core density ρ with temperatureas $$\rho =\frac{{\rho }_0}{1+\alpha (T_{\text{ceb}}-T_{\text{ref}})},$$(A.1)

where ρ₀ is the mean density yielded by solving the structure equations with the EoS of Seager et al. (2007). It is a function of the core mass M_core, the pressure P_ceb exerted by the gaseous envelope at the core-envelope boundary, and the ice mass fraction f_ice (see Mordasini et al. 2012b). The iron-silicate ratio is fixed at 2:1 in mass, inspired by Earth’s composition and condensation models of solar-composition gas (e.g. Santos et al. 2015). The other quantities are the thermal expansion coefficient α, the reference temperature T_ref, which we set to 300 K, and the temperature at the core-envelope boundary T_ceb.

The thermal expansion coefficient is found by using the ANEOS equation of state (Thompson 1990) and the Maxwell relations as $$\alpha =\frac{1}{\rho }{\left(\frac{\partial P}{\partial T}\right)}_{\rho }{\left(\frac{\partial \rho }{\partial P}\right)}_T,$$(A.2)

where the two derivatives are an output of the EoS. We considered the ANEOS data for water ice and dunite at higher pressures, and measured values at lower ones (Poirier 2000) for a set of pressure-temperature pairs representative of the Earth’s mantle, the Earth’s centre, Jupiter’s centre, and the centre of a 10 M_♃ super-Jupiter. Figure A.1 shows the α of these twomaterials, while the data is given in Table A.1. The plot shows that the measured (horizontal part) and ANEOS data can be approximated with a broken power law, so that we write $$\alpha =\min \left(a\times {\left(\frac{P_{\text{ref}}}{P_{\text{ceb}}}\right)}^b,{\alpha }_0\right).$$(A.3)

Fig. A.1 Thermal expansion coefficients of water and dunite as predicted by ANEOS (filled and empty circles), and the broken power law approximation used in our model (solid and dashed lines).

The parameters are a = 4 × 10⁻⁶ 1∕K, b = 0.45, and α₀ = 1 ×10⁻⁴ 1∕K for ice, and a = 2 × 10⁻⁶ 1∕K, b = 0.5, and α₀ = 1.5 ×10⁻⁵ 1∕K for dunite. The latter value is chosen to be somewhat lower than the typically used value of 2.0 × 10⁻⁵ 1∕K since we use it for the entire rocky core consisting of silicate and iron, which has a lower α (Poirier 2000). The reference pressure P_ref is 10¹³ dyn cm⁻² in both cases.

It is clear that our description represents only a simple approximation of the actual physical process, as it assumes for example a uniform expansion coefficient for the entire core. The comparisons with more complex models (see also Thomas & Madhusudhan 2016) and observational data presented below nevertheless indicate a relatively good match. A simple model as a first step also seems appropriate given our goal of studying the luminosities of young extrasolar planets, and not for example the detailed internal structure of solar system planets. In future work, we will still include the more accurate description of Alibert (2014).

References

All Tables

All Figures

	Fig. 1 Temporal evolution of the luminosity of the simplified solar system gas giants. The different contributions to the total luminosity are indicated with colours, the luminosity of today is shown as a black dot.
In the text

	Fig. 2 Temporal evolution of the radius of the simplified solar system gas giants. The radius of today is given in the respective colour of the planet as a cross. On this scale, the symbols for the radii of Uranus and Neptune partially overlap.
In the text

Fig. 3 Relative contributions to the total luminosity L over time for the water core simulation of GJ 436 b. The blue and black lines are from Baraffe (2015; priv. comm.), the red and brown lines are from this work and represent the following: solid blue/ brown line: energy release from the core, which is the sum of core cooling and core contraction; dashed blue/ brown line: core cooling; dash-dotted blue/ brown line: core contraction; dashed black/ red line: envelope cooling and contraction. A discussion can be found in Sect. 4.3.

In the text

	Fig. 4 Luminosity as a function of time for a 5 M⊕ planet with a 1% H/He envelope at a distance of 0.1 au from a 1 M⊙ star. The core is assumed to be isothermal. The different contributions are indicated in the plot. The solid black line shows the total luminosity as found by Lopez & Fortney (2014).
In the text

Fig. 5 Results from population synthesis calculations (blue circles) of combined planet formation and evolution that are used for the initial conditions. Orange dots give our final dataset. Left panel: envelope mass as function of total mass with the fit from this work (orange-red, Eq. (11)). For more massive planets the envelope mass as calculated from the relationship found by Thorngren et al. (2016; Eq. (12)) is shown for the 30 M⊕ planet with a green square. Right panel: post formation luminosity as a function of total planetary mass, together with the fit from Mordasini et al. (2017; dark-green, Eq. (14)) and from this work (orange-red, Eq. (13)).

In the text

Fig. 6 Temporal evolution of fundamental properties of the eight simulated planets with initial conditions given in Table 2. The change in intrinsic luminosity, entropy, temperature, and radius is shown. The colour code for the masses is given in the bottom left panel. The luminosities of Jupiter and Saturn are indicated as squares in the top left panel. The dots in the luminosity, effective temperature, and radius panel show the results from Baraffe et al. (2003). Solid lines show the evolution with a cloud-free atmosphere with solar metallicity from the petitCODE grid, dashed lines show it with an AMES-Cond atmosphere, dash-dotted lines show it with a solar metallicity, cloud-free atmosphere from the HELIOS grid, and dotted lines finally show the evolution calculated with a grey atmosphere.Since the three smallest masses are evolving outside of the AMES-Cond grid (see Fig. 7), these masses are not plotted for an AMES-Cond atmosphere. For the AMES-Cond, petitCODE and HELIOS atmospheres, the tracks are shown as long as the planet is evolving on the atmospheric grid. For better visibility, three resp. two evolution curves are not shown in the entropy resp. radius plot. The horizontal lines at 100 and 150 K in the Teff plot show the lower limit of the HELIOS and petitCODE grid, respectively.

In the text

Fig. 7 Evolution of the planets in the logg–Teff space, together with the coverage by the atmosphere grids given as rectangles. The colour code for the masses is given in the middle panel. Thick lines indicate cooling curves that are in the atmosphere grid and thin lines those that are outside of the grid. Top panel: evolution in the AMES-Cond grid. Middle panel: evolution in the petitCODE grid for the cloudy (fsed = 0.5) and solar metallicity atmosphere is shown. This is representative also for other types of petitCODE atmospheres. Bottom panel: evolution in the HELIOS grid for a cloud-free atmosphere with a [M/H] = 0.6 and is representative also for the solar -metallicity evolution paths.

In the text

Fig. 8 Temporal evolution of fundamental properties of the simulated planets with different metallicities and cloudy atmospheres. The colours represent different masses, where the colour code is given in the bottom left panel. The two squares in the top left panel indicate the luminosity of Jupiter and Saturn. Solid lines correspond to a cloud-free atmosphere with solar metallicity, dashed lines correspond to an atmosphere with clouds (fsed = 0.5) and solar metallicity (as in Fig. 6), dotted lines correspond to a cloud-free atmosphere with a metallicity [M/H] = 0.6, all of them using the petitCODE model. Dash-dotted lines indicate a cloud-free atmosphere with a metallicity [M/H] = 0.6 from the HELIOS grid. For better visibility, three resp. two planets are not shown in the entropy resp. radius plot. The horizontal lines at 100 and 150 K in the Teff plot show the lower limit of the HELIOS and petitCODE grid, respectively. The tracks are shown as long as the planet is evolving on the atmospheric grid.

In the text

Fig. 9 Influence of varied post formation luminosity (Lpf) on the intrinsic luminosity and temperature for four different masses. The colour code is given in the left panel, the line style is given in the right panel. All cooling curves are calculated with a cloud-free atmosphere with solar metallicity in the petitCODE grid. Solid lines show the evolution with the nominal Lpf of Table 2. Dashed lines show the evolution with a ten times bigger Lpf (hot scenario). Dash-dotted lines show the evolution with a ten times fainter Lpf (cold scenario).

In the text

	Fig. 10 Absolute magnitudes in JWST filter band F356W (JWST/NIRCam) for four different planetary masses and various atmosphere parameters and grids as well as a variation on post formation luminosity with a clear and solar metallicity petitCODE grid. The magnitudes are shown as long as they are in the atmosphere grid, and for the cloudy models (fsed = 1.0) as long as they are above 200 K.
In the text

	Fig. 11 Absolute magnitudes in JWST filter band F444W (JWST/NIRCam) for four different planetary masses and various atmosphere parameters and grids as well as a variation on post formation luminosity with a clear and solar metallicity petitCODE grid. The magnitudes are shown as long as they are in the atmosphere grid, and for the cloudy models (fsed = 1.0) as long as they are above 200 K.
In the text

Fig. 12 Spectra for cloud-free solar metallicity atmospheres from the petitCODE grid together with the theoretical blackbody for four planetary masses. The age is given in color, the x- and y-axis are the same for all the figures. The temperature of the blackbodies corresponds to the temperature of the planet at the age given in color. The grey dots show the background sensitivity limits for the JWST/NIRCam filters included in this work, the black dots those for JWST/MIRI. There are 2, 3, 5, and 7 spectra shown for the 5, 20, 100, and 318 M⊕ planet.

In the text

Fig. 13 Overview of all evolutionary models that were calculated for this work. As an example, the SPHERE/IRDIS H band magnitude is shown. Some of the lower mass planets might not be visible on the axis range chosen here. The petitCODE grid clear magnitudes for the nominal post formation luminosity that are shown as thick solid lines in the first column are repeated in the other panels as thin lines for comparison between clear versus cloudy cases. In the clear solar metallicity panel (first panel, second row), we also show the magnitudes corresponding to an evolution with a ten times higherpost formation luminosity (brighter, dotted) and with a ten times lower post formation luminosity (fainter, dashed) than the nominal case. This simulations were introduced in Sect. 6.3. The magnitudes are only shown as long as the planets evolve in the atmospheric grid, and for the cloudy models as long as they are above 200 K.

In the text

Fig. 14 Isochrones in the VISIR SiC filter for the HELIOS (solid) and petitCODE (dashed) solar metallicity grid. The different lines represent times in years, where the dashed isochronescorrespond to the same times as the adjacent solid isochrones. Masses from 5 to 636 M⊕ are covered. The tables are available for log(age yr−1) = 6–10 in steps of 0.1 dex while the planets are evolving within the grid of the atmospheric models.

In the text

	Fig. A.1 Thermal expansion coefficients of water and dunite as predicted by ANEOS (filled and empty circles), and the broken power law approximation used in our model (solid and dashed lines).
In the text