© The Authors 2024

1 Introduction

In recent decades, the discovery of strongly irradiated giant exoplanets (Mayor & Queloz 1995) has opened new avenues of exploration in planetary science. These exoplanets, named hot Jupiters, exhibit strikingly different atmospheres compared to the giant planets in our Solar System (e.g., Guillot et al. 1996; Guillot & Showman 2002), rendering them compelling subjects for in-depth investigation. Understanding their unique atmospheric characteristics and their impact on planetary observables has become crucial for advancing our knowledge of exoplane- tary systems. Furthermore, the discovery of directly imageable exoplanets (e.g. Chauvin et al. 2004, 2005; Marois et al. 2008, 2010; Lagrange et al. 2009) enables us to explore the early, still-warm stages of giant planets, providing valuable insights into their formation processes.

A better understanding of the gas giants of our Solar System has had meaningful consequences when it comes to studying exoplanets. On one hand, the difficulty in accurately evaluating the evolution of the ice giants raises questions regarding their internal structure (Scheibe et al. 2019). On the other, the results from the Juno mission highlight the ongoing difficulty in understanding the layout of Jupiter’s planetary core, combined with the limitations of current equations of state (e.g. Nettelmann et al. 2021; Idini & Stevenson 2022; Militzer et al. 2022; Howard et al. 2023). We still have far more data for our close neighbours than we do for exoplanets, which makes the former ideal objects of study to fine tune current modelling tools. However, the growing number of discovered exoplanets represents an ideal test bed for evaluating and exploring the wider uses of these tools and for performing statistical analyses.

In addressing these complexities, the integration of comprehensive grids of model atmospheres has played a pivotal role (e.g. Baraffe et al. 2003; Phillips et al. 2020; Marley et al. 2021). Each successive generation of models has been designed to refine and augment various physical aspects, whether through improved interior equations of state (EOS) and disequilibrium chemistry (Phillips et al. 2020) or the exploration of the effects of varying C/O ratios (Marley et al. 2021). As emphasised by Marley et al. (2021), the diversity of model grids is essential for their complementary nature and intercomparison, enhancing their overall efficacy. Furthering the development of exoplanet model grids requires the inclusion of stellar irradiation and improvement of the coherence between planet interior and atmosphere. The present work addresses both issues.

The modelling of exoplanet interiors often faces significant degeneracy, underscoring the importance of integrating interior models with atmospheric properties. This integration is exemplified in Bloot et al. (2023), where atmospheric metallicity is used to constrain interior designs. Additionally, Müller & Helled (2021) explore how atmospheric properties contribute to constraining planetary bulk composition. Their study of 51 Eri b showcases the use of observed luminosity, derived metallic- ity, and synthetic evolution tracks to estimate mass and bulk metallicity. The strength of such an approach lies in its coherent connection of various parameters. In our work, we propose to advance this approach further by integrating a more sophisticated atmospheric model capable of generating synthetic spectra. This advancement grants us the ability to derive atmospheric properties from observations while maintaining constraints imposed by the interior. Guzmán-Mesa et al. (2022) exemplify the strengths of this type of model when extracting the properties of GJ436 b, a sub-Neptune-like planet, through a coupled interior–atmosphere model with spectral retrieval. Our understand- ing of Sub-Neptunes is plagued by degeneracy, as highlighted by Valencia et al. (2013). Expanding upon this framework, we extend the approach to encompass Jupiter and super-Jupiter-like planets for both transit and direct imaging studies.

When irradiated, predicting the observed radii of hot Jupiters is distinctly challenging due to the inflated nature (Fortney et al. 2021). Accurately modelling these bodies necessitates the consideration of stronger intrinsic thermal fluxes than expected for the ages of the systems they are found in. Properly accounting for this phenomenon is crucial not only for 1D model grids but also for 3D general circulation models (GCMs). Recent studies have demonstrated the possibility of establishing a connection between the equilibrium temperature of a hot Jupiter and its intrinsic thermal flux, providing a plausible explanation for the observed radius (Thorngren & Fortney 2018; Thorngren et al. 2019). Currently, the two favoured hypotheses are ohmic dissipation (Batygin & Stevenson 2010) and advection of heat by atmospheric circulation (Tremblin et al. 2017). With the present work, we aim to address this issue by directly coupling atmospheric and interior modellings so as to evaluate the evolution of irradiated planets by resolving the thermal evolution equation across a grid of planets in equilibrium at different intrinsic temperatures.

In this work, we try to add to the models available by producing a self-consistent interior-atmosphere model that ties physical parameters together. We believe that this will allow better constraint of parameters such as mass, intrinsic temperature, metallicity, core mass, and radius from observations. In Section 2, we explore the atmosphere and interior models used and how we link the two into one that is self-consistent in order to build grids of models. We make literature comparisons to verify the veracity of the results. In Section 3, we evaluate how exoplanets are affected by stellar irradiation. In Section 4, we discuss the various planetary parameters and how they can lead to degeneracies. In Section 5, we explore how Bayesian inference tools can be used with our model to infer parameters that explain an observed spectrum. Finally, in Sections 6 and 7, we show how the model can be applied to infer planetary parameters from observations.

2 Model description

2.1 Atmosphere model

The atmosphere model used in this work is Exo–REM (Baudino et al. 2015, 2017; Charnay et al. 2018; Blain et al. 2021). Exo–REM is a 1D radiative-equilibrium model that calculates fluxes using the two-stream approximation assuming hemispheric closure. Radiative-convective equilibrium is solved assuming that the net flux (radiative + convective) is conserved. The conservation of flux over the grid of pressure levels is solved iteratively using a constrained linear inversion method. Thermodynamic quantities in Exo–REM are linked to one another via the ideal gas law. The vertical chemical profiles for the various species are calculated for a given temperature profile assuming some non-equilibrium chemistry between C-, O-, and N-bearing compounds. To do this, we use an analytical formulation based on a comparison of chemical time constants with vertical mixing time from Zahnle & Marley (2014). For this calculation, vertical mixing is parameterised by an eddy mixing coefficient, K_zz. To maximise the comparability of the models, we chose to use a constant value for the K_zz, as done by Phillips et al. (2020): we take log(K_zz) = 10 (with K_zz in cm²/s). Phillips et al. (2020) demonstrate that the eddy mixing coefficient has limited effects on evolution models. We include the option to use equilibrium chemistry. Condensed species that would otherwise have formed clouds are removed. Exo–REM is parameterised by a gravity at a pressure of 1 bar (g), an intrinsic temperature (T_int), an incoming stellar irradiation (T_irr), and a metallicity ratio (Z) expressed as a fraction of solar metallicity using values from Asplund et al. (2009). The intrinsic temperature is derived from an incoming flux at the bottom of the atmosphere model, which serves as a boundary condition, using Equation (1): $$T_{int}={\left(\frac{F_{int}}{{\sigma }_B}\right)}^{\frac{1}{4}}.$$(1)

