Evidence of Three Mechanisms Explaining the Radius Anomaly of Hot Jupiters

The radii of hot Jupiters are still not fully understood and all of the proposed explanations are based on the idea that these close-in giant planets possess hot interiors. We approach the radius anomaly problem by adopting a statistical approach. We infer the internal luminosity for the sample of hot Jupiters, study its effect on the interior structure, and put constraints on which mechanism is the dominant one. We develop a flexible and robust hierarchical Bayesian model that couples the interior structure of exoplanets to their observed properties. We apply the model to 314 hot Jupiters and infer the internal luminosity distribution for each planet and study at the population level ({\it i}) the mass-luminosity-radius distribution and as a function of equilibrium temperature the distributions of the ({\it ii}) heating efficiency, ({\it iii}) internal temperature, and the ({\it iv}) pressure of the radiative-convective-boundary (RCB). We find that hot Jupiters tend to have high internal luminosity leading to hot interiors. This has important consequences on the cooling rate and we find that the RCB is located at low pressures. Assuming that the ultimate source of the extra heating is the irradiation from the host star, we illustrate that the heating efficiency follows a Gaussian distribution, in agreement with previous results. We discuss our findings in the context of the proposed heating mechanisms and illustrate that ohmic dissipation, advection of potential temperature, and thermal tides are in agreement with certain trends inferred from our analysis and thus all three models can explain aspects of the observations. We provide new insights on the interior structure of hot Jupiters and show that with our current knowledge it is still challenging to firmly identify the universal mechanism driving the inflated radii.


Introduction
Two decades of observational and theoretical exploration have revealed that the anomalously large radii of close-in transiting giant planets holds firmly (e.g. Laughlin et al. 2011;Weiss et al. 2013). The radii of hot Jupiters are larger than what is predicted by standard interior structure models . Observations reveal that there is a strong correlation between the observed radii and the stellar incident flux (e.g. Enoch et al. 2012), with a threshold around ∼ 2 × 10 8 erg s −1 cm −2 , corresponding to an equilibrium temperature of about 1000 K (Demory & Seager 2011;Miller & Fortney 2011), below which the physical mechanism becomes inefficient. Sestovic et al. (2018) further demonstrated that the inflation extent is mass dependent, where the planets with the largest anomalous radii have masses less than ∼< 1 M J .
There has been a lot of investigations to explain the discrepancy between the observations and theoretical models. The proposed mechanisms can be divided into two categories: (i) slowing down cooling and contraction or (ii) depositing extra heat into the interior. Burrows et al. (2007) showed that slowing down the cooling and thus delaying contraction can be achieved by in-Send offprint requests to: Paula Sarkis, e-mail: paula.sarkis@space.unibe.ch creasing the atmospheric opacity. Another way to delay contraction is to reduce the heat transport efficiency due to compositional gradients (Chabrier & Baraffe 2007).
It is well established that heating up the interior of the planet increases its entropy and thus its radius (Arras & Bildsten 2006;Marleau & Cumming 2014). The source of heat is still not constrained and possible sources could be tidal dissipation of an eccentric orbit (e.g. Bodenheimer et al. 2001), advection of potential temperature, which is a consequence of the strong stellar irradiation (Tremblin et al. 2017;Sainsbury-Martinez et al. 2019), or dissipative processes powered by the stellar irradiation flux. The latter has received a lot of attention and the mechanism to transport fraction of the stellar incident flux into the interior is still an open question. One mechanism is atmospheric circulation, which leads to thermal dissipation of kinetic energy into the interior . Another mechanism is ohmic dissipation (Batygin & Stevenson 2010;Batygin et al. 2011;Perna et al. 2010a;Huang & Cumming 2012;Wu & Lithwick 2013;Ginzburg & Sari 2016), where the irradiation drives fast winds through the planets magnetic fields, giving rise to currents that dissipate ohmically in the interior. Other mechanisms are thermal tides (Arras & Socrates 2010) and the mechanical greenhouse (Youdin & Mitchell 2010). P. Sarkis: Inflated Hot Jupiters Some of these mechanisms come with a lot of approximations and uncertainties. For example, an important uncertain parameter in atmospheric circulation, ohmic dissipation, and advection of potential temperature is the wind speeds and the effect of magnetic drag in damping the winds (Perna et al. 2010a,b). Another uncertainty is how deep the wind zone extends, which is important to constrain the pressures at which the extra heat should be dissipated. Wu & Lithwick (2013) illustrate that if the wind zone is at shallow pressures, then a significantly larger heating efficiency is needed to achieve the same interior heating, compared to heating at deeper pressures. Komacek & Youdin (2017) argued that the extra heat should be deposited in the convective layers or at the radiative-convective-boundary (RCB), otherwise it will be re-radiated away. Huang & Cumming (2012) deposited the extra heat in the radiative layers and as a consequence showed that the RCB moves to deeper pressures. Fortney et al. (2007) showed that RCB is located at pressures of 1000 bar, where little is known about the wind speeds at such deep pressures. However, the Fortney et al. (2007) models were developed for irradiated planets and do not consider the high internal entropy that hot Jupiters are believed to possess.
All the mechanisms proposed have been tested and applied on single or a handful of planets.
It is yet to be demonstrated that these mechanisms can explain the radii of all the observed hot Jupiters. Within this context, in this paper we approach the radius inflation problem from a statistical point of view, similar to the approach of Thorngren & Fortney (2018) (hereafter TF18). We do not model any of the previously mentioned mechanisms but rely solely on the interior structure model and atmospheric model. We develop a hierarchical Bayesian model that allows us to couple the interior structure models to the observed physical properties of hot Jupiters while incorporating the measurement uncertainties. Our approach naturally accounts for non-Gaussian likelihoods. We first apply our model on the individual planets to infer the internal luminosity that reproduces the observed physical properties of hot Jupiters, namely radius, mass, and equilibrium temperature. Second, as a consequence of the high internal entropy, we find that the interior tends to be hot and show that the RCB moves to shallow pressures. Finally, we compare our findings to the proposed mechanisms and show that ohmic dissipation (Batygin & Stevenson 2010), advection of potential temperature (Tremblin et al. 2017), and thermal tides (Arras & Socrates 2010) can explain the anomalously large radii of hot Jupiters.
In a recent study, TF18 showed that the heating efficiency increases as a function of equilibrium temperature until a maximum of ∼ 2.5% is reached at around 1500 K, beyond which it decreases. The basic shape of (T eq ) provides evidence for ohmic dissipation. Building on the functional form of (T eq ), Thorngren et al. (2019) (hereafter T19) studied the effect of central heating on the interior structure of hot Jupiters and found that the internal temperature is much hotter than previous estimates, which pushes the RCB to lower pressures. Our approach is similar to TF18 but rather than modeling the extra heating as a function of , we do not assume explicitly a source for the extra heat. Instead, we consider the planet reached steady state and compute the internal luminosity given the planet mass, radius, and equilibrium temperature. The advantage of this approach is twofold: first, we can compare our results to heating mechanisms where the source of extra heat is not the stellar irradiation, and second, we self-consistently study the effect of high internal entropy on the interior structure of hot Jupiters, namely the internal temperature and pressure of the RCB. We note, however, that both approaches should lead to the same results. We also con-vert the internal luminosity to a heating efficiency and compare our results to TF18 in Section 6.3. We show that our results are qualitatively similar using a larger sample focused on FGK main-sequence stars and using an independent interior structure model.
The outline of this paper is as follows. Section 2 provides an overview of the sample selection criteria. In Section 3 we present the interior structure model used in this analysis. In Section 4 we outline the probabilistic framework used to link observations and theory and derive the basic equation which our method is based on (Equation (28)). We validate the statistical model by applying it on synthetic planetary data set generated using the Generation III Bern global model in Section 5. Readers interested in the results can safely skip to Section 6 where we present the results of our analysis. We discuss the results and the shortcomings of our approach in Section 7 and conclude in Section 8.

