© The Authors 2022

1 Introduction

Understanding the origin of giant planets and their characterisation are two key objectives of exoplanetary science. Determining their atmospheric and bulk compositions is particularly important: they can be compared with theoretical predictions and be used to constrain formation and evolution pathways (e.g. Guillot et al. 2006; Mousis et al. 2009; Johansen & Lambrechts 2017; Hasegawa et al. 2018; Ginzburg & Chiang 2020; Shibata et al. 2020; Helled et al. 2022). However, the planetary bulk composition cannot be directly measured and has to be inferred with theoretical models from mass-radius measurements and stellar ages (e.g. Fortney et al. 2007; Miller & Fortney 2011; Thorngren et al. 2016; Teske et al. 2019). The inferred composition depends on model assumptions such as the equations of state (Baraffe et al. 2008; Vazan et al. 2013), the atmospheric model (Poser et al. 2019) and the assumed atmospheric metallicity (Burrows et al. 2007; Müller et al. 2020a). This leads to an intrinsic uncertainty in the predicted bulk metallicity since there are several possible solutions for a given set of measured planetary parameters. Part of this uncertainty can be reduced with more accurate mass and radius measurements expected from ground-based facilities such as High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003), Near Infra-Red Planet Searcher (NIRPS; Bouchy et al. 2017) and Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2021), as well as the upcoming PLAnetary Transits and Oscillations of stars mission (Plato; Rauer et al. 2014) which will also provide an accurate determination of the stellar age.

Unlike the bulk composition, the atmospheric composition of giant exoplanets is measurable. In recent years, there have been detections of several chemical species in hot Jupiter atmospheres (e.g. Madhusudhan et al. 2014; Barstow et al. 2017; Sedaghati et al. 2017; Nikolov et al. 2018; Pinhas et al. 2019; Welbanks et al. 2019; Ben-Yami et al. 2020; Line et al. 2021). However, the atmospheric characterisation of warm Jupiters has so far been limited. The James Webb Space Telescope (JWST; Gardner et al. 2006; Barstow et al. 2015) and in the near future the Ariel mission (Tinetti et al. 2018) will change that with high-precision measurements of exoplanet atmospheres.

One reason atmospheric measurements are vital for the characterisation of giant exoplanets is that their inferred bulk compositions are degenerate with respect to the assumed atmospheric metallicity (Müller et al. 2020a; Müller & Helled 2021). Therefore, these measurements promise to further constrain the bulk composition and therefore the planetary interiors, and provide crucial information for linking the interiors and atmospheres of giant planets (Thorngren & Fortney 2019).

Clearly, the combination of accurate mass, radius, age and atmospheric composition measurements will launch a new era in exoplanet characterisation. In order to take full advantage of these exciting future data, detailed theoretical investigations are required. Such models can identify conditions and key targets for which the new measurements will yield the greatest gain in knowledge about planetary interiors.

In this work we investigate how the upcoming Plato and Ariel missions can improve our knowledge of the composition of giant exoplanets and suggest high-value targets for observations. The paper is structured as follows: in Sect. 2 we describe the thermal evolution models and how they are used to estimate the bulk metallicity of giant exoplanets. We then apply this framework to investigate under which conditions an accurate stellar age measurement can help constrain the heavy-element mass (Sect. 3). In Sect. 4 we investigate how the determination of the atmospheric composition with the Ariel mission can also independently be used as a further constraint. We determine the planets from the Ariel mission reference sample for which this is the case, and suggest several key targets for further investigation. We calculate the mass-metallicity relation for the planets in the Ariel sample and compare it to values from previous studies in Sect. 4.2. In Sect. 5 we investigate how collisions could contribute towards inflating giant planets, which could be important for future observations. We discuss our findings and the limitations of our models in Sect. 6. Finally, our conclusions are summarised in Sect. 7.

2 Methods

We considered a planet with a measured mass (M_obs), radius (R_obs) and host-star age (τ_obs). The age of the planet was assumed to be the same as the star. Since the protoplanetary disk dissipates within about a million years (Mamajek 2009) and the giant planet has to be formed before that, this is generally a valid simplification. The planet’s composition is unknown, but we can calculate evolution models for the measured mass with different bulk metallicities (Z), until we find the composition for which the prediction from the model matches the observed radius (R_model = R_obs), at time τ_obs. This would, in principle, yield a unique bulk metallicity for an observed exoplanet.

However, there are two complications: First, there are measurement uncertainties. Second, the theoretical models have to make a range of assumptions, for example, the equations of state for the different chemical elements, their distribution in the planetary interior, the atmospheric model and the treatment of incident stellar irradiation. These assumptions lead to additional uncertainties in the predictions (Burrows et al. 2007; Baraffe et al. 2008; Thorngren et al. 2016; Müller et al. 2020a). Therefore, the evolution models predict model-dependent probability distributions of the bulk metallicity p(Z | M_obs, R_obs, τ_obs) rather than a single value. These distributions are usually well approximated by a Gaussian p(Z | M_obs, R_obs, τ_obs) ~ 𝒩(µ_Z, σ_Z) with mean value µ_Z and standard deviation σ_Z.

In practice, µ_Z and σ_Z are inferred with a Monte Carlo approach by repeatedly drawing a mass (M), radius (R) and age (τ), from prior distributions. For the mass and radius priors, we assume normal distributions 𝒩(µ_M,R, σ_M,R) with mean values µ_M,R and standard deviations σ_M,R. The age prior is assumed to follow a uniform distribution 𝒰(τ_min, τ_max) with lower and upper limits (τ_min and τ_max) given by observations of the host stars. Each draw of a set of parameters (M, R, τ) is used to estimate the bulk metallicity as described above. If repeated sufficiently many times, the posterior distribution p(Z | M_obs, R_obs, τ_obs) can be estimated.