The expressions of T_irr and T_int are similar to those found in Guillot (2010).

Exo–REM can evaluate the gravity throughout the atmosphere using the hydrostatic balance equation (HSE). The gravity varies with pressure and altitude. It is hence possible to evaluate how gravity varies at the bottom of different atmospheres as a function of the model parameters, and most notably the stellar irradiation, as shown in Figure 1. Exo–REM is run from a pressure of 10⁻⁴ to 10⁴ bar. The list of molecules considered and the references for the line lists used are given in Table 1¹.

2.2 Interior model

The interior model used in this work is Exoris (Licari 2016; Mazevet et al. 2019, 2022), which solves the hydrostatic equations for a given set of ab initio equations of state. Two regions are considered: an H-He envelope and a core that can be divided into an outer ice and an inner rocky part. The envelope is considered adiabatic and consists of a mixture of hydrogen and helium at a given mass fraction denoted γ_He. The core of mass M_c is considered isothermal. The isentrope is computed by interpolation in an EOS dataset at a constant entropy that is fixed at the outer boundary for a given temperature and pressure. The hydrostatic equations are solved for a number of layers of fixed composition and mass and lead to a radius corresponding to the mass of the planet considered. The EOSs considered in this work for the envelope are the SCVH (Saumon et al. 1995), which is used to reproduce the results of previous work, and the reevaluated CMS ab initio EOS (Mazevet et al. 2022), which is used to produce our results. The water ab initio EOS (Mazevet et al. 2019) is used to describe the core, as is the MgSiO3 EOS from Mazevet et al. (2019).

For both H-He EOSs, the additive rule for mixing hydrogen and helium is used. Recent works, such as those of Chabrier & Debras (2021) and Howard & Guillot (2023), have shown that corrections to the ideal mixing should be taken into account. Most notably, Chabrier & Debras (2021) suggest that by adapting the hydrogen EOS, it is possible to take into account the interactions between hydrogen and helium. However, as pointed out in Mazevet et al. (2022), such an ad hoc procedure incorporates the intrinsic differences arising from either the ab initio results or the fitting procedure used to produce the EOS. As this introduces an uncontrolled parameter into the evaluation of the H-He EOS, we willingly chose not to adopt this approach.

Research on the internal structure of Jupiter has experienced a surge in recent years, marked by significant contributions from missions such as Juno. By integrating insights from these missions with formation models, there has been a notable increase in the complexity of interior models (Vazan et al. 2016; Venturini et al. 2016; Wahl et al. 2017; Debras & Chabrier 2019; Valletta & Helled 2019; Ormel et al. 2021; Debras et al. 2021; Miguel et al. 2022; Helled et al. 2022; Militzer et al. 2022; Idini & Stevenson 2022). In order to fit the higher-order gravitational moments (above J2), models proposing inhomogeneous interiors with dilute cores and/or superadiabaticity are being proposed. We however chose to remain with a two-layered approach: a homogeneous envelope and a solid core, which reduces the number of uncontrolled parameters. This approach appears to be validated by the latest results of Bloot et al. (2023). Through the comparison of a homogeneous and a non-homogeneous interior model applied to exoplanets, these authors show that the two are degenerate with the current observations. They also show that the differences in bulk metallicity are not substantial, and find little difference between the resultant estimations of the core mass in most cases. Given these findings, and the fact that the homogeneity of the interior is ill constrained with the current observations, we do not consider it in the remainder of this work. It should, however, be kept in mind that by using a two-layer model, we are formally representing only one example within a range of possible interior models and are neglecting the possibility that exoplanet cores could be diluted, as observed in the case of Jupiter.

Fig. 1 Pressure–temperature profiles for a 1 MJ planet across a range of stellar irradiation values. The red crosses indicate the point at which the interior model and atmosphere model are linked. The thicker zones delimited by black boxes correspond to the onset of convective dynamics in the atmosphere continuing down to the interior.

2.3 Joining models

We chose to link the atmosphere and the interior in a self-consistent way. In order to do this, we link on pressure (P_link), temperature (T), average molecular mass (µ), and gravity (𝑔).

The linkage pressure (P_link) is a function of the atmosphere’s input properties. We can consider two regimes in the atmosphere, radiative and convective, which can be identified in Figure 1 as the approximately isothermal regions for the former and the adiabatic regions for the latter. P_link approximately corresponds to the radiative–convective boundary. In the cases where there ismore than one radiative region, P_link is defined by the deepest radiative–convective boundary, at higher pressure. To correctly join the interior model to the atmosphere model, we need to ensure that we join uniquely within the convective region. As the interior model is calculated considering an isentropic profile, P_link is in practice either equal to or at a slightly higher pressure than the radiative–convective boundary pressure. Exo–REM determines the convective regions of a given profile by checking that the temperature gradient is greater than the adiabatic gradient. Once a convective region is determined, we can then use the uppermost temperature and pressure values within the given region to run Exoris. The deeper, most convective regions are depicted in Figure 1 by the black boxed zones with crosses marking the linkage pressure. As the irradiation temperature decreases below the approximate value of the intrinsic temperature, we see that the atmosphere is convective at lower pressures. Once Exoris has computed a radius for this given pressuretemperature point at the chosen planetary mass, the gravity can be evaluated. As discussed in Section 2.1, for Exo–REM, the gravity at 1 bar is an input parameter and the gravity is evaluated using the HSE throughout the atmosphere. Therefore, the gravity value at 1 bar needs to be adjusted to ensure it matches the gravity at Plink, the top pressure point of Exoris, located deeper in the atmosphere. Consequently, Exo–REM and Exoris are run iteratively until the gravity at Plink is matched.

The average molecular mass is linked in a direct way by adapting the mass ratio of helium: $$X_{H_2}+Y_{He}+Z_{Heavy}=1.$$(2)

Considering the interior mix as $$X_{h_2}+{\left(Y_{He}\right)}_{int}=1,$$(3) $${\left(Y_{He}\right)}_{int}={\left(Y_{He}+Z_{Heavy}\right)}_{atm,}$$(4)

which can be linked to the average molecular mass given by the atmosphere at the P_link using $$\frac{1}{<\mu >}={\displaystyle {\displaystyle \sum }_i}\frac{X_i}{{\mu }_i}.$$(5)

The hydrogen fraction is derived from the atmospheric model when considering the mixing ratio of various species. Hence, heavier elements found in the atmosphere are also accounted for in the interior envelope. The current formulation paves the way for future improvements by including, for example, the EOS of water as a proxy for heavier elements.

For a converging model, only a few iterations are required, although this varies across the parameter space. The convergence criterion is set for an error of below 0.2% for each parameter.

