Issue 
A&A
Volume 616, August 2018



Article Number  A76  
Number of page(s)  13  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201731454  
Published online  21 August 2018 
Investigating hotJupiter inflated radii with hierarchical Bayesian modelling
^{1}
Center for Space and Habitability, University of Bern,
Gesellschaftsstrasse 6,
Bern
3012,
Switzerland
email: marko.sestovic@csh.unibe.ch
^{2}
Astrophysics Group, Cavendish Laboratory,
J. J. Thomson Avenue,
Cambridge
CB3 0HE,
UK
^{3}
Observatoire de Genève, Universite de Genève,
51 chemin des Maillettes,
1290
Sauverny,
Switzerland
Received:
27
June
2017
Accepted:
27
March
2018
Context. As of today, hundreds of hot Jupiters have been found, yet the inflated radii of a large fraction of them remain unexplained. A number of mechanisms have been proposed to explain these anomalous radii, however most of these can only work under certain conditions and may not be sufficient to explain the most extreme cases. It is still unclear whether a single mechanism can sufficiently explain the entire distribution of radii, or whether a combination of these mechanisms is needed.
Aims. We seek to understand the relationship of radius with stellar irradiation and mass and to find the range of masses over which hot Jupiters are inflated. We also aim to find the intrinsic physical scatter in their radii, caused by unobservable parameters, and to constrain the fraction of hot Jupiters that exhibit inflation.
Methods. By constructing a hierarchical Bayesian model, we inferred the probabilistic relation between planet radius, mass, and incident flux for a sample of 286 gas giants. We separately incorporated the observational uncertainties of the data and the intrinsic physical scatter in the population. This allowed us to treat the intrinsic physical scatter in radii, due to latent parameters such as the heavy element fraction, as a parameter to be inferred.
Results. We find that the planetary mass plays a key role in the inflation extent and that planets in the range ~0.37−0.98 M_{J} show the most inflated radii. At higher masses, the radius response to incident flux begins to decrease. Below a threshold of 0.37 ± 0.03 M_{J} we find that giant exoplanets as a population are unable to maintain inflated radii ≿1.4 R_{J} but instead exhibit smaller sizes as the incident flux is increased beyond 10^{6} W m^{−2}. We also find that below 1 M_{J}, there is a cutoff point at high incident flux beyond which we find no more inflated planets, and that this cutoff point decreases as the mass decreases. At incident fluxes higher than ~1.6 × 10^{6} W m^{−2} and in a mass range 0.37−0.98 M_{J}, we find no evidence for a population of noninflated hot Jupiters. Our study sheds a fresh light on one of the key questions in the field and demonstrates the importance of populationlevel analysis to grasp the underlying properties of exoplanets.
Key words: planets and satellites: fundamental parameters / planets and satellites: atmospheres / methods: statistical
© ESO 2018
1 Introduction
Understanding the internal state and energy processes of giant planets is a key goal in exoplanet science. Gas giant interiors are generally modelled by assuming an adiabatic temperature profile and internal structure and calculating the evolution of the planet radius as it contracts (for a review see e.g. Fortney & Nettelmann 2010; Fortney et al. 2010; Baraffe et al. 2014).
Fortney et al. (2007) and Baraffe et al. (2008) modelled the evolution of planetary radii by including the effects of stellar irradiation and heavy element cores with stateoftheart internal equations of state. Such model predictions can be tested against the currently observed population of exoplanets with known radii, masses, and orbital distances. At low stellar irradiation, these models are thought to agree with observations, but the high irradiation regime still presents a challenge. At fluxes greater than ~ 2 × 10^{5}W m^{−2} (Miller & Fortney 2011; Demory & Seager 2011), a sizeable fraction of gas giants are found to be anomalously larger than predicted; these gas giants include hot Jupiters WASP17 b, WASP121 b, and Kepler435 b, which all have measured radii R > 1.8 R_{J} (Anderson et al. 2011; Almenara et al. 2015; Delrez et al. 2016).
There are several candidate explanations for this inflation process, including tidal dissipation (Bodenheimer et al. 2001, 2003; Arras & Socrates 2010; Jermyn et al. 2017), kinetic heating (Guillot & Showman 2002), enhanced atmospheric opacities (Burrows et al. 2007), double diffusive convection (Chabrier & Baraffe 2007; Kurokawa & Inutsuka 2015), vertical advection of potential temperature (Youdin & Mitchell 2010; Tremblin et al. 2017), and Ohmic heating through magnetohydrodynamic effects (Batygin & Stevenson 2010; Perna et al. 2010; Wu & Lithwick 2013; Ginzburg & Sari 2016).
In most cases to date, theoretical predictions are tested against individual planets or by producing radius predictions for a narrow set of parameters. However, the radii may depend on many parameters, including some that are nonobservable or poorly constrained (latent parameters) such as the core mass, system age, internal composition, and atmospheric opacity. The degeneracy between these parameters and the strength of the inflationary process limits the information we can obtain from a single planet. Thanks to the sheer count of exoplanets discovered so far, we are now in a position to study the issue of inflated radii using populationlevel analyses.
In such an attempt, the choice of forward model is central. The candidate model should not only reproduce a range of individual planet observations, but also predict the full distribution of the observable planet radii and their dependence on planet mass and stellar incident flux within reasonable assumptions on the latent parameters. Laughlin et al. (2011) already used the observed radii of 90 wellcharacterised transiting exoplanets to assess the incidence of Ohmic heating on the inflated radius anomaly (first introduced by Guillot et al. 2006). Using a populationlevel approach, one can also extract modelindependent conclusions from the data. For example, Enoch et al. (2012) used multivariate regression models to investigate the factors that affect the radii of 119 extrasolar gas planets such as planetary equilibrium temperature, planet mass, stellar metallicity, and tidal heating rate.
To that end, it is important to separate the real physical scatter caused by the latent parameters and the scatter caused by measurement uncertainties (see Wolfgang et al. 2016 in which this was addressed for superEarths and subNeptunes). For example, hotJupiter radii span from ~ 1.99 R_{J} (Kepler435 b; Delrez et al. 2016)to ~0.84 R_{J} (Kepler41 b; Santerne et al. 2011). We seek to understand how much of this spread can be accounted for by the measurement uncertainty in the radius and thus constrain the true range of radii that is driven by the currently unidentified physical process.
In this work, we infer the fluxmassradius (FMR) relation using a population of 286 transiting gas giants with measured masses and radii. We use a very similar approach to Wolfgang et al. (2016) and Shabram et al. (2016), by employing ahierarchical Bayesian model (HBM) in which we separate the effects of intrinsic scatter and measurement uncertainty (see also Kelly 2007; Hogg et al. 2010; Demory 2014; Rogers 2015). This also allows us to include the uncertainty in both the dependent parameter (planet radius) and the independent parameters (planet mass and incident flux), which cannot be carried out with basic leastsquares regression. We thus constrain the FMR relation as a distribution rather than a deterministic function, and characterise the physical scatter of the gas giant population.
It is as of yet unknown whether all hot Jupiters are inflated or if some may have radii that agree with noninflationary models (Fortney et al. 2007). Many of the candidate inflation mechanisms, for example, tidal heating for circularised orbits (Gu et al. 2003), would not apply for all planets. Similarly, Heng (2012) and Perna et al. (2012) showed that the depth and degree of Ohmic heating are sensitive to numerous atmospheric parameters such as opacity or thermal inversions and may be inhibited in some cases. Measuring the FMR distribution and its intrinsic scatter thus provides a way to determine whether such nonanomalous planets exist.
We present the planet data used in this study in Sect. 2 and present our parameterisation of the FMR distribution in Sect. 3. The methods for fitting the parameters and selecting our model are explained in Sect. 4, and the results of our fit for the FMR distribution are shown in Sect. 5. We test the agreement of our model with the data in Sect. 6, and we finally discuss the implications of our results in Sect. 7.
2 Data
As we require both the mass and radius of the planets, we include only transiting giants that have radial velocity (RV) or transit timing variation (TTV) mass measurements. Of those, nine planets have masses derived from TTVs, which are biased towards lower fluxes and masses, and have uncertainties that are larger on average than for the whole sample^{1}. We include these planets because they populate a relatively sparse region of the parameter space, however we note that they may have a different selection bias than the rest of the sample (see e.g. Wolfgang et al. (2016) in which such a bias was found for subNeptune planets). We also use the stellar temperature, stellar radius, and the semimajor axis of the planetary orbit to calculate a distribution for the incident flux on the planet (see Sect. 4).
The data was taken from the Exoplanet Orbit Database (EOD; Han et al. 2014) and supplemented by the NASA Exoplanet Archive (Akeson et al. 2013), both last accessed on 10 October 2016. We checked the archive and the source paper in cases in which the EOD data is incomplete. Planets that are missing required values in both databases were discarded. We also took information about the planet ages and metallicities from the exoplanet.eu database (Schneider et al. 2011). Where the raw data have asymmetric error bars, we averaged these data and took the uncertainty distribution to be Gaussian. This follows the treatment of Weiss et al. (2013) and Wolfgang et al. (2016), since we did not have access to the full likelihood distributions for all our parameters. While this may introduce a bias into our results, a short test (see Sect. 7.1) suggested the effect on our hyperparameters is well within our posterior distribution uncertainties.
Our main criterion for selecting our sample is the planet mass. We set 13 M_{J} as the upper cutoff to exclude planets where deuterium burning is thought to happen (Chabrier et al. 2005; Spiegel et al. 2011). Choosing the lower bound is related to the wider problem of whether there are any strong boundaries in gasgiant characteristics, which is one of the questions we seek to answer. Enoch et al. (2012) chose 0.1− 0.5 M_{J} to define Saturnmass planets and 0.5−2.0 M_{J} for Jupiterlike planets, finding that Saturnmass planets do not exhibit any inflation with temperature. Laughlin et al. (2011) also used 0.1 M_{J} as a lower bound for their sample.
We chose 0.1 M_{J} as an absolute lower bound with the understanding that this introduces planets into our model that do not have the same properties as the rest of our sample, i.e. Saturn and Neptunelike planets that may not exhibit inflation, as shown by Enoch et al. (2012). Finding mass boundaries between the hotJupiter planets that inflate and those that do not is part of our HBM implementation (see Sects. 3 and 4).
3 Forward model
In this section, we review the known trends and present our main parametrisation for the FMR distribution. We refer to the following as the baseline model. Our current understanding of giantplanet structure and evolution is that at low incident flux, thermal evolution models can reproduce the range of planet radii present in the exoplanet population (Miller & Fortney 2011; Thorngren et al. 2016). Although previous models have predicted that radii depend on stellar irradiation, particularly at fluxes of ~ 1 × 10^{5}W m^{−2} and greater (~0.1 AU), it is also clear that a large percentage of gas giants in shortperiod orbits are much more inflated than predicted (see Fortney et al. 2007, 2010; Baraffe et al. 2008, 2014; Fortney & Nettelmann 2010). Demory & Seager (2011) and Miller & Fortney (2011) showed that these inflated hot Jupiters start to occur at incident fluxes of greater than ~ 2 × 10^{5}W m^{−2}. Enoch et al. (2012) and Laughlin et al. (2011) also found that the radii of these hot gas giants can depend strongly on irradiation temperature. These trends can be seen in the top of Fig. 1.
Taking this into account when choosing the parametrisation of our empirical fluxmassradius relation (represented as the forward model for radius, R, given mass and flux, M and F), we modelled the effect of stellar irradiation by considering two regimes. For planets below a threshold F_{s} in incident flux, the radii are taken to be constant at size C, and independentof flux. Above F_{s} the radii increase proportionally with logF. This relation (Eq. (1)) qualitatively agrees with the distribution of radii and fluxes in Fig. 1 and is written as (1)
While planetary radii may depend on incident flux even at F < F_{s}, this dependence is very weak, and variations in age, core mass, and envelope composition are expected to play a significantly greater role; see for example Fortney et al. (2007); Baraffe et al. (2008); and Miller & Fortney (2011). Additionally, because of a lack of transiting planets with known masses found at large orbital distances (there are only 22 planets further out than 0.1 AU in our sample), combined with the observational uncertainty of the measurements, it is unlikely that we could extract meaningful constraints about the dependence of planetary radius with incident flux for gas giants at low irradiation.
Equation (1) provides a deterministic way of modelling the planetary radius based on the incident flux and model parameters, which we must infer (A, C, and F_{s}). However, owing to nonobservable parameters such as core mass, composition, and age, we expect to see a physical scatter in the radii at a given flux level. We include this uncertainty by attaching to the planet radius a Gaussian distribution of standard deviation σ_{R}, around a mean value μ_{R}, which is given the following same form as Eq. (1): (2) (3)
where ~N(μ, σ) means drawn from a normal distribution with mean μ and variance σ^{2}. Enoch et al. (2012) found that there were large differences in the radiustemperature relations between planets of different masses. We thus include the effect of planet mass in our model by splitting our total mass range into J separate bins, each with a different radiusflux relation of the same form as Eq. (2), but with independent model parameters (C, F_{s}, and A) as follows: (4)
where C(M) and F_{s} (M) follow the same form and represent the inner boundaries separating the mass bins. Unlike previous studies, we do not fix the mass boundaries to be constant, as their placement would be arbitrary and could bias the FMR relation. Instead, the mass bin boundaries are also model parameters that we infer.
Our final baseline model uses a mass and fluxdependent scatter, σ_{R} (M, F), with three possible values for three classes of planets: planets in the lowest mass bin (M < m_{1}), and inflated and weakly irradiated planets (F ≥ F_{s} and F < F_{s} respectively) shared for all higher mass bins (M ≥ m_{1}), written as: (5)
The baseline model, with 4 mass bins and 17 parameters, was chosen after a process of model selection, inwhich we considered several alternative parametrisations and compared their posterior results. We varied the number of mass bins, choosing between 3, 4, and 5 bins, as well as the form of σ_{R} and F_{s}. For the flux threshold, we compared with the case in which the threshold was the same for all mass bins. For the intrinsic scatter, we also set it todepend on just the mass instead, or on both mass and incident flux, and we tried models in which the intrinsic scatter was the same for all fluxes and masses (as a constant single parameter). We did this by considering which models converged to unimodal posteriors and with tests on our posterior predictive fit, which are described in Sect. 6.
Fig. 1 Top panel: distribution of the transit radii in our sample vs. the incident flux. The mass is colour coded. The radii begin to show a correlation with flux at ~ 2 × 10^{5} W m^{−2}, corresponding to where inflated hot Jupiters begin to appear (Demory & Seager 2011). Bottom panel: the transit radii plotted against mass, colour coded by incident flux. 
4 Hierarchical Bayesian model
Our aim is to infer the values of the model parameters^{2} of the distribution for R given F and M, or RF, M. We also refer to the parameters by shorthand, for example , or all of the parameters together as α = (A, C, F_{s}, m, σ_{R}). These are not necessarily vectors however.
The probabilistic model described in Sect. 3 defines the relationship between the true parameters R, M, and F of a planet. However when an observation is made, for example measuring the transit radius , the result is subject to noise. It can be taken as being drawn from a probability distribution that depends on the true radius R and the noise characteristics of the observation. This is represented as , where α_{obs} parametrises the noise of a particular observation. We take care to differentiate between the observed value and the true value R. In the most simple case, is drawn from a Gaussian distribution with mean R and a standard deviation, which is the observational uncertainty, σ_{obs}. As shorthand, we refer to the combined observational data as .
The distribution of our observables for a given planet is the distribution of the true parameters (R, M, F) convolved with the observational uncertainty distribution. Thus the true scatter of the radii, σ_{R}, is lower than the scatter of our observed data, and the posterior values for our parameters (R, M, F) are less spread out than the data (see Figs. 6 and 7). This is a general consequence in cases in which quantities drawn from the same distribution (the true values) are subject to an additional scatter; see Stein (1956) and Good (1965).
Treating this problem in a HBM framework allows us to separately incorporate both the measurement errors in all the parameters and the intrinsic scatter in the radius (which is a parameter to be inferred). Wolfgang et al. (2016), for example provided a brief overview of the advantages of HBMs for a very similar problem. By defining the probabilistic relations and placing priors on all the parameters, we can calculate the joint posterior distribution for R, F, M, and α, given the data. The forward model relation for R_{i}M_{i}, F_{i}, α has already been defined in Sect. 3; in this section, we define the following other probabilistic relations and priors (with their dependence relationships shown in Fig. 2):
where U(a, b) represents a uniform distribution with upper and lower bounds a and b, respectively, and SN(μ, σ, α) represents a skewnormal with the skewness parameter α (Azzalini 1985).
We use an uninformative prior for A_{j} (a Jeffreys prior for the slope of a line, see VanderPlas 2014) and essentially uninformative priors for σ_{R,j} and C_{j}. The cutoffs on σ_{R,j} and C_{j} are placed at values beyond which we can assume to have negligible likelihood, although this is arbitrary. An average cold radius of greater than 2 R_{J} or less than zero is clearly not supported by the data, and a scatter of less than 0.007 R_{J} is much less than expected based on possible variations in core mass and planet age. The posteriors are not found to have any significant probability density near these limits, so our choice is justified.
The mass bin boundaries are given a uniform prior that is constrained to prevent these bin boundaries from crossing each other. There are no planets withmass less than 0.1 M_{J} in our sample (see Sect. 2), thus we place a minimum cutoff for the lowest bound at 0.1 M_{J}. Among the model parametrisations that are tried, we find that in models with four bins (i.e. 3 boundaries), the highest boundary generally has a very wide distribution and would occasionally get stuck outside our mass range. However, we also find that when we fix the highest mass bin boundary to an arbitrary value larger than 2.0 M_{J}, we still obtain clearly distinct behaviour in each of the four mass bins (in terms of the radius response to inflation). In the interest of extracting the maximum amount of information from the data, we thus use four mass bins, but with the highest boundary (m_{3}) fixed at 2.5 M_{J}.
The prior for F_{s} is based on the results of Demory & Seager (2011); and Miller & Fortney (2011), which found that hotJupiter inflation is no longer obvious below an incident flux of ~ 2 × 10^{5}W m^{−2}. This constraint is only weakly informative and has a conservative 1σ interval between 2 × 10^{4} and 2 × 10^{6} W m^{−2}.
One of the main benefits of Bayesian techniques is the ability to set priors based on known physical constraints (for example a maximum density based on a fully iron giant planet). However in this case, we do not insert any physical priors. According to Fortney et al. (2007), radii of 0.5− 0.6 R_{J} are possible for low mass giants with heavy cores in a 50–50 rock/ice mixture and even smaller radii would be possible for pure rock and iron planets (down to less than 0.3 R_{J}), which is significantly lower than any observed radius in our sample. Thus, it is unlikely that adding a maximum density constraint would have affected our results. We also prefer not to insert constraints from interior models into our HBM to keep it purely observation driven, so we simply constrain the true mass, radii, and fluxes to be positive. Apart from the constraints on the mass bin boundaries, we do not find that our choice of priors significantly affects the posterior distributions for any of the other parameters.
The incident flux depends on the stellar effective temperature, orbital distance, and stellar radius as (6)
where σ is the Stefan–Boltzmann constant, a is the orbital distance (taken as the semimajor axis), and R_{⋆} and T_{⋆} are the stellar radius and temperature, respectively. For each planet, we calculate the distribution of flux based on those parameters, before solving the HBM. We find that the skewnormal distribution, described in Azzalini (1985), is a very good fit to the resulting asymmetric distributions of flux. We therefore fit a skewnormal distribution for each planet’s incident flux, and extract the bestfit parameters μ, σ, and the skewness α. This reduces our number of free parameters by nearly a factor of two, aiding in convergence time.
Our data come from a large range of programmes and sources, and we do not include the statistical biases introduced during observation. Observations currently favour planets with large radii and masses, however the nature of this bias depends on the instrument, and how close a parameter is to the detection threshold. The statistical biases in the flux parameter are also complicated; higher fluxes are favoured because of closer orbital distances, but disfavoured because of a smaller transit depth for larger stars. We must take this into consideration when forming our conclusions. There are also selection effects in the groundbased RV follow up that we are unable to model.
We findthe posterior distributions of our parameters by drawing samples from the joint posterior distribution , using a Markov chain Monte Carlo (MCMC) sampler, with a Metropolis–Hastings (M–H) algorithm (Hastings 1970; Chib & Greenberg 1995) implemented in the PyMC2 package in Python (Patil et al. 2010). Hamiltonian Monte Carlo samplers such the NoUTurn Sampler (Hoffman & Gelman 2011), implemented in PyMC3, cannot be implemented as a consequence of the discontinuous nature of our forward model. Other alternatives, such as the affine invariant ensemble sampler, implemented in emcee, are unsuitable owing to the nearly 900 dimensions in our model (Goodman & Weare 2010; ForemanMackey et al. 2013).
While increasing the number of hyperparameters can lead to a better fit to the data, it can also lead to overfitting. Increasing the model complexity gives our model more freedom, which can lead to multimodal posterior distributions that our M–H sampler cannot accurately sample, leading to unconverged posteriors. We fit multiple parametrisations and based our final choice firstly on whether it could converge and whether it had welldefined unimodal posteriors. To check for convergence, we run multiple MCMC chains and calculate the Gelman–Rubin convergence metric for our posterior samples, where values close to 1 indicate convergence and good mixing of the chains (Gelman & Rubin 1992).
To choose between converged models, we use two tools. Firstly, we calculate the Bayesian information criterion (BIC; Schwarz 1978), from (7)
where k is the number of hyperparameters, N is the number of samples (the number of planets), and the highest likelihood is defined as
where by marginalising out the second level of our model, we obtain a twolevel model, and Eq. (10) relies on the independence between observations. This is an assumption that we have to make as we cannot model the correlation between observations and between parameters with the publicly available information. The symbol refers to the model in question, and , , represent the model parameters set to their maximum likelihood values. A lower BIC indicates a better model, where differences of ~ 10 are requiredfor strong evidence in favour of a model (see Schwarz 1978). Because of the Gaussian approximation required for the BIC, which is not satisfied for all our posterior parameter distributions, the BIC only plays a minor role in our model selection. We only consider very large differences (more than 20 in the BIC) to be significant.
A better tool for model testing is to perform detailed statistical tests on our posterior predictive fit. We simulate the data that would be produced by our model and compare our model data to the real observed data. This essentially tests the question of if our fully marginalised model was correct, how likely would we be to produce the observed data. This model checking follows the procedure in Wolfgang et al. (2016) and is described in Sect. 6.
For thefinal baseline model, we decided to fix the highest mass boundary at 2.5 M_{J}. In models where the highest mass boundary was treated as an unknown parameter, its posterior distribution was very wide and poorly constrained, generally varying between 2 and 10 Jupiter masses. This is likely due in part to the sparsity of such high mass planets in the database. Fixing it in place greatly sped up convergence times, and allowed us to focus on the lower masses where there is greater variation.
Fig. 2 Our model shown as a Bayesian network. The forward model hyperparameters are indicated in light orange (the third level of the HBM), the true planet parameters in white (the second level of the HBM), and the fixed observed values in dark blue. The arrows denote conditional dependencies; the probability distribution of a child parameter is dependent on the values of its parent parameters (arrows go from parent to child). 
5 Results
5.1 Model posteriors
We ran 10 chains for 2 000 000 iterations each with a burnin of 1 000 000 iterations and a thin factor of 200. The Rubin–Gelman convergence metric was evaluated for the posterior distributions, and we achieved for all the model parameters, indicating a high level of convergence (Gelman & Rubin 1992).
The posterior distributions for our mass bin boundaries, m_{1} and m_{2} (see Fig. 4), and the fixed boundary at m_{3} = 2.5, split our sample into four mass regimes: planets with subSaturn masses (< 0.37 M_{J}), planets with subJupiter masses (0.37−0.98 M_{J}), and planets heavier thanJupiter (0.98−2.5 M_{J} and 2.5 + M_{J}). Above the inflation threshold, the radii have a clear dependence on incident flux for all but the subSaturn mass planets; the fluxradius relation for each mass bin is distinct in the sense that the joint posterior distributions for the hyperparameters that govern this relation, A, C, and F_{s}, have little overlap (see Fig. 3).
We plot the posteriors for C, A, and F_{s} in Fig. 3 for each mass bin and report the bestfit values (which we take to refer to the median sample values from here on) and error bars on the 16th and 84th percentiles for all the hyperparameters in Table 1. We find that radius inflation, Δ R = μ_{R} − C, governed primarily by A, is highly dependent on mass. The highest degree of inflation is found for planets with masses 0.37−0.98 M_{J}, and it decreases with mass for the 0.98−2.5 M_{J} and 2.5 + M_{J} regimes. From the bestfit values of our parameters, the radius dependence at F > F_{s} is (11)
For the lowest mass bin, M < 0.37 M_{J}, the posteriors are much broader and the scatter is much larger, as shown in Fig. 5. In Figs. 6 and 7, we see that the radiusflux distribution is more complicated than for the other mass bins, featuring a change in the behaviour of the radii at ~ 10^{6} W m^{−2}. A number of planets appear not to follow the trend of increasing radii with flux^{3}, and we also note a lack of inflated planets at fluxes higher than ~ 10^{6} W m^{−2} in this massrange. Our chosen model, which assumes that the planets in a mass range behave uniformly with regards to their radius relationship with flux, does not seem to be a good fit for the leastmassive planets in our sample. As a consequence of its large uncertainties however the model fit below 0.37 M_{J} may not be reliable. The result is still informative as it places a lower bound on the range of masses over which our model works. For a further discussion see Sect. 7.2.
We find that the mean low irradiation radius C is highest for the class of planets more massive than Jupiter (0.98−2.50 M_{J}), and decreases slightly for both heavier and lighter planets in line with what is expected from model predictions (Fortney et al. 2007; Baraffe et al. 2008). However, it would also be affected by any trends in the heavy element masses of the planets (see Miller & Fortney 2011; Thorngren et al. 2016).
The inflation threshold is only well constrained in 0.37−2.5 M_{J}, and the inflation mechanism requires higher fluxes to be activated when mass is increased, as F_{s} is about 1.3 times greater for the 0.98−2.50 M_{J} than for the 0.37−0.98 M_{J} gas giants. For 2.5 + M_{J} the distribution of F_{s} provides a poor constraint with an extended tail at low incident fluxes in Fig. 3. This is likely due to the lack of heavy planets at low irradiation in our sample. Similarly for 0.1−0.37 M_{J} there is also a poor constraint, with a long tail in Fig. 3, likely as a result of the poor fit of our model in this mass range.
Figure 5 shows the physical scatter posteriors, σ_{R}. Low mass planets (M < 0.37 M_{J}) have by far the largest scatter, at 0.21 R_{J}. For masses above 0.37 M_{J}, we find the scatter is slightly higher in the inflated planets above the flux threshold than for weakly irradiated planets (0.12 R_{J} vs. 0.10 R_{J}). Most of the proposed mechanisms for inflation would introduce their own latent parameters, which control the degree of inflation, for example opacity, tidal heating rate, and Ohmic heating depth. Variations in these parameters may add to the scatter already present from heavy element content (see e.g. Heng 2012).
Fig. 3 Posterior distributions of the A, C, and log F_{s} parameters for each mass regime. 
Bestfit values for the hyperparameters with error bars for the 16th and 84th percentiles.
Fig. 4 Posterior distribution of the two variable mass bin boundaries, m_{1} and m_{2}. We kept m_{3} fixed at 2.5 M_{J}. 
Fig. 5 Posterior distributions for the intrinsic physical scatter of the 3 regimes considered (see Eq. (5)). 
5.2 Marginalised posterior distribution
The functional form of the FMR relations, taking the bestfit parameter values, is shown in Fig. 6, with scatter limits as dotted lines. The points represent the raw observational dataset and its observational uncertainties.
Figure 7 represents a more complete picture of our posterior FMR relation. The points are the posterior true parameters and their uncertainties, which can be seen to be different from the observed values. The shaded regions show the posterior predictive distribution for the true radii and represent the parameter space within which we expect the central 68% of planets to be contained for a given incident flux and planet mass. The boundaries of the region represent the 1σ contour lines of the distribution as a function of F for specific masses^{4}. is our forward model marginalised over the posterior distribution of the hyperparameters, (12)
where is the posterior distribution of our hyperparameters. The dashed line represents the median point of the distribution in Eq. (12), and the dotted line is the central 95% coverage interval. In the Bayesian sense, Eq. (12) represents the most accurate statement of our knowledge of the FMR relation, having inferred the posterior distributions of all our model parameters using the available data, as it includes the uncertainties we have in those parameters. However, the bestfit distribution, , with from Table 1 is often a good enough approximation.
The values in Fig. 6 are more scattered than the true R_{i} in Fig. 7. This is expected (see in Sect. 4) since the scatter from observational noise is removed from the posterior distribution for R_{i}. While the form of model we chose affects these posteriors, they tend to always shrink towards a mean line. We justify this result in the model checking in Sect. 6. The posterior distribution should be seen as the distribution of true planet parameters most likely to reproduce the data when observational noise is applied to the values.
Fig. 6 Relations of our model plotted for parameters set at their bestfit values. The solid line is μ_{R} (F, M), and the dotted lines are μ_{R}(F, M) ± σ_{R}(F, M); in the bestfit model, 68% of true values of the planets should lie within. The points represent the data values of the planet parameters of our sample with observational uncertainty error bars. 
Fig. 7 Marginalised posterior distribution of radii given flux, , plotted for different planet masses; see Eq. (12). The shaded region is the central 68% interval within which 68% of planet true radii should lie for a given incident flux and planet mass. It represents the 1σ points of ; the central dashed line is the median and the dotted lines are the 95% interval. The red points represent the posterior true values of the planet parameters of our sample with posterior uncertainty error bars. 
6 Model checking
We performed detailed tests to justify our model and checked its consistency with the data using a posteriorpredictive test. Our aim was to test the ability of our model to reproduce datasets similar to the observed dataset and to quantify the discrepancy. Owing to the similar nature of our investigation, we followed the procedure of Wolfgang et al. (2016), which provides a more indepth explanation. For a detailed overview of hierarchical model checking, see also Bayarri & Castellanos (2007).
We calculate the posterior predictive distribution of model , defined as the probability of observing a new planet with parameters , given our original observed dataset,
where θ refers to the combined model parameters (R, M, F, α). Equation (14) is to be read as the expectation value of , for values of θ drawn from their posterior distribution, . This followsfrom Eq. (13) by the definition of the expectation value (mean in this case), and is of course trivial to calculate with the samples produced by our MCMC. From here on, we drop the notation and treat it as implicit.
To perform the posteriorpredictive test, we draw points from the posterior predictive distribution in Eq. (13) to create a mock dataset, , of the same size N as our planet sample. By drawing thousands of such mock datasets, we can compare these statistically to our observed dataset.
The comparison in this case involves quantifying certain aspects of the distributions in some derived statistics. We use the same two statistics as Wolfgang et al. (2016): i.e. f_{1σ}, the fraction of the simulated radii of a given dataset that fall within the 68% coverage interval (the shaded region in Fig. 7, described in Sect. 5.2), and f_{μ}, the fraction of data points that have radius 1σ error bars that cross the median of the posterior distribution (the dashed line in Fig. 7). While we use the posterior true distributions to aid in performing this comparison, it is important to remember that we are comparing the observed data, i.e. , produced from the true values (e.g. R), based on observational uncertainties.
The two statistics test two key attributes of our model. The value f_{1σ} is a measure of the scatter of data produced by our model (controlled mainly by σ_{R}), while f_{μ} is a measure of how centred the data is about the median line. The value f_{μ} is thus a proxy for the shape of the dataset distribution, and checks whether the normal distribution chosen in Eq. (2) is justified. As an example, if the scatter of our model proved to be accurate (in reproducing f_{1σ}), but f_{μ} was too small, it could suggest that we should have used a distribution with heavier tails and less clustering near the mean.
We must also address how to draw a posterior predictive dataset, x_{new}, from Eq. (13), where we also follow the treatment of Wolfgang et al. (2016). There are two methods, each tied to one level of our HBM, where the second level represents the true values (R, F, M) and the third level represents the hyperparameters α (see Fig. 2).
For the first method (drawing from the second level), we draw each mock dataset from the posterior true value distributions of the real planet sample. Thus for theith planet of our mock dataset, we draw x_{new,i} from
represents the published observational uncertainty distributions of that particular planet. Performing the test on the datasets produced above checks the posterior distributions of our true parameters for the real planet population (i.e. the points in Fig. 7).
The second method generates an entirely new population of hypothetical planet true values from the marginalised model and is written as
From that hypothetical planet sample, we generate a new dataset as in the first method, i.e. (20)
In this case, we draw the experimental uncertainties from distributions consistent with our data, conserving the trends between the parameter uncertainties and parameter values^{5}. Our hyperparameters do not parametrise the distributions of F and M for our population, thus we must first draw M and F from the same distributions as our population to use them in Eqs. (19) and (20). Performing the model checking statistics on the hypothetical datasets produced by this method is thus a test of our posteriors for the third level of our HBM, i.e. the hyperparameters.
We calculate f_{1σ} and f_{μ} for a large number of datasets produced by both the above methods, and also for the single real dataset. If our model is consistent with the data, we expect that the statistics calculated from the mock datasets and from the real data agree with each other. We define agreement as meaning that the statistic calculated for the real dataset is within the 1σ interval of the statistics calculated from the mock datasets. In other words, the real dataset must have the same properties as a typical mock dataset.
Figure 8 shows the results of our model checking for the second level of our HBM, and Fig. 9 shows the results for the third level. We measure modeldata consistency by computing the percentile of the generated mock datasets that the real datasets fall into (measured by the statistic of choice).
We find that the observed f_{μ} fall within the typical values produced by our model for datasets produced by both methods, coming in at the 50th and 82nd percentile for the second and third levels, respectively. The observed f_{σ} falls in the 73rd percentile of data generated from the second level of our HBM, indicating that the spread of our posterior true values (R_{i}, M_{i}, F_{i}) is consistent with the data, and justifying the shrinking effect seen in the points of Fig. 7. For the data generated by the third level of our HBM, f_{μ} falls in the 82nd percentile, while f_{σ} falls in the 97th percentile. In other words, only ~6% of the generated mock data were more extreme than our observed data. The result suggests that the form of our model may have room for improvement, as we overestimate the scatter of the data. However, it does not negatively affect our conclusions in Sect. 7.3.
Fig. 8 Distribution of f_{1σ} (left) and f_{μ} (right) for 10 000 generated datasets from the second level of our HBM (the true posteriors, see Eq. (16)). The comparison to the values for the real dataset is represented by the vertical line. The real data falls within the 73rd and 50th percentiles of the generated distributions for f_{1σ} and f_{μ}, respectively. 
Fig. 9 Distribution of f_{1σ} (left) and f_{μ} (right) for 50 000 generated datasets from the third level of our HBM (the hyperparameters, see Eqs. (19) and (20)). The comparison to the values for the real dataset are represented by the vertical line. The real data falls within the 97th and 82rd percentiles of the generated distributions for f_{1σ} and f_{μ}, respectively. 
7 Discussion
7.1 Fluxmassradius relation
Our results show a strong trend of decreasing radius inflation with mass; the most inflated planets are in the mass range 0.37−0.98 M_{J} (see Fig. 1). This effect is also expected from models without additional inflation mechanisms, in which mass increases past a few Jupiter masses lead a to decreasing radius response to incident flux (Fortney et al. 2007; Baraffe et al. 2008). Below the inflationthreshold, we also find that planet mass has a small effect and that the largest cold giant radii have masses of 0.98−2.50 M_{J}.
Thorngren et al. (2016) have found that the heavy element content of noninflated gas giants is strongly and positively correlated with planet mass. Unless there is a different formation mechanism behind the hottest Jupiters, or in cases of significant evaporation, this trend should apply in our sample as well. This would lead to a further reduction in the radii as planet mass increases and would agree with the observed trends for inflated hot Jupiters. As more detailed predictions from models with inflation become available in the future (with predictions given for wide ranges of parameters such as in Fortney et al. 2007), we will similarly be able to probe the core mass distribution of inflated hot Jupiters.
The incidence of using binning in the mass parameter must also be examined. To do this we look at the radius residuals, Δ R_{i} = R_{i} − μ_{R}(F_{i}, M_{i}), and their relation with mass. We fit a simple linear regression to the relationship, with an additional scatter parameter, to determine whether a massradius gradient of zero is excluded to more than 1σ. We find that within the 0.10−0.37 M_{J} bin there is a strong positive correlation of radius with mass, and a negative correlation for the 2.5 + M_{J} planets (seeFig. 10), and for both a gradient of zero is more than 1σ from the mean of the distribution. This additional dependence with mass may be the reason for the large scatter in the lowest mass bin, M < 0.37. We find no evidence of a massradius correlation for the middle two mass bins. Future models could take this into account.
Extending this analysis to other parameters, we also find a weak negative correlation of radius with metallicity for all mass bins (see Fig. 10), albeit with a very large scatter. This could be caused by a correlation between heavy element content and stellar metallicity (Guillot et al. 2006; Miller & Fortney 2011). However, we find no correlation of the radius with host star age. Our results agree broadly with Enoch et al. (2012).
We also performed a short test to determine the firstorder bias introduced by averaging our observational error bars where they are asymmetric (see Sect. 2). The asymmetry in the uncertainties is most significant for the radii, where it is biased towards larger upper error bars. To determine the bias, we fit a simple threelevel hierarchical model to the residual observed radii, , treating thepopulation as being drawn from a normal distribution with mean μ and variance σ^{2}, and having the same measurement uncertainties as the data. We then take two cases: one where we average the upper and low uncertainties as before and one where we allow them to vary. For the latter, we use a discontinuous split normal distribution (21)
to simulate asymmetric uncertainty distributions with upper and lower error bars σ_{u} and σ_{l} respectively. Comparing the results between the two cases, we find that the difference in our posteriors for σ and μ are 0.05% and 0.1%, respectively. This is well within the error bars on the posteriors in this simple HBM and within the errors bars we present for our hyperparameters in Sect. 5, which are all of the order of 1− 10%.
Fig. 10 Radius residuals, ΔR_{i} = R_{i} − μ(F_{i}, M_{i}), plotted against host star metallicity (top), and planet mass for the lowest mass bin (middle) and the highest mass bin (bottom). The shaded grey region is the marginalised 68% coverage interval of our linear regression, and the shaded blue region represents the 1σ limits of themean line. 
7.2 Lowmass inflation cutoff
A key finding we present is the inference of a boundary at ~ 0.37 M_{J}, above which we see the most inflated gas giants, with a clear positive dependence on incident flux. Below we find a lack of highly inflated radii with R > 1.5 R_{J} and a more complicated trend (see Figs. 6 and 7).
From visual inspection of the data in the lowest mass bin, we see the beginning of a trend of inflating radii as a function of flux ~ 10^{5} W m^{−2}, particularly for planets closer to 0.37 M_{J}. However this stops at ~ 10^{6} W m^{−2}, and there are no inflated planets beyond this threshold. Instead we see a reversal, with decreasing radii and denser planets, although we have a small sample size in that regime (five planets). Giants at the top end of the mass range, which would otherwise be inflated, also show decreasing radii (e.g. HD 149026 b has a mass of 0.37± 0.03 M_{J} and evolutionary models find it to be very dense, with a predicted core of 80− 110 M_{⊕}; Burrows et al. 2007). Previous studies about hotJupiter inflation (e.g. Laughlin et al. 2011; Enoch et al. 2012) generally use manually chosen minimum mass cutoffs for their samples. However, we recommend that future studies focus on planets heavier than ~ 0.37 M_{J}, to avoid contaminating the samples with planets where other processes may start to play a role alongside inflation.
In Fig. 11 we plot the incident flux of all the planets in our sample that have radii greater than 1 R_{J} (and are thus potentially inflated) and find that there is a maximum incident flux that depends on the planet mass and decreases as we decrease themass.
To explain the reversal of the inflation trend and the lack of such inflated low mass giants at F > 10^{6} W m^{−2}, an obvious cause could be bulk evaporation or Rochelobe overflow leading to significantly denser planets if the outer gaseous layers are stripped (see e.g. Kurokawa & Nakamoto 2014). This could also lead to such evaporated giants leaving remnants too low in mass to enter our sample (M < 0.1 M_{J}). Rochelobe overflow and asymmetric evaporation, such as described in Gu et al. (2003) and Baraffe et al. (2005), could also stop planetary migration beyond a certain point, or even lead to outward migration, which may explain the paucity of planets in this mass range closer to their star (relative to the other mass ranges).
We must also consider the observational biases that may be present. There are few data points available in this regime, however our conclusion relies not only on the trend of overdense planets closer to the star, but also on the clear lack of R > 1.0 R_{J} gas giants at incident fluxes greater than ~10^{6} W m^{−2}. Below this threshold, there are numerous such large giants, and we would expect that as orbital distance is decreased and incident flux increased, the probability of detecting similar objects should increase. Similarly, larger radii should also be easier to detect, thus if observational surveys are able to find planets with R < 1.0 R_{J} in the mass ranges that we examine, we would also expect these surveys to be able to find more inflated planets (if they exist). Thus it is unlikely that our conclusions are impacted by an observational bias.
Fig. 11 Massflux distribution for all gas giants that have a radius greater than 1 R_{J}. The blue dashed lines show our bestfit mass boundaries. As the mass decreases below 1 M_{J}, the maximum incident flux at which we find inflated planets also decreases (black line). 
7.3 Lack of noninflated hot Jupiters
Of the proposed mechanisms for radius inflation, many can only work under specific circumstances. For example, tidal dissipation requires nonzero eccentricity, unless the eccentricity can be excited (Arras & Socrates 2010) and Ohmic heating requires heat deposition deep enough near the radiativeconvective boundary; see e.g. Heng (2012); Perna et al. (2012); Ginzburg & Sari (2016). Thus, not only may we expect some planets not to be inflated, but the varying degrees of inflation caused by each mechanism should leave an imprint on the distribution of radii and its scatter. As an example, we make a simple comparison of our posterior FMR distribution with planetary models that do not include inflation mechanisms (Fortney et al. 2007) to determine what fraction of radii could be explained without resorting to inflation.
To compare our observationbased results with theory, we must not only account for the mass and incident flux, but also for planet age and the unknown core mass. We assume that planet ages are distributed uniformly in the range 0.316−8.0 Gyr (in rough agreement with the age distribution of our planet sample). The core mass, or heavy element fraction, is not an observable quantity. However Thorngren et al. (2016) showed that there exists a strong correlation between core mass and planet mass at low insolation, and inferred the relation through a regression fit. For the purposes of this work, we assume that the inferred core mass (M_{C}) to planet mass (M) relation of Thorngren et al. (2016) holds equally well for highly irradiated planets. Their relation is denoted as a probability distribution ^{6}, for the highest likelihood point of their model parameters .
Fortney et al. (2007) gave a detailed table of predictions of planet radii as a function of mass, incident flux (orbital distance), core mass, and planet age. However, using the Thorngren et al. (2016) relation, planets heavier than ~1 M_{J} can have core masses greater than the domain of predictions in Fortney et al. (2007). Thus, we can only carry out the following analysis for the lowest inflated mass bin, 0.37−0.98 M_{J}. There is also a systematic error in this comparison, as the Fortney et al. (2007) models considered the heavy elements to be contained fully in a rockice core, while the Thorngren et al. (2016) models have the heavy elements also distributed in the envelope. Models with all the heavy elements contained in the core tend to produce larger radii than models with heavy elements mixed throughout the envelope, so our model radii are overestimated (see Baraffe et al. 2008; Thorngren et al. 2016).
We produce a distribution of model radii as a function of incident flux, p_{M}(RF), fully marginalised over the relevant distributions of core mass and planet age (t) and marginalised over the width of the mass bin, m_{i} − m_{i+1}, (22)
where f(M, M_{C}, t, F) refers to the model predictions of Fortney et al. (2007), which are deterministic, δ is the delta distribution, and U(t) is the uniform distribution for the planet age. Thus we calculate the model prediction for the distribution of radius as a function of incident flux, similar to our posterior FMR distribution. We compare the model prediction against our observationdriven posterior FMR distribution in Fig. 12.
This is intended as a rough and qualitative exploration of how we could compare marginalised theoretical models to uncertainty decoupled data. We see that below the flux threshold, our theorypredicted p_{M} (RF) overestimates the radius, but falls within the 1σ interval. This could be because we underestimated the core masses or planet age, but may also be caused by the previously mentioned increase in radii, when all the heavy elements are placed into a pure rockice core, compared to being distributed in an envelope.
With that in mind, we look at the predictions above ~ 10^{6} W m^{−2}. The posterior FMR relation and the model prediction drastically diverge, with negligible overlap beyond F ~ 1.6 × 10^{6} W m^{−2}. Only one real planet, HATS9 b (Brahm et al. 2015), has a posterior density that has significant overlap with p_{M} (RF), and we cannot rule out the fact that it may be an inflated planet with an unusually large core. Such overdense planets have been found in the lowflux regime, for example see Kepler539 b, with total mass 1.00 M_{J} and heavy element mass 0.49 M_{J} (Thorngren et al. 2016). Furthermore, what little overlap exists is only significant above the mean line of p_{M} (RF) (solid blue line). Below the mean line, where we would expect to find denser planets with large cores, there is negligible overlap, even with the 95% interval of our the FMR relation of our HBM.
While a more detailed study would be warranted, our results suggest that classical planets that do not exhibit radius inflation do not exist above fluxes of ~10^{6} W m^{−2} for this mass range. Even if some of the observed planets could be explained as being young noninflated gas giants with light cores, if a population of such young lowdensity planets were to exist, we would also expect to find similar planets with mediummass or heavy cores and older ages (which would fall below the mean line of p_{M} (RF)).
Although we found a slight discrepancy between our hyperparametermarginalised HBM posteriors and the data in Sect. 6, this discrepancy suggests that our HBM is overestimating the intrinsic scatter. With a lower scatter, the overlap between the observed data and the model prediction would be decreased. Although we found that we were overestimating the radii at low irradiation, increasing the input core masses to counter this would only serve to further drive p_{M} (RF) away from our FMR distribution at high fluxes.
If reinflation of hot Jupiters is impossible (see Wu & Lithwick 2013; Ginzburg & Sari 2016), this would mean that all of the current population of highly irradiated 0.37−0.98 M_{J} giants must have migrated to their present locations very early in their lifetimes. This could however be avoided with tidal reinflation (e.g. Jermyn et al. 2017) followed by delayed contraction from Ohmic heating or another mechanism.
Furthermore, the findings of Heng (2012) showed that variations in atmospheric scattering and absorption should produce a large scatter in the degree of Ohmic dissipation. In certain cases (such as atmospheres with temperature inversions and high optical opacities, or planets with too weak a magnetic field) we may even expect no significant Ohmic dissipation at depth (Heng 2012; Perna et al. 2012). This should be reflected in the radius distribution of hot Jupiters. Thus if Ohmic heating is the dominant inflation process, the lack of noninflated hot Jupiters must be explained by either a lack of high optical opacities and temperature inversions in this temperature regime, or would require a different inflation mechanism for the cases where Ohmic heating might fail.
We also briefly consider the observational biases that may be present in our data. While there is a detection bias towards finding larger radii and masses that we cannot fully incorporate into our model, since we use various sources of data for which we may not have welldocumented completeness functions, we note that the noninflated planets that we are comparing with have roughly the same radii as gas giants found below the inflation threshold. Thus surveys that can detect the giants below the flux threshold should detect closerin gas giants of the same radii as well. Planets closer to their star should also have a higher chance of transiting, and thus being detected, so we should be biased towards finding planets in the high incident flux regime.
A more thorough treatment of the observational biases would be the next step, as well as a comparison with newer models that include inflation mechanisms. Within the framework of such a hierarchical model, we could also extract a more quantitative measure of the percentage of gas giants that are noninflated using a mixture model (see McLachlan & Peel 2000; Celeux 2007). Another approach, applied recently by Thorngren & Fortney (2018), is to infer the heating required without specifying an underlying mechanism. These authors find that the functional form of the heating rate ϵ(F) required to explain the radii most closely matches the profile predicted by Ohmic heating. Given our findings that all 0.37−0.98 M_{J} hot Jupiters beyond ~10^{6} W m^{−2} show inflation, it would be informative to further explore the scatter of ϵ(F) at a particular flux within their framework.
Fig. 12 Our marginalised posterior radiusflux distribution for , with the central 68% interval shaded in dark grey, and the 95% interval in light grey. Plotted against a fluxradius relation p_{M} (RF) predicted by a noninflationary planet model, with the central 68% interval shaded in blue. 
8 Conclusions
We have constructed a theoryindependent HBM to constrain the probabilistic relation between the planet radius, mass, and incident fluxfor 286 gas giants. We thus find the intrinsic scatter of the gas giant population, decoupled from the observational uncertainties in both the radius parameter, and the mass and incident flux parameters. The posterior distribution we inferred may be used to test and constrain theoretical models for gas giant radius evolution, especially involving inflation mechanisms. Constraining the intrinsic scatter allows us to determine the real range of radii that need to be reproduced by theory at each regime of incident flux and mass.
Our keyresults are summarised below:

We find that planetary mass plays a significant role in the degree of radius inflation, with the most inflated hot Jupiters coming from 0.37−0.98 M_{J} range, and that the response of radius to incident flux decreases as mass is increased. This could be caused by either a trend of increasing core masses with increasing total mass, similar to that found for gas giants below the inflation threshold (see Thorngren et al. 2016), or by the increase in surface gravity with mass, which would make them harder to inflate.

Below 0.37 Jupiter masses, there is an abrupt lack of heavily inflated radii, and we find that radii begin to decrease with incident fluxes near ~10^{6} W m^{−2}. Such atrend of decreasing radii is not present at higher masses. We also note that there is a cutoff point in incident flux beyond which hot Jupiters no longer exist and that for inflated planets below 1.0 M_{J} this cutoff point decreases with decreasing mass.

We use our inferred posterior distribution to show that there is no evidence for noninflated hot Jupiters at fluxes greater than F ~ 1.6 × 10^{6} W m^{−2} and masses of 0.37−0.98 M_{J}. In this case we define noninflated to mean that their radii could be reproduced by an evolutionary model that only considers incident flux energy deposition in the photosphere with no additional inflationary effects (Fortney et al. 2007).
To extend this study, a mixture model could be a useful tool to provide a more quantitative measure of the fraction of gas giants that do not exhibit any inflation, and could be incorporated within a Bayesian hierarchical model (see McLachlan & Peel 2000; Celeux 2007). With increasingly detailed inflated model predictions given for wide ranges of latent parameters (such as those in Fortney et al. (2007) for core mass and age), we will begin to extract information about hotJupiter inflation from the entire observed distribution and its scatter. The forward model of our HBM, p(RM, F, α), may be replaced by a theoretically driven model for model selection purposes, and/or to constrain latent parameters in such models. Using a mixture model, we could also find the fraction of objects that are consistent with a particular model, and in which regions of mass, incident flux, and other parameters this may be the case. With the increasingly large number of discovered exoplanets, such statisticsbased studies will continue to become more important.
Acknowledgements
We thank the referee, Mathieu Havel, for a comprehensive review that improved our paper. We thank Kevin Heng for valuable discussions and ideas. M.S acknowledges support from the Swiss National Science Foundation (PP00P2163967). B.O.D. acknowledges support from the Swiss National Science Foundation in the form of a SNSF Professorship (PP00P2163967). This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. Calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.
References
 Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989 [NASA ADS] [CrossRef] [Google Scholar]
 Almenara, J. M., Damiani, C., Bouchy, F., et al. 2015, A&A, 575, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anderson, D. R., Smith, A. M. S., Lanotte, A. A., et al. 2011, MNRAS, 416, 2108 [Google Scholar]
 Arras, P., & Socrates, A. 2010, ApJ, 714, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Azzalini, A. 1985, Scand. J. Stat., 12, 171 [Google Scholar]
 Baraffe, I., Chabrier, G., Barman, T. S., et al. 2005, A&A, 436, L47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baraffe, I., Chabrier, G., Fortney, J., & Sotin, C. 2014, Protostars and Planets VI (Tucson: University of Arizona Press), 763 [Google Scholar]
 Batygin, K., & Stevenson, D. J. 2010, ApJ, 714, L238 [NASA ADS] [CrossRef] [Google Scholar]
 Bayarri, M. J., & Castellanos, M. E. 2007, Statist. Sci., 22, 322 [CrossRef] [Google Scholar]
 Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466 [NASA ADS] [CrossRef] [Google Scholar]
 Bodenheimer, P., Laughlin, G., & Lin, D. N. C. 2003, ApJ, 592, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Bonomo, A. S., Sozzetti, A., Lovis, C., et al. 2014, A&A, 572, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brahm, R., Jordán, A., Hartman, J. D., et al. 2015, AJ, 150, 33 [Google Scholar]
 Burrows, A., Hubeny, I., Budaj, J., & Hubbard, W. B. 2007, ApJ, 661, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Butler, R. P., Wright, J. T., Marcy, G. W., et al. 2006, ApJ, 646, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Celeux, G. 2007, Mixture Models for Classification, eds. R. Decker, & H. J. Lenz (Berlin, Heidelberg: Springer), 3 [Google Scholar]
 Chabrier, G., & Baraffe, I. 2007, ApJ, 661, L81 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. H. 2005, ArXiv eprints [arXiv:astroph/0509798] [Google Scholar]
 Chib, S., & Greenberg, E. 1995, Am. Stat., 49, 327 [Google Scholar]
 Delrez, L., Santerne, A., Almenara, J.M., et al. 2016, MNRAS, 458, 4025 [NASA ADS] [CrossRef] [Google Scholar]
 Demory, B.O. 2014, ApJ, 789, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Demory, B.O., & Seager, S. 2011, ApJS, 197, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Enoch, B., Collier Cameron, A., & Horne, K. 2012, A&A, 540, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 Fortney, J. J., & Nettelmann, N. 2010, Space Sci. Rev., 152, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Baraffe, I., & Militzer, B. 2010, Giant Planet Interior Structure and Thermal Evolution, ed. S. Seager, 397 [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Ginzburg, S., & Sari, R. 2016, ApJ, 819, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Good, I. J. 1965, The Estimation of Probabilities: An Essay on Modern Bayesian Methods, Research Monograph No. 30 (Cambridge, USA: MIT Press) [Google Scholar]
 Goodman, J., & Weare, J. 2010, Comm. App. Math. Com. Sci., 5, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Gu, P.G., Lin, D. N. C., & Bodenheimer, P. H. 2003, ApJ, 588, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Guillot, T., & Showman, A. P. 2002, A&A, 385, 156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guillot, T., Santos, N. C., Pont, F., et al. 2006, A&A, 453, L21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Han, E., Wang, S. X., Wright, J. T., et al. 2014, PASP, 126, 827 [NASA ADS] [CrossRef] [Google Scholar]
 Hastings, W. 1970, Biometrika, 57, 97 [Google Scholar]
 Heng, K. 2012, ApJ, 748, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Hoffman, M. D., & Gelman, A. 2011, ArXiv eprints [arXiv:1111.4246] [Google Scholar]
 Hogg, D. W., Myers, A. D., & Bovy, J. 2010, ApJ, 725, 2166 [NASA ADS] [CrossRef] [Google Scholar]
 Jermyn, A. S., Tout, C. A., & Ogilvie, G. I. 2017, MNRAS, 469, 1768 [NASA ADS] [CrossRef] [Google Scholar]
 Kelly, B. C. 2007, ApJ, 665, 1489 [NASA ADS] [CrossRef] [Google Scholar]
 Kurokawa, H.,& Nakamoto, T. 2014, ApJ, 783, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Kurokawa, H.,& Inutsuka, S.I. 2015, ApJ, 815, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Laughlin, G., Crismani, M., & Adams, F. C. 2011, ApJ, 729, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Maxted, P. F. L., Anderson, D. R., Collier Cameron, A., et al. 2016, A&A, 591, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McLachlan, G. J. & Peel, D. 2000, Finite mixture models, (New York, Chichester: John Wiley & Sons, Inc.), 419 [Google Scholar]
 Miller, N., & Fortney, J. J. 2011, ApJ, 736, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Patil, A., Huard, D., & Fonnesbeck, C. 2010, J. Stat. Softw., 35, 1 [CrossRef] [EDP Sciences] [Google Scholar]
 Perna, R., Menou, K., & Rauscher, E. 2010, ApJ, 724, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Perna, R., Heng, K., & Pont, F. 2012, ApJ, 751, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, L. A. 2015, ApJ, 801, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Santerne, A., Bonomo, A. S., Hébrard, G., et al. 2011, A&A, 536, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, G. 1978, Ann. Statist., 6, 461 [Google Scholar]
 Shabram, M., Demory, B.O., Cisewski, J., Ford, E. B., & Rogers, L. 2016, ApJ, 820, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Spiegel, D. S., Burrows, A., & Milsom, J. A. 2011, ApJ, 727, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, C. 1956, Proc. 3rd Berkeley Sympos. Math. Stat. Prob. (Berkeley, CA: University of California Press), 197 [Google Scholar]
 Thorngren, D. P., & Fortney, J. J. 2018, AJ, 155, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Thorngren, D. P., Fortney, J. J., MurrayClay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Torres, G., Winn, J. N., & Holman, M. J. 2008, ApJ, 677, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Chabrier, G., Mayne, N. J., et al. 2017, ApJ, 841, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Van Eylen, V., Albrecht, S., Gandolfi, D., et al. 2016, AJ, 152, 143 [NASA ADS] [CrossRef] [Google Scholar]
 VanderPlas, J. 2014, ArXiv eprints [arXiv:1411.5018] [Google Scholar]
 Weiss, L. M., Marcy, G. W., Rowe, J. F., et al. 2013, ApJ, 768, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y., & Lithwick, Y. 2013, ApJ, 763, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Mitchell, J. L. 2010, ApJ, 721, 1113 [NASA ADS] [CrossRef] [Google Scholar]
The five planets are HD 149026 b (Butler et al. 2006; Torres et al. 2008), K239 b (Van Eylen et al. 2016), Kepler101 b (Bonomo et al. 2014), Kepler41 b (Santerne et al. 2011; Butler et al. 2006), and WASP126 b (Maxted et al. 2016).
The masses of each plot are fixed and chosen to be in the middle of each mass bin (away from the boundaries). Near the boundaries, the scatter would be larger (and the shaded region wider) as it transitions from one bin to the next. This is because marginalising over the hyperparameters, including the mass bin boundaries, means that a particular mass could lie in either of the two adjacent mass bins with some probability. In essence, the RF relations plotted above vary smoothly with mass (except at the 2.5 M_{J} boundary), although only cross sections are shown in Fig. 7.
The values σ_{M,obs} and σ_{F,obs} are strongly correlated with M and F, respectively. We thus fit a linear regressive model (including a dispersion term) for log σ_{M,obs} ~ N(f_{1}(log M), σ) and log σ_{F,obs} ~ N(f_{2}(log F), σ), and we draw σ_{M,obs} and σ_{F,obs} from the resulting highest likelihood distribution. We also find that σ_{R,obs} is correlated with F. To assign σ_{R,obs}, we split the range of fluxes into four bins, and draw σ_{R,obs} from a random planet that is in the same F bin.
The fitted model in Thorngren et al. (2016) has the form M_{C} ∕ M_{⊕} = k × 57.9 × M^{0.61}, where the multiplicative scatter, k is drawn from log_{10}k ~ N(μ = 0, σ = log_{10}1.83).
All Tables
Bestfit values for the hyperparameters with error bars for the 16th and 84th percentiles.
All Figures
Fig. 1 Top panel: distribution of the transit radii in our sample vs. the incident flux. The mass is colour coded. The radii begin to show a correlation with flux at ~ 2 × 10^{5} W m^{−2}, corresponding to where inflated hot Jupiters begin to appear (Demory & Seager 2011). Bottom panel: the transit radii plotted against mass, colour coded by incident flux. 

In the text 
Fig. 2 Our model shown as a Bayesian network. The forward model hyperparameters are indicated in light orange (the third level of the HBM), the true planet parameters in white (the second level of the HBM), and the fixed observed values in dark blue. The arrows denote conditional dependencies; the probability distribution of a child parameter is dependent on the values of its parent parameters (arrows go from parent to child). 

In the text 
Fig. 3 Posterior distributions of the A, C, and log F_{s} parameters for each mass regime. 

In the text 
Fig. 4 Posterior distribution of the two variable mass bin boundaries, m_{1} and m_{2}. We kept m_{3} fixed at 2.5 M_{J}. 

In the text 
Fig. 5 Posterior distributions for the intrinsic physical scatter of the 3 regimes considered (see Eq. (5)). 

In the text 
Fig. 6 Relations of our model plotted for parameters set at their bestfit values. The solid line is μ_{R} (F, M), and the dotted lines are μ_{R}(F, M) ± σ_{R}(F, M); in the bestfit model, 68% of true values of the planets should lie within. The points represent the data values of the planet parameters of our sample with observational uncertainty error bars. 

In the text 
Fig. 7 Marginalised posterior distribution of radii given flux, , plotted for different planet masses; see Eq. (12). The shaded region is the central 68% interval within which 68% of planet true radii should lie for a given incident flux and planet mass. It represents the 1σ points of ; the central dashed line is the median and the dotted lines are the 95% interval. The red points represent the posterior true values of the planet parameters of our sample with posterior uncertainty error bars. 

In the text 
Fig. 8 Distribution of f_{1σ} (left) and f_{μ} (right) for 10 000 generated datasets from the second level of our HBM (the true posteriors, see Eq. (16)). The comparison to the values for the real dataset is represented by the vertical line. The real data falls within the 73rd and 50th percentiles of the generated distributions for f_{1σ} and f_{μ}, respectively. 

In the text 
Fig. 9 Distribution of f_{1σ} (left) and f_{μ} (right) for 50 000 generated datasets from the third level of our HBM (the hyperparameters, see Eqs. (19) and (20)). The comparison to the values for the real dataset are represented by the vertical line. The real data falls within the 97th and 82rd percentiles of the generated distributions for f_{1σ} and f_{μ}, respectively. 

In the text 
Fig. 10 Radius residuals, ΔR_{i} = R_{i} − μ(F_{i}, M_{i}), plotted against host star metallicity (top), and planet mass for the lowest mass bin (middle) and the highest mass bin (bottom). The shaded grey region is the marginalised 68% coverage interval of our linear regression, and the shaded blue region represents the 1σ limits of themean line. 

In the text 
Fig. 11 Massflux distribution for all gas giants that have a radius greater than 1 R_{J}. The blue dashed lines show our bestfit mass boundaries. As the mass decreases below 1 M_{J}, the maximum incident flux at which we find inflated planets also decreases (black line). 

In the text 
Fig. 12 Our marginalised posterior radiusflux distribution for , with the central 68% interval shaded in dark grey, and the 95% interval in light grey. Plotted against a fluxradius relation p_{M} (RF) predicted by a noninflationary planet model, with the central 68% interval shaded in blue. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.