For an easier comparison between different cases, we often use the relative uncertainty r_x ≡ µ_x/x instead of the uncertainty σ_x. The symbols we used throughout this work are listed and described in Table 1.

To calculate the planetary evolution we used the python package planetsynth (Müller & Helled 2021). Instead of solving the planetary structure and evolution equations (e.g. Kippenhahn et al. 2012) directly, planetsynth generates synthetic cooling tracks by interpolating on a large grid of pre-calculated thermal evolution models. For our purposes, planetsynth calculated the planetary radius R(t | M, Z, Z_atm, F_*) as a function of time t and the planetary parameters mass M, bulk metallicity Z, atmospheric metallicity Z_atm and incident stellar irradiation F_*.

Underlying the synthetic cooling tracks are thermal evolution models calculated with a modified version (Müller et al. 2020b) of the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2022). As stated above, the output of these simulations depends (among other parameters) on the equation of state and the atmospheric model. Müller & Helled (2021) used the Chabrier et al. (2019, hereafter CMS) hydrogen-helium equation of state together with a 50-50 rock-water mixture equation of state for the heavy elements based on the Quotidian Equation of State (QEoS; More et al. 1988; Vazan et al. 2013). The atmospheric model includes the effect of the heavy elements on the radiative opacity (Freedman et al. 2014) as well as the incident stellar irradiation. A limitation of the Freedman et al. (2014) opacity is that the heavy elements are assumed to be present in solar ratios, which is unlikely to be the case for giant planets in general. However, radiative opacities that consider different elemental ratios and can be used in thermal evolution models are currently unavailable. While our fiducial models do not include the effect of grains, we briefly investigate their effect in Appendix C. More details on the evolution models and their equation of state, atmospheric model, and initial conditions can be found in Müller & Helled (2021).

3 The importance of measuring the planetary age and the connection to the plato mission

In this section we explore how the accurate determination of stellar ages can further constrain the composition of giant exoplanets. In order to guide the target selection for the Plato mission, we also study how this depends on the stellar age and the irradiation that a planet receives.

We considered a synthetic planet with R = 1 R_J, M = 1 M_J and solar atmospheric metallicity that is being irradiated by two different stellar fluxes of either F_* = 10⁵ (low irradiation) or F_* = 10⁸ erg cm⁻² s⁻¹ (high irradiation), and three stellar ages of 0.1, 3 and 8 Gyr. We then considered three different relative stellar age uncertainties of r_τ = 50, 25 and 10%. For the mass and radius priors, we used a base-line standard deviation corresponding to an observational uncertainty of 1%. Overall, these combinations of parameters resulted in 18 synthetic planets for which the metallicity was estimated with 1000 000 samples (for each planet) as described in Sect. 2.

The resulting probability densities for the bulk metallicities are shown in Fig. 1. For very young planets (0.1 Gyr), reducing r_τ = 50% to r_τ = 10% reduces the relative uncertainty in the metallicity significantly. This depends on whether the irradiation is weak or strong, with the weakly irradiated planets benefiting more from the lower age uncertainty. For intermediate-age and old planets (3 and 8 Gyr) the opposite is true: an accurate stellar age is more important for highly irradiated planets.

Independent of the irradiation is the fact that older planets have (in absolute terms) more accurately determined compositions compared to younger ones. Also, the age uncertainty has less influence on r_Z for older planets, with r_τ = 10% causing σ_z to be smaller by ~ 10–30% compared to the r_τ = 50% case. This is because the intrinsic luminosity decreases as planets age and therefore their radii only change slowly with time. At early times, the difference in radius for two different times can be much more dramatic.

Figure 2 shows the inferred bulk metallicities and their uncertainties from the models with 50 and 10% stellar age uncertainties. This is depicted as a function of age for high (left panel) and low stellar irradiation (right panel). In general, the bulk metallicity for old, weakly irradiated giant planets is better constrained. As discussed above, this only modestly depends on the stellar age uncertainty. However, in particular young and intermediate-aged planets benefit greatly from a low r_τ. Compared to a 50% uncertainty in stellar age, a 10% uncertainty decreases the uncertainty in the bulk-metallicity by up to a factor of two. It is therefore clear that accurate stellar age measurements are highly valuable for the characterisation of giant exoplanets.

Fig. 1 Distributions of the inferred metallicities for a synthetic planet (R = 1 RJ and M = 1 MJ) at an age of 0.1 (top row), 3 (middle row), and 8 Gyr (bottom row). The panels in each row show the probability densities for decreasing stellar age uncertainties: rτ = 50% (left column), 25% (middle column), and 10% (right column). The lines show the probability densities for low stellar irradiation, and the filled curves for high irradiation (F* = 105 and 108 erg cm−2 s−1, respectively) Normal and bold font weights correspond to the low and high irradiation models.

Fig. 2 Bulk metallicities and their 1σ uncertainties for the synthetic planet with high irradiation (left panel) and low irradiation (right panel). The grey circles and coloured hexagons are for relative age uncertainties of rτ = 50% and rτ = 10%, respectively.

4 The importance of atmospheric measurements and the connection to the Ariel mission