To complete the analysis presented in this work, grids of individual models are constructed along the various physical dimensions (M, T_int, T_irr, Z, M_c). Individual models within the grid have different output values (e.g. radius and spectrum) along these dimensions but are built without altering these input parameters. Table 2 gives the layout of the grid, including the grid limits and the step used.

2.4 Model comparison

An overview of the features available in the models used here can be found in Table 3. Marley et al. (2021) is the first of a series of papers, and Karalidi et al. (2021) is the second paper in this series and adds disequilibrium chemistry. Clouds can be found in the latter Morley et al. (2024) model. The models are for the non-irradiated case; they are developed for emission spectra and do not provide transit spectra.

Using the public evolution grid from Phillips et al. (2020), as well as that from Marley et al. (2021), we can compare the results given by our models for the non-irradiated case. Figure 2 illustrates a comparative analysis of planetary models within the 1–5 M_J range. We conducted this comparison using both the ab initio EOS (Chabrier et al. 2019) and the SCVH EOS (Saumon et al. 1995), aiming to provide a comprehensive description of the variations between these models.

It is noteworthy that when directly comparing the calculations (non-interpolated), the disparities between the models and our findings remain modest, and below 5% in most of the evaluated range. Interestingly, our results consistently exhibit slightly larger radii at 200–400 K, regardless of the EOS used. Above 400K, our outcomes align closely with those of Phillips et al. (2020). However, when considering lower masses, a noticeable divergence emerges between our findings and those of Marley et al. (2021) as temperatures increase. While differences in atmosphere models could contribute to the observed variations, it is important to acknowledge their relatively minor impact compared to the substantial discrepancies arising from distinct EOS selections. For the remainder of this work, we use the ab initio EOS for effective temperatures above 400 K.

The generation of complete grids along each dimension presents a substantial computational and temporal burden. To address this, we propose leveraging multidimensional interpolation tools such as the linear barycentric interpolation function from the Python SciPy interpolation toolbox (LinearNDInter- polator). Careful consideration is essential when using interpolations within a grid, as inaccuracies at specific points or grid boundaries can disproportionately influence the resulting values.

In our verification process, we assess the reliability of the interpolation tool by comparing our interpolated results with those of Phillips et al. (2020), using a dataset computed at masses distinct from those within our target grid. Remarkably, the interpolation yields highly comparable results and the rightmost graphs of Figure 2 are comparable to the equivalent graphs on the left, which use actual grid points.

The above comparison shows the overall validity of our model compared to currently available models, this for the two different EOSs tested here. While this comparison is made for the non-irradiated case, there is no difference in the methodology used when adding stellar irradiation. The capacity to interpolate within the grid allows for more flexible use and reduced overall computation time when evaluating planetary parameters as done throughout the present work.

3 Stellar irradiation

Considering stellar irradiation, we can define the effective temperature of a planet as $$T_{eff}^4=\left(1-A_{Bond}\right)\cdot T_{irr}^4+T_{int}^4.$$(6)

Planets cool over time by evacuating internal energy into space and reducing their entropy (S). This is summarised by the following equation, $${\displaystyle {\int }_{M_{core}}^{M_p}T\frac{\text{\hspace{0.17em}d}S}{\text{\hspace{0.17em}d}t}}\text{\hspace{0.17em}d}m=-4\pi \sigma R_p^2\left(T_{eff}^4-\left(1-A_{Bond}\right)\cdot T_{irr}^4\right).$$(7)

By assuming an adiabatic (i.e. isentropic) envelope and by neglecting the contribution of the upper atmosphere (i.e. the radiative region) in the total entropy, the thermal evolution is given by $$\frac{\text{d}t}{\text{\hspace{0.17em}d}S}=\frac{{\displaystyle {\int }_{r=0}^{r=R_p}T}(r)\rho (r)r^2\text{\hspace{0.17em}d}r}{\sigma R_p^2{T}_{int}^4}.$$(8)

For the remainder of this work, the Bond albedo (A_Bond) derived from Exo–REM is approximately zero, assuming an M star for the light source (BT-Settle spectrum at 3500 K from Allard et al. 2012) and a cloud-free atmosphere. As such, for irradiated planets, the equilibrium temperature is approximately equivalent to the irradiation temperature.

In the context of a given pressure–temperature profile for a planet, distinct radiative and convective regions become discernible. The location of the radiative-convective boundary, which demarcates these regions, hinges upon the equilibrium between stellar irradiation and the internal flux associated with planetary cooling. The result of this is shown in Figure 1, where an array of profiles is depicted, each characterised by a different irradiation temperature. Heightened irradiation extends the reach of the radiative domain to higher pressures within the atmosphere, a trend consistently revealed in Figure 3 within the (T_int;T_irr) domain. This same figure shows how the radiative- convective equilibrium pressure (RCEP) evolves as a planet cools. For non-irradiated planets, there is little to no evolution of the RCEP with intrinsic temperature. Above approximately 400 K in T_irr, there is significant burying of the RCEP for decreasing intrinsic temperatures, equivalent to a planet ageing. This result helps explain why non-irradiated models are valid for fixed linkage pressures between atmosphere and interior and why irradiated models need variability in this pressure. As the RCEP is driven to lower strata due to irradiation, the validity of certain assumptions regarding atmosphere models is put into question. Notably when the atmosphere model relies on the ideal gas law, significant deviation can arise at pressures above 10³ bars. While we do not treat this in the present work, it remains important to ensure that the density does not deviate too greatly between the interior model and the atmosphere model. This can be the case for very old irradiated planets with T_int < 150 K, which Figure 4 shows to be a sparsely populated region considering current observations. Furthermore, models with fixed gravity throughout the atmosphere can suffer from deviations of up to 10% in gravity between the top and bottom of the atmosphere (1bar to RCEP). The adaptive linking implemented in this work accounts for this variability, as discussed in Section 2.3.

As the heat penetrates further into the atmospheric layers, it notably elevates the entropy of the selected adiabatic gradient for the internal structure. Consequently, this upshifted gradient engenders a seemingly decelerated cooling rate for the planet. This is where a requirement emerges to establish an intricate relationship between age, mass, and radius within an irradiated evolution model for Jupiter-like exoplanets.

The self-coherent linking of the gravity between the atmosphere model and the interior model requires that we solve Equation (7) on an irregular grid in entropy. To complete this, Equation (7) can be rewritten as Equation (8) by considering a negligible contribution from the isothermal region of the core, where the entropy is variable. This is shown in Mordasini et al. (2012); the luminosity of the core is negligible during the evolution phase of a planet when compared to the internal luminosity. Mordasini et al. (2012) defines the internal luminosity as $$T_{int}^4=\frac{L_{int}}{4\pi \sigma R^2}.$$(9)

