Abstract
The inflated radii of hot Jupiters have been explored by various theoretical mechanisms. By connecting planetary thermal evolution models with the observed properties of hot Jupiters using hierarchical Bayesian models, a theoretical parameter called the heating efficiency has been introduced to describe the heating of the interiors of these planets. Previous studies have shown that the marginal distribution of this heating-efficiency parameter has a single-peak distribution along the planetary equilibrium temperature (Teq). Since the observed properties of hot Jupiters are the foundation of these Bayesian inference models, there must be a corresponding feature in the observed data that leads to the inferred single-peak distribution of the heating efficiency. This study aims to find the underlying cause of the single-peak heating-efficiency distribution without relying on specific theoretical models. By analyzing the relationships between different observed physical properties, we obtain a similar single-peak distribution of the radius expansion efficiency of hot Jupiters along Teq, which can be explained by the correlation with the stellar effective temperature. However, a detailed investigation suggests that this single-peak distribution is actually the result of straightforward physical processes. Specifically, the increase in heating efficiency can be attributed to the increase in incident stellar flux, while the decrease in heating efficiency can be attributed to the rise in the gravitational binding energy associated with the increase in planetary mass.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Being the first type of exoplanet discovered, hot Jupiters have garnered significant attention, due to their unusual orbital and physical properties (M. Mayor & D. Queloz 1995; D. Charbonneau et al. 2000; G. W. Henry et al. 2000). Observational data from a variety of detection methods, including spectroscopy, have helped us to characterize these highly irradiated objects. Such information can help us to understand how hot Jupiters form and evolve, and to gain a better understanding of Jovian planets in general (H. A. Knutson et al. 2007; J. J. Fortney et al. 2008; J. J. Fortney & N. Nettelmann 2010; P. Mollière et al. 2015; D. P. Thorngren et al. 2016).
One puzzling aspect of hot Jupiters is that they have unusually large radii compared to the gaseous planets in our solar system (A. Burrows et al. 2000; D. R. Anderson et al. 2011; B. Enoch et al. 2012; J. M. Almenara et al. 2015; L. Delrez et al. 2016). Although the abnormal radii of hot Jupiters can be intuitively explained by the high levels of radiation they receive from their host stars, theoretical models that take the incoming stellar flux into account still predict radii that are smaller than observed (J. J. Fortney et al. 2007; I. Baraffe et al. 2008; B.-O. Demory & S. Seager 2011; N. Miller & J. J. Fortney 2011; B. Enoch et al. 2012; L. M. Weiss et al. 2013). A number of mechanisms have been proposed to explain the radius anomaly of hot Jupiters. Some require additional energy to heat up the interiors of hot Jupiters, which can increase their entropy and thus their radii (P. Arras & L. Bildsten 2006; G.-D. Marleau & A. Cumming 2014). For example, the power of tidal dissipation of an eccentric orbit (P. Bodenheimer et al. 2001; P. Arras & A. Socrates 2010; A. S. Jermyn et al. 2017), advection of potential temperature due to strong stellar irradiation (A. N. Youdin & J. L. Mitchell 2010; P. Tremblin et al. 2017; F. Sainsbury-Martinez et al. 2019), atmospheric circulation that creates the thermal dissipation of kinetic energy (T. Guillot & A. P. Showman 2002; A. P. Showman & T. Guillot 2002), and ohmic dissipation of the irradiation energy through a planet's magnetic fields (K. Batygin & D. J. Stevenson 2010; R. Perna et al. 2010; K. Batygin et al. 2011; X. Huang & A. Cumming 2012; Y. Wu & Y. Lithwick 2013; S. Ginzburg & R. Sari 2016), etc. Some mechanisms slow down the normal cooling of a planet. For example, the effect of enhanced atmospheric opacities (A. Burrows et al. 2007), double diffusive convection (G. Chabrier & I. Baraffe 2007; H. Kurokawa & S. I. Inutsuka 2015), and the mechanical greenhouse (A. N. Youdin & J. L. Mitchell 2010), etc.
Because hundreds of hot Jupiters have been discovered, statistical methods can be used to study the mechanisms that are responsible for their large radii. Hierarchical Bayesian models are particular helpful in this field, as they can separate the natural variation of a physical quantity (called “inherent scatter” in this paper) from the observational biases that affect the measurement of that quantity (called “observational bias”; B. C. Kelly 2007; L. A. Rogers 2015; A. Wolfgang et al. 2016). Using a hierarchical Bayesian model to analyze the statistical distribution of planetary radii, M. Sestovic et al. (2018) showed that at a specific irradiation flux, the planetary mass plays a critical role in regulating the magnitude of the planet's radius enlargement.
Departing from purely data-driven statistical models, D. P. Thorngren & J. J. Fortney (2018) and P. Sarkis et al. (2021) integrated the thermal evolution of hot Jupiters into their hierarchical Bayesian model. They showed that the marginal distribution of the heating efficiency, a parameter that describes how efficiently a planet's interior is heated, exhibits a peaked shape with respect to the planetary Teq. The heating efficiency increases as Teq increases, up to a maximum at Teq of 1500–1800 K, and then decreases as Teq continues to increase. While the single-peak distribution of the heating efficiency can be partially explained by various heating mechanisms of the planetary interior, a precise relationship between these mechanisms and the heating-efficiency distribution has not yet been established.
When incorporating planetary thermal evolution into hierarchical Bayesian models, a substantial number of hyperparameters are also introduced into the statistical model. This can obscure the features that are directly related to the observed physical quantities, due to the hierarchical dependence of hundreds of parameters. Because the single-peak distribution of the heating efficiency is inferred from the hyperparameters that are on top of directly observed physical quantities (D. P. Thorngren & J. J. Fortney 2018; P. Sarkis et al. 2021), there must be a corresponding feature in the directly observed quantities that can explain the peaked marginal distribution of the heating efficiency. In this work, we search through the directly observed physical quantities of hot Jupiters to find the cause of the peaked shape of the marginal distribution of the heating efficiency. We aim to identify and characterize the factors that determine the radius expansion of hot Jupiters using a data-driven approach.
This article is structured as follows. Section 2 describes and validates our two-level, two-parameter hierarchical model. Section 3 demonstrates that the radius expansion of a hot Jupiter is sensitive to both the effective temperature of its host star and the amount of incoming stellar flux the planet receives. Furthermore, the strength of this dependence on stellar effective temperature can reproduce the single-peak distribution of the heating efficiency along the Teq of hot Jupiters. Section 4 analyzes the dependence of relevant physical properties, particularly focusing on the stellar effective temperature. It finds that the single-peak distribution observed in our model can be attributed to fundamental physical mechanisms, namely the combined influence of incoming stellar flux and planetary mass. Section 5 provides a concise discussion and summary of the work.
2. Modeling and Validation
2.1. Two-level, Two-parameter Hierarchical Model
In the context of observed data related to a specific physical quantity of a population of objects, at least two sources of variation are anticipated. The first source of variation is the natural dispersion of the physical quantity itself. For a quantity that can be influenced by numerous factors, such as the radius of a hot Jupiter, it is reasonable to model this variation as a normal distribution, following the central limit theorem. The second source of variation is observational bias, which can also be modeled as a normal distribution using maximum likelihood estimation.
Hierarchical Bayesian models are effective at distinguishing between the intrinsic variability of a physical quantity and its observational bias. As a result, they are commonly used to accurately derive the radius relationships of exoplanets (L. A. Rogers 2015; A. Wolfgang et al. 2016; J. Chen & D. Kipping 2017; M. Sestovic et al. 2018).
To effectively distinguish the inherent scattering of planetary radii from observational bias, we utilize a straightforward two-level, two-parameter hierarchical Bayesian model to examine how different physical factors affect the radii of hot Jupiters. This two-parameter approach allows us to focus on the relationship between planetary radius and one specific physical variable at a time. This method effectively marginalizes the influences of all other physical factors related to the variable under investigation. The model can be outlined as follows:
The first level of the hierarchical model is defined by Equations (1) and (2). This level includes two variables, x and y, each associated with the observational standard deviations 
and 
, respectively. These two variables have been recorded as the observed sets xobs and yobs. The second level of the hierarchical model is described by Equation (3), which illustrates a power-law relationship between the variables x and y. We employ a power-law function to model the relationships between various pairs of physical quantities, as it effectively captures a broad range of dependencies with notable precision. Additionally, its concise formulation makes it well suited to examining the radius relationships of hot Jupiters (G. Laughlin et al. 2011).
When applying this model to analyze the radius relationships of hot Jupiters, y represents the mean planetary radius corresponding to a specific value of the physical quantity x. The variable x can encompass a variety of physical quantities, such as the irradiation flux, planetary mass, stellar effective temperature, and others. It is important to note that, in reality, the Gaussian errors should be modeled independently for each planet and each observation. The hierarchical model outlined in Equations (1)–(3) simplifies the planetary radius relationship by assuming that the radii of the entire hot-Jupiter population can be approximated as a single Gaussian distribution. This simplification allows for effective inference of the general characteristics of hot-Jupiter populations (M. Sestovic et al. 2018; D. P. Thorngren & J. J. Fortney 2018; P. Sarkis et al. 2021).
Based on the aforementioned two-level, two-parameter model, the corresponding hierarchical inference model is as follows:
where Equation (12) provides the posterior distribution of the model. The subsequent two subsections will illustrate the retrieval capabilities of this hierarchical model, as well as the simplifications employed when using it to investigate the radii of hot Jupiters.
2.2. Grid Test of Hierarchical Inference
One question concerning our hierarchical inference model is whether it can effectively differentiate between the two levels of dispersion. To investigate this, we conducted a grid test utilizing simulated data generated by Equations (1)–(3), employing various combinations of σy and 
. The values for log(C), γ (y, x), and 
for all test runs were fixed at 0.3, 0.5, and 0.1, respectively. For each value of x, an average of 18 simulated observations of y were generated. The simulated data for the grid test are presented in Figure 1. For each run, we use the parallel-tempering Markov Chain Monte Carlo (MCMC) code Nii-C (S. Jin et al. 2024) to sample the posterior distributions specified by Equation (12). In total, 5 million iterations were sampled for each run, with the first 1 million being discarded as initial burn-in. We examined the corner plot of each run to verify that convergence was achieved. A typical corner plot illustrating the mass–radius relationship fitting is displayed in Figure 2, while similar corner plots from other runs have been omitted for simplicity. In Section 4.1, we also compute the Gelman–Rubin criterion for the separate-binning group (Table 2), which further confirms that our parallel-tempering MCMC has reached convergence. We summarize the fitted mean values of each parameter in Table 1. The power-law relationship, along with the 1σ regions of this relationship derived from the fitted values of log(C), γ (y, x), and σy, are presented in Figure 1.
Figure 1. A grid test was performed to evaluate the fitting of the two-level, two-parameter hierarchical model, using simulated data generated by various combinations of σy and 
. In each panel, the blue dots represent the simulated data points, while the red line shows the power-law relationship given by the fitted values of log(C) and γ (y, x). The red dotted lines denote the 1σ regions of the power-law relationship, computed from the fitted σy.
Download figure:
Standard image High-resolution image
Figure 2. The left panel displays the mass–radius distribution of hot Jupiters, with mass ranging from 0.5 to 5.0 MJupiter, along with the power-law relationship fitted using the two-level, two-parameter hierarchical models. The blue dots represent the hot Jupiters. The red line indicates the fitted power-law relationship, and the red dashed lines depict the 1σ regions of this relationship. The fitted power-law index, γ (Rplanet, Mplanet), is precisely equal to 0.00
. The right panel shows the corner plot derived from parallel-tempering MCMC sampling.
Download figure:
Standard image High-resolution imageTable 1. Grid Test for the Effectiveness of Our Hierarchical Inference Model in Differentiating Between Intrinsic and Observational Scatters of y
| RUN | σy(Model) |
 (Model) | σy(Fitted) |
 (Fitted) | γ(Fitted) | log(C)(Fitted) |
 (Fitted) |