Atmospheric measurements of giant planets can significantly improve our understanding of the planetary bulk composition. First, the atmospheric metallicity can be assumed to be a lower-bound for the bulk (Z_atm ≤ Z). For a homogeneously mixed planet, the atmospheric composition would be the same as the bulk (Z_atm = Z). However, from interior models of Jupiter and Saturn that match observational constraints we expect deep interiors with a higher metallicity than the envelope (e.g. Wahl et al. 2017; Debras & Chabrier 2019; Miguel et al. 2022). Therefore, the lower-bound assumption is valid as long as the location where the heavy elements are measured is not enriched compared to the rest of the outer envelope. Second, the inferred metallicity from mass-radius measurements is degenerate with respect to the assumed atmospheric composition (Müller et al. 2020a). At a given age, the theoretically predicted radius can be significantly different depending on the atmospheric metallicity, because giant planets with enriched atmospheres cool more slowly (Burrows et al. 2007). An atmospheric measurement is therefore an additional and important constraint for the evolution modelling of giant planets.

In this section we demonstrate the importance of atmospheric metallicity measurements and suggest key targets for future observations, focusing on the planets in the Ariel mission reference sample (Edwards & Tinetti 2022). The sample was reduced to include cool-to-warm giant planets with measured masses, radii and stellar ages. We define cool-to-warm giant planets to have a mass between 1–10 M_J with an incident stellar irradiation smaller than F_* ≤ 2 × 10⁸ erg cm⁻² s⁻¹, which is the conventional irradiation flux cutoff for inflated hot Jupiters (Miller & Fortney 2011). They were excluded since their inflation mechanism is still not well understood (Fortney & Nettelmann 2010; Weiss et al. 2013; Baraffe et al. 2014), and therefore their bulk metallicities cannot be estimated properly.

We identify 27 planets in the Ariel mission reference sample that match these criteria. There are no cool giants (with effective temperatures similar to Jupiter) among these planets, and their semi-major axes are all below 0.3 au. The lowest irradiation flux is received by TOI-1899b with a value of F_* ≃ 4 × 10⁶ erg cm⁻² s⁻¹. For comparison, this is still about two orders of magnitude higher than the flux that Jupiter receives (F_* ≃ 5 × 10⁴ erg cm⁻² s⁻¹). The planets are diverse in terms of their masses, and cover the entire parameter space. There are a few planets with similar masses but rather different radii, suggesting a variety of ages or compositions. Further details on the 27 planets are given in Appendix A. There, the planets are listed in Fig. A.2, along with their semi-major axis and their incident stellar irradiation. We also show the measured masses and radii of these 27 planets in Fig. A.1.

Fig. 3 Probability densities of the inferred metallicities for 11 planets from our sample. The colours are the distributions for different atmospheric metallicities (see the legend). The depicted planets show at least moderate separation between the peaks in the different Gaussian distributions (see text for details).

4.1 Inferred bulk composition

We inferred the composition of each planet for four different assumed atmospheric metallicities: Z_atm = 0.5, 1, 4, 8Z_⊙, where Z_⊙ is the solar metallicity. For each Z_atm we inferred the metallicity from 1 000 000 samples (see Sect. 2 for details) constructed from the observed masses, radii, stellar ages and their associated uncertainties¹. This resulted in estimates for the bulk-metallicity probability distributions for each planet, allowing us to investigate whether they are overlapping or significantly different. In the latter case, a measurement of Z_atm would significantly constrain the bulk composition.

Figure 3 shows the probability densities for the inferred bulk metallicities of 11 planets in our sample for each atmospheric metallicity. These planets and their parameters are listed in Table 2. They show a moderate to clear separation of the Gaussian peaks, and are therefore promising targets for atmospheric characterisation by Ariel. For the other 16 planets from the sample, the probability densities for different atmospheric metallicities were moderately to strongly overlapping. As a result, given the current observational uncertainties, an atmospheric measurement would not help constrain their bulk compositions. However, we stress that this may change in the future as planetary parameters are measured more precisely.

There are two planets whose metallicities could not be inferred. First, HAT-P 20b with R = 0.87 ± 0.03 R_J, M = 7.26 ± 0.18 M_J, 2.9 ≤ τ [Gyr] ≤ 12.4 and F_* = 1.99 × 10⁸ erg cm⁻² s⁻¹ (Bakos et al. 2011). This planet is known to be extremely dense and metal-rich, and could not be simulated with our models due to equation of state limitations. Previous estimates suggest a very high total heavy-element mass of M_Z ≃ 660 M_⊕ (Thorngren et al. 2016). We note that a more recent measurement (Esposito et al. 2017) inferred R = 1.03 ± 0.05 R_J at a similar mass, in disagreement with the radius given by Bakos et al. (2011). This alternative is consistent with a solar or few times solar metallicity.

Second, TOI-1899 b has R = 1.15 ± 0.05 R_J, M = 0.66 ± 0.07 M_J, 2.7 ≤ τ [Gyr] ≤ 11.8, and F_* = 3.95 × 10⁶ erg cm⁻² s⁻¹ (Cañas et al. 2020). For its age, the measured radius is inflated by ~15% compared to predictions from the evolution model, despite the irradiation flux being two orders of magnitude below the hot Jupiter cutoff. Giant exoplanets with very high or very low densities such as HAT-P 20b and TOI-1899 b demonstrate the need for more detailed evolution simulations, as well as additional observations to confirm and further constrain their masses and radii.

There are a few general trends to be noted. First, a higher atmospheric metallicity increases the inferred bulk metallicity. This is because a metal-enriched atmosphere is more opaque, which delays the cooling and leads to a larger radius at a given age (Burrows et al. 2007; Müller et al. 2020b). Second, for some planets only the low Z_atm = 0.5, 1 Z_⊙ and high Z_atm = 8 Z_⊙ cases were clearly separated, with the intermediate Z_atm = 4 Z_⊙ overlapping the two distributions by various degrees.