We solve Equation (8) using a Runge Kutta 4 numerical method along an interpolated regular grid in T_int for a given set of planetary parameters. We use the specific entropies of young planets under the assumption of a hot start from Marley et al. (2007) as initial entropies. A rigorous exploration of the impact of the initial entropy should be undertaken for the case of young giant planets whose ages are estimated to be below the Kelvin- Helmholtz timescale. Using equation 11 from Guillot & Gautier (2015), we find a timescale of ≈700 Myr for 51 Eridani b using the best retrieval from Section 7.

In their study, Thorngren et al. (2019) established a connection between stellar irradiation and intrinsic temperature to explain the radius inflation of hot Jupiters by building upon their earlier research (Thorngren & Fortney 2018). This earlier work established a link between the heating efficiency of incoming stellar radiation and its role in our understanding of the observed radius of a group of transiting planets. These authors used the Thorngren et al. (2016) mass–metallicity relationship to further increase the Bayesian heating efficiency fit. They propose a Gaussian fit between the model-based intrinsic temperature and the measured equilibrium temperature and radius. We chose to derive, via interpolation, intrinsic temperature values for planets within the NASA Exoplanet archive, taking into account the mass–metallicity trend identified in Thorngren et al. (2016).

Figure 4 shows how the intrinsic temperature varies with observables. The intrinsic temperature values for each planet are obtained by interpolating from the grid of models. We identify a link between the intrinsic temperature and the equilibrium temperature. Figure 4a shows that increased irradiation leads to an increased intrinsic temperature following a quasi-linear trend given by Equation (10). $$T_{int}=0.16\times T_{irr}+\mathrm{44.6.}$$(10)

From a modeling perspective, this can be interpreted as a warming of the internal adiabat. Higher-mass planets exhibit a broader range of intrinsic temperatures, as shown in Figure 4b, likely due to their slower cooling rates and longer timescales for evacuating intrinsic energy compared to lower-mass planets. Figure 4c shows that the intrinsic temperature that we derive is linked to stellar age. We would expect that a planet with a lower intrinsic temperature would be older. We observe that planets with higher irradiation temperatures maintain higher intrinsic temperatures. Finally, Figure 4d shows that the intrinsic temperature increases for more inflated planets. There is a potential plateau that appears around the 400 K region, where the radius is more sensitive to the intrinsic temperature. The stand-out samples far above the 400 K intrinsic temperature region should be studied further.

Comparing our results with those of Thorngren et al. (2019), we observe no significant change in the trend of the intrinsic flux with increasing irradiation. In addition, the derived values for the intrinsic temperature are generally lower than those found by these authors. On the lower end of the irradiation temperature range, there is relatively good agreement between our results and theirs. However, as shown in Figure 4a, the majority of planets have a predicted intrinsic temperature of between 200 and 400 K, which is as low as half of that shown in the fit by Thorngren et al. (2019). This result corresponds to a heat transfer of ≈ 0.09% to the interior for a HD 2O9458b-like planet (Henry et al. 2000), which is in agreement with the results of Guillot & Showman (2002) for the flux transfer required to explain the radius of HD 209458b, taking into account energy dissipation in the interior. Compared to the efficiency of heat transfer required in Thorngren et al. (2019), this makes the observed inflated radii easier to reproduce by physical processes such as ohmic dissipation (Batygin & Stevenson 2010) or advection of heat by the atmospheric circulation (Tremblin et al. 2017). We believe that Equation (10) should serve as an approximate temperature to be used when running global climate models for hot Jupiters, as is done for instance in Komacek et al. (2022).

Fig. 2 Comparison between results from this work and available work using similar physical properties. Top row (2): (a) Comparison between the radius obtained in this work using the ab initio EOS of Mazevet et al. (2022) and the results of Phillips et al. (2020), who use the ab initio EOS of Chabrier et al. (2019). (b) Comparison between the radius obtained using the developed interior-atmosphere model using the SCVH EOS with the results of Marley et al. (2021), who also use the SCVH EOS (Saumon et al. 1995). (c) Our interpolated ab initio EOS grid compared to Phillips et al. (2020). Bottom row (2a, b, and c): Corresponding percentage differences.

Fig. 3 Schematical representations of the radiative-convective boundary (in pressure height) for varying intrinsic and irradiation temperatures for 1 MJ planet.

Fig. 4 Intrinsic temperatures for planets between 0.5 and 10 Jupiter masses. (a) Intrinsic temperature with equilibrium temperature showing a linear increasing trend, with stellar age indicated qualitatively by point size. (b) Intrinsic temperature with mass, showing a greater spread at higher mass values, stellar age indicated qualitatively by point size. (c) Intrinsic temperature with stellar age, showing planets cooling with age, irradiated planets maintain a higher intrinsic temperature. Interpolated trends are shown for 1000, 1600, and 1900 K irradiation temperatures. (d) Intrinsic temperature with radius showing the contribution to the radius from the irradiation. The list of planets used can be found in the Appendix A.1.

4 Breaking degeneracies

To assess planetary parameters based on observable data, a commonly employed method is a statistical approach, such as Bayesian inference. Nevertheless, when faced with limited observables, it becomes intriguing to explore whether it is feasible to resolve parameter degeneracies using physical models. In this context, we focus on two specific parameters: metallicity and core size. Our objective is to investigate whether discrepancies in these parameters, within the (T_int; T_irr) space at a fixed mass, result in discernible variations in planetary radii.

4.1 Metallicity

By linking the molecular mass of the atmosphere in the convective region with the helium fraction of the interior, we can adjust the interior structure to take into account varying metallic- ities. We can then compare the (T_int;T_irr) space between different metallicities and uncover regions where the irradiation, age, and radius are sufficient variables to distinguish the planet’s metallicity.

Figure 5a illustrates that a planet with ten times the solar metallicity and a very high intrinsic temperature may exhibit a larger radius than an equivalent planet with solar metallicity. In contrast, the right panel shows that a planet with 30 times the solar metallicity consistently maintains a smaller radius across the range of intrinsic and irradiation temperatures. The increase and subsequent decrease in the radius with metallicity is better shown in the lower panel of Figure 5. As we increase the temperature, the radius inflates for all cases, as expected. When increasing metallicity at fixed temperature, the radius also increases, but beyond a certain metallicity (≈0.5) the radius decreases again. We therefore observe two regions, one where an increasing metallicity leads to an increase in absorption and hence the thermal effects dominate, and a second region where the molecular mass dominates and gravity forces the radius to contract. We therefore observe a ‘C’ shape in the (Z;T_int) space. Thus, for a fixed radius, higher metallicities necessitate a lower intrinsic temperature because the radiative–convective boundary occurs at a lower pressure due to the increased opacity and greenhouse effect. However, beyond a certain metallicity threshold, the effect of the mean molecular mass on the atmospheric scale height dominates the effect of the opacity. Hence, higher intrinsic temperatures are required to maintain the radius and counteract the effects of gravity. The threshold is temperature dependent, 0.2 dex for 200 K and 1 dex for 800 K.