|---|---|---|---|---|---|---|---|
| A1 | 0.5 | 0.5 | 0.527 | 0.487 | 0.50 | 0.30 | 0.097 |
| A2 | 0.5 | 1.0 | 0.533 | 0.970 | 0.51 | 0.28 | 0.098 |
| A3 | 0.5 | 2.0 | 0.726 | 1.935 | 0.48 | 0.38 | 0.097 |
| A3+ | 0.5 | 2.0 | 0.577 | 1.974 | 0.50 | 0.28 | 0.096 |
| B1 | 1.0 | 0.5 | 1.051 | 0.477 | 0.49 | 0.33 | 0.099 |
| B2 | 1.0 | 1.0 | 1.051 | 0.961 | 0.49 | 0.35 | 0.097 |
| B3 | 1.0 | 2.0 | 1.094 | 1.954 | 0.48 | 0.36 | 0.096 |
| C1 | 2.0 | 0.5 | 2.014 | 0.482 | 0.52 | 0.21 | 0.096 |
| C2 | 2.0 | 1.0 | 2.014 | 0.970 | 0.51 | 0.27 | 0.097 |
| C3 | 2.0 | 2.0 | 1.954 | 1.954 | 0.52 | 0.21 | 0.096 |
Download table as: ASCIITypeset image
We find that the ability to differentiate between the two levels of dispersion depends on the relative magnitudes of the intrinsic scatter and observational bias of y. This is demonstrated by the fitted values of σy and 
for each run, as shown in Table 1. The only run that did not yield accurate values for σy and 
was A3. In this run, the simulated data were generated with a σy of 0.5 and a 
of 2.0, indicating a significant observational bias coupled with a small intrinsic scatter. However, the fitting result could be significantly improved by increasing the number of observations. In an additional A3+ run, we generated an average of 48 simulated observations of y for each value of x. In this scenario, the fitted values of σy and 
were 0.577 and 1.974, respectively, which are in close proximity to the original model.
When we apply this two-level, two-parameter hierarchical model to the hot-Jupiter population, σy and 
represent the standard deviations of the mean and observed planetary radii (σR and 
) for a given value of an independent physical quantity. We find that the fitted values of σR and 
depend on the specific independent physical quantity, such as the incoming stellar flux, planetary mass, stellar effective temperature, and others. In all instances, the value of σR is three to five times greater than that of 
, suggesting that our two-level, two-parameter model is appropriate for retrieving the intrinsic radii of hot Jupiters.
2.3. Marginalization over Planetary Mass
The fundamental assumption of our two-level, two-parameter hierarchical Bayesian inference model is that the effects of all other physical quantities can be disregarded when analyzing the relationship between the radii of the hot-Jupiter population and a specific physical quantity. This is indeed a bold assumption. Previous work investigating the radius relationships of hot Jupiters from a data-driven perspective has utilized both planetary mass and incoming flux as independent variables (M. Sestovic et al. 2018). In the work of M. Sestovic et al. (2018), the authors primarily concentrate on the impact of the incoming stellar flux, and as a result, they fit the radius relationships of hot Jupiters separately across different mass ranges (mass bins). In this study, we utilize all the hot Jupiters with masses between 0.5 and 5 MJupiter when analyzing the radius relationship for a given physical quantity, rather than binning them into different mass ranges. Our approach assumes that the mass of a hot Jupiter has a negligible influence on its radius expansion when evaluated at the population level, a notion that is supported by the distribution of hot-Jupiter radii as a function of planetary mass.
Figure 2 shows the mass–radius relationship for all hot Jupiters with mass ranging from 0.5 to 5.0 MJupiter. The prototype function used to fit this relationship is the power-law function given by Equation (3), which is suitable for fitting a diverse range of functional forms with a relatively high accuracy. The fitted power-law index, γ (Rplanet, Mplanet), is equal to 0.00, with a standard deviation of 0.01 for the mass–radius relationship. This indicates that the mean radius of hot Jupiters is not correlated with planetary mass. Such a fitted flat relationship reproduces the mass–radius power index for Jovian worlds identified by J. Chen & D. Kipping (2017), which reported a value of −0.04 ± 0.02, showing that the radius of a Jupiter-mass exoplanet is nearly degenerate with respect to planetary mass. The weak correlation between the population-wide mean radius and the planetary masses of hot Jupiters lends support to our two-parameter model.
It is important to note that the weak dependence of planetary radius on mass is statistically valid only when considering the entire population of hot Jupiters collectively. In fact, the differences in gravitational potential energy associated with varying planetary masses do influence the extent of the radius expansion observed in hot Jupiters. If we partition the entire hot-Jupiter population into smaller groups, as we do in Sections 3.3 and 4, we may inadvertently introduce bias into the results. This, in turn, could lead to incorrect interpretations if the mass distributions of the various groups are not adequately considered. In Section 4, we will demonstrate that the increase in planetary mass for hot Jupiters at ultrahigh Teq is likely the primary factor contributing to the previously inferred decrease in heating efficiency at higher Teq.
2.4. Data Binning in Incoming Flux
In addition to planetary mass, several other factors also affect the radius of a hot Jupiter. The rationale behind using a simple two-parameter inference model is to circumvent the mutual dependence among different parameters. However, caution is necessary when fitting a two-parameter relationship model, as a third parameter may implicitly affect the two-parameter relationship being examined. Therefore, for a specific physical parameter that is known to influence the radius of a hot Jupiter, it is necessary to group hot Jupiters based on that parameter. We can then fit two-parameter relationships between other physical parameters for the grouped subpopulations of hot Jupiters that share similar values of the specific physical parameters known to affect their radii. One such specific parameter that can influence the radius of a hot Jupiter is the strength of the incoming stellar flux (M. Sestovic et al. 2018; D. P. Thorngren & J. J. Fortney 2018; P. Sarkis et al. 2021). Consequently, we categorize the hot Jupiters into different groups based on the amount of flux they receive, when analyzing the two-parameter relationships between planetary radius and other physical quantities. We then fit the power-law indices separately for the hot Jupiters in these different groups. This specific binning strategy employed in this work can be referred to as flux binning.
To effectively compare the strength of the incoming flux received by hot Jupiters orbiting around different types of stars at varying orbital semimajor axes, we employ a normalized incoming stellar flux F, defined as follows:
where a is the orbital semimajor axis, σ is the Stefan–Boltzmann constant, and T* and R* represent the temperature and radius of the host star, respectively. Additionally, F is standardized to 1 at a distance of 0.1 au from a solar-like star with an effective temperature of 5777 K, which corresponds to 100 times the flux received by the Earth.
3. Single-peak Distribution: Initial Observation
3.1. Dependence on Flux Affected by Tstar
Our normalized incoming stellar flux provides only a rough estimation of the total radiation that a hot Jupiter receives, and it does not capture the variations in the types of radiation experienced by each planet. Specifically, different types of stars emit different spectral profiles, resulting in the intensity of their radiation being allocated differently across various spectral bands. In this section, we examine the correlation between the planetary radius and the incoming stellar flux for hot Jupiters orbiting different types of stars. Our goal is to determine whether the spectral types of the host stars have a significant impact on the efficiency of the radius expansion in hot Jupiters.
Figure 3 shows the relationships between planetary radius and incoming flux for hot Jupiters orbiting different types of stars, as fitted by our two-level, two-parameter hierarchical relationship model. In this context, x denotes the incoming stellar flux, while y denotes the planetary radius, based on our two-level, two-parameter model. For host stars with effective temperatures ranging from 6500 to 7500 K, 5200 to 6000 K, and 3500 to 5200 K, the fitted power-law index γ (Rplanet, flux) is 0.13
, 0.10
, and 0.08
, respectively. It appears that as the stellar effective temperature increases, the correlation between the extent of the radius expansion in hot Jupiters and the normalized incoming flux becomes stronger. At first glance, it seems that the stellar effective temperature also influences the radius expansion of hot Jupiters. A more thorough investigation of this topic will be undertaken in the following sections.
Figure 3. The flux–radius distributions of hot Jupiters (blue points) orbiting stars with varying effective temperatures, along with the fitted power-law relationships (red lines) derived from the two-level, two-parameter hierarchical models, are displayed. The red dashed lines indicate the 1σ regions of the relationships. The fitted power-law indices, γ (Rplanet, flux), for hot Jupiters orbiting host stars with effective temperatures in the ranges of 6500–7500 K, 5200–6000 K, and 3500–5200 K are 0.13
, 0.10
, and 0.08
, respectively. The incident flux has been normalized to 100 times the amount received by the Earth.
Download figure:
Standard image High-resolution image3.2. Dependence on Tstar Affected by Flux
The results presented in Section 3.1 highlight the inherent difficulties in accurately evaluating how the radii of hot Jupiters depend on specific physical parameters. Although we have made attempts to establish the relationship between the radii of hot Jupiters and the intensity of the incoming stellar radiation, it appears that the effective temperature of the host star also has an indirect impact on this relationship, as shown in Figure 3. In this section, we explore the influence of stellar radiative properties on the radii of hot Jupiters from a different angle, emphasizing the impact of stellar effective temperature. We categorize the hot Jupiters into three groups based on the amount of incoming flux they have received and individually fit the relationship between planetary radius and stellar effective temperature for each group, using the flux-binning strategy outlined in Section 2.4.
The amount of incident flux received by a planet is normalized using Equation (13). The hot Jupiters are categorized into three groups, according to the levels of normalized incoming stellar flux they receive. The first group receives flux values ranging from 1 to 10, the second from 10 to 20, and the third from 20 to 50. For each of these groups, we fit the relationship between the stellar effective temperature and the planetary radius using our two-level, two-parameter hierarchical model, where x represents the stellar effective temperature and y represents the planetary radius, in this context.
Since the normalized incoming flux can be converted to planetary equilibrium temperatures Teq, which provide a clearer physical interpretation, we will refer to the three groups with normalized fluxes of 1–10, 20–20, and 20–50 as the subpopulations of hot Jupiters with Teq in the range of 883–1565 K, 1565–1857 K, and 1869–2311 K, respectively. Figure 4 shows the three groups of hot Jupiters, along with the fitted curves representing the relationships between planetary radii and stellar effective temperature. The mean values of the power-law index γ (Rplanet, Tstar) fitted using our hierarchical power-law relationship are 0.45, 0.97, and 0.43, with corresponding standard deviations of 0.08, 0.15, and 0.16 for each group. The fitted values of γ (Rplanet, Tstar) for these three groups indicate that the γ (Rplanet, Tstar) does vary monotonically with the amount of incoming stellar flux. In the second group, where the planetary Teq ranges from 1565 to 1857 K, the positive correlation between the extent of the radius expansion and the stellar effective temperature is the strongest. In Section 3.1, Figure 3 demonstrates that the radius expansion of hot Jupiters becomes more significant as stellar effective temperatures rise, showing a correlation between radius expansion efficiency and stellar effective temperature. Here, Figure 4 indicates that for hot Jupiters subjected to extremely intense stellar radiation, which leads to ultrahigh planetary Teq, the correlation between their radius expansion efficiency and stellar spectral type is relatively weak. This implies that at extremely high Teq, there could be a fundamental physical mechanism that mitigates the effects of radius expansion.
Figure 4. The distribution of stellar effective temperature vs. planetary radius for hot Jupiters (blue dots) with Teq in the ranges of 883–1565 K, 1565–1857 K, and 1869–2311 K, respectively. The red lines and red dashed lines show the fitted power-law relationships along with the 1σ regions of those relationships. The fitted power-law indices, γ (Rplanet, Tstar), for the three groups of hot Jupiters are 0.45 ± 0.08, 0.97 ± 0.15, and 0.43 ± 0.16, respectively.
Download figure:
Standard image High-resolution image3.3. Single-peak Distribution
The three groups of hot Jupiters discussed in Section 3.2 represent a coarse grid of planetary Teq. The three fitted values of γ (Rplanet, Tstar) indicate that the strength of the power-law relationship between planetary radius and stellar effective temperature changes according to the subpopulation-level Teq of the hot Jupiters within a group. To more precisely assess how the dependence on stellar effective temperature varies with the subpopulation-level Teq for a group of hot Jupiters, we conduct a larger set of similar fittings using a finer grid. We continuously and overlappingly bin the hot Jupiters based on the amount of normalized flux they have received, with bins covering the ranges of 1–11, 2–12, 3–13, ..., 30-40, 31–41, and 32–50, respectively. For all 32 groups, we apply the two-level, two-parameter hierarchical power-law relationship to fit the correlation between stellar effective temperature and planetary radius. In the fitting results, we focus on the key parameter γ (Rplanet, Tstar), which reflects the strength of the power-law relationship.
Figure 5 plots the γ (Rplanet, Tstar) values for all 32 groups. As in Section 3.2, we have also converted the incoming stellar flux into the Teq of hot Jupiters, which is represented on the x-axis of the figure, to provide a clearer physical interpretation. The values of the power-law index γ (Rplanet, Tstar), which quantifies the relationship between the radius of hot Jupiters and the stellar effective temperature, exhibit a single-peak distribution when plotted against the mean values of Teq of a group of hot Jupiters. Specifically, the γ (Rplanet, Tstar) initially increases with rising Teq, reaches a peak at approximately 1800 K, and then subsequently decreases. The top panel of Figure 5 displays the single-peak curves of heating efficiency as derived by D. P. Thorngren & J. J. Fortney (2018) and P. Sarkis et al. (2021). Furthermore, the shape of the single-peak distribution of γ (Rplanet, Tstar) and the position of the peak at Teq show a resemblance to the distribution of the radius anomaly of hot Jupiters as a function of Teq (see Figure 2 in G. Laughlin et al. 2011). The magnitude of the heating efficiency, the radius anomaly, and the power-law index γ (Rplanet, Tstar) are not directly comparable, since they represent different physical quantities. However, because all these quantities are associated with the heating process of hot Jupiters due to stellar radiation, it is reasonable to compare the characteristics of these curves. The single-peak distribution of γ (Rplanet, Tstar) with respect to Teq shows a trend that is similar to the single-peak distribution of heating efficiency in relation to Teq, though there is a minor discrepancy in the position of the peak. To quantitatively compare our single-peak distribution of γ (Rplanet, Tstar) with the unimodal distributions of heating efficiency provided by P. Sarkis et al. (2021) and D. P. Thorngren & J. J. Fortney (2018), we select and normalize the regions between 800 and 2425 K from these three distributions, then calculate the locations of the 25th, 50th, and 75th percentiles for each. For the heating-efficiency distributions, the three quantiles correspond to Teq values of ∼1485, 1786, and 2062 K in the results obtained by P. Sarkis et al. (2021), while in the findings from D. P. Thorngren & J. J. Fortney (2018), the quantiles are around 1404, 1614, and 1855 K. For our single-peak distribution of γ (Rplanet, Tstar), the three quantiles are located at approximately 1259, 1615, and 1821 K.
Figure 5. For different groups of hot Jupiters with Teq within different ranges, the bottom panel shows the values of γ (Rplanet, Tstar) for the power-law relationship between the effective temperature of their host star and their planetary radii, fitted using our two-level, two-parameter hierarchical model. Each blue point represents the mean Teq for a group of hot Jupiters along with the fitted mean value of γ (Rplanet, Tstar) for that group. The error bars along the γ (Rplanet, Tstar) axis represent the fitted 1σ locations, while the error bars along the Teq axis indicate the distribution range of Teq within each group. In the middle panel, each blue point and error bar show the mean effective temperature of the host stars for that group and the corresponding distribution range. The top panel displays the heating-efficiency curves for planets with varying Teq, as derived by D. P. Thorngren & J. J. Fortney (2018) and P. Sarkis et al. (2021).
Download figure:
Standard image High-resolution imageAt first glance, the single-peak distribution of γ (Rplanet, Tstar), derived from the combination of planetary Teq and stellar effective temperature, appears to replicate the single-peak shape of the distributions of the heating-efficiency parameter obtained from planetary thermal evolution models. However, the following section will show that its interpretation may be misleading, due to the influence of several interdependent physical parameters, supported by a more comprehensive analysis.
4. Single-peak Distribution: Reanalysis
In this section, we re-evaluate the relationship presented in Section 3.3 by employing a rigorous binning strategy and incorporating the influence of additional physical quantities. This approach aims to uncover the underlying factors contributing to the observed single-peak distribution.
4.1. Separate Binning
We group the hot Jupiters into separate bins to fit the power-law index γ (Rplanet, Tstar) using the hierarchical relationship. Given the limited total number of confirmed hot Jupiters, we establish a total of eight bins based on the Teq of hot Jupiters, to ensure that each bin contains a sufficient number of planets. The hot Jupiters are divided into eight bins with Teq ranges of 500–1000 K, 1000–1200 K, 1200–1400 K, 1400–1600 K, 1600–1800 K, 1800–2000 K, 2000–2200 K, and 2200–2700 K, respectively. The bins in the middle are spaced 200 K apart, while the bins at the ends are spaced 500 K apart, due to the small number of planets in those ranges. Although the separate-binning strategy leads to a coarser grid along the Teq axis, it avoids the complications associated with the moving average that can occur when using overlapping bins. Additionally, this approach produces bins with fixed centers and widths.
To ensure that our sampling yielded convergent results, we conducted eight independent sampling processes using different initial random seeds to fit each bin of hot Jupiters. For each of the eight bins, we gathered the means and standard deviations of γ (Rplanet, Tstar) from the eight independent runs. Subsequently, the Gelman–Rubin criteria (A. Gelman & D. B. Rubin 1992) were utilized to evaluate the convergence of the Markov Chains in each bin. The results are presented in Table 2, which demonstrates that all our MCMC samplings produced convergent results, as all calculated Gelman–Rubin criteria are below 1.03. This indicates that the MCMC results satisfy the convergence criteria for rigorous diagnostics (A. Gelman & D. B. Rubin 1992; S. P. Brooks & A. Gelman 1998).
Table 2. The γ (Rplanet, Tstar) Values (Mean ± Standard Deviation) Obtained from Eight Independent Sampling Procedures Within Each of the Separate Bins, Along with the Corresponding Gelman–Rubin Criteria
| Teq | run1 | run2 | run3 | run4 | run5 | run6 | run7 | run8 | Rc Value |
|---|---|---|---|---|---|---|---|---|---|
| (K) | |||||||||
| 500–1000 | −0.030 ± 0.140 | 0.058 ± 0.148 | 0.030 ± 0.169 | 0.006 ± 0.143 | 0.038 ± 0.143 | 0.036 ± 0.139 | 0.034 ± 0.159 | 0.030 ± 0.145 | 1.016 |
| 1000–1200 | 0.171 ± 0.182 | 0.218 ± 0.165 | 0.182 ± 0.176 | 0.172 ± 0.148 | 0.194 ± 0.150 | 0.210 ± 0.155 | 0.152 ± 0.117 | 0.185 ± 0.142 | 1.009 |
| 1200–1400 | 0.390 ± 0.103 | 0.385 ± 0.132 | 0.402 ± 0.120 | 0.385 ± 0.110 | 0.395 ± 0.119 | 0.383 ± 0.128 | 0.371 ± 0.122 | 0.400 ± 0.126 | 1.004 |
| 1400–1600 | 0.152 ± 0.134 | 0.164 ± 0.139 | 0.172 ± 0.150 | 0.137 ± 0.135 | 0.154 ± 0.154 | 0.115 ± 0.152 | 0.159 ± 0.125 | 0.143 ± 0.138 | 1.008 |
| 1600–1800 | 0.902 ± 0.157 | 0.877 ± 0.187 | 0.893 ± 0.177 | 0.886 ± 0.128 | 0.885 ± 0.183 | 0.886 ± 0.157 | 0.885 ± 0.179 | 0.919 ± 0.192 | 1.003 |
| 1800–2000 | 1.049 ± 0.266 | 1.167 ± 0.252 | 1.243 ± 0.258 | 1.104 ± 0.229 | 1.124 ± 0.212 | 1.116 ± 0.221 | 1.120 ± 0.238 | 1.111 ± 0.241 | 1.027 |
| 2000–2200 | 0.496 ± 0.164 | 0.476 ± 0.187 | 0.414 ± 0.205 | 0.447 ± 0.196 | 0.486 ± 0.196 | 0.448 ± 0.170 | 0.464 ± 0.197 | 0.467 ± 0.193 | 1.009 |
| 2200–2700 | −0.109 ± 0.173 | −0.108 ± 0.184 | −0.075 ± 0.209 | −0.100 ± 0.189 | −0.143 ± 0.169 | −0.128 ± 0.172 | −0.108 ± 0.172 | −0.096 ± 0.199 | 1.006 |
Download table as: ASCIITypeset image
Figure 6 displays the fitted values of γ (Rplanet, Tstar) for these eight groups. Similar to Figure 5, a single-peak distribution is observed, with the highest γ (Rplanet, Tstar) value recorded for hot Jupiters with Teq ranging from 1800 to 2000 K. Figure 6 also plots the 1σ ranges of the fitted γ (Rplanet, Tstar) for each group. The results show that the hot Jupiters with Teq ranging from 1600 to 2000 K exhibit significantly higher γ (Rplanet, Tstar) values compared to other groups, even when accounting for the 1σ fluctuations.
Figure 6. The fitted values of γ (Rplanet, Tstar) for eight groups of hot Jupiters created using a separate-binning approach. The red dots represent the mean values of Teq along with the fitted mean values of γ (Rplanet, Tstar) for each group. The gray rectangles indicate the distribution ranges of Teq for each group of hot Jupiters in the horizontal direction, while in the vertical direction, they reflect the fitted 1σ ranges of γ (Rplanet, Tstar). The vertical labels denote the number of hot Jupiters in each group.
Download figure:
Standard image High-resolution image4.2. The Rising Part: Incident Flux
Sections 3.3 and 4.1 have shown that a similar single-peak distribution of the planetary radius expansion efficiency along Teq can be observed due to the influence of stellar effective temperature. However, it is important to verify the validity of this finding, as the radius of a hot Jupiter is influenced by numerous factors. An implicit interdependent parameter that may be overlooked is the amount of incident flux received by a planet, as the strength of the stellar flux is correlated with the stellar effective temperature. Consequently, we investigate whether the single-peak distribution we have identified is influenced by an implicit dependence on incident flux. Specifically, for the eight subpopulations of hot Jupiters in Figure 6 across different Teq bins, how do the relationships, at the subpopulation level, between incident flux and stellar effective temperature vary?
The left panel of Figure 7 presents the Teq distributions of the eight subpopulations of hot Jupiters plotted against the effective temperatures of their host stars. We employ the same two-level, two-parameter hierarchical model to fit the power-law relationship between planetary Teq and stellar effective temperature for these eight bins of hot Jupiters, and we collect the means and standard deviations of the power-law index γ (Teq, Tstar). Table 3 presents the means and standard deviations of γ (Teq, Tstar) obtained from the eight bins. The means for these eight bins, listed from the lowest Teq to the highest, are −0.086 ± 0.120, −0.053 ± 0.061, −0.029 ± 0.054, 0.068 ± 0.043, 0.121 ±0.053, 0.093 ± 0.055, 0.110 ± 0.064, and 0.189 ± 0.092, respectively. It is noteworthy that these values exhibit a monotonically increasing trend across the entire range of Teq, reaching the highest value in the 2200–2700 K group. Since different Teq correspond to varying amounts of incident flux, our result indicates that the relationships between incident flux and stellar effective temperature change across the eight subpopulations of hot Jupiters in different Teq bins. In general, for planet subpopulations in higher-Teq bins, the incident flux they receive becomes increasingly dependent on the effective temperature of their host stars. As a consequence, the planetary radii will also become more influenced by the stellar temperature. This implies that the rising part of the single-peak distribution observed in Sections 3.3 and 4.1 is primarily driven by the correlation between incident flux and stellar effective temperature.
Figure 7. The left panel displays the Teq of the hot Jupiters in the eight separate bins plotted against the effective temperatures of their host stars. The right panel shows the distribution of planetary masses for hot Jupiters across each of the eight bins, indicating that the masses of hot Jupiters with ultrahigh Teq are significantly larger.
Download figure:
Standard image High-resolution imageTable 3. The Fitted Means and Standard Deviations of Three Important γ in the Eight Bins
| 500–1000 K | 1000–1200 K | 1200–1400 K | 1400–1600 K | 1600–1800 K | 1800–2000 K | 2000–2200 K | 2200–2700 K | |
|---|---|---|---|---|---|---|---|---|
| γ (Teq, Tstar) | −0.086 ± 0.120 | −0.053 ± 0.061 | −0.029 ± 0.054 | 0.068 ± 0.043 | 0.121 ± 0.053 | 0.093 ± 0.055 | 0.110 ± 0.064 | 0.189 ± 0.092 |
| γ (Rplanet, Mplanet) | −0.011 ± 0.033 | 0.019 ± 0.026 | −0.017 ± 0.018 | −0.009 ± 0.018 | −0.129 ± 0.018 | −0.094 ± 0.035 | −0.024 ± 0.044 | −0.135 ± 0.042 |
| γ (Rplanet, Tstar) | −0.030 ± 0.140 | 0.171 ± 0.182 | 0.390 ± 0.103 | 0.152 ± 0.134 | 0.902 ± 0.157 | 1.049 ± 0.266 | 0.496 ± 0.164 | −0.109 ± 0.173 |
Download table as: ASCIITypeset image
This conclusion is easier to grasp from a physical standpoint. First, the intensity of the incoming flux received by a hot Jupiter is correlated with the effective temperature of its host star. Second, a higher planetary Teq typically indicates that the planet is closer to its host star, making its structure more sensitive to variations in the incoming stellar flux. The combination of these two factors leads to the observed rising part of the single-peak distribution of Teq, which spans from approximately 1500 to 1800 K, where the radius expansion of hot Jupiters is more efficient. But why is the radius expansion efficiency of hot Jupiters observed to decrease at even higher Teq, above 2000 K?
4.3. The Falling Part: Planetary Mass
The reason for employing a two-parameter model in our fit is due to the finite number of available hot Jupiters. We do not wish to further subdivide our sample by introducing additional physical dependencies, as our goal is to perform a statistical analysis based on a population with the largest possible number of cases. However, when we split the hot Jupiters into smaller bins in Sections 3.3 and 4.1, the number of effective samples is further reduced. Therefore, we must be cautious not to introduce additional bias by relying on an even smaller set of samples. Another important physical quantity of hot Jupiters is the planetary mass. While Section 2.3 demonstrates a weak correlation between the population-wide mean radius and planetary mass, this correlation may not be applicable to smaller sample sizes. The right panel of Figure 7 shows the distribution of planetary mass for hot Jupiters across the eight separate bins. It indicates that for hot Jupiters with Teq above 2000 K, the distribution of planetary mass does not follow the same trend as that observed for hot Jupiters with lower Teq; specifically, there is a noticeable absence of planets with lower masses. Indeed, when we calculate the mean planetary masses for hot Jupiters in the eight bins, the values obtained are 1.628, 1.073, 1.091, 1.590, 1.364, 1.411, 1.771, and 2.456 (MJupiter), respectively. The large planetary masses at ultrahigh Teq can dampen the radius expansion of hot Jupiters, as planetary mass is a crucial factor in determining the extent of a planet's radius enlargement at a particular irradiation flux (M. Sestovic et al. 2018).
To investigate how the mass distribution of the limited number of hot Jupiters in the eight bins may influence their planetary radius relationship, we applied the same two-level, two-parameter hierarchical Bayesian model to fit the mass–radius relationship for the eight separate bins. We employed the same approach as in Section 2.3, albeit with smaller sample sizes. Table 3 presents the means and standard deviations of the fitted power-law index γ (Rplanet, Mplanet). The means of the fitted γ values for the eight bins, ordered from the lowest Teq to the highest, are −0.011 ± 0.033, 0.019 ± 0.026, −0.017 ± 0.018, −0.009 ± 0.018, −0.129 ± 0.018, −0.094 ± 0.035, −0.024 ± 0.044, and −0.135 ± 0.042, respectively. Clearly, as the number of planets in each bin decreases, the correlation between planetary mass and radius becomes stronger. This is especially pronounced for the bin with Teq in the range of 2200–2700 K, where the fitted γ is −0.135, indicating a significant reduction in radius at larger planetary masses. Furthermore, the masses of the hot Jupiters in the two bins with planetary Teq greater than 2000 K are significantly larger compared to those in the other bins. These factors, when combined, lead to the substantial decrease in radius expansion efficiency observed in Figures 5 and 6.
The analysis in this section elucidates the reasons for the decrease in the radius expansion efficiency of hot Jupiters at high Teq, based on the mass distribution of the currently detected hot Jupiters across different Teq bins. However, this does imply that hot Jupiters are inherently more massive at higher Teq, rather than being solely the result of an observational bias—because hotter planets are predominantly found around hotter stars with effective temperatures exceeding 6000 K, and it is known that radial velocity measurements become increasingly difficult as stars exceed the Kraft break at approximately 6250 K (R. P. Kraft 1967).
5. Discussion and Summary
The primary aim of this work is to uncover the underlying cause of the single-peak distribution of the heating efficiency of hot Jupiters, as inferred by previous studies. We fit the relationships between different pairs of observed physical quantities using two-layer, two-parameter hierarchical Bayesian models. Our analysis is solely based on the fitting of observed data and does not depend on any theoretical model. We find that a parameter characterizing the planetary radius expansion efficiency, specifically the power-law index γ (Rplanet, Tstar) of the radius relationship, exhibits a single-peak shape along the planetary Teq. This feature resembles previously derived heating-efficiency parameters found in complex inference models that incorporate planetary thermal evolution (D. P. Thorngren & J. J. Fortney 2018; P. Sarkis et al. 2021). Through a comprehensive analysis, we find that this single-peak distribution can be explained by straightforward physical processes. The rising portion of the single-peak distribution can be attributed to the increase in incident stellar flux, while the the falling part is associated with the increase in gravitational binding energy related to larger planetary masses.
An important implication of this work is the necessity for caution when making statistical inferences about exoplanet populations, since many planetary physical quantities are interdependent, and there exist subpopulations of exoplanets that exhibit distinct distributional properties. For instance, our initial finding indicates that the single-peak distribution of the radius expansion efficiency of hot Jupiters can be explained by the strength of the correlation with the effective temperature of their host stars. However, further analysis reveals a gradual strengthening of the correlation between incident stellar flux and stellar effective temperature as the planetary Teq increases. Additionally, the planetary masses of the subpopulation of hot Jupiters with Teq above 2000 K are considerably larger. Because the incident flux and planetary masses are more fundamental quantities that can influence the structure of a planet, they are likely to be the primary factors contributing to the observed single-peak distribution of the heating efficiency.
Although the effects of incident flux and planetary mass offer a more fundamental physical explanation for the radius expansion efficiency of hot Jupiters, it is still possible that the stellar effective temperature may also play a role, as indicated by the distribution of γ (Rplanet, Tstar) along the planetary Teq in Figures 6 and 7. On the one hand, different stellar effective temperatures correspond to different stellar spectral types, which can affect the specific details of the radiative transfer process occurring in planetary atmospheres. Theoretical studies indicate that the dissociation and recombination of hydrogen in ultrahot Jupiters influence the atmospheric structure and the heat transport efficiency from the dayside to the nightside (T. J. Bell & N. B. Cowan 2018; T. D. Komacek & X. Tan 2018; X. Tan & T. D. Komacek 2019; A. Roth et al. 2021). Consequently, the atmospheric structure and dynamics of ultrahot Jupiters are highly sensitive to both the planetary Teq and the type of host star (X. Tan et al. 2024). On the other hand, observations have shown that the spectra of ultrahot Jupiters appear featureless throughout the 1.1–1.7 μm region (J. Arcangeli et al. 2018; L. Kreidberg et al. 2018; M. Mansfield et al. 2018). This is because the majority of molecules are thermally dissociated, and alkaline elements are ionized within the dayside photospheres of ultrahot Jupiters, resulting in continuum opacity due to dissociated hydrogen (T. J. Bell et al. 2017; J. D. Lothringer et al. 2018; V. Parmentier et al. 2018). From this perspective, the radiative transfer process in the atmospheres of ultrahot Jupiters will be less sensitive to differences in host star types.
It is also noteworthy that subpopulations of hot Jupiters exhibit distinct distribution characteristics. For instance, the right panel of Figure 7 shows that the mass distribution of hot Jupiters in subpopulations with Teq above 2000 K differs from that of subpopulations with lower Teq. The masses of the detected hot Jupiters with ultrahigh Teq are significantly larger. First, detecting planets through the radial velocity method is extremely challenging around hot stars, which introduces observational biases. Second, gaseous planets with high Teq orbiting hotter stars may have experienced significant atmospheric escape processes, which can subsequently influence their final physical characteristics (I. Baraffe et al. 2004). Ignoring the differences in the distributional characteristics of the subpopulations of hot Jupiters in statistical inference could result in erroneous conclusions.
Since the results presented in this work are derived exclusively from statistical analyses of observational data pertaining to hot-Jupiter populations, our findings can serve as a useful reference for future theoretical studies of hot Jupiters.
Acknowledgments
We thank the anonymous referees for their constructive comments, which have significantly improved this paper. This work is supported by the National Natural Science Foundation of China (Nos. 11973094 and 12103003), the Youth Innovation Promotion Association CAS (No. 2020319), the incubation program for recruited talents (No. 2023GFXK153), and the doctoral start-up funds from Anhui Normal University. S.J., Y.C., and Z.G. acknowledge support from the Zijin exploration program provided by the High School Affiliated to Nanjing Normal University and the Nanjing Branch of the Chinese Academy of Sciences.
Software: Nii-C (S. Jin et al. 2024), gnuplot (T. Williams & C. Kelley 2022).