There are, however, planets for which the peaks of the distributions are narrow and clearly separated (except for sub-solar atmospheric metallicities). This is generally the case for low-metallicity planets, such as TOI-1259 A b. One reason for this is that the assumption Z ≥ Z_atm is enforced when the bulk metallicity is inferred, resulting in a strict lower bound and therefore a strong constraint. This separates the low Z_atm cases, which have their peaks at low Z, from the higher Z_atm ones. Notable exceptions are planets that have very low observational uncertainties; for example HATS-72 b. For most planets, there is also a clear difference between the lowest and highest atmospheric metallicities.

From Fig. 3 it is clear that the uncertainty over the entire range of atmospheric metallicities is much larger than for any individual case. To quantify this uncertainty, we assume that each Z_atm is equally likely and calculate the full uncertainty as follows: let the mean values of the bulk metallicity for Z_atm = 0.5, 8 Z_⊙ be ${\mu }_{Z_{0.5,8}}$ with the corresponding standard deviations ${\mu }_{Z_{0.5,8}}$. We calculate the lower bound as $Z_l={\mu }_{Z_{0.5}}-{\sigma }_{Z_{0.5}}$ and the upper $Z_u={\mu }_{Z_{\text{8}}}-{\sigma }_{Z_{\text{8}}}$. The full uncertainty is then given by σ_Z = Z_u − Z_l. The corresponding mean metallicity is then simply the mean within that range.

We now compare this full uncertainty to an individual case where the atmospheric composition is known. Here, for illustration purposes and without loss of generality, we assume that each planet has its atmospheric metallicity determined as Z_atm = 3 Z_⊙. Of course in reality, the measured atmospheric metallicity will have an uncertainty that depends on the individual target and on the stellar properties. It should be noted that our conclusions are independent of this choice, because for each atmospheric metallicity the individual ${\sigma }_{Z_i}$ are nearly identical. The value was chosen to facilitate the comparison of the uncertainties since the mean bulk metallicities line up at around 3 Z_⊙.

The comparison between the cases with the full uncertainty σ_Z (grey circles) against the ${\sigma }_{Z_3}$ (purple hexagons) is shown in Fig. 4. Without information on the atmospheric composition, the relative bulk-metallicity uncertainty (r_Z) is mostly between 60 and 90%. This improves drastically if Z_atm is known, with r_Z typically smaller than 30%. In other words, the uncertainty is reduced by a factor of four to eight. In the most extreme case (HATS-72 b), the uncertainty is reduced by a factor of 16. It is therefore clear that knowledge of the atmospheric composition can significantly reduce the uncertainty in the bulk composition.

Fig. 4 Top: inferred bulk metallicities and their uncertainties for the 11 key targets from the Ariel reference mission sample. Bottom: Relative bulk-metallicity uncertainties in percent for the same planets. The grey dots and error bars represent the metallicities and uncertainties when the atmospheric metallicity is unknown between 0.5 and 8 Z⊙. The purple hexagons are the results for when the atmospheric metallicity is known to be three times solar. See text for details.

4.2 Mass-metallicity relation

We used Bayesian linear regression to infer the mass-metallicity relation for the warm gas giants from the Ariel mission reference sample (see Sect. 4 and Appendix A). The linear regression was performed in log-space, which resulted in a power-law relation between the planetary mass M[M_J] and the heavy-element enrichment of the planet Z_pl/Z_st, where Z_st = 0.014 × 10^[Fe/H] is the metallicity of the host star. The log-likelihood function for a linear model mx + b is $$\begin{array}{l}\ln p\left(y_o,{\mu }_x|x,m,b,{\sigma }_x,{\sigma }_y,f\right)=\hfill \\ -\frac{1}{2}{\displaystyle \sum _n=\left[\frac{{\left(y_{o,n}-m{x}_n-b\right)}^2}{{\sigma }_{y,n}^2}+\ln \left(2\pi s_n^2\right)\right],}\hfill \end{array}$$(1)

where x_n are the independent variables, y_o,n the dependent variables, and $s_n^2={\sigma }_{y,n}^2+f^2{\left(m{x}_n+b\right)}^2$ with σ_y,n being the uncertainty on y_o,n. The parameter f accounts for the possibility that the uncertainties σ_y,n are underestimated by a (constant) fraction of y_o,n. The prior for x is the (normal) distribution from observations: p(x) ∝ ∏_n 𝒩(µ_x,_n, σ_x,n). We used uniform non-informative priors p(m, b, f) for the model parameters. The posterior probability function (up to some constant) is then $$\begin{array}{l}p\left(x,m,b,f|y_o,{\mu }_x,{\sigma }_x,{\sigma }_y\right)\propto \hfill \\ p\left(m,b,f\right)p\left(x\right)p\left(y_o,{\mu }_x|x,m,b,{\sigma }_x,{\sigma }_y,f\right).\hfill \end{array}$$(2)

As shown in Fig. 3, the assumed atmospheric metallicity can have a large influence on the inferred bulk composition. To account for this in our fitting procedure, we created 1000 bootstrap samples of the planets with a random atmospheric metallicity Z_atm [Z_⊙] drawn from the set {0.1, 0.5, 1, 2, 4, 8}. Since we have no information about the atmospheric composition of these planets, we assumed that they were equally likely. For each bootstrap sample, we used a Markov chain Monte Carlo (MCMC) method to estimate the posterior distribution for the fitting parameters. The combined bootstrap and MCMC samples were then used to calculate the best fit parameters. The best fit line, as well as the one-σ error contour are shown in Fig. 5. For readability purposes, the error bars show the uncertainty of the inferred heavy-element enrichment for a solar metallicity atmosphere. Figure 6 shows the two-dimensional projection of the posterior probability distributions of the power-law Z_pl/Z_st = β × M[M_J]^α lit parameters α ≡ m and β≡ 10^b.