Fig. 5 Radius difference as a percentage between a 1 MJ planet with (a) 10 times or (b) 30 times solar metallicity and an equivalent planet with 1 times solar metallicity. (c) Iso-radii lines in the metallicity–intrinsic temperature frame.

4.2 Core size

Exploring the properties of cores within gaseous planets is a topic of significant interest, and fundamental questions persist regarding their distribution within these planets and their contribution to both mass and radius. The layout of the core and even the presence of one, is a topic that has received significant attention where formation mechanisms are considered (Pollack et al. 1996; Mordasini et al. 2012). However, when studying giant exoplanets, the primary question that arises pertains to whether variations in core characteristics have a discernible impact on observable parameters, such as the planet’s radius. To explore this topic, we introduced an additional dimension to our analysis by considering different core sizes, from 1 to 40 times Earth mass (M_E).

In Figure 6, we observe that at low intrinsic temperatures, the size of the core appears to be a degenerate parameter, showing minimal to negligible impact on the radius of a 1 M_J planet. However, as we examine considerably younger and warmer planets with substantial radial extent, we notice an increasing influence of the core on the planet’s radius. In fact, if a larger proportion of the planet’s mass resides within the core, its ability to expand with rising temperature becomes constrained. This observation is supported by the slight tilt in the iso-radii differences.

Interestingly, we observe a weaker effect on the radius with increasing irradiation temperature compared to increasing intrinsic temperature. Consequently, the core of young planets has the potential to become a quantifiable parameter, so long as they have a sufficiently high intrinsic temperature. Exploration of this realm of the parameter space is poised to offer valuable insights into alternative interior designs, such as those featuring diluted cores.

5 Grid interpolation and retrieval

As discussed Section 2.4 and shown in Figure 2c, it is possible to obtain the modelled radius of a planet at any point within the grid. Hence the radius is a continuous parameter. In order to efficiently develop grid retrieval tools, it is equally interesting to be able to interpolate the spectra (transit or emission) of a simulated planet with a given set of input parameters, defined as follows: $$Planet=\left({\text{M}}_p,{\text{T}}_{int},{\text{T}}_{irr},\text{Z},{\text{M}}_c\right).$$(11)

It is not possible to do this simply using the LinearNDIn- terpolator interpolator as for the radius. The issue is that this will lead to leakage between different wavelengths and the resulting spectrum will be heavily influenced by the distance between the old and new wavelengths, as well as the other planetary parameters. One solution is to create an array of interpolators, following: $$f_{spectra}(Planet)=\left[f_{{\lambda }_i}(Plane)\mid {\lambda }_i\in \Lambda \right].$$(12)

For each wavelength of the generated spectra (Λ), we create a unique interpolation function that takes the planetary parameters as input. It is no longer possible to leak information from one wavelength into another, as the unique interpolators are independent of one another. The drawback of this method is that it is essential to re-interpolate each point on the grid to a desired wavelength range before obtaining a new spectrum for a given planet.

Using Markov chain Monte Carlo (MCMC) sampling on a continuous grid with the emcee Python library, we aim to ascertain the optimal planetary parameters that can elucidate an observed spectrum. To evaluate this approach, we introduce Gaussian noise (σ = 0.15× < Flux >) to a spectrum extracted from the grid, and therefore not interpolated.

The selected spectrum pertains to an artificial planet for which we fixed the parameters to 1 Jupiter mass, featuring intrinsic and irradiation temperatures of 800 and 100 K, respectively, alongside solar metallicity and a 10 M_E core. The resolution used is chosen so as to prepare for the test case in Section 7 (R ≈ 33). By considering a H-band spectrum and irradiation temperature as observable factors, our assessment involves sampling the spectrum while treating the mass, intrinsic temperature, metallicity, and core mass as variables. We consider that the planet is placed at 10 parsecs from the observer. We adjust the flux level to account for this and the model radius using: $$F_{observed}=F_{emitted}\cdot \frac{R_p^2}{D_{Earth-Planet}^2}.$$(13)

The general expression of the likelihood function that we use in the MCMC is given as follows: $$\begin{array}{l}{\chi }^2={\displaystyle \sum _j^{n_{\text{insr\hspace{0.17em}}}}\frac{1}{n_j}}{\displaystyle \sum _i^{n_j}\frac{{\left(f_{{\lambda }_{ji}}(\text{\hspace{0.17em}Planet\hspace{0.17em}})\right)-{\text{\hspace{0.17em}spectra\hspace{0.17em}}}_{{\lambda }_{ji}})}^2}{{\sigma }_{spectr{a}_{{\lambda }_{ji}}}^2}}\\ +{\displaystyle \sum _k^{n_{obs}}{w}_k}\frac{(\text{\hspace{0.17em}Planet\hspace{0.17em}}{-ob{s}_k)}^2{}_{{\sigma }_{ob{s}_k}}}{{\sigma }_{ob{s}_k}}.\end{array}$$(14)

We consider two distinct terms, the first consisting of the double sum. This takes into account the fact that a spectrum can be a mosaic of spectra from multiple instruments. In this case, we propose to allow each instrument to be treated as an equally weighted observable. The second term is the simple sum that takes into account one or more physical characteristics of the planet, such as a mass measurement using radial velocity or a radius using the average transit depth. The proposed wk term is given to allow greater importance to be given to certain observed characteristics.

For the example studied here, Equation (14) is reduced to the first term. We use uniform priors, except at the grid limits where we impose an infinite barrier. The priors are chosen as follows: $${\text{M}}_p\in [0.5;3]\left|{\text{T}}_{int}\in [500;1100]\right|Z\in \left[{10}^{-2};10\right]\mid {\text{M}}_c\in [2;35].$$(15)

The outcome of this evaluation is illustrated in Appendix C.1, which presents a corner plot of the MCMC-based interpolated grid retrieval. The radius distribution is added by interpolating the four other parameters (M_p,T_int,Z,M_c) for each walker and iteration. Figure 7 shows the retrieved interpolated spectrum compared to the degraded and original spectra.

The results indicate a good capacity to retrieve a spectrum using interpolators and MCMC sampling; they reveal a relatively accurate retrieval of the intrinsic temperature, metallicity, and mass. The core mass is retrieved correctly, albeit with a relatively broad posterior distribution. This indicates that future core-mass estimations will at best be given with around a 5 M_E uncertainty for a spectrum with similar characteristics to the one shown in Figure 7.

Fig. 6 Radius difference as a percentage between a 20 ME core and 5 ME core.

Fig. 7 Spectra of the overall best fit with 100 sample best fits overplotted on an actual grid spectrum and a downgraded grid spectrum in H-band.

Fig. 8 Posterior distributions of retrieval of the WASP-39 b JWST NIRSpec-G395H spectrum: mass and radius (blue), mass and spectrum (red).