Sample Selection
For the purpose of our study, we required that all the planets have measured masses and radii. Sestovic et al. (2018) showed that the radii of planets with masses less than 0.37 M J do not show a clear dependence on the stellar incident flux. Photoevaporation plays an important role in the evolution of such low-mass closein planets (Owen & Jackson 2012;Jin et al. 2014). Baraffe et al. (2004) also showed that these planets are subject to undergo Roche-lobe overflow. We therefore restrict our analysis to planets with masses 0.37 < M p < 13 M J with semi-major axis a < 0.1 au. In our study, we make no attempt to correct for selection effects where it is still challenging to detect "mediuminflated" hot Jupiters around F stars using ground-based surveys (see the discussion in Section 7.5).  suggested that giant planets around stars leaving the main-sequence experience a high level of irradiation that could ultimately increase their radii. However, other studies argued that ohmic heating cannot re-inflate planets after they have already cooled (Wu & Lithwick 2013;Ginzburg & Sari 2016). A handful of re-inflated planets have been discovered around giant stars (Grunblatt et al. 2016(Grunblatt et al. , 2017Hartman et al. 2016). Since different mechanisms can be at play around evolved stars, we exclude such planets and only consider hot Jupiters around solar-like stars. Specifically, we consider stars with stellar temperature T * = 4000 − 7000 K and surface gravity log g = 4 − 4.9.
The data was taken from the Transiting Extrasolar Planet Catalogue (TEPCat 1 ; Southworth 2011), last accessed on November 2018. The aforementioned constraints on the planet mass, semi-major axis, and stellar temperature and surface gravity, lead to a final sample consisting of 314 hot Jupiters. The equilibrium temperature (T eq ) values in the literature are often not homogeneous, where different teams use different assumptions for the albedo and heat redistribution. To mitigate this, we compute the equilibrium temperature for all the planets assuming a circular orbit, zero albedo, and full heat redistribution from the day-side to the night-side (Guillot 2010) where T * and R * are the stellar temperature and radius, respectively, and a is the semi-major axis. Figure 1 1. Equilibrium temperature-radius diagram (left panel) and mass-radius diagram (right panel) color-coded by entropy for the 314 hot Jupiters selected for our analysis. The solid black and the red dashed lines compare the radii computed using our model completo21 and TF18, respectively. Both models are for a 1 M J planet with a pure H/He composition with Y = 0.27 at 5 Gyr and without accounting for inflation. The entropy was computed using the observed physical properties and an assumed heavy-element fraction of 0.2. Planets with large radii tend to have high internal entropy, with a weak dependence on planetary mass. and mass-radius (right panel) diagrams color coded by the entropy 2 . The entropy was calculated for all the planets given the observed physical properties of each system and assuming the fraction of heavy element is 20% the planet mass. Note that this value was chosen arbitrarily and for the rest of the results presented in this paper, we use the mass-heavy-element mass relation (Thorngren et al. 2016, see also Section 3.2). The solid black line is the radius at 5 Gyr computed using the interior structure model (see Section 3) for a 1 M J planet with a pure H/He composition and the He mass fraction set to Y = 0.27. The dashed red line is the same model computed by TF18. These models do not account for inflation and the radii can be considered as an upper limit for radii expected in the absence of inflation mechanisms. The radii of most of the planets with T eq > 1000 K are larger than the predicted values found with standard planet evolution models (e.g. . It is also evident that larger internal entropy leads to larger radii as noted by previ-2 When comparing to other work, it is crucial to use the same entropy zero-point or to correct for this. See Footnote 2 of Mordasini et al. (2017). ous work (e.g. Arras & Bildsten 2006;Spiegel & Burrows 2013;Marleau & Cumming 2014), with a weaker dependence on planetary mass. Planets with the largest radii have high equilibrium temperatures, masses below 1 M J , and high entropy in their deep convective interior. There is thus a compelling evidence from observations that the proximity to the star, planet mass, and the incident stellar flux play a major role in keeping hot Jupiters at high entropy.

Interior Structure Model
The primary way to gain insights into the interior structure of exoplanets is typically derived from theoretical structure models by matching the observed mass and radius. Such models are often used to constrain the planet bulk composition. Given the age of the host star and the mass of the planet, the amount of heavy elements is determined by matching the observed radius with the radius predicted from structure models. This has been applied to warm Jupiters (e.g. Thorngren et al. 2016), sub-Neptunes (e.g. Valencia et al. 2013), and super-Earths (e.g. Dorn et al. 2019) but is challenging to apply for hot Jupiters because the radii are inflated.
The aim of our study is to characterize the interior structure of hot Jupiters within a probabilistic framework. This allows us to gain insights into the physical properties governing the interior. We are specifically interested in inferring the internal luminosity of the planets based on the observed mass, radius, and equilibrium temperature. This is in turn will provide constraints on the heating efficiency, internal temperature, and the pressure at the radiative-convective-boundary (RCB). The standard interior structure model is briefly outlined in Section 3.1 and we discuss in Section 3.2 our approach to account for heat dissipation. The main model assumptions and limitations are addressed in Section 3.3.

Standard Model
The planetary evolution model completo21 was presented in Mordasini et al. (2012) and several modifications have been introduced since such as photoevaportation (Jin et al. 2014;Jin & Mordasini 2018) and coupling the interior to a non-gray atmospheric model (Linder et al. 2019;Marleau et al. 2019). In the following sections, we provide a brief description of the code relevant to our work and discuss in Section 3.3 the limitations of the model.
The internal structure of a gas giant planet is modeled using the 1D equations below. Equation (2) defines the conservation of mass. We assume that the planet is in hydrostatic equilibrium (Equation 3) and that the luminosity is constant with radius (Equation 4). Mordasini et al. (2012) showed that the latter assumption does not significantly affect the evolution of the planet when the heating occurs deep, as we assume (see below). Finally, Equation (5) is the energy transport equation describing the transport of energy either via radiation or convection.
In the above equations, r is the planetary radius as measured from the center, m the total mass inside r, ρ density, P pressure, T temperature, l planet internal luminosity, G the gravitational constant, and ∇ is the temperature gradient which depends on the process energy is transported. We use the classical SCvH EOS of hydrogen and helium (Saumon et al. 1995) with a He mass fraction Y = 0.27. Our model does not include a central core and all the heavy elements are homogeneously mixed in the gaseous envelope, see Section 3.3.1 for a discussion on the distribution of heavy elements. We model the heavy elements as water and adopt the widely used EOS of water ANEOS (Thompson 1990;Mordasini 2020). H/He and water are mixed according to the additive volume law (Baraffe et al. 2008). The transit radius is defined at P = 20 mbar.

Atmospheric Model
The atmospheric boundary conditions control the cooling rate of irradiated giant planets. The evolution of the planet and its final structure are thus sensitive to the upper boundary conditions . Jin et al. (2014) calibrated the semi-gray model of Guillot (2010) against the fully non-gray atmospheric models of Fortney et al. (2008) in order to determine the value of γ, the ratio of the optical to the infrared opacity. They used a nominal value of T int = 200 K. For our study, hot Jupiters are thought to be inflated due to dissipation or advection of heat into the interior, which thus leads to T int > 200 K. Hence, using the tabulated values of Jin et al. (2014) will lead to different PT structures and therefore alter significantly the interior structure of the planet. Indeed, we find that for T eq = 1500 K, T int = 500 K, and log g = 3, the relative change in the radius between using the improved version of the semi-gray model and using a non-gray model is around ∼ 7%, where the semi-gray model tend to lead to larger radii. It is essential then to have realistic atmospheric boundary conditions by using wavelength dependent radiative transfer atmospheric models.
Following a similar approach to Linder et al. (2019), we compute a grid of fully non-gray atmospheric models calculated using the petitCODE (Mollière et al. 2015. We included the following line absorbers CH 4 , H 2 O, CO 2 , HCN, CO, H 2 , H 2 S, NH 3 , OH, C 2 H 2 , PH 3 , Na, K, TiO, VO, and SiO, and the following pseudo-continuum absorbers H 2 -H 2 Collision Induced Absorption, H 2 -He Collision Induced Absorption, H − boundfree, H − free-free, H 2 Rayleigh scattering, and He Rayleigh scattering. The reference for these opacities can be found in Mollière et al. (2019). These grids are then used to relate the planet atmospheric temperature and pressure to the planet internal structure. The atmospheric grid was calculated assuming solar composition and covering a range of 2.5-4.5 in log g, 500-2700 K in equilibrium temperature, and 100-1000 K in internal temperature. The equilibrium temperature and surface gravity were chosen to cover the range of all the hot Jupiters selected in our sample.
The coupling between the atmosphere and the interior is done in the interior adiabat, following the first convective layer below the RCB. Details are given in Marleau et al. (2019). For a given log g, equilibrium temperature, and internal temperature, the corresponding pressure and temperature were used as boundary conditions to calculate the inward interior structure. The outward structure was calculated using the petitCODE structure and assuming hydrostatic equilibrium (Equation 3) between the pressure at the coupling point and 20 mbar, i.e. the pressure at which the transit radius is defined. We verify that coupling at a high fixed pressure, P = 1000 bar, or following the RCB layer does not significantly affect the transit radius with relative change around ∼ 0.3%.
The atmospheric PT structures assume constant log g. In fact, log g changes slightly in the radiative layers. Assuming that the change in log g in the radiative layers during the planet evolution is around ∼ 0.05, then the change in entropy is only around ∼ 0.05 kB/baryon for an internal temperature (T int ) of 700 K and an equilibrium temperature (T eq ) of 2500 K. It would take a change of 0.5 in log g to have a significant change in entropy (around 0.5 kB/baryon for T int =700 K and T eq =2500 K). We confirm that the change in entropy is negligible across the entire grid except for models with T eq > 2500 K, T int > 700 K, and log g < 3.5. In our sample, only WASP-12 b has T eq = 2580 K and log g = 3.0 (Collins et al. 2017) where the change in entropy is between 0.06 -0.08 kB/baryon. The radius of only one planet in our sample could be slightly underestimated, and therefore a constant log g in the PT structures is not a strong assumption.

Heat Dissipation
It is well established that, compared to cold Jupiter-like planets, the high internal entropy of a hot Jupiter increases its radius (Spiegel & Burrows 2013;Marleau & Cumming 2014). For example, the planet interior can gain entropy through ohmic or tidal heating. In this work, we do not attempt to model a mechanism to transport heat into the interior. We assume the planet is in steady state and thus do not calculate the planetary thermal evolution. We use the planet mass, radius, and equilibrium temperature (technically the stellar luminosity and the semimajor axis) from observations along with the mass-heavy-element-mass relation from Thorngren et al. (2016), to quantify the present internal luminosity L int of the planet. At steady state, L int is identical to the extra heating power deposited and thus where is the fraction of stellar irradiation transported into the interior, i.e. the heating efficiency, σ is the Stefan-Boltzmann constant, R p the planetary radius, and F is the flux the planet receives at the substellar point as a function of the stellar temperature T * , stellar radius R * , and the semi-major axis a (Guillot 2010). We assume that the heat dissipated is absorbed at τ = 2/3 and deposited at the center of the planet. Komacek & Youdin (2017) showed that heating at any depths larger than 10 4 bar yields nearly similar radii. However see the discussion relevant to this assumption in Section 3.3.2. Our definition agrees well with TF18, where they also deposit the extra heat at the center.
We note that 1D models without extra heating do not transfer energy into the interior on their own. The main effect of irradiation is that it decreases the cooling rate and thus the contraction rate of irradiated giant planets (Burrows et al. 2000). Planets with higher T eq will have a larger radius compared to an identical planet with lower T eq but still, not as large as the observed radii. The black line in Figure 1 shows the radius for a 1 M J planet with a pure H/He envelope at 5 Gyr at different T eq . All planets have R p < 1.25 R J . The difference in the radius between the highly and least irradiated planets is 0.14 R J . As such, our definition of is valid where all the extra energy is transported into the interior via a physical mechanism and it is not due to the 1D irradiated models transporting energy at high T eq .

Distribution of Heavy Elements
The distribution of heavy elements in the interior of exoplanets is still an open question. Some models assume for simplicity that all the heavy elements are in the core (Mordasini et al. 2012). For warm Jupiters, Thorngren et al. (2016) set an upper limit of 10 M ⊕ of heavy elements in the core and the rest is mixed homogeneously in the envelope. Current models developed to explain the anomalously large radii of hot Jupiters mix all the heavy elements in the envelope and do not include a central core (e.g. TF18; Komacek & Youdin 2017).
From the Juno mission, we now know that Jupiter has a diluted core (Wahl et al. 2017) based on the measurements of Jupiter's low-order gravitational moments (Folkner et al. 2017), yet these findings are challenging to explain from standard formation models (Muller et al. 2020). Even though the interior structures are highly affected by the chosen equation of state, the prediction of an enriched envelope still holds (Wahl et al. 2017). Planet formation models based on core accretion and that include the effect of envelope enrichment, also suggest that gas giant planets can be formed, notably at an accelerated rate (Venturini et al. 2016). Envelope enrichment compared to the Sun has also been observed for all of our four giant planets (Guillot & Gautier 2014).
In this work, all the heavy elements are mixed homogeneously in the convective part of the interior and are made up entirely of water. A central core is therefore not included. We compare the effect of the distribution of the heavy elements in the core versus in the envelope on the transit radius of the planet and hence on the heating efficiency . We find that for HD209458 b, 42 M ⊕ distributed in the core or in the envelope do not change significantly the radius when we account for heating in the interior. The absolute relative change in the radius is less than 2% for ranging between 0 − 5%. These results are also in agreement with Thorngren et al. (2016), which reached the same conclusion without accounting for heat dissipation. The median relative uncertainties on the radii measurements from observations in our sample is 4.3%, thus the distribution of the heavy elements has little effect on the inference of L int and therefore . We also show in Section 6 that the uncertainty on the heating efficiency is mainly dominated by the amount of heavy elements in the planet rather than their distribution within the planet.

Depth of Internal Heating
In our model, we assume that the heat is deposited in the interior of the planet. However, the pressures at which heat is deposited is still not constrained. Within the context of ohmic dissipation (Batygin & Stevenson 2010), the depth of internal heating is mainly dominated by the electrical conductivity profile and by the depth of the wind zone. The layers that contribute the most are the layers close to the RCB. At lower pressures heat is reradiated, whereas at higher pressures ohmic heating is not efficient due to the high conductivity there (Batygin & Stevenson 2010;Batygin et al. 2011). Huang & Cumming (2012) deposit the extra heat in the radiative layers and do not include ohmic heating below pressures of 10 bar. Under these assumptions, the RCB moves to higher pressures. Wu & Lithwick (2013) showed that heat deposited at deep layers requires significantly less heating efficiency in comparison to depositing the extra heat at shallow pressures. For planetary parameters similar to TrEs-4 b and using the same heating efficiency, the model of Batygin & Stevenson (2010) yields a planetary radius of 1.9 R J , while under a similar model Wu & Lithwick (2013) yields 1.6 R J . Differences in the radial profiles of the conductivity and wind might explain this difference. This however shows the difficulty in comparing models under the same heating mechanism but using different assumptions.
Komacek & Youdin (2017) studied systematically the effect of varying the depth of heating on the radius and found that heat deposited in the convective layers can explain the radii of hot Jupiters. Modest heating at pressures larger than 100 bar is enough, on condition that the heating is applied at an early age while the interior at such pressures is still convective. Heating at any pressure deeper than 10 4 bar leads to similar radii.
All the results we show are based on the assumption that heat is deposited in the deep interior. Therefore, the heating efficien-cies we compute could be underestimated. This potentially has also an effect on the interior structure of hot Jupiters, where we show that the RCB moves to lower pressures.

Statistical Model
Our goal is to estimate the internal luminosity and heating efficiency for the individual planets and for the population of hot Jupiters, while accounting for the uncertainties on the observed parameters. In this section, we describe the method used to infer the distribution of the internal luminosity and thus the heating efficiency for each planet, by establishing a probabilistic framework to link the observed planetary radius to the predicted one from the theoretical model described in Section 3. We start by describing how the internal luminosity for each individual planet is computed in Section 4.1. We will refer to this step as the lower level of the hierarchical model. In Section 4.2, we then combine the individual posterior samplings to study the global distribution of the full population. This will be referred to as the upper level of the hierarchical model. For each planet n (n = 1, 2, . . . , N), the planetary radius R p,n depends in our model on the planetary mass M p,n , the fraction of heavy elements Z p,n , the planet internal luminosity L int,n , and the stellar incident flux F p,n , which further depends on the stellar luminosity L * ,n and on the semi-major axis a n . In what follows, all the quantities refer to the individual hot Jupiter's physical parameters. In this framework, we define ω n , the parameters that determine the planetary radius for each individual hot Jupiter ω n ≡ (M p,n , Z p,n , L int,n , L * ,n , a n ) and thus the predicted radius from the theoretical models R t,n is a deterministic function of ω n , where R t,n = f (ω n ). R t,n is determined using the internal structure model described in Section 3. Given the observed planetary mass, semi-major axis, and stellar luminosity, and using the mass-heavy-element mass relation from Thorngren et al. (2016), we aim to infer the distribution of L int,n that reproduces the observed radius. We thus intend to answer the question: What is the internal luminosity of the planet given the observable parameters and our assumption on the fraction of heavy elements? Therefore, we define the likelihood function, the probability to observe the data given a specific set of model parameters, as P( D n |ω n ) = P(R p,n |M p,n , Z p,n , L int,n , L * ,n , a n ).
Finally, the posterior probability function, the probability of the parameters ω n given the data D n , is ∝ P(R p,n |M p,n , Z p,n , L int,n , L * ,n , a n ) × P(M p,n , Z p,n , L int,n , L * ,n , a n ) ∝ P(R p,n |M p,n , Z p,n , L int,n , L * ,n , a n ) × P(Z p,n |M p,n )P(M p,n )P(L int,n )P(L * ,n )P(a n ). (12) In the last line in Equation (12) we assume that L int,n , L * ,n , and a n are independent of each other and that Z p,n depends on M p,n following the mass-heavy-element mass relation (Thorngren et al. 2016). This inference allows us to account for data uncertainties. The semi-major axis is known precisely from observations and hence we fix the value to the observed one. We then marginalize over M p,n , Z p,n , and L * ,n to infer the distribution of the internal luminosity. We assume that the distribution of each of the observed parameter is a Gaussian distribution centered on the true quantity with a scatter given by the measurement uncertainties. Following the standard statistical notation, we can write where α, β, and σ Z are the values taken from the mass-heavyelement mass relation established by Thorngren et al. (2016). We use α = 57.9/317.828, β = 0.61, and σ Z = 10 1.82 /317.828 where 1M J = 317.828 M ⊕ and M p is in Jovian mass M J . Here, y | µ, σ ∼ N(µ, σ) implies that y is drawn from a normal distribution N with mean µ and standard deviation σ. U denotes that is sampled from a uniform distribution. We perform the inference twice each time using a different prior for the internal luminosity where we set τ 0 = (a, b) = (0, 5). LU and U implies that L int is drawn from a log-uniform and uniform distribution, respectively, and L int is in Jovian luminosity L J . Note that in our analysis, we do not sample , we sample L int and at each step in the Markov Chain Monte Carlo (MCMC) compute using which was obtained by combining Equations (6) and (7) and the relation between the stellar luminosity and flux. We further set a uniform prior on over the range 0 − 5% (Equation (17)). In Equation (18a), L int,n is sampled from a log-uniform distribution LU. We choose this prior because the internal luminosity covers a wide range of values and little is known about the true underlying distribution. This prior however does not lead to a uniform distribution in (see Section 4.1.1 and the right panel of Figure 3 for details), we therefore also consider a prior distribution uniform in linear space (Equations (18b)). The distribution of is uniform under this prior. In Section 4.1.1 we show in detail how the choice of prior on the internal luminosity affects the prior on and we discuss its effect on the inference. Finally, we can use the structure models to compute the internal temperature T int . As discussed in Section 3.1.1, the atmospheric models were computed for T int between 100 and 1000 K. We therefore set an upper limit of T int < 1000 K in order to avoid extrapolation.
The statistical model described in Equations (13)-(18b) and setting T int < 1000 K contain all the relevant distributions to evaluate Equation (12). All the results shown in Section 6, were 6 P. Sarkis: Inflated Hot Jupiters The posterior distributions inferred for HD209458 b using our model (Equation (12)). The gray dashed lines show the observed value for the relevant parameters. The effect of using different prior distribution leads to different posterior distributions for L int , , and Z p . The inferred posterior distributions for the other parameters (L * , M p , and R p ) are almost identical for both priors since they are constrained well from observations. produced by running MCMC using emcee (Foreman-Mackey et al. 2013). For each planet, we ran MCMC with 50 walkers each with 1000 steps and discard the first half as burn-in. At each iteration we compute the heating efficiency using Equation (19). Using 25,000 samples we then marginalize over the nuisance parameters and infer the posterior distribution of L int,n and of . The average acceptance ratio was around ∼ 0.5 for almost all the planets in the sample. As a by-product of this analysis, we also keep track of the PT profiles and thus infer the distribution of the pressure at the RCB and the planet internal temperature T int . This is useful to gain insights on the interior structure of hot Jupiters and we present the analysis in Section 6.4.

Choice of Prior on the Internal Luminosity
In the lower level of the hierarchical model (Section 4.1), we use non-informative uniform distributions in log and linear space as prior for the internal luminosity. It is worth studying the effect of the prior distribution on the final results. Figure 2 shows the marginalized distributions for HD209458 b using the two different priors. The luminosity distribution is shown in log-scale for  both distributions for illustrative purposes. Red shows the samples using a log-uniform distribution while blue using a uniform distribution in linear space. Note the strong correlation between the fraction of heavy elements Z p and the internal luminosity with a Pearson correlation coefficient ρ > 0.9. The observed parameters (R p , M p , and L * ) are reproduced in both cases and the distributions look almost identical. But the distributions of L int , the main parameter of interest, are different leading thus to different distributions in heating efficiency . We are in a regime where the data size is small and the choice of the prior distribution is important and dominates the inference. Note that Figure 2 shows the radius distribution even though we do not sample this parameter. This is useful to validate the model and to check that it predicts the observed data. Such plots are referred to as posterior predictive plots and we will apply them in Section 6.1 to validate the model for each planet.
Ideally, we would want to learn about the internal luminosity of the planet by relying entirely on the observed parameters while the choice of the prior should have minor effects on the posterior inference. Even though both distributions are noninformative, the data is not enough that the prior dominates. To put it in another way, more data is needed to be able to infer the distribution of L int independently of the choice of prior. Unfortunately, the physical parameters that can be observed for exoplanets in general and transiting planets specifically are very limited. One promising avenue might be inferring precisely the internal temperature, which was for the first time recently estimated for WASP-121 b (Sing et al. 2019) with T int = 500 K. In our results for WASP-121 b, the T int distributions look similar using both priors and therefore it is not possible to put tighter constraints on L int . Another promising approach is to put tighter constraints on the planet mass-heavy-element mass relation, which translates to tighter constraints on L int due to the large degeneracy between L int and Z p . This can be achieved by increasing the number of confirmed transiting warm Jupiters, i.e. giant planets with T eq < 1000 K. Such relatively cool planets are not inflated (Demory & Seager 2011). This allows to infer the fraction of heavy elements for such planets and re-calibrate the relation between the planet mass and fraction of heavy elements, similar to what was done by Thorngren et al. (2016) but with a larger sample.
It is important to explicitly mention that given the setup of the statistical model, the prior distributions for the individual planets are not the same because of the imposed upper limit of = 5%, which further depends on the observed parameters (Equation (19)). This can be understood by looking at the bottom line in Equation (12) 3 , where it is clear that each planet has different L * , M p , a, and Z p distributions due to differences in the observed physical properties. We confirm this by sampling the prior probability density function (PDF), i.e. by running the model on an empty data set D n for two different planets EPIC-211418729 b and WASP-48 b. By not sampling D n in Equation (12), we effectively sample the prior PDF. The left panel of Figure 3 illustrates this concept where we show that the internal luminosity prior distributions are different under the linear-uniform prior for both planets. Note though the log scale for better visualization. Even though we imposed a uniform distribution between 10 0 − 10 5 L J , L int larger than 10 2.5 L J for EPIC-211418729 b are not sampled and thus are ruled out. This cutoff in the distribution at high L int values is a consequence of the upper limit imposed on and the low stellar luminosity which translates to low T eq . With an equilibrium temperature roughly of T eq = 700 K, a heating efficiency of 5% for EPIC-211418729 b is equivalent to a maximum L int = 10 2.5 L J . On the other hand, WASP-48 b with T eq = 2000 K (i.e. high L * ), an upper limit of 5% on the heating efficiency is equivalent to a maximum of L int ∼ 10 5 L J . Note that for WASP-48 b low L int values are not ruled out but are less probable. To summarize, even if the initial prior imposed on L int is U(10 a , 10 b ) with a = 0 and b = 5, the actual prior distributions for the individual planets are different with different a and b values. This is a consequence of the additional prior on ( < 5%). Planets with low T eq , their distributions are truncated at high L int values (with b < 5). While this is not the case for planets with high T eq (with b = 5). The importance of a and b is relevant for the discussion in Section 4.2.
It is also worth studying the consequence of using different L int priors (U and LU) on the heating efficiency prior PDF since the relationship between the two parameters is deterministic following Equation (19). We follow the same procedure described in the previous paragraph, i.e. we run the model on an empty data set for EPIC-211418729 b. The right panel of Figure 3 shows samples from the prior distribution on for EPIC-211418729 b using the linear-uniform and log-uniform cases. It is evident that a log-uniform prior distribution on L int does not lead to a uniform prior on and the inference is biased towards small values. Whereas this is not the case when assuming a linear-uniform prior on L int . We want to stress that this holds for almost all of the planets in our sample and not only for EPIC-211418729 b, which was chosen arbitrarily.
From a statistical point of view, a log-uniform prior distribution is favored because of the large range of values and it is therefore easier to explore the entire parameter space in log space. However, this prior leads to biases in the distribution. To mitigate this, in the following section (Section 4.2) we develop a flexible hierarchical Bayesian model that accounts for the choice of prior. We study the population distributions under both priors in Section 6 and show that the inference at the population level is independent on the choice of prior.

Upper Level of the Hierarchical Model: Population Level
Posterior Samplings

General Framework
In Section 4.1, we inferred the distributions of L int , , T int , and pressure at the RCB (P RCB ) for each planet individually. In this Section, we derive the equations needed to study the general distribution of the (i) internal luminosity as a function of planet radius, (ii) heating efficiency, (iii) internal temperature, and (iv) pressure at the RCB as a function of T eq . The distributions (i), (iii), and (iv) provide insights into the interior structure of hot Jupiters while (ii) gives insights into the efficiency of transporting energy into the interior, similar to the work of TF18. Distributions (i) and (iv) are modeled using a 4 th degree polynomial The set of parameters describing the population is referred to as hyperparameter and defined as τ = {a 0 , a 1 , a 2 , a 3 , a 4 }. x is the planet radius R p for (i) and equilibrium temperature T eq for (iv).
There are many benefits of using polynomial regression compared to other parametric and non-parametric approaches. One important factor is that these models are flexible and can take a variety of shapes and curvatures to fit the data, making the results thus less model dependent compared to parametric models. Another important factor is that polynomial regression is similar to fitting a linear model and thus is computationally inexpensive and very fast to compute, unlike non-parametric models such as Gaussian process. A disadvantage to this approach is the curse of dimensionality, where the number of model parameters grows much faster than linearly with the growth of degree of the polynomial. In our case, we use univariate polynomial regression with degree 4 and thus the total number of model parameters is 5. Distribution (ii) is modeled using both a 4 th degree polynomial and a Gaussian function where the hyperparameters τ = max , T eq0 , s are the amplitude, the temperature at max , and the width of the Gaussian function, respectively. Finally, distribution (iii) is modeled using a Gaussian function with the hyperparameter τ = T int,max , T eq0 , s .

Derivation
In what follows, we derive the key equation which the inference is based on (Equation (28)) but first provide the motivation and simple description of the method.
We aim to use the single distributions we inferred in the lower level of the hierarchical model to infer the set of population parameters τ, which we will refer to as hyperparameters. The general form of the full posterior distribution in the hierarchical framework is In this equation N is the total number of planets, p(τ) is the prior probability distribution on the hyperparameters, p(ω n ) and p( D n | ω n ) are the prior and likelihood distributions for the individual planets, respectively. The population posterior distribution is a strong function of the prior imposed at the lower level of the hierarchical model. This can be understood if we assume that p(ω n ) is the same for all planets. Using this assumption, It is crucial therefore to make sure that the distribution we infer for the population has physical origins rather than is an output of the choice of prior. Hence, in order to account for the prior distribution imposed at the lower level of the hierarchical model, we apply the importance sampling algorithm. We follow closely the pioneering work established by Hogg et al. (2010) (see also the Appendix of Price-Whelan et al. 2018). This method has been used by Foreman-Mackey et al. (2014) to infer the occurrence rate of planets as a function of period and radius and by Rogers (2015) to infer the radius at which the composition transition from rocky super-Earth to volatile-rich sub-Neptunes. Briefly, we re-weight the individual posterior samples by the ratio of the value of the hyperparameters τ evaluated given the new hyperprior distribution to the old prior on which the individual sampling is based on evaluated at the old default τ 0 values. We derive below the marginal likelihood distribution.
For each n of N planets, we obtain K posterior samples of the parameters that determine the planetary radius θ n = (M p,n , Z p,n , L * ,n , a n ) and L int,n . Following similar notation to Section 4.1 and for brevity, we define the full set of parameters as We use the individual posterior samplings to compute the likelihood of the hierarchical model. For a single planet, the likelihood given the hyperparamters τ is where in the last equation we applied Bayes' theorem on the posterior distribution p(ω n | D n , τ 0 ), which is the posterior distribution for a single planet computed using Equation (12). The set of parameters from which the previous inference was generated is denoted by τ 0 . For example, as described in the previous section, the parameters describing the distribution of L int are τ 0 = (a, b) = (0, 5). We can then apply the Monte Carlo integral approximation to estimate the marginalized likelihood distribution over θ Essentially, we are assuming that all the probability integrals can be approximated as sums over samples. In case of infinite samples, this approximation becomes exact. Having derived the marginalized likelihood distribution for a single planet (Equation 27), the full marginal likelihood is then the product of the individual likelihoods We can then choose a prior probability distributions for the hyperparameter τ and the posterior probability distribution is Inside the sum, the numerator is the new probability distribution that we want to infer given a new set of hyperparameters τ, while the denominator is the value of the default prior on which the single posterior samples is based at the previously assumed values of τ 0 . We then re-weight the y nk posterior samples by the ratio. This approach of using the posterior samples from the lower level of the hierarchical model like data in the upper level has been first addressed by Hogg et al. (2010) (see also Foreman-Mackey et al. 2014, and TF18). Ideally, the inference of τ and ω n for all the planets should be done simultaneously, however this is computationally very expensive as it involves solving 4N + m integrals, where N is the number of planets and m is the number of hyperparameters in our model. Equation (28) is the main equation we use to infer the general distributions of (i), (ii), (iii), and (iv) defined at the beginning of this Section. We use Kernel Density Estimation (KDE) to estimate the probability density function (PDF) of each of the previously inferred distributions to compute p(y nk | τ), where we discuss below the functional forms. Note that even though we define a flat distribution for the internal luminosity, Equations (18a) and (18b), and set τ 0 = (a, b) = (0, 5), this is not strictly the case because additionally we truncate the heating efficiency 0 < < 5 % and require 100 < T int < 1000 K. Also, as noted in Section 4.1.1, each planet has a different prior probability distribution, leading thus to different values of τ 0 for each planet (for an example see left panel of Figure 3). Therefore to evaluate p(y nk | τ 0 ), we sample Equation (12) for each planet on an empty data set similar to what was done in Section 4.1.1 and then estimate the PDF using KDE.
For each of the four distributions, we define the general form y nk = g(x nk ), specifically y nk = where R p, nk and T eq, nk are the samples of the individual posterior distributions for the planetary radius and equilibrium temperature, respectively. The latter was computed at each iteration in the MCMC at the lower level of the hierarchical model and the values were stored.

Computational Details
We summarize below the computational procedure. First, at each iteration in the MCMC we sample the hyperparameters τ and evaluate the function y nk = g(x nk ) using the sampled values of τ. Second, we compute p(y nk | τ) and p(y nk | τ 0 ) using the pre-computed KDE functions. Finally, we evaluate the loglikelihood of Equation (28) lnp where in the last equation we compute the log of the sum of exponentials (log-sum-exp). In practice, this is numerically more stable compared to evaluating Equation (35).
For all the results presented below, we use emcee to sample from the posterior probability distribution (Equation (30)). The functional forms of g(x nk ) are either a 4 th degree polynomial or a Gaussian function or both. These are specified in Section 6. In what follows, we draw K = 2000 random samples from the single posterior samples when evaluating the mass-luminosityradius relation. For the other relations we set K = 1 and use the observed T eq values. This is possible since the equilibrium temperature is often well constrained from observations. We verified that accounting for the uncertainties does not effect the results. We adopt 44 walkers and run the sampler for 4000 iterations where the first half are discarded as burn-in and retain only every 20 th sample in the chain to produce independent samples. We monitored convergence by computing the acceptance ratio and by visually inspecting the trace plots and corner plots. Note that for each relation, we execute this procedure twice, each time using the samples drawn under the different prior, log-uniform LU and uniform U. By running this process twice, in Section 5 and Section 6 we show that the results are not biased by the choice of prior, unlike the lower level of the hierarchical model. All the data and results are available at https://tinyurl.com/bloated-hjs-results and the source code can be found at https://github.com/psarkis/bloatedHJs.

Model Validation using Planet Population Synthesis
To validate the statistical method, we applied the hierarchical model to a synthetic catalog based on planet population synthesis. The true distribution under which the synthetic dataset was generated is known. Applying thus our hierarchical model on this dataset allows us to evaluate the quality of the fit and to check whether the statistical model gives an accurate representation of the real distribution based on the observed data.

Generating Synthetic Catalog
The data set was generated using the Generation III Bern model of planetary formation and evolution (Emsenhuber et al. 2020). Inflation was accounted for by including a parameterized bloating model with a small addition during the formation phase compared to Equation (6) defined as where τ mp is the optical depth in the disk midplane from the star to the planet. This relation takes into account that at early times the disk is optically thick and the planet is at large semimajor axis, therefore bloating is inefficient. At later times, the planet migrates inwards, the disk dissipates, and the heating becomes relevant. Mol Lous & Miguel (2020) showed that migration can affect the inflation and radius of the planet only when high fraction of energy is deposited into the interior ( > 5%) but has no effect for smaller values.
For the heating efficiency , we use the Gaussian relation (Equation 21). Specifically, we use the values we infer using the log-U and presented in Table 1. For more details check Section 6.3. Our model also assumes the heating efficiency is constant in time and the stellar mass was fixed to 1 M . Using the same assumptions discussed in Section 3, the heavy elements are distributed homogeneously in the envelope and we use the fully non-gray atmospheric models of the petitCODE.
We perform the same cut on the synthetic data, i.e. we select only planets with 0.37 < M p < 13 M J and semi-major axis a < 0.1 au. Since the population synthesis did not produce hot Jupiters with T eq > 2250 K, we manually moved the planets inwards by 0.04 au after the formation epoch. This however does not have an effect on the inference. The population synthesis consists of 30000 single embryo per disk systems (population NG73) out of which 174 hot Jupiters made it into the synthetic sample.
One of the main advantages of the statistical model is the ability to account for uncertainties on the parameters. We generate synthetic uncertainties by calculating the relative uncertainty for M p , R p , T * , and R * based on the observed data and then taking the median of the computed values. The median of the relative uncertainty for M p and R p is 7% and 4%, respectively. Whereas the median of the relative uncertainty based on the observed data for T * and R * is 1% and 4%, respectively. These parameters were then used to calculate the uncertainty on L * .

Performing Statistical Inference on the Synthetic Catalog
The Bern planet population synthesis model is based on the coreaccretion model. As such, the model self-consistently computes the accretion of gas and solids onto the protoplanets, which we keep track of. We find that the mass of heavy elements is lower in the synthetic planets than inferred by Thorngren et al. (2016). We therefore refrain from using this relation in the lower level of the hierarchical model and replace Equation (14) by where Z ps,n is the value from the population synthesis with a standard deviation of 0.05, which is equivalent to a relative uncertainty of 5%. With this only modification to the original lower level of the hierarchical model described in Section 4.1, we apply the method to infer the distribution of L int and for each of the synthetic planets. We compare the marginalized posterior distributions of the parameters to the simulated values and confirm that we were able to reproduce M p , R p , L * , and thus T eq for all the synthetic planets. We repeated this procedure twice each time assuming L int follows a log-uniform LU distribution or a linearuniform U distribution. The individual posterior distribution for most of the synthetic planets are flat, which highlights the need for a hierarchical model that combines the individual distributions to extract useful information at the population level. This is one of the main advantages of using hierarchical Bayesian model (Loredo & Hendry 2019).
We therefore use the marginalized distribution for each planet to infer the heating efficiency -equilibrium temperature (HEET) for the synthetic population following the method described in Section 4.2. We model the HEET distribution with both a Gaussian function and a 4 th degree polynomial. The former function is used to test the ability of our hierarchical model to retrieve the input parameters of the Gaussian function and the latter function to test whether our model can indeed predict a Gaussian-like pattern.

Results Using Synthetic Data
With this procedure, we end up with four posterior distributions, which are shown in 4 th degree polynomial and also in agreement with the inferred Gaussian distribution. The median posterior using the linearuniform prior underestimates slightly the heating efficiency at the 68% (1σ) level but the true model is contained within the 90% (2σ) credible interval. This test shows that the statistical framework is able to retrieve the true distribution. The light gray distribution in the right panel is the inference done assuming log-uniform distribution without correcting for the choice of prior at the lower level. This shows the importance of understanding the prior at the lower level and highlights the need to re-weight the distributions.

Results Using Real Data
We now apply the model described in Section 4.1, i.e. the lower level of the hierarchical model, to infer the distribution of L int , , T int , and P RCB for each of the detected planets. In Section 6.1, we present diagnostic tools to validate the lower level of the hierarchical model. We then use the inferred posterior distributions to study the mass-luminosity-radius (MLR), T int -T eq , P RCB -T eq , and heating efficiency -equilibrium temperature (HEET) distributions for the population of hot Jupiters following the model introduced in Section 4.2. In Section 6.2, we show that by properly correcting for the choice of prior, the MLR distribution at the population level is prior independent. We hence present the rest of the results under the uniform in linear space prior in Sections 6.3 -6.4. For completeness, we show the results using both priors in Appendix A.

Posterior Predictive Checks
For each system, we infer the distribution of the internal luminosity that reproduces the observed radius, mass, and stel-lar luminosity while fixing the semi-major axis to the observed value. We visually inspect each system to double check that the marginalized posterior distributions of the observed parameters, M p , R p , L * , and thus T eq , are reproduced. Such plots are important to check that the model is a good fit and is thus capable of generating data that resemble the observed data. There are in total 17 systems where the observed mass and/or radius was not reproduced and thus we exclude these systems from the data set and do not include them in the analysis presented below. For most of the planets the radii are not possible from theoretical models as they are at the edge of the computed grid for a given planet mass, stellar luminosity, and semi-major axis. The observed radii tend to be larger than what is possible from the theoretical grid and most of these planets have masses M p > 2.5 M J . Note that for three systems the stellar luminosity and therefore the equilibrium temperature was not reproduced (HAT-P-20, Qatar-2, and WASP-43). We decide however to keep these systems since the difference in the equilibrium temperature is on the order of ∼ 30 K and hence the change in the internal luminosity is almost insignificant.

Mass-Luminosity-Radius (MLR) distribution
We divide the samples into four mass ranges, similar to the mass bins estimated by Sestovic et al. (2018) but further divide their second mass bin into two: the sub-Jupiter planets (0.37 − 0.7 M J and 0.7 − 0.98 M J ) and the massive-Jupiter planets (0.98 − 2.5 M J and > 2.5 M J ). The number of planets in each group is 86, 59, 119, and 33 planets, respectively. To infer the MLR distribution, we run the model (Equation (28) or equivalently Equation (36)) for each mass bin by specifying the functional form of g p (x) as a 4 th degree polynomial using Equation (20). As such, x is credible interval (shaded area) assuming a uniform prior in log (blue) and linear (red) space. Using either prior leads to almost identical results. The internal luminosity is high with the largest planets having a luminosity ∼ four orders of magnitude larger than Jupiter.
the planet radius R p in Equation (20) and the hyperparameter τ = {a 0 , a 1 , a 2 , a 3 , a 4 }. At each iteration in the MCMC, we compute following Equation (19), where the semi-major axis is fixed to the observed value and L * and R p are drawn from the individual marginalized posterior distributions. We further impose an additional log-normal prior on ∼ LN(−1, 1) for the planets with an equilibrium temperature less than 1000 K. This reflects our be-liefs that planets with low equilibrium temperatures are not inflated (Demory & Seager 2011), and thus should be small. We tested several prior probability distributions on and verify that our results are not affected by the choice prior. We repeat the full procedure twice each time drawing samples from the lower level of the hierarchical model under the different priors at the lower level (LU and U) and assign uniform uninformative priors on the hyperparameters. In Table A.1 and Table A.2 in Appendix A we give the 68% credible interval values assuming linear-uniform and log-uniform priors and provide the chains online 4 . Figure 5 shows the posterior distribution inferred for all mass bins under the two priors, uniform in log (red) and linear (blue) space. Notice that the lower right panel has a different scale to better visualize the results. Each data point is represented by a small line at the bottom of the plot at the corresponding radius. Such plots are called rug plots and are used to visualize the distribution of the data. The posterior distributions under both priors are almost identical and indistinguishable inline with the conclusion reached in Section 5 by validating the hierarchical model on synthetic data. There are few differences between both models, such as at small radii for the least massive planets and at large radii for the most massive ones. These differences are mainly dominated by the small number of planets in these regions. This highlights the importance of re-weighting the samples by dividing by the prior used to do the sampling at the lower level of the hierarchical model. For the rest of the paper, we show the results under the prior uniform in linear space, but confirm that the choice of prior at the lower level of the hierarchical model does not affect the main results and conclusions.
The basic shape of the MLR relation is similar across all mass bins, where as expected larger planets have higher internal luminosity with a plateau around 1.6 R J beyond which the luminosity is almost constant. The small drop towards high radii has little statistical significance and likely reflects the choice of a fourth-order polynomial. The inferred internal luminosity for most of the planets is several orders of magnitude larger than Jupiter, reaching even up to four orders of magnitude. We also find that the internal luminosity is mass dependent, with the most massive planets having the highest internal luminosity.
A noticeable feature is that the sub-Jupiter planets with masses 0.37 − 0.98 M J and radii less than 1 R J have an internal luminosity larger than Jupiter. At first glance, one might expect such planets to have an internal luminosity smaller than Jupiter's. We note however that the planets that have an equilibrium temperature less than 1000 K, indeed tend to have L int ∼ 3 L J and not more. A higher luminosity is expected because, even with T eq < 1000 K, these planets are still much closer than Jupiter, which reduces the cooling rate and thus leads to higher internal luminosity. As for the planets that have equilibrium temperature larger than 1000 K, they tend to have higher fraction of heavy elements distributed in the envelope. There are only two sub-Jupiter planets in our sample that have radii less than 0.7 R J , K2-60 and WASP-86, both of which require large fraction of heavy elements, 0.64 and 0.8, respectively, ruling out values less than 0.5. The high fraction of heavy elements explains the high luminosity values and the small number of planets with radii less than 1 R J is why the distribution is poorly constrained in this regime.

Heating Efficiency Equilibrium Temperature (HEET) distribution
Similar to the previous section, we also apply the model defined in Section 4.2 to study the HEET relation using both function forms: g p a 4 th degree polynomial (Equation (20)) and g g a Gaussian function (Equation (21)) with τ = max , T eq0 , s . The former is a flexible function that allows us to constrain the general shape of the relation by relying entirely on the data as motivated in the previous section, while the latter allows us to compare our results to TF18 and to theoretical predictions. Following the same methodology applied to the MLR relation, we further impose for the g p model the LN(−1, 1) prior on the heating efficiency for planets with equilibrium temperatures less than 1000 K. Note that the individual distributions are flat, similar to the distributions of the synthetic planets and useful information can only be extracted by combining the individual distributions.
In Table A.3 in Appendix A we give the 68% credible interval values assuming LU and U priors using the polynomial model. The Gaussian models are shown in Table 1 and the MCMC chains are available online 5 . The true distribution that was used to generate the synthetic data in Section 5 are the values we obtained using the log-U prior and shown in Table 1.
The left panel of Figure 6 shows that the posterior distributions are similar under both functional forms, with the polynomial function leading slightly to higher efficiencies. Using an independent interior structure model and a larger sample focused on FGK main-sequence stars, our results are qualitatively consistent with TF18. We confirm the Gaussian pattern holds independent of the choice of prior (see Figure A.1 in Appendix A). This pattern was predicted by ohmic dissipation first based on simulations (e.g. Menou 2012) and then later supported by TF18. Our analysis provides further evidence of the Gaussian-like distribution.

Comparison to TF18
To compare our results to TF18, we report the median and the 68% credible interval of TF18 in Table 1. We also show the posterior distributions in the right panel of Figure 6. The heating efficiency increases until a maximum is reached at T eq0 , beyond which the efficiency decreases. Our result regarding the maximum heating efficiency agrees well within 1σ with TF18, where we determine max ∼ 2.50 %, compared to ∼ 2.37 %. In our model, the peak occurs at ∼ 1860 K, while TF18 estimate the transition at ∼ 1566 K. This discrepancy can be attributed either to differences in the statistical framework or in the interior structure model. We will address both next.
While TF18 used a non-parametric Gaussian Process (GP) approach to model the HEET distribution, they found consistent results with the Gaussian function. In our study, instead of modeling the HEET distribution with a non-parametric GP model, we use a flexible 4 th degree polynomial that we stress is very fast to compute 6 and find consistent results with the Gaussian function. To test whether this discrepancy could be due to the There is a good agreement between the Gaussian and poly models, which shows that indeed the HEET distribution follows a Gaussian function. Our results are in agreement with the findings of TF18 although the peak in our models is shifted to higher equilibrium temperatures. statistical framework, we ran our statistical model using the individual distributions inferred by the analysis of TF18, which were shared with us. Note that using their data, there is no need to re-weight the distributions. See Section 6.3.2 for a detailed explanation. We confirm we were able to recover their posterior distribution using both a Gaussian function and a 4 th degree polynomial. There is a very good agreement at the 1σ level, except for T eq < 1000 K where the results are slightly different. The amplitudes are in agreement at the 1σ level even though we find tighter credible intervals at the 1σ level but very good agreement at 2σ. With this we conclude that the differences are not due to the statistical framework.
We now study the differences in the interior structure model by comparing the solid black and dashed red models in the left panel of Figure 1 computed using our model completo21 and by TF18 7 , respectively. Both of these models are for a 1 M J planet with a pure H/He envelope without accounting for inflation. In our case the planets are 5 Gyr old. Using our structure model, R p ranges between 1.12 − 1.22 R J for T eq between 800 − 2500 K. In comparison, R p is between 1.11 − 1.31 R J using the TF18 models for the same T eq interval. While the the radii at low T eq are almost identical, the differences at high T eq are up to ∼ 0.1 R J . We notice that both models lead to different radii starting at T eq > 1500 K. This difference could explain the higher heating efficiency we infer at T eq > 2000 K. Since the planets in our model have smaller radii starting at 1500 K, then more energy needs to be transported into the interior to reproduce the observed radius, leading to higher values compared to TF18. Note that this is a simple case scenario where the models are for planets made entirely of pure H/He. While this scenario explains the trend, more tests are needed to compare the radii at different T eq for different fraction of heavy elements since the details of the EOS for the heavy element could in principle also be a source of discrepancy between the models.
The discrepancy in the radii could be caused by differences in the atmospheric modelling. We next compare the atmospheric models of petitCODE and Fortney et al. (2007), which was used by TF18 and did not include TiO and VO. The previously computed atmospheric grid using petitCODE include TiO and VO (see Section 3.1.1). We therefore calculate the PT structure for a typical hot Jupiter at solar compostion, T int = 100 K, T eq = 2000 K, and log g = 3.27 without accounting for these absorbers. We find that in the absence of TiO and VO no inversion was formed with similar profiles using both atmospheric models. We then calculate the entropy of both structures using our EOS. We find that the entropy in the convective layers at pressure of 10 4 bar is 7.55 kB/baryon using petitCODE compared to an entropy of < 7.65 kB/baryon at pressure of 1000 bar using the models of Fortney et al. (2007). Note that at this pressure the structure is still not convective and thus the entropy is smaller than 7.65 kB/baryon in the convective layers and most likely the difference is < 0.1 kB/baryon between both models. Note that when including TiO/VO the entropy is ∼ 7.5 kB/baryon. A higher entropy leads to larger radii (e.g. Spiegel & Burrows 2013;Marleau & Cumming 2014) and as such we conclude that the difference between both models presented in Figure 1 could be due to differences in the opacities high up in the atmosphere, namely TiO/VO, which then have a larger effect on the deep atmosphere due to the (anti) greenhouse effect.
A more systematic comparison between both atmospheric models for different T eq and log g is required to further quantify the discrepancies, which is beyond the scope of this paper. We note that both studies do not account for systematic differences in the structure models. Such comparisons will therefore allow similar future studies to account for the systematic differences and thus infer more reliable credible intervals.

Are the Results of TF18 Prior Dependent?
In short, no.
In our study, we sample L int and then compute using Equation (19). We imposed two different prior distributions on L int because we do not have a priori knowledge which distribution best represent the population. Within a statistical framework, a log-uniform distribution is preferred in order to explore the entire parameter space. However, as we showed in the right panel of Figure 3, this leads to biases giving more weight to lower heating efficiency. Whereas, a linear-uniform prior distribution on L int leads to approximately a uniform prior distribution on .
As discussed in Section 4.2.2, the choice of prior distribution is important as the posterior distribution scales to the number of planets N. Thus the need to re-weight the distributions at the upper level. Another way to approach this study is to perform a full hierarchical Bayesian modeling where the inference on both the individual planets and population is made simultaneously (Wolfgang & Lopez 2015;Wolfgang et al. 2016).
In the study of TF18, the setup is different. They sample and impose a uniform prior between 0 − 5%. There are no additional conditions that truncate the distribution, which itself is flat non-informative prior. Therefore, there is no need to reweight the distributions.
In general, it is always a good practice to sample the prior PDF distribution (Hogg & Foreman-Mackey 2018). This step is important to check whether MCMC samples correctly the specified prior distributions.

Is the Decrease in Efficiency at High T eq Real?
In order to check whether the decrease in the heating efficiency at high T eq is real or not, we re-ran the lower-level model for all the planets with T eq > 2000 K. We assumed a linear-uniform prior for L int to ensure the results are not biased towards small values. Additionally, we did not put constraints on T int . This is important for the highly irradiated planets where T int ∼ 1000 K translates to < 5%.
The main goal of this exercise is to check whether the interior structure model allows for high values for all the planets with T eq > 2000 K. The cutoff was chosen to be close to the peak of the Gaussian function as inferred previously in Section 6.3 (see also Table 1). This allows us to compare the distribution for planets with T eq close to 2000 K to the highly irradiated ones. If the structure model allows for high values for the mostly irradiated planets, then we do not have enough evidence that the decrease in the heating efficiency is real. Otherwise, there is evidence that the decrease is real. Figure 7 shows the posterior distributions of the heating efficiency . Planets with T eq < 2250 K are in blue and planets with T eq > 2250 K are shown in red. As can be seen, all but one of the The colors indicate planets with T eq < 2250 K and T eq > 2250 K in blue and red, respectively. Six out of the seven planets shown in red favor small heating efficiency values with the most probable value close to ∼ 1%. This provides evidence that the interior structure model disfavor high values and thus the decrease seen in the HEET distribution is real given our structure model. planets with T eq > 2250 K disfavor high values with the most probable value around ∼ 1%. The only planet where high values are likely is the massive hot Jupiter WASP-18 b (10.52 M J ; Maxted et al. 2013). At this mass, the radius is a weak function of L int and as it is difficult to inflate massive planets (Sestovic et al. 2018). We therefore consider WASP-18 b an exceptional case, especially that all the planets with T eq > 2000 K have This analysis illustrates that hot Jupiters with T eq > 2250 K require low heating efficiencies to reproduce their radii using our interior structure model, which supports the Gaussian-like pattern and the decrease at 2000 K. With 20 planets having T eq > 2000 K out of which only 7 planets have T eq > 2250 K, future ultra-short hot Jupiters discoveries are essential to further confirm or refute this trend.

Distributions of Internal Temperature and Pressure at the RCB
Having inferred the population level distributions of the internal luminosity distribution and the heating efficiency, it is interesting to study the effect of energy dissipation on the interior structure of the planet. In particular, we show that as a consequence of transporting energy into the interior, hot Jupiters have very hot interiors which in turn pushes the RCB to low pressures. Our findings are in agreement with Thorngren et al. (2019) (hereafter T19), where they used the HEET relation presented in TF18 to compute T int and then generate PT atmospheric models for a range of T eq and surface gravities to locate the P RCB . poly Fig. 8. T int -T eq and P RCB -T eq diagrams in the upper and lower panel, respectively. The dark and light shaded regions present the 68% and 95% credible intervals. Although the analytical approach overestimates the internal temperature at T eq between 1000 − 1800 K, there is a good agreement at T eq < 1000 K and T eq > 1800 K. Due to the increase in the internal temperature with equilibrium temperature, the P RCB moves to lower pressures with increasing T eq , reaching up to ∼ 3 bar for the most irradiated planets.
As mentioned in Section 4.1, we keep track of the PT profiles, and thus we can infer the distribution of the internal temperature and the pressure of the RCB for each planet. We again apply the model defined in Section 4.2 to study the distributions of T int and P RCB as a function of T eq . We model the distributions, T int -T eq and P RCB -T eq , as a Gaussian function and 4 th degree polynomial, respectively. At steady state, where the last equation was obtained by replacing L int = 4πR 2 p σT 4 int in Equation (6) and combining Equations (1) and (7). We use the samples from our previous analysis using the Gaussian model (see Section 6.3) to compute T int using Equation (40) and compare the results to the hierarchical Bayesian approach. We refer to the former method as the analyti-cal approach. For all the models, we assign uniform distributions on all the hyperparameters. Figure 8 shows the inferred posterior distribution for the internal temperature (upper panel) and pressure at the RCB (lower panel) as a function of the equilibrium temperature. The analytical approach leads similar results to the Bayesian approach at the lowest and highest equilibrium temperatures. However, T int is overestimated at the 2σ level for T eq between 1000 and 1800 K. This difference could be because we did not account for intrinsic scatter in the model, which we leave for future work.
For both models, almost all hot Jupiters have T int larger than 200 K, while, for comparison, the internal temperature of Jupiter is 100 K (Li et al. 2012;Guillot & Gautier 2014). This is expected given the observed inflated radii. WASP-121 b is the only exoplanet to date whose internal temperature was constrained from observations of Mg and Fe in the transmission spectrum, with a reported value of 500 K (Sing et al. 2019). With an equilibrium temperature of T eq = 2358 ± 52 K (Delrez et al. 2016), we infer T int ∼ 800 K and by inspecting the individual posterior distribution of WASP-121 b, we rule out values below 500 K. This is the first hint from observations that hot Jupiters possess hot interiors, which is associated with a high internal entropy.
Another notable parameter to study is the pressure of the RCB as this partly controls the planetary cooling rate (Arras & Bildsten 2006;Spiegel & Burrows 2013). It is known that high equilibrium temperature pushes the RCB deeper into the planet (e.g. Fortney et al. 2007), however high internal temperature pushes the RCB to lower pressures. Therefore, the location of the RCB is not known beforehand for planets with high equilibrium and internal temperatures. The lower panel of Figure 8 shows that the RCB is situated at low pressures or at shallow depths for high T eq . The effect of the high internal temperature is thus dominant. The planets receiving high stellar irradiation tend to have hot interiors, typically around ∼ 800 K, which pushes the RCB to low pressures, reaching ∼ 3 bar for the most extreme cases.
Our results agree well with T19. While we report a maximum T int of 800 K at T eq ∼ 2500 K, T19 finds the maximum T int of 700 K at T eq ∼ 1800 K. The difference is mainly due to the differences in the (T eq ) distribution (see Section 6.3). We estimate the RCB to be at 100 bar and 4 bar for T eq = 1000 K and 2000 K, respectively, in agreement with the findings of T19. Qualitatively, both models show the same pattern where the hot interior of hot Jupiters drive the RCB to lower pressures.
We provide the 68% credible interval values for the Gaussian model under both priors for the T int -T eq distribution in Table A.4. The values for the P RCB -T eq distribution are shown in Table A.5, also under both priors. For both distributions the chains are available online 8 .

Discussion
Building on the assumption that hot Jupiters are inflated because of a process leading to high internal luminosity, we infer for each planet the internal luminosity distribution that reproduces the radius given the planet mass and equilibrium temperature from observations and using the mass-heavy-element relation (Thorngren et al. 2016) as a prior for the fraction of heavy elements. We then combine the individual distributions to constrain the population mass-luminosity-radius (MLR) distribution. Assuming that the source of extra heat in the interior is the irradiation by the host star (e.g. tides or magnetic fields), we then compute the fraction of the incident flux deposited in the interior and study the heating-efficiency-equilibrium-temperature (HEET) distribution for the full population. Finally, as a byproduct of our structure model, we can also gain insights into the interior structure of the planets by inferring the distributions of the internal temperature and the pressure at the RCB.
In what follows, in Section 7.1 we discuss the consequences of the hot interior hot Jupiters possess on the internal structure. Then we discuss our results within the context of the competing heating mechanisms, mainly ohmic dissipation in Section 7.2 and advection of potential temperature in Section 7.3. In Section 7.4, we give a general comparison with analytical relations and discuss the limitations and caveats of our results in Section 7.5.

Insights into the Interior Structure of Hot Jupiters
We have shown that hot Jupiters have hot interiors, with an internal temperature as high as 800 K. This has important conse-quences on the location of the RCB, which in turn is important for the heating mechanism. Komacek & Youdin (2017) showed that heat dissipated in the convective layers suppresses cooling and thus enables the planet to maintain a large radius. Heat deposited in the radiative layer, however, does not significantly inhibit cooling. Most it is re-radiated away leading therefore to small radii. The location of the RCB is hence important to constrain the minimum depth at which the heat should be deposited and thus the efficiency of the heating mechanism. We find that the RCB is around 100 bar for planets with equilibrium temperatures of about 1000 K, and can reach 3 bar for the highly irradiated planets, which is significantly lower than previous estimates of 1000 bar without accounting for a bloating mechanism (Fortney et al. 2007). Our results are in agreement with T19 based on coupling the heating efficiency relation (TF18) to a planetary interior structure model.
Mechanisms based on transporting heat into the deep interior, such as atmospheric circulation , ohmic dissipation (Batygin & Stevenson 2010), or advection of potential temperature (Tremblin et al. 2017) rely on the existence of winds in the interior. While the extra heat must be deposited in the convective layer in order to inflate the planet (Komacek & Youdin 2017), the actual wind speeds are not constrained from Global Circulation Models (GCMs) due to inaccurate coupling between the atmosphere and deep interior. Recently, Carone et al. (2019) showed that through a better treatment of the lower boundary condition, i.e. by accounting for a hot interior, shallow zonal winds are present at 100 bar. With new estimates and better understanding of the internal temperature and pressure at the RCB, the depth of the wind zone and wind speeds can be constrained from GCM models, which in turn will be key inputs to further study the efficiency of the proposed mechanisms.

Comparison to Ohmic Dissipation
The general idea of ohmic dissipation is that equilibrium temperatures larger than 1000 K lead to thermally ionized atmospheres that couples to the magnetic field and in the presence of strong winds produces currents, which then dissipate thermally in the deep interior (Batygin & Stevenson 2010;Batygin et al. 2011). However, in the high equilibrium temperature regime and therefore high atmospheric ionization fraction, ions slow down the winds due to Lorentz force, which in turn decrease the efficiency of ohmic dissipation (Perna et al. 2010a,b). Scaling law relations based on ohmic dissipation showed that indeed the heating efficiency increases with equilibrium temperature until a maximum is reached beyond which the efficiency decreases (Menou 2012), which was also confirmed by TF18 and now in our study. The scaling laws also suggest that the location of the peak depends on the strength of the magnetic field. Therefore, studying the functional form of the HEET distribution provides insights within the context of ohmic dissipation.
Based on our analysis, we find that the HEET distribution can be modeled by a Gaussian function, in agreement with TF18 and with the theoretical predictions. We find however that the location of the peak is at 1860 K, which is higher compared to the work of TF18 that reported the peak around 1566 K (see Table 1). Menou (2012) showed that the transition is a function of the strength of the magnetic field (see his Figure 4) where stronger magnetic fields push the peak to higher equilibrium temperatures (the peak is at ∼ 1800 K for a 30 G field). Ginzburg & Sari (2016) estimate the transition around ∼ 1500 K based on analytical models and Rogers & Komacek (2014) at ∼ 1500 − 1600 K based on magnetohydrodynamic simulations. Yadav & Thorngren (2017) estimate the surface magnetic field strength of hot Jupiters using the energy flux scaling law from Christensen et al. (2009) and account for the extra heat injected using the heating efficiency relation presented by TF18. They found magnetic field strengths around 50 − 100 G for the most inflated hot Jupiters. There are no theoretical atmospheric circulation models with such strong magnetic fields, which might hence change the location of the peak. The transition is still not well constrained and might depend on the field strength but the Gaussian distribution is robust and most importantly is prior independent. Future observations of magnetic field strengths could potentially provide a better overview but for now they remain unconstrained from an observational point of view (for a current review see Griessmeier 2017; Lazio 2018).

Comparison to Advection of Potential Temperature
Another source of heat could be the movement of high-entropy fluid parcels deeper into the atmosphere, a process known as advection of potential temperature. Within this context, Tremblin et al. (2017) suggested, using two-dimensional (2D) circulation model, that this process leads to a hot interior that can naturally explain the radius anomaly of hot Jupiters. This was further supported recently by 3D GCM simulations (Sainsbury-Martinez et al. 2019). The 2D models show that a stronger stellar incident flux leads to hotter interior adiabat (see their Figure 5). We compare our results based on the 1D model to the 2D models by selecting four planets from our sample that matches their simulation parameters, i.e. log g = 2.97 ± 0.15 with the corresponding equilibrium temperatures. We do not include the model with the lowest equilibrium temperature (∼ 500 K) as it does not match any of the selected systems in our sample. The planets we selected as a function of increasing stellar incident flux are HAT-P-17 b, Corot-4 b, HD209458 b, and HATS-35 b. We then compare the temperatures at 100 bar (T 100 ) using the PT profiles based on the 2D models to the ones based on our 1D model presented in Section 3. The results are illustrated in Figure 9, where the derived temperatures at 100 bar are shown in blue circles and the red dashed line shows the fit to the data. The gray dashed line shows the 1:1 relation on which the points would lie if their model and our data derived from observations would predict identical temperatures. We find that roughly the results agree well with a slope of 1.25, deviating from the 1:1 relation. We note however that these values are model dependent and any change in the treatment of the atmospheric model, e.g. including clouds and new opacity sources, will change these values. The temperatures estimated from the average PT profiles using the 2D circulation models are larger than the values predicted by our model, varying from 6% up to 15% for the most irradiated planets. This is expected since the 2D models tend to overestimate the radii compared to the observed ones (Tremblin et al. 2017). Our results concerning the adiabatic profile are also in agreement, where the 2D and 3D atmospheric circulation models suggest a hot adiabat starting at ∼10 bar, significantly at lower pressures compared to standard irradiated models (e.g. Fortney et al. 2007). This is in agreement with our findings and conclusions that future GCM models should account for the extra heat in the interior of inflated hot Jupiters and in-line with the work of Carone et al. (2019). We note however that convection is not included in the models of Tremblin et al. (2017). In this context, the RCB should not be interpreted as a Radiative-Convective-Boundary but rather as a proxy for the Radiative-Circulation-Boundary. As such, the energy flux is downwards and not upwards, which in turn leads to a hotter adiabat.

General Comparison to Previous Studies
It is useful and informative to compare the results of our model with analytical relations. We consider the analytical approximations of the internal luminosity based on ohmic dissipation (L Huang ; Equation (14) In the above equations, B φ0 is the toroidal component of the magnetic field at a reference pressure of 10 bar, σ t is the electrical conductivity in the dissipation region, and P is the orbital period.
To compute L Huang , we fix σ t to the nominal value 10 6 s −1 and consider two different cases for B φ0 . In the first case, we fix B φ0 to 10 G and in the second case, we compute the mean magnetic field strength at the surface of the dynamo based on the scaling law of Christensen et al. (2009) where M, L, and R are the mass, luminosity, and radius of the planet normalized to solar units. We assume B φ0 = B dyn . Note that using this relation, B dyn ranges roughly between 30 and 480 G for our sample, in agreement with the previous estimates of Yadav & Thorngren (2017). We refer to these cases as L Huang, Bfixed and L Huang, Bvar , respectively. It is straightforward then to calculate L Socrates , L Huang, Bfixed and L Huang, Bvar for each hot Jupiter in our sample using the relevant physical properties.
We also examine our results within the context of advection of high-entropy material based on models of Tremblin et al. (2017). Our aim is to compare the internal luminosity of the planets that this mechanism predicts to the internal luminosities derived in Section 4.1. Tremblin et al. (2017) computed 2D PT profiles only for four planets with different T eq values. We thus need to estimate the internal luminosity of all the planets based on the model of advection of potential temperature, for which we follow the procedure described next. We first compute the entropy using the SCvH EOS (Saumon et al. 1995) and T 100 , which was derived from the 2D PT profiles based on four fiducial planets with different T eq (see Section 7.3 for more details). Second, we fit a relation between the equilibrium temperatures of the four planets and their estimated entropy. Finally, to convert the entropy into an internal luminosity, we use the entropymass-luminosity relation from an updated version of the population synthesis of Mordasini (2018). The second step allows us to compute the entropy for all the selected hot Jupiters in our sample using the observed T eq . Having calculated the entropy and knowing M p from observations, the last step allows us to compute the internal luminosity of the planets. With this procedure, we therefore calculate the internal luminosity of the planets predicted by this mechanism based on these fits and based on T eq and M p from observations. We consider three cases for comparison by assuming the planets are composed of H/He and setting the fraction of heavy elements to 0%, 10%, and 20%. We refer to these models as L Tremblin 0 , L Tremblin 10 , and L Tremblin 20 , respectively. We point out that the values should be taken with caution as there are strong approximations involved in this approach.
Finally, to compare our results to TF18, we use the analytical (T eq ) (Equation (34) in their paper) to compute and then estimate L TF18 using Equation (6). Figure 10 compares our results to the various studies where the dashed lines are the 1:10, 1:1, and 1:0.1 relations. The predicted luminosities based on the analytical solution of thermal tides as suggested by Socrates (2013) and the advection of potential temperature (Tremblin et al. 2017) are on the same order of magnitude as the ones we derive in this work based on observations. The advection of potential temperature (Tremblin et al. 2017) predicts high luminosity values for the least luminous planets in our sample. This is expected since their model tend to overestimate the radii compared to observations, even for planets with incident flux below the threshold of inflation (stellar incident flux of ∼ 2 × 10 8 erg s −1 cm −2 or T eq ≈ 1000 K).
The relation of Huang & Cumming (2012) based on ohmic dissipation leads to small internal luminosity values. Note that this relation is an order-of-magnitude estimation of the total ohmic power. We therefore caution that these results do not provide evidence against ohmic dissipation, but rather that this relation underestimates the ohmic power. Based on our results and the work of TF18, there is compelling evidence from the HEET relation that ohmic dissipation can explain the radii of hot Jupiters. The ohmic power values estimated by Batygin & Stevenson (2010) and Menou (2012) are up to three orders of magnitude higher than the values predicted by Huang & Cumming (2012) and thus on the same order of magnitude estimated in this work. Moreover, the small internal luminosity values using the relation of Huang & Cumming (2012) could also explain the findings of , where it was shown that the relation did not lead to re-inflation of hot Jupiters.
For our models with L int < 10 2 L J , the model of TF18 predicts smaller values of L int . This difference is a direct consequence of the discrepancy in as shown in the right panel of Figure 6, where as discussed in Section 6.3 we predict higher heating efficiencies for the least and the most irradiated planets.
Converting the luminosity values to a heating efficiency using Equation (6), the models of Socrates (2013) and Tremblin et al. (2017) do not lead to a decrease in the heating efficiency at the highest equilibrium temperatures. The former predicts a continuous increase as was shown by TF18 with values as high as 20-25% and the latter seems to increase moderately up to 30%, 10%, and 2% for Z p = 0, 0.1, and 0.2, respectively. This is expected given the steeper increase in the luminosity values above 10 4 L J for both models. These are the highly inflated and highly irradiated hot Jupiters (R p > 1.4 R J and T eq > 1900 K). Note that the peak in the HEET distribution in our model occurs close to 1900 K (see Section 6.3 and Table 1), beyond which decreases for higher T eq . This explains why the models of Socrates (2013) and Tremblin et al. (2017) do not predict a Gaussian function, i.e. why does not decrease at high T eq . We stress that these models can nevertheless explain the observed radii of most of the hot Jupiters and can be the dominant mechanisms responsible for inflation even in the absence of the Gaussian function. It could be thus that these mechanisms are too efficient in inflating hot Jupiters at temperatures above than 1900 K. Thermal tides have received less attention within the context of the radius anomaly problem and thus more work is needed to understand the physical regime where this mechanism is efficient.
In summary, we provide evidence that thermal tides and advection of potential temperature can reproduce the large observed radii of most of the hot Jupiters based on the internal luminosity predicted using these models. Moreover, the HEET distribution suggests that ohmic dissipation can also explain the radii of the close-in giant planets (see Section 7.2). We therefore conclude that all of these three mechanisms can explain the inflation of hot Jupiters. This is in line with our main goals where we stress that these mechanisms were tested on only a handful of exoplanets and a statistical approach is necessary to confirm or refute these mechanisms for the entire population.

Limitations and Caveats
There are important caveats and limitations related to this work that should be explicitly mentioned.
Our results and conclusions are based on a simple 1D interior structure model. Hot Jupiters however are tidally locked, which gives rise to a temperature gradient between the day-side and the night-side. The RCB at the night-side might thus be at lower pressures compared to the day-side leading to uneven cooling. As a consequence of that, Spiegel & Burrows (2013) showed using a 1+1D model that the net effect of incorporating night-side cooling leads to higher cooling rates compared to the default 1D models. 2D circulation models also showed that the location of the RCB differs from the day-side to the night-side, which further enhances the cooling rate (Rauscher & Showman 2014) and thus requires even higher efficiency to explain the radii of hot Jupiters. This is especially important for the highly irradiated planets as it was shown that the day-side-night-side temperature differences increases with stellar irradiation (Komacek & Showman 2016;).
In addition, we assume that the heat is deposited in the interior of the planet and we do not account for dissipation in the intermediate layers. A better treatment would be to deposit the heat over a range of depths similar to Ginzburg & Sari (2016) or Komacek & Youdin (2017). Moreover, even though we showed that the Gaussian profile of the HEET distribution is in agreement with ohmic dissipation there are few shortcomings to this. A key component for ohmic dissipation is the electrical conductivity σ, where the ohmic power is proportional to 1/σ (Batygin & Stevenson 2010). The electrical conductivity increases dramatically in the interior leading to efficient heating only at lower densities and thus at lower pressures. However, the layers that contribute to the inflation are not at the surface where the conductivity is maximum but rather at deeper layers (between 100 and 1000 bar; Batygin et al. 2011). Wu & Lithwick (2013) confirmed these results by showing that heat deposited at 100 bar requires significantly less heating efficiency in comparison to 10 or 3 bar (0.3% compared to 3% and 200%, respectively, see their Figure 3). It is therefore unclear whether the Gaussian functional form holds for energy dissipated at lower pressures. The depth of the heating has also direct consequences on the interior structure. For example, Huang & Cumming (2012) included ohmic heating only in the radiative layers deeper than 10 bar and showed that as a consequence of that the RCB moves to deeper pressures. However, their model cannot reproduce the radii of massive planets. Understanding the location of the RCB is crucial as it regulates the planetary cooling rate and thus the contraction rate (Arras & Bildsten 2006;Marleau & Cumming 2014). Future developments of state-of-the-art GCM models that solve the complete equations without approximations and that couple the upper atmosphere with the deep convective layers will provide a complete picture of the underlying physical processes.
Finally, in this work we did not account for observational biases. A large number of the hot Jupiters discovered to date are discovered using ground based telescopes, such as WASP (Pollacco et al. 2006) and the HATNet and HATSouth (Bakos 2018) exoplanet surveys. There is a lack of hot Jupiters with radii smaller than ∼ 1.4 R J around early-and mid-F stars. This is because detecting such planets is still challenging from the ground as the transit depths are shallow and less than 0.5%. Heng (2012) showed that if ohmic dissipation can explain the anomalously large radii of hot Jupiters, then this naturally leads to scatter in the radii at a given stellar incident flux due to variations in the opacity, albedo, cloud/hazes properties, and the magnetic fields strength. It is therefore still not quite clear whether the lack of "medium-inflated" hot Jupiters around F stars is due to observational biases or variations in the efficiency of the heating mechanism. The NASA Transiting Exoplanet Survey Satellite mission (TESS; Ricker et al. 2015) will discover such planets if they exist and will help to better constrain the efficiency of the heating mechanisms either by the lack or existence of such planets. Subsequently high precision follow-up observations with the CHaracterising ExOPlanet Satellite (CHEOPS; Broeg et al. 2013) will help to get very accurate radii.

Conclusion
In this work, we developed a flexible and robust hierarchical Bayesian model to couple the observed physical parameters of hot Jupiters to an interior structure model. The model accounts for observational uncertainties and for the scatter in the relation between planet mass and heavy-element fraction. We validated the statistical method by applying it to synthetic planets based on planet population synthesis and showed that we are able to retrieve the true distribution. We then applied this method to quantify the internal luminosity needed to explain the radii of a sample of 314 hot Jupiters. We tested this model under two different priors (assuming a log-uniform and a linear-uniform distributions for L int ) and showed that the population level distributions are prior independent ( Figure 5). This provides useful and robust constraints on the interior structure of hot Jupiters. We find that such planets tend to have hotter interiors compared to previous assumptions, and as a result, the RCB is located at low pressures, in agreement with recent work by Thorngren et al. (2019) (Figure 8).
Assuming the planet has reached steady state and assuming that the additional source of heat is the stellar irradiation, we compute the heating efficiency , defined as the fraction of stellar irradiation deposited into the interior of the planet that is needed to explain the observed inflated radii. We find that the heating-efficiency-equilibrium-temperature relation is described by a Gaussian function (Figure 6), in agreement with previous work by TF18, however, the peak is not consistent in both studies. We found that the models of TF18 predict larger radii than our models for T eq > 1500 K, which we attribute due to differences in the atmospheric modelling. The Gaussian-like pattern is more importantly in agreement with theoretical predictions based on the ohmic dissipation model (Menou 2012). We also show that thermal tides (Arras & Socrates 2010;Socrates 2013) and advection of potential temperature (Tremblin et al. 2017) can explain the observations of most of the planets in our sample and thus are possible mechanisms responsible for the anomalously large radii of hot Jupiters (Figure 10).
To conclude, we provide new insights into the interior of hot Jupiters by coupling observations to theoretical models within a powerful statistical framework. With a better understanding of the interior, we highlight the importance of accounting for the extra heat flux in the interior in 3D GCM models, which will further improve our understanding of wind speeds and hence on the efficiency of the heating mechanisms.
The future of hot Jupiters is exciting and bright. Simulations of the exptected TESS yield (Barclay et al. 2018) predict that TESS will discover more than 250 hot Jupiters suitable for RV follow-up (R p > 1 R J ) with orbital periods < 10 days orbiting bright stars (V < 14 mag), almost doubling the number of hot Jupiters discovered. The mission already detected few hot Jupiters (e.g. Kossakowski et al. 2019;Wang et al. 2019) with many yet to be discovered. Furthermore, CHEOPS (Broeg et al. 2013) is capable of detecting the phase curves of hot Jupiters, which provide information on the day-night temperature contrast. CHEOPS will therefore play a major role in providing clues into the efficiency of energy transport in hot Jupiter atmospheres (e.g. HD189733 b; Knutson et al. 2007). With a better understanding of the interior structure of hot Jupiters thanks to the development of flexible and computationally efficient statistical tools, we will be able to provide further constraints on the radius inflation conundrum.
The source code for the hierarchical model is open source and available at https://github.com/psarkis/bloatedHJs under the MIT open source software license. Part of the code is still being added and available upon request.
The posterior samples at the population level are also available online at https://tinyurl.com/bloated-hjs-results.
Software importance sampling technique and multilevel modeling and Daniel Thorngren for sharing the posterior distributions of the individual planets and the data used in Figure 1. P.S. thanks Saavi and Gabriele for exchanging priority time on the aida server. equilibrium temperature (HEET), T int -T eq , and P RCB -T eq distributions are also similar using both priors. Figure A.1 and Figure A.2 show the HEET and both the T int -T eq and P RCB -T eq distributions, respectively. Tables A.1 and A.2 present the 68% credible interval values for the model parameters for the MLR distribution assuming linear-uniform and log-uniform priors. Similarly, Table A.3 for the HEET distribution using a 4 th degree polynomial, Table A.4 for the T int -T eq distribution using a Gaussian function, and finally Table A.5 for the P RCB -T eq distribution using a polynomial function. and the log-uniform (middle) priors using a Gaussian and 4 th degree polynomial. The shaded region shows the 68% credible interval. There is a good agreement between both models using the same prior. To better compare the same model using different priors, the right panel shows the Gaussian models using log (red) and linear (blue) uniform priors. Table A.1. 68% credible interval values of the parameters for the mass-luminosity-radius (MLR) distribution for the linear-U case modelled as a 4 th degree polynomial g p (x) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 where x = R p .  Table A.2. 68% credible interval values of the parameters for the mass-luminosity-radius (MLR) distribution for the log-U case modelled as a 4 th degree polynomial g p (x) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 where x = R p .  T int -T eq and P RCB -T eq diagrams in the upper and lower panel, respectively. The shaded regions show the 95% credible interval. Both distributions are similar at the 95% level using the linear-uniform (red) and the log-uniform (blue) priors. Table A.3. 68% credible interval values of the parameters for the heating-efficiency-equilibrium temperature (HEET) distribution for the linear-U and log-U cases using the 4 th degree polynomial model g p (x) = a 0 +a 1 x+a 2 x 2 +a 3 x 3 +a 4 x 4 , where x = T eq /1000.   A.4. 68% credible interval values of the parameters for the T int -T eq distribution for the linear-U and log-U cases using the Gaussian function Equation (21), where x = T eq and T int is in K.  Table A.5. 95% credible interval values of the parameters for the P RCB -T eq distribution for the linear-U and log-U cases under the polynomial function Equation (20), where x = T eq /1000 and g p /100 in bar.