The best lit parameters are α = −0.37 ± 0.14 and β = 7.85 + 1.50. We calculated Kendall’s tau correlation and the corresponding p-value from the bootstrap sample, which yielded a correlation of −0.32 (p-value: 0.035), where the p-values of the bootstrap samples were calculated as described in Efron & Tibshirani (1994). Since we were using bootstrap samples, the correlation and the p-value accounted for the uncertainties in the independent (planetary mass) as well as the dependent (inferred metallicities) variables.

For comparison, Thorngren et al. (2016) inferred α_T16 = −0.45 ± 0.09, β_T16 = 9.70 ± 1.28 and a correlation of −0.44 (p-value: 1.3 × 10⁻⁵) for the planets in their sample. They also calculated Kendall’s tau; however they did not account for the uncertainties and used the mean values instead. While similar, our best lit for the Ariel sample yields a lower metal-enrichment at a given planetary mass, drops off more slowly (flatter slope) and is more uncertain. However, the Thorngren et al. (2016) best lit is mostly consistent within one σ of our results.

The lower metal-enrichment trend can be explained by the different hydrogen-helium equations of state that were used. Here, the evolution models used the CMS hydrogen-helium equation of state, while Thorngren et al. (2016) used the Saumon et al. (1995, hereafter SCvH) one. It has been shown previously that CMS generally leads to smaller planetary radii, and therefore a lower heavy-element content will be inferred for a given observed planet (Müller et al. 2020a). On the other hand, the different slope and the higher uncertainties of the fit parameters are likely the results of our models considering non-solar atmospheres, and having fewer (and a few different) planets for the fit. These are also the likely causes for the lower and less statistically significant correlation.

As discussed in Thorngren et al. (2016), these results are inconsistent with planetary population synthesis models (e.g. Mordasini 2014). Hasegawa et al. (2018) suggested that a power law with roughly Z_pl/Z_st ∝ M[M_J]^−0.4 emerges if the metal enrichment is a tracer for the final stages of giant planet formation, and the heavy elements are mostly accreted from gapped planetesimal disks. This interesting topic should be investigated further by developing more comprehensive formation models and their comparison with future accurate data.

Fig. 5 Heavy-element enrichment as a function of the planetary mass. The black dots with error bars show the inferred metallicity of the potential Ariel targets. The solid purple line shows the best-fit line, the purple-shaded region shows the 1σ error contour around the best fit, and the dashed grey line shows the best fit from Thorngren et al. (2016).

Fig. 6 Corner plot showing the one- and two-dimensional projections of the posterior probability distributions of the fit parameters α and β.

5 Radius inflation due to collisions

Our models assumed that the planet is isolated (except for the irradiation). Therefore, we did not account for events such as giant impacts, which could affect the subsequent evolution and therefore the observed radius. While no direct observations exist, theoretical studies suggest that giant impacts are common during the early evolution of planetary systems (e.g. Liu et al. 2015). Collisions have also been used to explain, for example, Jupiter’s apparently dilute core (Liu et al. 2019), the high metallicity of some exoplanets (e.g. Ogihara et al. 2021) or their eccentricities (Frelikh et al. 2019). Typically, evolution models do not include the effects such impacts could have on the predicted sizes of exoplanets.

Here, we investigated how the deposited energy from collisions inflates the radius of a giant planet, and the timescale associated with that inflation. Our approach was as follows: We started with a somewhat typical warm giant (M = 1M_J, Z = 0.1, F_* = 10⁷ erg cm⁻² s⁻¹), and let it evolve for 2 Gyr (since most observed warm giants are commonly at least this old). To simulate an impact, we injected a certain amount of energy into the giant planet on a very short timescale, and ignored the mass deposition for simplicity. The specific energy of the impact (erg g⁻¹) was deposited uniformly between a pressure of P_top = 10⁶ Ba and P_bot, which we varied from 10⁹–10¹³ Ba or used the pressure at the centre. This accounts for the fact that the impactor would likely be ablated and ripped apart by the dense envelope of the target. Afterwards, the planet evolved and cooled usual.

Guided by the simulations of the giant impacts for Jupiter of Liu et al. (2019), we considered impact energies of E_i = 10³⁸, 10³⁹, 10⁴⁰, 10⁴¹, and 10⁴² erg. The upper limit roughly corresponds to the difference in energy of Jupiter before and after the collision with a 10 M_⊕ planetary embryo in Liu et al. (2019), and is similar to the impactors kinetic energy. If the impact velocity were roughly at the escape speed of Jupiter today, this would correspond to the kinetic energy of the impact with approximate masses of M_i ≈ 10⁻³, 10⁻², 10⁻¹, 1, and 10 M_⊕. Clearly such impacts would also alter the planetary interior structure, deposit mass and could also cause envelope mass-loss. We also did not consider that a very energetic impact could create composition gradients and inhibit convection, leading to a non-adiabatic evolution. Instead, we intentionally kept our model simple in order to focus on studying the effect of the post-impact radius inflation that is solely associated with the energy of the impact.

Figure 7 shows the radius evolution after the impact for all the cases we considered. Each row shows a different impact energy, and the colours correspond to the various P_bot values. For impact energies beyond 10³⁹ erg, using small P_bot values resulted in a non-hydrostatic configuration, and therefore the evolution could not be simulated further. It is unlikely that energetic impacts are associated with impactors that only penetrate to such low pressures (see e.g. the 1 M_⊕ impact of Liu et al. 2019). Depending on the impact energy and P_bot, shortly after the impact the giant planet can be significantly inflated, up to almost twice of its original size (for the highest energy we considered). However, the post-impact cooling is relatively rapid in most cases. We find that for the radius to remain inflated beyond a few times 10³–10⁴ yr, either P_bot has to be small or the impact must be highly energetic, requiring very massive impactors that are unlikely to be present at ages of Gyrs. We do, however, find that impacts could affect the early evolution of young giant planets. This could affect the interpretation of directly imaged planets and we hope to address this in future research.