6 Case study: WASP-39 b, a transiting hot-Saturn

WASP-39 b is an inflated Saturn-like planet orbiting a G-type star with a period of P~ 4d (Faedi et al. 2011). Radial velocity and transit photometry measurements as described by Faedi et al. (2011) indicate a mass of 0.28 ± 0.03 M_J and a radius of 1.27 ± 0.04M_J, suggesting an inflated planet. Ahrer et al. (2023) detected the presence of CO₂ from its strong 4.4 µm band from JWST NIRSpec-G395H observations, as part of the Early Release Science (ERS) program. This detection suggests a relatively high metallicity (~10×solar metallicity). We propose to use this JWST transit spectrum in an exemplary application of the model. Two approaches are adopted: in one instance, we assume no spectrum and uniquely mass and radius measurements, and hence Equation (14) is reduced to the second term; while in another, we use the JWST ERS spectrum with the radial-velocity mass measurement. In this case, Equation (14) includes the first term. The objective is to understand the impact that spectral information has on the constraint of parameters such as the intrinsic temperature derived for example in Figure 4. We do not include a radius measurement as it is given within the transit depth at the JWST NIRSpec-G395H wavelengths. Radius measurements from, for example, Faedi et al. (2011) correspond to the radius in the visible and near-infrared (R and Z bands).

Figure 8 shows the posterior distributions for the two scenarios. The retrieved values are given in Table 4. We see that without the spectrum, the posterior distributions show a highly degenerated ensemble of potential parameters. Most notably, there are no constraints on the metallicity or the core mass. Possible values for the intrinsic temperature considering a Gaussian posterior with 1 σ variation are between 150 and 320 K. For the case with the spectrum, we ran the retrieval based on our grid of models using a chemistry at equilibrium. Stronger constraints are obtained for the metallicity and core mass. A relatively high metallicity of about ten times solar is necessary to adequately fit the CO₂ band as shown in the best fit of Figure 9. A high T_int of ${297.39}_{-16.9}^{+8.95}\text{\hspace{0.17em}K}$ is required to reproduce the inflated radius, given the equally high metallicity, but also to convert CH₄ into CO and CO₂ leading to the observed CO₂ spectral signature.

The core proposed is relatively large, namely of between 20 and 30 times Earth’s mass, and therefore represents up to 30% of the total planetary mass. To match the higher atmospheric metallicity, the helium mass fraction in the interior is approximately 0.50. As such, considering a solar-like helium mass fraction of 0.246 (Lodders 2019), the mass of heavy elements (core mass + excess helium) is approximately 45% of the total mass of the planet. Bloot et al. (2023) finds a core mass of 0.88 times Earth’s mass using a homogeneous model and 0.90 times Earth’s mass using a non-homogeneous core model. JWST data were not used in either case. Furthermore, Thorngren & Fortney (2019) use a bulk metallicity of 22% of the planet’s mass to derive an atmospheric metallicity of 40.51 times solar, again without using JWST data. We now know the atmospheric metallicity to be lower, and within the framework of their approach this could change how well mixed the planet is. In all cases, we find a core mass that is substantially higher with the JWST data. Without the JWST data, the result is very degenerate and would include both possibilities. Figure 9 shows the fit over-plotted on the NIRSpec-PRISM spectrum from Rustamkulov et al. (2023). The CO₂ signature is shown to be well fitted by the model. Discrepancies arise at shorter wavelengths due to the absence of clouds in the current model. The absence of photo-chemistry and SO₂ in the model explains the missing absorption band around 4 µm. Nevertheless, this has no impact on the overall fit. However, we are able to reproduce the sodium signature at around 0.5 µm. We conclude that adding spectroscopic measurements enables us to break the degeneracies between metallicity, intrinsic temperature, and core mass. The fact that the metallicity and core values found for the spectrum-constrained case are outside the 68% region of the unconstrained posterior distribution indicates that the spectrum provides essential information that is not included with mass and radius alone. The metallicity derived here using the transit spectrum is in agreement with the approximately ten times solar value found by the extensive range of models applied to this planet and spectrum (Powell et al. 2024; Constantinou & Madhusudhan 2024; Carone et al. 2023; Feinstein et al. 2023; JWST-ERS Team et al. 2023; Alderson et al. 2023). Using the same spectrum and the retrieval code TauREx from Al-Refaie et al. (2021) - and used for example in Panek et al. (2023) -, we find a metallicity ranging from 11.99 to 12.85 times solar. This result is also in agreement with the derived atmospheric metallicity found using the developed model. Alderson et al. (2023) find a wide range of intrinsic temperatures using three atmosphere models: Atmo (100K), PHOENIX (400K), and PICASO+Virga (100K). This could be due to the lack of coupling between the atmosphere model and the interior model, to which particular attention was paid in the present study.

Fig. 9 Fit of the WASP-39 b JWST NIRSpec PRISM spectrum (Rustamkulov et al. 2023) with contributions from the species considered. This fit was found for an intrinsic temperature of ${297.39}_{-16.9}^{+8.95}\text{\hspace{0.17em}K}$ and a metallicity of ${1.16}_{-0.09}^{+0.06}\text{\hspace{0.17em}dex}$.

7 Case study: 51 Eridani b, a young remote massive-Jupiter

51 Eridani b (hereafter 51 Eri b) is a massive-Jupiter planet orbiting the young (20 ± 6 Myr) F0 type star 51 Eridani. The high uncertainty on its age makes it difficult to interpret the photometric data in order to derive the planet’s physical parameters. 51 Eridani b was dentified in 2015 in high-contrast imaging on the basis of a CH₄ absorption in its near-infrared spectrum (Macintosh et al. 2015) using the Gemini Planet Imager (GPI, Macintosh et al. 2014). In a subsequent paper, Rajan et al. (2017) used the same GPI Integral Field Spectrograph (IFS) data for 51 Eri b, finding an effective temperature ranging between 605 and 737 K, solar metallicity, and a surface gravity of log(g) = 3.5–4.0. Samland et al. (2017) show a supersolar metallicity at $[\text{Fe}/\text{H}]={1.03}_{-0.11}^{+0.10}\text{dex}$, a ${759.51}_{-21.97}^{+21.34}\text{\hspace{0.17em}K}$ effective temperature, and a surface gravity of $\log (\text{g})={4.26}_{-0.25}^{+0.24}\text{\hspace{0.17em}K}$. More recently, Brown-Sevilla et al. (2023) used photometric data obtained with the Spectro-Polarimetic High-contrast imager for Exoplanets REsearch (SPHERE, Beuzit et al. 2019) at the VLT to constrain its gravity to log(g)=4 05 ± 0 37 and its radius to 0 93 ± 0.04 R_J. Unlike the previous case study described in Sect. 6, only spectrophotometric data are available. We can then reduce Equation (14) down to the first term. This equally means that we have fewer observational constraints to infer the structure of the exoplanet.