6 Discussion

The characterisation of giant exoplanets relies on evolution simulations that require assumptions and simplifications. The limitations of the evolution models used in this work are discussed in detail in Müller & Helled (2021). Here, we briefly summarise some of the most influential assumptions and parameters. Important factors are the equations of state (Baraffe et al. 2008; Vazan et al. 2013), atmospheric model (Burrows et al. 2007; Poser et al. 2019; Müller et al. 2020a), the distribution of heavy elements in the planetary interior (Baraffe et al. 2008; Thorngren et al. 2016), and whether the interiors are fully convective (Leconte & Chabrier 2012; Kurokawa & Inutsuka 2015). These parameters are not fully constrained and are therefore limitations of evolution models in general. In most cases, they can influence the predicted radius on the percent level, which is larger than some observational uncertainties today.

An important uncertainty is whether grains exist in the atmospheres of exoplanets, which could slow the planetary cooling due to an enhancement of the opacity (Movshovitz et al. 2010; Vazan et al. 2013; Mordasini et al. 2014). While our evolution models do not include the grain opacity, we briefly demonstrate how the cooling could be affected in Appendix C. If grains cause the atmosphere to be significantly more opaque, the inferred metallicities would be shifted towards higher values.

An additional possible mechanism for changing the radius of a giant exoplanet is a giant impact, which we have investigated in Sect. 5. It should be noted that collisions may be unique to any given planet, and that the possible characteristics of potential collisions include a very large parameter space (mass of the impactor, impact timing, impact velocity, etc.), which is difficult to generalise. While the aim of this work is to investigate exoplanets on a population level, it is nonetheless valuable to estimate the degree to which collisions leave a planet inflated. Our results suggest that only when collisions are extremely energetic (>10⁴¹ erg) a significant inflation lasting longer than a few thousand years is expected. Unless giant impacts are very frequent or extremely energetic, which is unlikely for older planetary systems (Liu et al. 2015; Frelikh et al. 2019), collisions are unlikely to significantly inflate giant exoplanets over long timescales.

For planets younger than a few tens of megayears, additional uncertainties are the unknown initial thermal state of the planets (e.g. Baraffe et al. 2003; Marley et al. 2007; Berardo & Cumming 2017) and the fact that the evolution of the host star is not considered. Since in this work we focus on older planets, these uncertainties are relatively unimportant.

Another important source of uncertainty is linked to the commonly used hydrogen-helium equations of state (e.g. Saumon et al. 1995; CMS) that do not include the non-ideal interactions of the two elements. Recently, Chabrier & Debras (2021) showed the importance of this effect in the low-temperature and high-density conditions expected in giant planets (and brown dwarfs). They provided a new effective hydrogen table, which can be used to calculate hydrogen-helium mixtures including these interactions. It is, however, unclear to what extent this would influence the mass-radius relation of giant planets. While the implementation of this new equation of state into the evolution models is beyond the scope of this work, we calculated hydrostatic models of pure hydrogen-helium planets in order to quantify this effect. We find that for planets with masses between 0.1 and 10 M_J, the inclusion of the non-ideal effects from Chabrier & Debras (2021) yielded a 1–2% larger radius compared to the CMS equation of state. Therefore, this uncertainty is comparable in magnitude to the other model assumptions.

In this work we did not focus on the characterisation of giant exoplanets but rather on the conditions for which more precise or additional measurements can lead to lower uncertainties in the bulk composition. Our calculations depend on the above assumptions: Changing, for example, the equation of state would lead to a systematic shift in the bulk metallicities. It would, however, not change the uncertainty of the result for a given observation and theoretical model. Therefore, this does not affect our main conclusions.

In Sect. 3 we investigated how the uncertainty in the bulk composition depends on the uncertainty of the stellar age measurement. In order to isolate this variable, uncertainties in the mass and radius had to be a fixed value. We set this value to 1%. The reason behind this choice was that larger uncertainties, in particular for the radius, would quickly be the main contributors to the bulk-metallicity uncertainty. Currently, there is not a large number of warm giant exoplanets with similar measurement errors (see Figs. B.1 and B.2). However, this is expected to change in the near future, with more precise mass-radius measurements from ground missions as well as the Plato mission. We also set the mass and radius of the synthetic planets that we have investigated to R = 1 R_J and M = 1 M_J. While this choice is biased towards a Jupiter analogue, it is a reasonable representation of a somewhat typical warm gas giant (see Fig. A.1). We performed additional calculations with different masses and radii, and confirmed that they lead to the same conclusions.

In Sect. 4 we found that there are currently only 27 warm giants with measured masses and radii that are part of the Ariel mission reference sample. Detections of hot Jupiters are much more numerous, and the number of observed warm giants is comparably small. This is clearly not ideal for statistical inferences, such as constructing a mass-metallicity relation. Hopefully, in the near future new observations will increase the number of observed cool and warm giants substantially, which will make statistical inferences more reliable.

We note that in our models the atmospheric metallicity is used as a single free parameter, which is clearly a simplification. In reality, the metallicity is made up of different chemical elements present in various ratios, influenced by the composition of the host star and the formation history of the planet (e.g. Turrini et al. 2018). In addition, exo-atmospheric observations determine the abundance of specific compounds at relatively high regions in the atmosphere and not the overall atmospheric and planetary bulk metallicity. It takes an additional, model-dependent step to infer the overall atmospheric metallicity (e.g. Heng 2018). Many of the current atmospheric composition measurements are limited to a narrow wavelength, and cannot measure the abundance of all the chemical species that carry the most common heavy elements (Zellem et al. 2019). The Ariel spacecraft is expected to have a broad spectral and molecular coverage, which should provide much better estimates of the overall atmospheric metallicity (Tinetti et al. 2018). We suggest that detailed investigations of the connection between the compositions of the outer atmosphere, the deeper atmosphere and the bulk interior should be investigated in future work.