Using the same observation as in Brown-Sevilla et al. (2023) and Rajan et al. (2017), we extracted a YJH spectrum using SPHERE/IFS (Claudi et al. 2008) and the H spectrum using two GPI/IFS epochs. The observing conditions of the observations can be found in Table D.1. The extraction of a spectrum can be complicated and the uncertainties associated with this process must be kept in mind when subsequently running retrievals.

The pre-processing steps consist of corrections to the sky, the flat field, and the background, as well as bad pixel interpolation and wavelength calibration. These steps use the standard pipeline of the High Contrast Data Center, formerly known as the Sphere Data Center (Delorme et al. 2017). Regarding GPI, every raw frame was dark subtracted, flat-fielded, cleaned of correlated detector noise, bad pixel corrected, and flexure corrected using the standard GPI Data Reduction Pipeline (DRP; version 1.6) documented in Perrin et al. (2014, 2016). The resulting set of frames was combined into a spectral datacube using Pyklip (Wang et al. 2015), which allowed a more accurate wavelength and PSF recalibration in addition to frame centring.

As a post-processing algorithm, we consistently use PACO ASDI (Flasseur et al. 2020a) for all datasets. PACO ASDI is an ADI-based algorithm that models the noise at a local scale in small patches. It relies on the learning of local temporal and spectral correlations of the noise to provide a reliable and statistically grounded signal-to-noise-ratio map and photometry confidence intervals. A full and detailed explanation of the algorithm can be found in Flasseur et al. (2018, 2020a,b) and its implementation on the Data Center and spectral optimization in Chomez et al. (2023).

After measuring the contrast spectrum of the planet with PACO, that is the ratio of the flux coming from the planet to the stellar PSF flux, we obtained a physical spectrum by multiplying this contrast by a synthetic spectrum of the star based on BT-NextGen AGSS09 atmospheric models (Allard et al. 2012). The stellar spectrum is normalised by 2MASS (Skrutskie et al. 2006) H magnitude. In addition to the H-band uncertainty, the final spectrum error budget also accounts for systematic uncertainties on the star PSF intensity (Wang et al. 2014), as well as for random uncertainties related to the PSF variations over the exposure and to PACO’s fit. The two GPI spectra were finally averaged to get the H-band spectrum used in the present section. A similar procedure was applied to SPHERE data. Error bars for both SPHERE and GPI are displayed at 1 σ.

Appendix E.1 shows the corner plot of the spectral retrievals using MCMC, the best fit is shown in Figure 10, and the retrieved values are given in Table 5. While this fit shows similar spectral features and an overall flux amplitude that is close to the observed spectrum, there are notable differences below 1.1 µm. These could be due to the presence of clouds, which are not included in the current version of the model. Charnay et al. (2018) show how clouds have an effect on emission spectra and most notably in the Y, J, and H bands. Appendix E.1 indicates a multi-nomial distribution with three mass solutions found, one of which is further degenerate in intrinsic temperature and metal-licity. All solutions are given in Table 5. In order to break the degeneracy, we include the age estimations found by solving Equation (8) in the corner plot. The three proposed ages are to be compared to the age estimate of the star of 20 ± 6 Myr proposed by Macintosh et al. (2015). We identify the higher-mass case $\left({3.13}_{0.04}^{0.05}{\text{M}}_J\right)$ as being the most compatible configuration. This configuration has a corresponding age of 19.5 ± 1.8 Myr and is hence in strong agreement with the age of the system. It is however important to mention that the planet and the star need not be coeval. However, Mordasini et al. (2017) show that a 5 M_J would take no longer than 10 Myr to finish its mass accretion with observations indicating less than 5 Myr. As such, when comparing the next closest age solutions (solutions 2 and 3) with the best match (solution 4), we see that the age of 19.5 Myr is a better match than 7–8 Myr.

As discussed in Section 3, the initial entropies used to evaluate the age are under the assumption of a hot start. The MCMC posterior distribution shows a core mass of ${31.86}_{-0.18}^{+0.32}{\text{M}}_E$ for the best matching age case. This result, within the models approximations, indicates the presence of a core and suggests a formation with core-accretion processes combined with the hot-start scenario. Such scenarios have been studied for example by Berardo & Cumming (2017) and Mordasini et al. (2017). Mordasini et al. (2017) show that it is possible for planets to undergo core accretion and have similar luminosities to those expected in both hot-start and warm-start scenarios. We explored the cold-start hypothesis for initial entropies, but all solutions found were below the estimated system age.

The radius of 1.25 R_J found here is inconsistent with the ‘New nominal’ radius found by Brown-Sevilla et al. (2023) of 0.93 ± 0.04R_J. The radius is a parameter deduced from the input parameters as defined by Equation (12). The model flux is scaled according to Equation (13). As such, the radius deduced through the model is directly linked to the flux scaling as well as the input parameters. This is not the case with purely atmospheric retrievals, where the radius becomes a free input parameter.

Fig. 10 Best fit of the observed spectrum of 51 Eridani b with 100 best fits.

8 Discussion

The two examples show the advantages of a coherent model such as the one proposed here. The internal links between the mass, radius, metallicity, core size, and intrinsic temperature allow physically coherent results to be derived. Trends in parameters, such as intrinsic temperature derived from mass, radius, and effective temperature, need to be treated with utmost caution. The case of WASP-39b exemplifies this; albeit a planet with a mass below those studied in Figure 4, it represents an interesting case where the metallicity and chemical composition derived from the spectrum drives the model towards a higher intrinsic temperature and helps explain the apparent inflated radius. WASP-39 b showed that deriving the intrinsic temperature as well as the internal structure from the mass and radius can lead to a wide range of possible chemical and thermal structures. Repeating this over the range of observed exoplanets can lead to increased uncertainties on overall trends. Using spectral data, it is possible to constrain the thermal and structural properties in a way that should pave the way for a re-evaluation of planetary statistics. Thorngren et al. (2019) use metallicity trends from Thorngren et al. (2016), which provides a constraint that corresponds to having a spectrum. However, the uncertainty related to this metallicity trend is large. This leads us to believe that it is best to derive the metallicity from the spectrum on a planet-wise basis. This will become possible with the range of planets observed with JWST (Greene et al. 2019; Beichman et al. 2019) and subsequently the dedicated Ariel mission (Tinetti et al. 2020).

The fits obtained for directly imaged planets need to rely on physical models that are self-consistent in order to correctly determine physical characteristics. Commonly used methods for decoupling the radius from the model will potentially lead to erroneous conclusions due to the contribution of nonphysical radii in adjusting the overall amplitude of the spectrum. The model proposed here has the potential to provide the most physically coherent parameters for directly imaged planets, even without prior knowledge of their age.

Four essential elements are missing from the current model, namely the C/O ratio, the H₂O EOS, clouds, and a diffuse core. The C/O ratio is famously hard to derive; it would be interesting to evaluate our capacity to constrain this parameter within a physically coherent model. The H₂O EOS is essential for exploring higher-metallicity planets, and especially Neptune and sub-Neptune-like planets, and might also have an impact on the distribution of heavy elements in exoplanet interiors should for instance metals be over-represented by the core mass. Clouds represent an important feature within the spectrum, changing the flux intensity between shorter and longer wavelengths and changing the thermal evolution of a planet. As we have shown, not including this parameter with a self-consistent model such as the one presented here can lead to physically incoherent results for other parameters. While it may be possible to account for differences in a parametric way, a self-coherent model would be in the spirit of the current work. Finally, incorporating diffused cores and examining constraints on the extent of the core in the context of an increasing amount of exoplanet data may provide better insights into this process. As pointed out in Section 1, diffuse cores provide better explanations for the Juno data than non-diffuse cores. While the impact on the radius is limited for older planets and those with lower intrinsic temperatures, warmer, younger planets might provide greater insight into the interior design. This addition to the model might also help us to better understand the distribution of heavier metals found within the core and the envelope, leading to a better understanding of the core mass.

9 Conclusions

We present a combined atmosphere and interior model, HADES (Heat Atmosphere Density Evolution Model), which will serve as a powerful tool for proposing potential structural configurations of Jupiter-like planets. We linked the atmosphere and interior in a self-consistent manner by converging the physical parameters at the interface. We subsequently reevaluated the intrinsic temperatures corresponding to the internal fluxes of a range of exoplanets. This led us to the conclusion that there is reduced need for ohmic dissipation or other mechanisms that lead to higher intrinsic temperatures to explain the observed radius. We show that, by using a grid of models, we can use Bayesian inference to derive potential parameters for observed exoplanets. The case of WASP-39 b demonstrates the importance of using spectra to constrain parameters and hence reinforces the importance of atmospheric surveys such as the one anticipated as part of the Ariel mission. The case of 51 Eri b showcases the use of a physically coherent model to correctly scale the flux with radius even if the fit is of lower quality. Our application to 51 Eri b also demonstrates the ability of the interior-atmosphere model to infer formation processes when a system’s age is available. Further work is required to complete the model, most notably by adding clouds.

Acknowledgements

We thank the anonymous reviewer for their valuable comments on the manuscript. We also thank Raphaël Galicher and Lucas Tein-turier (LESIA, Observatoire de Paris, France) for their time saving advice on MCMC usage and expert opinions on results. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COBREX; grant agreement n° 885593). This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This work was supported by CNES, focused on AIRS on Ariel, and by the Programme National de Planétologie (PNP) of CNRS/INSU, co-funded by CNES. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF - Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zürich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). This research has made use of 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.

Appendix A List of planets used

Fig. A.1 List of planets used to derive trends.

Appendix B List of planets used

Appendix C Test MCMC result

Fig. C.1 Corner plot of posterior distributions explaining downgraded grid test spectrum. Red cross-hairs indicate location of initial parameters of test spectrum.

Appendix D 51 Eridani observing logs

Appendix E 51 ׃ Eridani b MCMC result

Fig. E.1 Corner plot with posterior distributions of MCMC retrieval of the spectra of 51 Eridani b. Radius, log(g) and log(age) are derived parameters from the model. The red colored distribution indicates the most compatible parameters with the systems age.

References

All Tables

All Figures

	Fig. 1 Pressure–temperature profiles for a 1 MJ planet across a range of stellar irradiation values. The red crosses indicate the point at which the interior model and atmosphere model are linked. The thicker zones delimited by black boxes correspond to the onset of convective dynamics in the atmosphere continuing down to the interior.
In the text

Fig. 2 Comparison between results from this work and available work using similar physical properties. Top row (2): (a) Comparison between the radius obtained in this work using the ab initio EOS of Mazevet et al. (2022) and the results of Phillips et al. (2020), who use the ab initio EOS of Chabrier et al. (2019). (b) Comparison between the radius obtained using the developed interior-atmosphere model using the SCVH EOS with the results of Marley et al. (2021), who also use the SCVH EOS (Saumon et al. 1995). (c) Our interpolated ab initio EOS grid compared to Phillips et al. (2020). Bottom row (2a, b, and c): Corresponding percentage differences.

In the text

	Fig. 3 Schematical representations of the radiative-convective boundary (in pressure height) for varying intrinsic and irradiation temperatures for 1 MJ planet.
In the text

Fig. 4 Intrinsic temperatures for planets between 0.5 and 10 Jupiter masses. (a) Intrinsic temperature with equilibrium temperature showing a linear increasing trend, with stellar age indicated qualitatively by point size. (b) Intrinsic temperature with mass, showing a greater spread at higher mass values, stellar age indicated qualitatively by point size. (c) Intrinsic temperature with stellar age, showing planets cooling with age, irradiated planets maintain a higher intrinsic temperature. Interpolated trends are shown for 1000, 1600, and 1900 K irradiation temperatures. (d) Intrinsic temperature with radius showing the contribution to the radius from the irradiation. The list of planets used can be found in the Appendix A.1.

In the text

	Fig. 5 Radius difference as a percentage between a 1 MJ planet with (a) 10 times or (b) 30 times solar metallicity and an equivalent planet with 1 times solar metallicity. (c) Iso-radii lines in the metallicity–intrinsic temperature frame.
In the text

	Fig. 6 Radius difference as a percentage between a 20 ME core and 5 ME core.
In the text

	Fig. 7 Spectra of the overall best fit with 100 sample best fits overplotted on an actual grid spectrum and a downgraded grid spectrum in H-band.
In the text

	Fig. 8 Posterior distributions of retrieval of the WASP-39 b JWST NIRSpec-G395H spectrum: mass and radius (blue), mass and spectrum (red).
In the text

	Fig. 9 Fit of the WASP-39 b JWST NIRSpec PRISM spectrum (Rustamkulov et al. 2023) with contributions from the species considered. This fit was found for an intrinsic temperature of ${297.39}_{-16.9}^{+8.95}\text{\hspace{0.17em}K}$ and a metallicity of ${1.16}_{-0.09}^{+0.06}\text{\hspace{0.17em}dex}$.
In the text

	Fig. 10 Best fit of the observed spectrum of 51 Eridani b with 100 best fits.
In the text

	Fig. A.1 List of planets used to derive trends.
In the text

In the text

	Fig. C.1 Corner plot of posterior distributions explaining downgraded grid test spectrum. Red cross-hairs indicate location of initial parameters of test spectrum.
In the text

	Fig. E.1 Corner plot with posterior distributions of MCMC retrieval of the spectra of 51 Eridani b. Radius, log(g) and log(age) are derived parameters from the model. The red colored distribution indicates the most compatible parameters with the systems age.
In the text