The results from Sects. 3 and 4 suggest that the bulk-metallicity uncertainty due to not knowing the atmospheric composition is larger than the one related to the stellar age uncertainty. This highlights the great power of atmospheric characterisation in determining the composition of giant exoplanets, and also shows that there is a synergy between more accurate age and atmospheric measurements. While our focus in this work was on the planets from the Ariel mission reference sample, our findings are also valid for observations with JWST.

Fig. 7 Radius evolution after the impact. The rows show the evolution for different impact energies. The colours correspond to the upper pressure boundary (lower boundary in terms of radial coordinates) of where the energy is deposited (see text for details). The dashed grey line shows the radius evolution without an impact. Some coloured lines for specific combinations of the impact energy and pressure depth are missing because they resulted in non-hydrostatic configurations.

7 Summary and conclusions

In this paper we have studied in detail how the accurate determination of stellar ages from Plato and atmospheric measurements from Ariel will improve the characterisation of giant exoplanets. Our main results can be summarised as follows:

1. The accurate determination of stellar ages from the Plato mission will significantly constrain the bulk composition of young to intermediate-aged planets (below a few gigayears) and reduce the uncertainty by up to a factor of two. This is true for both weakly and highly irradiated warm giants.

2. The Ariel mission reference sample includes 27 warm gas giants with parameters that would allow an accurate determination of their bulk compositions.

3. Our Bayesian linear regression resulted in the following power-law relation between the planetary mass and metal enrichment for these 27 warm giants: Z_pl/Z_st = (7.85 ± 1.50) × M[M_J]^{−0.37±0.14}. This suggests a lower metal enrichment at a given planetary mass compared to previous studies.

4. We find that the observed radius of TOI-1899 b is ~15% larger than values from evolution models. We therefore suggest that it must be inflated by a mechanism that is not considered in standard evolution models.

5. Measuring the atmospheric metallicity of giant planets can reduce the uncertainty in the bulk composition by a factor of about four to eight. The largest improvements are expected for planets that have low radius, mass, and age measurement uncertainties (in this order).

6. We identify 11 high-priority warm giants from the Ariel mission reference sample for atmospheric characterisation (see Fig. 3 and Table 2). For these planets, measuring the atmospheric metallicity will break the degeneracy in their inferred bulk compositions.

7. We find that impacts are unlikely to cause old giant planets to be significantly inflated. However, they could influence the early evolution of young giant planets, which could affect the interpretation of directly imaged planets.

Our study shows the important role of evolution models in interpreting the upcoming accurate measurements from space missions dedicated to giant exoplanet detection and characterisation. It is clear that there is a great synergy between atmospheric measurements and a more accurate determination of planetary masses, radii, and stellar ages, and that their fusion will considerably improve our understanding of giant exoplanets.

Acknowledgements

We thank the anonymous referee for thoroughly reading our paper and providing useful comments and suggestions. We gratefully acknowledge support from the UZH Forschungskredit grant FK-21-132, the SNSF grant 200020_188460 and the National Centre for Competence in Research ‘PlanetS’ supported by SNSF. This research used data from the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Extensive use was also made of the MESA code (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2022) and the python packages NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Jupyter (Kluyver et al. 2016), emcee (Foreman-Mackey et al. 2013) and planetsynth (Müller & Helled 2021).

Appendix A Potential Ariel targets

As described in Sect. 4, we selected the warm giants (F_* ≤ 2 × 10⁸ erg cm⁻² s⁻¹) with measured masses, radii and stellar ages from the Ariel mission reference sample. We also limited their masses in the range of 0.1 – 10 M_J. These planets are shown in the mass-radius plane in Fig. A.1 coloured according to their incident stellar irradiation. It can be seen that the planets occupy the full mass range, and for any given mass there are planets with different radii.

Fig. A.1 Mass-radius relation of the planets in the sample. The colours indicate the stellar irradiation received by the planets.

We also list the warm giants in Fig. A.2 together with their semi-major axis. It is clear that all the warm giants in the Ariel mission reference sample have tight orbits around their host star. None of them are cool giants.

Fig. A.2 Planets in the sample as a function of their semi-major axis. The size of a dot corresponds to the planet’s radius, and the colour shows its incident stellar irradiation.

Appendix B Observational uncertainties of giant planets

As discussed throughout this work, the measurement uncertainties of the planetary and stellar parameters have a large influence on the inferred bulk-metallicity uncertainty. In Fig. B.1, we show the distributions of the relative measurements uncertainties (σ_x/x) for the planetary masses and radii, as well as the stellar ages for all giant exoplanets with masses of 0.1 ≤ M [M_J] ≤ 10. This population includes hot Jupiters, which are currently unavailable for characterisation due to their inflated radii. The relative uncertainties peak at around 5% for the mass, 3% for the radius and 50% for the stellar age.

We show the same distributions for the potential Ariel targets (see Fig. A.2), which does not include hot Jupiters. Compared to Fig. B.1, these planets have higher observational uncertainties for the mass, similar for the radius and lower for the stellar age.

Fig. B.1 Histograms of the relative measurement uncertainties in mass (left), radius (centre), and age (right) of observed giant exoplanets with 0.1 < M < 10MJ. The data are taken from the NASA Exoplanet Archive.

Fig. B.2 Histograms for the relative uncertainties in mass (left), radius (centre), and age (right) of the potential Ariel targets (see Fig. A.2).

Appendix C The influence of grain-enriched atmospheres

While the opacity used in this work accounts for the influence of metals, we did not consider an opacity enhancement due to grains. Grain-rich atmospheres are more opaque, effectively trapping heat from the interior for longer, particularly at young ages (Vazan et al. 2013). Currently it is still unclear whether grains are present in exoplanet atmospheres over large timescales, or settle quickly to deeper regions (Movshovitz et al. 2010; Mordasini 2014).

Here, we investigated the effect of grains on the cooling of giant exoplanets. We calculated the evolution of a M = 1M_J, Z = 0.1 exoplanet with a modified version of MESA that includes a grain opacity (Müller et al. 2020a; Müller & Helled 2021). We used the analytical fit from Valencia et al. (2013) to the interstellar medium (ISM) grain opacity from Ferguson et al. (2005). ISM grains are very small, and therefore strongly enhance the opacity. However, the full ISM grain opacity is likely about two orders of magnitude too large for exoplanets (Movshovitz et al. 2010; Mordasini 2014). Therefore, we calculated the evolution several times using the full, 10%, 1% and no grain opacities. Since the grain opacity is temperature dependent, we also used two different stellar irradiation fluxes of 10⁴ and 10⁷ erg cm⁻² s⁻¹, which are typical fluxes for cold and warm Jupiters.

In Fig. C.1 we show the radius evolution for all cases. With both stellar fluxes, the full ISM grain opacity yields significantly larger radii, even at an age of a few gigayears. When the grain opacity is reduced by a factor of ten, the inferred radii are still significantly larger for the low-flux case, while for the high-flux case only a small inflation remains. Once the grain opacity is reduced to 1% of its original value, the radius of the low-flux planet is very similar to the case without grain opacity.

A strongly enhanced opacity would lead to larger radii at a few gigayears, and therefore change the inferred bulk metallicity. As noted above, in exoplanets the grain opacity is likely to be ≲ 1% of the ISM grain opacity, and therefore the effect is expected to be small. There are other possible way to enhance the opacity, such as thick cloud layers (Poser et al. 2019). However, these cloud models have many free parameters that are not well constrained. It is therefore important to develop realistic atmospheric models that include grains and clouds in order to take full advantage of future observations.

Fig. C.1 Radius evolution for a M = 1 MJ, Z = 0.1 planet with different grain opacities (see the legend) and stellar irradiation fluxes of 104 (left) and 107 erg cm−2 s−1 (right).

References

All Tables

All Figures

Fig. 1 Distributions of the inferred metallicities for a synthetic planet (R = 1 RJ and M = 1 MJ) at an age of 0.1 (top row), 3 (middle row), and 8 Gyr (bottom row). The panels in each row show the probability densities for decreasing stellar age uncertainties: rτ = 50% (left column), 25% (middle column), and 10% (right column). The lines show the probability densities for low stellar irradiation, and the filled curves for high irradiation (F* = 105 and 108 erg cm−2 s−1, respectively) Normal and bold font weights correspond to the low and high irradiation models.

In the text

	Fig. 2 Bulk metallicities and their 1σ uncertainties for the synthetic planet with high irradiation (left panel) and low irradiation (right panel). The grey circles and coloured hexagons are for relative age uncertainties of rτ = 50% and rτ = 10%, respectively.
In the text

	Fig. 3 Probability densities of the inferred metallicities for 11 planets from our sample. The colours are the distributions for different atmospheric metallicities (see the legend). The depicted planets show at least moderate separation between the peaks in the different Gaussian distributions (see text for details).
In the text

Fig. 4 Top: inferred bulk metallicities and their uncertainties for the 11 key targets from the Ariel reference mission sample. Bottom: Relative bulk-metallicity uncertainties in percent for the same planets. The grey dots and error bars represent the metallicities and uncertainties when the atmospheric metallicity is unknown between 0.5 and 8 Z⊙. The purple hexagons are the results for when the atmospheric metallicity is known to be three times solar. See text for details.

In the text

	Fig. 5 Heavy-element enrichment as a function of the planetary mass. The black dots with error bars show the inferred metallicity of the potential Ariel targets. The solid purple line shows the best-fit line, the purple-shaded region shows the 1σ error contour around the best fit, and the dashed grey line shows the best fit from Thorngren et al. (2016).
In the text

	Fig. 6 Corner plot showing the one- and two-dimensional projections of the posterior probability distributions of the fit parameters α and β.
In the text

Fig. 7 Radius evolution after the impact. The rows show the evolution for different impact energies. The colours correspond to the upper pressure boundary (lower boundary in terms of radial coordinates) of where the energy is deposited (see text for details). The dashed grey line shows the radius evolution without an impact. Some coloured lines for specific combinations of the impact energy and pressure depth are missing because they resulted in non-hydrostatic configurations.

In the text

	Fig. A.1 Mass-radius relation of the planets in the sample. The colours indicate the stellar irradiation received by the planets.
In the text

	Fig. A.2 Planets in the sample as a function of their semi-major axis. The size of a dot corresponds to the planet’s radius, and the colour shows its incident stellar irradiation.
In the text

	Fig. B.1 Histograms of the relative measurement uncertainties in mass (left), radius (centre), and age (right) of observed giant exoplanets with 0.1 < M < 10MJ. The data are taken from the NASA Exoplanet Archive.
In the text

	Fig. B.2 Histograms for the relative uncertainties in mass (left), radius (centre), and age (right) of the potential Ariel targets (see Fig. A.2).
In the text

	Fig. C.1 Radius evolution for a M = 1 MJ, Z = 0.1 planet with different grain opacities (see the legend) and stellar irradiation fluxes of 104 (left) and 107 erg cm−2 s−1 (right).
In the text