Separating Super-puffs versus Hot Jupiters among Young Puffy Planets

We select confirmed planets younger than 1 Gyr with radii larger than 6R_⊕ and incident bolometric flux above 3.17 × 10⁶ erg s⁻¹ cm⁻² (∼2.3 × Earth insolation flux) observed by TESS (G. R. Ricker et al. 2015), Kepler (W. J. Borucki 2016), and K2 (S. B. Howell et al. 2014) from the NASA Exoplanet Archive (2023a, 2023b, 2024), hereafter "Archive." We further remove Kepler-1976 and Kepler-302 from our sample. Their default ages are drawn from the nominal value in the Q1–Q8 KOI table (C. J. Burke et al. 2014). The age of Kepler-302 has since been updated to 6.61 Gyr (T. D. Morton et al. 2016). Recently, L. G. Bouma et al. (2024) revisited the ages of planetary systems in the Kepler field younger than 4 Gyr and found no detection of stellar rotation for Kepler-1976 or Kepler-302, suggesting that their ages are likely beyond a few Gyr, too old for our selection criteria.

Of the 19 planet-hosting stars in our sample, three have age estimates below 100 Myr: V1298 Tau, HIP 67522, and TOI-837. These youngest stars are the ones for which it is most important to have a reliable age estimate. These three stars all belong to young clusters with known ages. V1298 Tau, part of the stellar association Group 29, also has lithium content and rotational period measurements that confirm and further constrain the stellar age (A. Suárez Mascareño et al. 2021). For HIP 67522, part of the Scorpius-Centaurus OB association, the effective temperature and luminosity of the star are fit to stellar isochrones PARSEC 1.2 s (A. Bressan et al. 2012) and BHAC15 (I. Baraffe et al. 2015) to determine the stellar age and mass (A. C. Rizzuto et al. 2020). The age estimate for TOI-837, part of the open cluster IC 2602, is taken to be an absolute range encompassing the results from age estimates obtained through different stellar isochrone fittings and lithium aging techniques (see Table 3 of L. G. Bouma et al. 2020). Given that the quoted age estimate for each of the above systems is verified against multiple independent techniques, we consider the ages of these stars to be reliable and continue with our analysis.

Some planets have multiple parameter sets available, so we prioritize data sets with planetary mass measurements and otherwise choose the default parameter from the Archive. We further require that all planets studied have the necessary information available to obtain incident bolometric flux values, such as the bolometric luminosity of the host star and the planet's semimajor axis. If these two parameters are not available, then we estimate the stellar luminosity from the stellar effective temperature and radius using the Stefan–Boltzmann law. Likewise, we calculate the planet's semimajor axis from the orbital period and stellar mass using Kepler's third law. We obtain an uncertainty on all flux values by propagating the errors on the measured values.

Based on our selection criteria, we retrieve 20 planets with known mass measurements and two planets either without a mass measurement or only with upper limits on their mass, as illustrated in Figure 1. We label hot Jupiters as those with incident bolometric flux ≳1.72 × 10⁸ erg s⁻¹ cm⁻², equivalent to solar insolation flux at 10 days. We additionally label warm Jupiters as those with incident bolometric flux that ranges between 1.72 × 10⁸ and 3.17 × 10⁶ erg s⁻¹ cm⁻², equivalent to solar insolation flux at 200 days. The final selection of planets and their host stars are shown in Tables 1 and 2.

Zoom In Zoom Out Reset image size

With the target list ready, we now use thermal evolution models from D. P. Thorngren et al. (2016) to constrain the masses of the planets given their radius, incident bolometric flux, and age. The thermal evolution of a planetary interior is constructed by solving the structure equations (hydrostatic equilibrium, mass conservation, and energy conservation, in that order):

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial m}=-\displaystyle \frac{{Gm}}{4\pi {r}^{4}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial r}{\partial m}=\displaystyle \frac{1}{4\pi {r}^{2}\rho },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial L}{\partial m}=-T\displaystyle \frac{\partial S}{\partial T},\end{eqnarray} \tag{ 3 }$$

where P is the pressure, m is the mass, r is the radius, ρ is the density, L is the luminosity, T is the temperature, S is the entropy, and G is the gravitational constant. We therefore solve for the radius–mass relationship for a given metal mass, gas-to-metal mass ratio, incident flux, and age. We make the distinction between core mass and metal mass since not all heavy elements go into the planet's core. While giant planets require a ∼10 M_⊕ core to nucleate by runaway gas accretion (e.g., J. B. Pollack et al. 1996; R. R. Rafikov 2006; E. J. Lee et al. 2014), concomitant solid accretion and/or post-formation pollution by solid infall can increase the total metal mass content. Furthermore, observational data indicate that sub-Neptunes have rocky cores with masses that can go up to ∼20 M_⊕ (e.g., J. F. Otegi et al. 2020). Therefore, if the planet is composed of less than or equal to 20 M_⊕ of heavy elements, all metals go into the core and we consider a H/He envelope with solar metallicity. For planets with over 20 M_⊕ of heavy elements, we put 20 M_⊕ into the core and uniformly mix the remaining metals into the planet's H/He atmosphere using additive volumes. ⁵ The final grid covers metal masses between 0.3 M_⊕ and 200 M_⊕ in 60 logarithmic bins, gas-to-metal mass ratio between 10⁻⁶ and 10³ in 33 logarithmic bins, incident fluxes between 10⁴ and 10¹⁰ erg s⁻¹ cm⁻² in 19 logarithmic bins, and ages from 10 Myr to 10 Gyr in 100 logarithmic bins.

Since many of our selected planets are under intense incident flux (see Table 1), we account for extra heating as described in D. P. Thorngren & J. J. Fortney (2018), where some of the irradiation from the star is being converted into heating the planet. The effects of this extra heating are illustrated in Figure 2. When the incident flux is high, the extra heating effectively puffs up the planet, resulting in a larger radius for the same mass. This effect is more pronounced for smaller metal masses. Since these planets have less heavy elements, the H/He envelope is more susceptible to expansion due to weaker gravitational potential.

Zoom In Zoom Out Reset image size

Table 1. The Mass (M_p ), Radius (R_p ), Incident Flux (f_⋆), Semimajor Axis (a), and Period (P) for the Target Planets in Our Study

Planet	Mp	Rp	f⋆	a	P	Additional References
	(M ⊕)	(R⊕)	(erg s−1 cm−2)	(au)	(days)
V1298 Tau b	${13.10}_{-5.30}^{+5.30}$	${9.95}_{-0.35}^{+0.37}$	${5.01}_{-0.52}^{+0.53}\times {10}^{7}$	0.172	24.14	J. H. Livingston et al. (2024, in preparation)
V1298 Tau d	${6.00}_{-0.70}^{+0.70}$	${6.34}_{-0.30}^{+0.30}$	${1.22}_{-0.13}^{+0.13}\times {10}^{8}$	0.110	12.40	⋯
V1298 Tau e	${15.30}_{-4.20}^{+4.20}$	${9.50}_{-0.49}^{+0.51}$	${2.07}_{-0.22}^{+0.22}\times {10}^{7}$	0.267	46.77	⋯
HATS-36 b	${1022.14}_{-19.71}^{+19.71}$	${13.84}_{-0.48}^{+0.48}$	${7.89}_{-0.88}^{+0.79}\times {10}^{8}$	0.054	4.18	⋯
HIP 67522 b	<20	${10.07}_{-0.47}^{+0.47}$	${4.28}_{-0.42}^{+0.41}\times {10}^{8}$	0.076	6.96	P. C. Thao et al. (2024)
TOI-837 b	${120.62}_{-19.41}^{+18.46}$	${8.63}_{-1.01}^{+1.01}$	${2.48}_{-0.97}^{+0.97}\times {10}^{8}$	0.083	8.32	O. Barragán et al. (2024)
TOI-1268 b	${96.40}_{-8.30}^{+8.20}$	${9.10}_{-0.60}^{+0.60}$	${1.63}_{-0.77}^{+0.78}\times {10}^{8}$	0.071	8.16	⋯
TOI-1431 b	${991.62}_{-57.21}^{+57.21}$	${16.70}_{-0.56}^{+0.56}$	${7.82}_{-1.09}^{+1.56}\times {10}^{9}$	0.046	2.65	⋯
TOI-2046 b	${731.01}_{-88.99}^{+88.99}$	${16.14}_{-1.23}^{+1.23}$	${3.90}_{-1.55}^{+1.55}\times {10}^{9}$	0.027	1.50	⋯
TOI-4087 b	${232.01}_{-44.50}^{+44.50}$	${13.05}_{-0.27}^{+0.28}$	${1.06}_{-0.06}^{+0.07}\times {10}^{9}$	0.045	3.18	⋯
TOI-2152 A b	${899.45}_{-117.60}^{+120.77}$	${14.36}_{-0.52}^{+0.56}$	${2.48}_{-0.40}^{+0.40}\times {10}^{9}$	0.051	3.38	⋯
TOI-201 b	${133.49}_{-9.53}^{+15.89}$	${11.30}_{-0.17}^{+0.14}$	${4.09}_{-0.83}^{+0.57}\times {10}^{7}$	0.300	52.98	⋯
TOI-622 b	${96.30}_{-22.88}^{+21.93}$	${9.24}_{-0.33}^{+0.31}$	${8.41}_{-1.41}^{+1.25}\times {10}^{8}$	0.071	6.40	⋯
KOI-1783.01	${71.00}_{-9.20}^{+11.20}$	${8.86}_{-0.24}^{+0.25}$	${8.76}_{-0.38}^{+0.38}\times {10}^{6}$	0.513	134.46	⋯
Kepler-51 b	${6.60}_{-2.70}^{+2.70}$	${6.89}_{-0.14}^{+0.14}$	${1.61}_{-0.07}^{+0.07}\times {10}^{7}$	0.251	45.15	K. Masuda et al. (2024) a
Kepler-51 c	${6.07}_{-0.42}^{+0.42}$	${8.98}_{-2.84}^{+2.84}$	${6.92}_{-0.28}^{+0.30}\times {10}^{6}$	0.384	85.31	⋯
KOI-351 g	⋯	${7.72}_{-0.21}^{+0.21}$	${5.21}_{-0.16}^{+0.16}\times {10}^{6}$	0.710	210.61	L. M. Weiss et al. (2024)††
Kepler- 76 b	${638.84}_{-111.24}^{+117.60}$	${13.11}_{-0.56}^{+0.59}$	${5.32}_{-0.79}^{+0.80}\times {10}^{9}$	0.027	1.54	T. A. Berger et al. (2018)††
Kepler-289 c	††† ${157.00}_{-7.00}^{+7.00}$	${11.21}_{-0.47}^{+0.50}$	${6.26}_{-0.81}^{+0.80}\times {10}^{6}$	0.510	125.85	T. A. Berger et al. (2018)††
Kepler-43 b	${998.00}_{-28.00}^{+27.70}$	${12.43}_{-0.54}^{+0.57}$	${1.42}_{-0.22}^{+0.21}\times {10}^{9}$	0.046	3.02	T. A. Berger et al. (2018)††
Kepler-74 b	${192.00}_{-30.00}^{+28.00}$	${12.82}_{-0.61}^{+0.65}$	${4.98}_{-0.26}^{+0.28}\times {10}^{8}$	0.078	7.34	T. A. Berger et al. (2018)††
Kepler-539 b	${308.30}_{-92.17}^{+92.17}$	${8.30}_{-0.41}^{+0.36}$	${4.95}_{-0.29}^{+0.28}\times {10}^{6}$	0.499	125.63	T. A. Berger et al. (2018)††

Notes. All planetary parameters are obtained from the NASA Exoplanet Archive (2023a, 2023b, 2024), with individual references specified in Table 2, with the exception of f_⋆ , which is calculated. Period values are rounded to 3 decimal places. The extra reference is for the planet mass, except for ^††, where it is instead for planet radius. ^†††Updated mass from M. Greklek-McKeon et al. (2023).

^a The values tabulated here are slightly different but within the 1σ error of what is reported in K. Masuda et al. (2024) because our analyses use the values that were communicated to us by K. Masuda and J. E. Libby-Roberts in private prior to their publication. Since the numbers agree well within 1σ, we keep our original values of Kepler-51 planets.

Download table as:  ASCIITypeset image

Table 2. The Stellar Properties, Including Metallicity [M/H], Age (t_*), Luminosity (L_*), Effective Temperature (T_*), Radius (R_*), and Mass (M_*), for All the Planet-hosting Stars in Our Study

Star	[M/H]	t*	L*	T*	R*	M*	References
	(dex)	(Myr)	(log10(L⊙))	(K)	(R⊙)	(M⊙)
V1298 Tau	${0.10}_{-0.15}^{+0.15}$	${25}_{-5}^{+5}$	${0.018}_{-0.043}^{+0.044}$	${5050}_{-100}^{+100}$	${1.35}_{-0.03}^{+0.03}$	${1.16}_{-0.06}^{+0.06}$	J. Sikora et al. (2023)
HATS-36	${0.28}_{-0.04}^{+0.04}$	${620}_{-550}^{+550}$	${0.215}_{-0.048}^{+0.043}$	${6149}_{-76}^{+76}$	${1.16}_{-0.04}^{+0.04}$	${1.22}_{-0.03}^{+0.03}$	D. Bayliss et al. (2018)
HIP 67522	0.00	${17}_{-2}^{+2}$	${0.243}_{-0.023}^{+0.022}$	${5675}_{-75}^{+75}$	${1.38}_{-0.06}^{+0.06}$	${1.22}_{-0.05}^{+0.05}$	A. C. Rizzuto et al. (2020)
TOI-837	$-{0.07}_{-0.04}^{+0.04}$	${35}_{-5}^{+11}$	${0.084}_{-0.164}^{+0.164}$	${6047}_{-162}^{+162}$	${1.02}_{-0.08}^{+0.08}$	${1.12}_{-0.06}^{+0.06}$	L. G. Bouma et al. (2020)
TOI-1268	${0.36}_{-0.06}^{+0.06}$	${245}_{-135}^{+135}$	$-{0.234}_{-0.193}^{+0.193}$	${5300}_{-100}^{+100}$	${0.92}_{-0.06}^{+0.06}$	${0.96}_{-0.04}^{+0.04}$	J. Šubjak et al. (2022)
TOI-1431	${0.09}_{-0.03}^{+0.03}$	${290}_{-190}^{+320}$	${1.068}_{-0.047}^{+0.078}$	${7690}_{-250}^{+400}$	${1.92}_{-0.07}^{+0.07}$	${1.90}_{-0.08}^{+0.10}$	B. C. Addison et al. (2021)
TOI-2046	$-{0.06}_{-0.15}^{+0.15}$	${450}_{-21}^{+430}$	${0.290}_{-0.092}^{+0.092}$	${6250}_{-140}^{+140}$	${1.21}_{-0.07}^{+0.07}$	${1.13}_{-0.19}^{+0.19}$	P. Kabath et al. (2022)
TOI-4087	${0.24}_{-0.08}^{+0.08}$	${800}_{-600}^{+1200}$	${0.176}_{-0.020}^{+0.025}$	${6060}_{-67}^{+74}$	${1.11}_{-0.02}^{+0.02}$	${1.18}_{-0.04}^{+0.04}$	S. W. Yee et al. (2023)
TOI-2152 A	${0.28}_{-0.08}^{+0.07}$	${830}_{-580}^{+1100}$	${0.653}_{-0.067}^{+0.069}$	${6630}_{-290}^{+300}$	${1.61}_{-0.05}^{+0.06}$	${1.52}_{-0.10}^{+0.09}$	J. E. Rodriguez et al. (2023)
TOI-201	${0.24}_{-0.04}^{+0.04}$	${870}_{-490}^{+460}$	${0.415}_{-0.017}^{+0.016}$	${6394}_{-75}^{+75}$	${1.32}_{-0.01}^{+0.01}$	${1.32}_{-0.03}^{+0.03}$	M. J. Hobson et al. (2021)
TOI-622	${0.09}_{-0.07}^{+0.07}$	${900}_{-200}^{+200}$	${0.474}_{-0.010}^{+0.010}$	${6400}_{-100}^{+100}$	${1.42}_{-0.05}^{+0.05}$	${1.31}_{-0.08}^{+0.08}$	A. Psaridi et al. (2023)
KOI-1783	${0.11}_{-0.04}^{+0.04}$	${200}_{-1300}^{+1300}$	0.213	${5922}_{-60}^{+60}$	${1.14}_{-0.03}^{+0.03}$	${1.08}_{-0.03}^{+0.04}$	S. Vissapragada et al. (2020)
Kepler-51	${0.05}_{-0.04}^{+0.04}$	${500}_{-250}^{+250}$	−0.142	${6018}_{-107}^{+107}$	${0.94}_{-0.50}^{+0.50}$	${1.04}_{-0.12}^{+0.12}$	J. E. Libby-Roberts et al. (2020)
KOI-351	$-{0.12}_{-0.18}^{+0.18}$	${530}_{-880}^{+880}$	0.269	${6080}_{-170}^{+260}$	${1.20}_{-0.10}^{+0.10}$	${1.20}_{-0.10}^{+0.10}$	J. Cabrera et al. (2014)
Kepler-76	$-{0.10}_{-0.20}^{+0.20}$	${620}_{-480}^{+480}$	${0.451}_{-0.022}^{+0.015}$	${6409}_{-95}^{+95}$	${1.32}_{-0.08}^{+0.08}$	${1.20}_{-0.20}^{+0.20}$	L. J. Esteves et al. (2015)
Kepler-289	${0.05}_{-0.04}^{+0.04}$	${650}_{-440}^{+440}$	${0.061}_{-0.023}^{+0.022}$	${5990}_{-38}^{+38}$	${1.00}_{-0.02}^{+0.02}$	${1.08}_{-0.02}^{+0.02}$	J. R. Schmitt et al. (2014)
Kepler-43	${0.33}_{-0.11}^{+0.11}$	${690}_{-540}^{+540}$	${0.326}_{-0.048}^{+0.043}$	${6041}_{-143}^{+143}$	${1.42}_{-0.07}^{+0.07}$	${1.32}_{-0.09}^{+0.09}$	A. S. Bonomo et al. (2017)
Kepler-74	${0.42}_{-0.11}^{+0.11}$	${800}_{-500}^{+900}$	0.332	${6000}_{-100}^{+100}$	${1.12}_{-0.04}^{+0.04}$	${1.18}_{-0.04}^{+0.04}$	A. S. Bonomo et al. (2017)
Kepler-539	$-{0.01}_{-0.07}^{+0.07}$	${910}_{-840}^{+840}$	$-{0.060}_{-0.023}^{+0.022}$	${5820}_{-80}^{+80}$	${0.95}_{-0.02}^{+0.02}$	${1.05}_{-0.04}^{+0.04}$	L. Mancini et al. (2016)

Note. All stellar parameters are obtained from the NASA Exoplanet Archive (2023a, 2023b, 2024). The metallicity for V1298 Tau was obtained from A. Suárez Mascareño et al. (2021) since it was not available from the data source given in the table. The quoted metallicity for HIP 67522 is an assumed solar value and not a measurement.

Download table as:  ASCIITypeset image

Using the RegularGridInterpolator function from scipy.interpolate to interpolate between the grid points, we obtain the range of plausible total mass of a planet given their measured radius, incident flux, and age. ⁶ For each of the 22 planets in our sample, we step through each of the 60 log-spaced metal mass grid points and use the COBYLA method of scipy.optimize.minimize to solve for the gas-to-metal mass ratio that produces a radius that matches the observed value within the tolerance of 10⁻⁵ R_⊕. Our models find two different solutions for each metal mass, one at low total mass (usually planetary) and another at high total mass (usually brown dwarf or stellar; see the right panel of Figure 2). This degeneracy breaks and always favors planetary mass solution for planets with mass measurements (see the left panel of Figure 2). Kepler-74 b is a special case where both families of solutions are within the planetary mass regime, so we consider both in our analysis.

We now describe how we rule out the high-mass family of solutions for planets with no mass measurements. For KOI-351 g we can rule out the high-mass solution since it lies beyond 100 M_J, well within the stellar mass regime, implying a density that is unphysically large. For HIP 67522 b, we check the high-mass solution against predicted mass–radius curves at different ages for brown dwarfs from I. Baraffe et al. (2003), as shown in Figure 3. Compared to the expected brown dwarf cooling curve, the high-mass solution for HIP 67522 b is way too dense for its estimated age, so we can safely rule such solutions out. This lack of high-mass solutions is consistent with the empirical paucity of short-period (≲200 days) brown dwarf companions to main-sequence stars (e.g., D. Grether & C. H. Lineweaver 2006; G. Duchêne & A. Kraus 2013; F. Kiefer et al. 2021).

Zoom In Zoom Out Reset image size

We now place further mass constraints on the remaining low-mass family of solutions. We first consider photoevaporative mass loss to place a lower limit on the mass. We integrate the following equation over the system age to determine the amount of mass lost:

$$\begin{eqnarray}&&\dot{M}=-\eta \displaystyle \frac{\pi {R}_{p}^{3}{f}_{{\rm{XUV}}}}{{{GKM}}_{p}}\end{eqnarray} \tag{ 4 }$$

where η is the efficiency parameter for which we adopt the empirical fit to numerical simulations computed in Appendix A of A. Caldiroli et al. (2022), f_XUV is the XUV flux on the planet from the host star, R_p is the planet radius, M_p is the planet mass, and K is the reduction factor which accounts for stellar tidal forces (N. V. Erkaev et al. 2007), given by

$$\begin{eqnarray}&&K=1-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{R}_{{\rm{Hill}}}}{{R}_{p}}\right)}^{-1}+\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{R}_{{\rm{Hill}}}}{{R}_{p}}\right)}^{-3},\end{eqnarray} \tag{ 5 }$$

where ${R}_{\mathrm{Hill}}=a{\left(\tfrac{{M}_{p}}{3{M}_{\star }}\right)}^{1/3}$ is the Hill radius, a is the orbital distance of the planet, and M_⋆ is the mass of the host star. While the form of Equation (4) resembles that of energy-limited approximation (see, e.g., N. V. Erkaev et al. 2007; J. Sanz-Forcada et al. 2011), η is a nontrivially varying parameter with respect to f_XUV and the gravitational potential of the planet that fully captures the non-energy-limited regime (see A. Caldiroli et al. 2022 for details).

2.4.1. Comparing and Choosing a Source for the XUV Flux Evolution

To calculate mass loss, we need the time evolution of the XUV luminosity of the planet's host star over the lifetime of the system. We consider two stellar evolution grids, which provide the bolometric luminosities throughout the star's lifetime (C. P. Johnstone et al. (2021, hereafter J21), which provides grids for stars ranging from 0.1 to 1.2 M_⊙, and the MESA Isochrones and Stellar Tracks (MIST; A. Dotter 2016; J. Choi et al. 2016), which provides grids for stars with masses between 0.1 and 300 M_⊙ and metallicities from –4 to 0.5 dex, where zero corresponds to solar metallicity. We ultimately choose MIST and describe why below.

First, we take the bolometric luminosity for a given stellar mass from the corresponding stellar evolution track. We note a slight discrepancy between the present-day measured values for the luminosity and the values from the stellar evolution grid at the age of the star. To account for this, we scale the grid luminosity to the measured value at system age. We then use the empirical relation to obtain the time-dependent X-ray luminosity:

$$\begin{eqnarray}{L}_{{\rm{X}}}=\left\{\begin{array}{cc}{L}_{{\rm{sat}}} & t\lt 100\,{\rm{Myr}}\\ {L}_{{\rm{sat}}}{\left(\displaystyle \frac{t}{100{\rm{Myr}}}\right)}^{-1.42} & t\geqslant 100\,{\rm{Myr}}\end{array}\right.,\end{eqnarray} \tag{ 6 }$$

where the value of the luminosity during the saturation interval is L_sat = 10^−3.6 L_*(t) (see, e.g., O. Vilhu & F. M. Walter 1987; N. J. Wright et al. 2011; G. W. King & P. J. Wheatley 2021), L_*(t) is the time-dependent bolometric luminosity over the lifetime of the star, and 100 Myr corresponds to a typical saturation timescale. The power-law index describing the time dependence of X-ray luminosity is taken to be –1.42, which corresponds to the median value obtained by L. Tu et al. (2015), who use rotational evolution models to predict how stellar X-ray luminosity decays with time. Typically, this value is found to be between –1.5 and –1.2 (see, e.g., I. Ribas et al. 2005; A. P. Jackson et al. 2012; M. W. Claire et al. 2012). The EUV flux evaluated at the stellar surface, F_EUV, is then obtained from the following ratio:

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{{\rm{EUV}}}}{{F}_{{\rm{X}}}}=\beta {{F}_{{\rm{X}}}}^{\gamma },\end{eqnarray} \tag{ 7 }$$

where F_X is the X-ray flux evaluated at the surface of the star, γ is the power-law index, and β is the scaling factor.

We compare the X-ray to EUV surface flux scaling equations from J21 to the equations presented in G. W. King & P. J. Wheatley (2021, hereafter KW21). Comparing the best-fit parameters from J21 and KW21, we find that the slopes and the general trends of the EUV flux as a function of X-ray flux are very different. The best-fit parameters for the ratio of hard (100–360 Å) and soft (360–920 Å) EUV to X-ray surface flux obtained by KW21 are the power-law indices ${\gamma }_{\mathrm{hard}}=-{0.35}_{-0.15}^{+0.07}$ and ${\gamma }_{\mathrm{soft}}=-{0.76}_{-0.04}^{+0.16}$, with corresponding β_hard = 176 and β_soft = 3040. ⁷ The best-fit parameters presented by J21 are γ_hard = −0.319 and γ_soft = −0.373 with β_hard = 110 and β_soft = 34.4. ⁸ KW21 generally find a steeper slope than J21, particularly in the soft EUV. KW21 also find that the hard EUV and soft EUV have a different slope (factor of ∼2 difference), while J21 find that the two have similar slopes. In an attempt to better understand the equations presented in each paper and to get a better understanding of which of the two is more reliable, we attempt to reproduce the results of J21 and KW21 using the source data.

For KW21, we use TIMED/SEE data (T. N. Woods et al. 2005) as described in KW21 to obtain the results shown in their Figure 2. ⁹ Using scipy.optimize.curvefit, we fit Equation (7) to the TIMED/SEE data and obtain γ_hard = −0.4039157 ± 0.0000006 and γ_soft = −0.7551290 ± 0.0000009, with corresponding β_hard = 270 ± 4 and β_soft = 4000 ± 1000. The left panel of Figure 4 shows the resulting TIMED/SEE data and fit we obtain compared to the results from KW21. We find that our γ_hard and γ_soft agree with theirs within their quoted error, while our β's differ by factors of order unity. We perform two independent mass-loss calculations using each set (ours versus KW21) of parameters and find that they produce almost the same final result. Therefore, given that our final results are not significantly affected by the small difference in the parameters, we perform all calculations using the corrected best-fit parameters from KW21: the power-law indices ${\gamma }_{\mathrm{hard}}=-{0.35}_{-0.15}^{+0.07}$ and ${\gamma }_{\mathrm{soft}}=-{0.76}_{-0.04}^{+0.16}$, and corresponding β_hard = 176 and β_soft = 3040.

Zoom In Zoom Out Reset image size

We perform a similar exercise on the J21 data. Their EUV fluxes are also obtained from solar data, using the Flare Irradiance Spectrum Model (FISM; P. C. Chamberlin et al. 2007). Since FISM is not publicly available, we use the updated version of this model, FISM2 (P. C. Chamberlin et al. 2020). FISM2 is an improved, more reliable version of the original FISM, so we consider this a valid comparison to verify whether the parameters obtained by J21 match up with solar data. In their paper, J21 find that their models agree with the solar data from FISM. However, we find that FISM2 does not agree with the scaling relations presented in J21, particularly in the soft EUVs, where the EUV surface flux is almost an order of magnitude larger at low X-ray surface fluxes. The slopes of the FISM2 data are steeper (factor of 2 larger) than what is given by J21. The results of this exercise are shown in the right panel of Figure 4. We see that for both hard and soft EUV flux, the equations presented in J21 do not agree with the FISM2 data. The discrepancy between the data and the fit is significantly more severe for J21 than for KW21.

We conclude, given the discrepancy between the J21 equations and the FISM2 data that we cannot fully explain that the J21 grid and scaling equations are not reliable enough for our purposes. Since our final mass ranges are robust to the difference in the fit parameters between our analysis and the results from KW21, we use the scaling equations from KW21 for all mass-loss calculations. We also find that we are limited by the range of stellar masses for which stellar evolution tracks are available from J21 since some of our stars have masses larger than 1.2 M_⊙. Therefore, we use MIST instead of the J21 stellar evolution tracks. We perform all calculations using the time-varying bolometric luminosities from MIST (equivalent evolutionary points tracks for 0.4× the critical spin, rounded to the closest stellar mass and metallicity values for which a grid is available—where the stellar metallicity measurement is not available, we use solar metallicity), scaled to match the observed luminosities where necessary, Equation (6) to obtain the X-ray luminosity, and Equation (7) with the best-fit parameters from KW21 to obtain the EUV surface fluxes. We then calculate the incident X-ray and EUV fluxes on the planet from the star and add them to obtain the incident high-energy flux, f_XUV.

2.4.2. Calculating Mass Loss

For each metal mass, we solve for the gas-to-metal mass ratio given the radius, incident flux, and age of the planet, as described in Section 2.2. We consider this gas-to-metal mass ratio as the initial ratio and numerically integrate Equation (4) using the trapezoid method, going forward in time starting at 5 Myr up to the age of the system. We take the time steps from the MIST grid, which gives 454 log-spaced times from 0.032 yr to 1.806 Gyr, and slice the array to correspond to times 5 Myr < t < t_*. We have experimented with twice as coarse and twice as fine time resolutions and found our final result to converge at the adopted MIST grid time step. For the initial condition, we recalculate the radius of the planet given this gas-to-metal mass ratio at 5 Myr. At each iteration, we update the radius, mass, and incident flux, then the efficiency η and the reduction factor K. The radius is updated using the radius grid as described in Section 2.2. By doing so, we implicitly neglect the adiabatic expansion of the envelope, which can underestimate the amount of mass loss (e.g., D. P. Thorngren et al. 2023).

If we find that the entire envelope is lost, we end the calculation and consider the solution unstable to mass loss. Otherwise, we iterate until t = t_*, the stellar age, and check if the final radius after mass loss is within the error bounds of the observed radius. If it is, we consider the solution stable against mass loss, and otherwise, we rule it out. The top panel of Figure 5 shows the result of this test for V1298 Tau b, where all solutions with a final radius below the lower limit set by the 1σ uncertainty on the measured radius are ruled out. We perform this calculation for each of the 60 metal mass and gas-to-metal mass ratio solutions obtained. This allows us to narrow down the mass range, constraining the lower limit since solutions with small metal masses are particularly susceptible to atmospheric escape due to weak surface gravity.

Zoom In Zoom Out Reset image size

2.4.3. "Backward" Mass-loss and Formation Constraints

Finally, we check our model-inferred metals-to-gas ratio against empirical mass–metallicity trends from D. P. Thorngren et al. (2016), which presents a study of the bulk compositions of transiting giant planets with masses of 20 M_⊕ < M_p < 20 M_J. We compare our results to the relationship they obtain between the planet's heavy element mass M_met and total mass M_p , given by

$$\begin{eqnarray}&&{M}_{{\rm{met}}}=(57.9\pm 7.03){M}_{p}^{(0.61\pm 0.08)}.\end{eqnarray} \tag{ 8 }$$

They further analyzed the correlation between the planet's bulk heavy element enrichment and the host star's metallicity, finding

$$\begin{eqnarray}&&\displaystyle \frac{{Z}_{p}}{{Z}_{\star }}=(9.7\pm 1.28){M}_{p}^{(-0.45\pm 0.09)},\end{eqnarray} \tag{ 9 }$$

where Z_p = M_met/M_p is the planet metallicity and Z_⋆ = 0.0134 × 10^[M/H] is the stellar metallicity. For each metal mass and gas-to-metal mass ratio solution, we verify whether the planet lies within the 1σ and 2σ region from the best fit. The intrinsic spread is given by 10σ = 1.82 ± 0.09, so our 1σ region corresponds to the area within a factor of 1.82 from the best-fit line.

We consider these constraints only for solutions with a total mass of M_p > 20 M_⊕ since D. P. Thorngren et al. (2016) studied giant planets, and their results do not extend to lower mass planets. An example is shown in Figure 6, which illustrates how the check on planet and stellar metallicity places an upper limit on the plausible planet and heavy element mass.

Zoom In Zoom Out Reset image size

Our default analysis produces final mass estimates for 19 of the 20 planets with observed masses. Figure 7 shows our final mass ranges given all the constraints described in Section 2, where we also show the results of extending the metallicity constraint out to 2σ. When comparing the planetary and stellar metallicities to trends from D. P. Thorngren et al. (2016), as described in Section 2.5, we find that for only five planets, the measured mass agrees with our metallicity constraint within the 1σ range. While the 1σ limit is a good reference, it is not unreasonable to find planets beyond that. Figures 7 and 11 of D. P. Thorngren et al. (2016) show the best-fit lines for the planetary and stellar metallicity trends, respectively, compared with the planets used to obtain the trends. Although most of the planets used to obtain the trends lie within the 1σ region, a significant number have metallicities outside the 1σ region. In each case, ∼40%–50% of the planets studied by D. P. Thorngren et al. (2016) lie outside the 1σ region. Therefore, we can safely extend our metallicity constraint out to 2σ without invalidating our results, within which we find six planets with solutions consistent with their measured masses.

Zoom In Zoom Out Reset image size

We now discuss planets that fail even the 2σ metallicity constraint. Kepler-76 b does not have any masses that survive our constraints, either by mass-loss or by metallicity constraints; for example, the mass solutions within 1σ agreement with the measured mass lie beyond the 2σ metallicity constraints. Similarly, for Kepler-289 c and TOI-201 b, all the solutions that agree with the measured mass fall outside the 2σ range of the metallicity trends and are therefore ruled out (unlike Kepler-76 b, these planets have solutions that survive our mass-loss constraints). In the case of HATS-36 b, a subset of the solutions that are consistent with its measured mass lie just at the edge of the 2σ metallicity constraint. We note that while the results from D. P. Thorngren et al. (2016) provide a good foundation for general metallicity trends obeyed by short-period gas giants, we cannot entirely disregard the possibility of outliers.

The planets in the Kepler-51 and V1298 Tau systems are too puffy and require special formation constraints, where the planets likely accreted dust-free envelope farther out and migrated into their present-day orbits (in contrast to our fiducial in situ dusty accretion). The Kepler-51 and V1298 Tau planets are discussed in more detail in Sections 4.1 and 4.2, respectively.

For Kepler-43 b, the bounds of our interior structure grid (Section 2.2) do not allow for any mass solutions in agreement with the measured mass. Therefore, for this planet, we extend our analysis to include solutions that fall within the 1σ bounds of the measured radius to obtain a final mass range that encompasses the measured mass. Kepler-74 b has two sets of solutions since both the high-mass and low-mass family of solutions fall below the brown dwarf regime.

Having shown that our methods can recover a mass estimate for 19 of the 20 planets with mass measurements, given their measured radius, age, and incident flux, we now move on to estimating the plausible mass range of the two planets without mass measurements: KOI-351 g and HIP 67522 b. The final mass range obtained from our methods is shown in Figure 7, where the error bars without markers show the planets without mass measurements. We find that our method produces a solution that survives all constraints for both planets, straddling the boundary between low-mass and giant planets. However, HIP 67522 b has an upper limit of 20 M_⊕ (P. C. Thao et al. 2024), which agrees with the lower end of our mass range, making it a puffed-up Neptune mass planet. We discuss this planet in more detail in Section 4.

Kepler-51 is a well-known system harboring multiple super-puffs (K. Masuda 2014; J. E. Libby-Roberts et al. 2020). While the revised masses of member planets generally increase their bulk density (K. Masuda et al. 2024), they remain puffier than expected for their size and age (see Figure 7). E. J. Lee & E. Chiang (2016) hypothesized that the super-puffs likely formed beyond the water ice line where the rate of gas accretion under dust-free accretion greatly accelerates so that low-mass cores can reach the gas-to-core mass ratio (GCR) well beyond 10%. These planets can also appear spuriously large due to high-altitude hazes (L. Wang & F. Dai 2019; P. Gao & X. Zhang 2020), planetary rings (A. L. Piro & S. Vissapragada 2020), or tidal heating (S. Millholland 2019). Of these alternate scenarios, tidal heating is likely not the main cause of puffing up Kepler-51 planets as they are too far from the star. Neither high-altitude hazes nor planetary rings are mutually exclusive with initially high GCR (albeit lowering the initial GCR by some amount), so we proceed with the hypothesis of dust-free gas accretion followed by inward migration.

First, over a range of initial GCR and metal mass, we deduce the orbital period at which a planet of given metal mass would have been to reach a given GCR (since the metal mass relevant for Kepler-51 planets is always less than 20 M_⊕, their metal mass is equivalent to core mass and their gas-to-metal mass ratio is identical to GCR). We adopt the scaling relationship for dust-free accretion reported in (E. J. Lee et al. 2022, their Equation (7)),

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{GCR} & = & 0.27{\left(\displaystyle \frac{t}{0.1\,\mathrm{Myr}}\right)}^{0.4}{\left(\displaystyle \frac{{M}_{\mathrm{metals}}}{20\,{M}_{\oplus }}\right)}^{1.7}\\ & & \times {\left(\displaystyle \frac{{T}_{d}}{1000\,{\rm{K}}}\right)}^{-1.5}\mathrm{Exp}\left(\displaystyle \frac{t}{2.2\,\mathrm{Myr}}{\left(\displaystyle \frac{{M}_{\mathrm{metals}}}{20\,{M}_{\oplus }}\right)}^{4.2}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

where the normalization is taken from visual inspection of Figure 3 of E. J. Lee et al. (2022), t is the accretion time varied over 1–10 Myr, the exponential term mimics runaway accretion (E. J. Lee 2019, see also S. Ginzburg & E. Chiang 2019), and T_d = 1000 K(a/0.1 au)^−3/7 is the disk midplane temperature for an irradiated disk (E. I. Chiang & P. Goldreich 1997). The dependence on disk gas surface density is insignificantly weak (see also E. J. Lee & E. Chiang 2016), so we ignore it. We note that the steep dependence of GCR on T_d arises from the temperature dependence of gaseous opacity, which is benchmarked to R. S. Freedman et al. (2014) and E. J. Lee et al. (2018), the latter of which is based on the calculations outlined in J. W. Ferguson et al. (2005). Our chosen opacities are apt to test the hypothesis of dust-free gas accretion; however, any opacity that drops at colder temperatures—generally expected for gaseous opacities from the freeze-out of internal degrees of freedom—would produce similar results.

Next, over the same grid of initial GCR and metal mass, we compute atmospheric mass loss as described in Section 2.4.2 at the present location of each planet over the system age. We evaluate the final mass and radius of the planet using our interior structure model (Section 2.2) and compare them to the measured mass and radius. Figure 8 summarizes our calculation, where we illustrate in colors the degree of agreement between our model and observations in units of 1σ measurement error. For ease of visualization, we sum the mass and radius deviation; i.e., each color represents |∣measured mass − model mass|/σ_mass + |measured radius −model radius|/σ_radius. Highlighted in thin white lines are the GCR-metal mass pairs that are consistent with both the present-day mass and radius within 1σ after mass loss. We find a strictly negative correlation between the plausible initial GCR and metal mass for Kepler-51 c, as expected for a planet that undergoes a minimal level of mass loss (i.e., to explain the same total mass and radius, we need lower GCR for higher metal mass). Kepler-51 b is closer to the star, so the higher irradiation flux contributes to extra heating, puffing up the planet, especially for lower metal mass, so we require slightly lower initial (and final; see the right column of Figure 8) GCR for lower metal mass down to ∼3.5 M_⊕.

Zoom In Zoom Out Reset image size

Among the GCR-metal mass pairs that lead to final radii and masses that agree with the measured values within 1σ error, we locate the most likely formation location for each planet assuming that the two planets are within 5% of 2:1 orbital period ratio, their present-day near-resonant configuration. We find four sets of solutions, highlighted in colored boxes in Figure 8. Specifically, these sets of solutions are:

• 1.  
  Kepler-51 b

  • (a)  
    GCR_init = [0.64, 0.68, 0.73]
  • (b)  
    GCR_fin = [0.52, 0.51, 0.44]
  • (c)  
    M_metal = [4.39, 3.76, 2.92] M_⊕
  • (d)  
    a_init = [2.53, 4.22, 9.12] au
• 2.  
  Kepler-51 c

  • (a)  
    GCR_init = [0.65, 0.83, 1.08]
  • (b)  
    GCR_fin = [0.56, 0.72, 0.91]
  • (c)  
    M_metal = [3.74, 3.58, 3.11] M_⊕
  • (d)  
    a_init = [3.96, 6.61, 14.3] au

where GCR_init is the initial GCR, GCR_fin is the final GCR, and a_init is the formation location. While Kepler-51 d is not in our sample because it is too far from the host star, it is currently near 3:2 resonance with planet c, so its corresponding a_init = [5.19, 8.66, 18.8] au.

All of our solutions indicate formation beyond the typical water ice line of ∼1 au and therefore subsequent large-scale migration to bring both planets to their present locations at ∼0.25 au (planet b) and ∼0.38 au (planet c), consistent with the hypothesis of E. J. Lee & E. Chiang (2016). We note that our calculation is limited to two-layer planets with an inner metallic core and H/He envelope on top with atmospheric metallicity that of solar. One consequence of formation beyond the water ice line is the deposition of water during and after formation. Our solutions for GCR_init and a_init will likely change if we allow for a significant amount of water deposition (P. Gao et al. 2024, in preparation; E. Vlahos et al. 2024).

The estimated age of Kepler-51 is 500 ± 250 Myr. Over the full 1σ error in the age, we find that while the exact resonant pair solutions for Kepler-51 b and c can vary, their a_init are always beyond 1 au with minimal change to the initial and final GCR (see the Appendix), demonstrating the robustness of our qualitative result.

We apply the same analysis as in Section 4.1 for V1298 Tau planets in our sample. We limit the gas accretion time to 1 Myr because any longer duration leads to the formation location of planet d being inside its present location, which we find to be unlikely. Physically, shorter gas accretion time implies a later assembly of planetary cores.

Similar to Kepler-51 b (and unlike Kepler-51 c), a positive GCR-metal mass scaling that provides a set of solutions that agree with measured mass and radius within 1σ error is observed for V1298 Tau d (see Figure 9). Owing to the extra puffing of planets by stellar irradiation (which drives up the rate of mass loss), the required initial GCR is low ∼0.09. Even for planets b and e, which are slightly farther out from the host star, we see a generally positive correlation between the initial GCR and metal mass that matches the measured radius and mass. V1298 Tau is a much younger system compared to Kepler-51, and both planets b and e have shorter orbital periods as compared to Kepler-51 b. Consequently, all V1298 Tau planets are subject to higher irradiation flux, increasing the extra heating that enlarges the planets significantly toward lower metal mass, requiring lower GCR to explain the same radius (see dashed lines in Figure 2).

It has been argued that planets b and e are near 2:1 resonance (A. D. Feinstein et al. 2022), although not necessarily formally in resonance (D. Turrini et al. 2023). While the quoted orbital periods of planet d suggest it may be close to 2:1 resonance with planet b, whether this planetary system is in a resonant chain or not remains to be confirmed. Our inferred initial formation location of planet d is close to its present-day location, both of which are within ∼20 days and close to the typical spin periods of pre-main-sequence stars (e.g., J. Irwin et al. 2008). We surmise that planet d formed close to the inner edge of the disk set by the corotation radius with respect to the host stellar spin (P. Ghosh & F. K. Lamb 1979; A. Koenigl 1991), whereas the outer planets formed much farther away and underwent a large-scale inward migration. Following the same analysis as Section 4.1, the sets of solutions for planets b and e under the assumption that their initial locations were close to 2:1 resonance are:

• 1.  
  V1298 Tau b

  • (a)  
    GCR_init = [0.72, 0.72, 0.92]
  • (b)  
    GCR_fin = [0.71, 0.71, 0.91]
  • (c)  
    M_metal = [6.84, 6.84, 8.11] M_⊕
  • (d)  
    a_init = [1.90, 1.90, 1.76] au
• 2.  
  V1298 Tau e

  • (a)  
    GCR_init = [0.81, 1.04, 0.92]
  • (b)  
    GCR_fin = [0.81, 1.03, 0.91]
  • (c)  
    M_metal = [6.21, 7.16, 6.84] M_⊕
  • (d)  
    a_init = [2.97, 2.97, 2.78] au

The level of mass loss is more muted for V1298 Tau planets compared to Kepler-51 planets (in spite of higher irradiation flux) because the former are more massive. Nevertheless, like Kepler-51, we find that a similar formation scenario—formation beyond ∼1 au followed by large-scale inward migration—can explain the current mass and radii of V1298 Tau b and e. At face value, such a scenario may be in tension with the reported substellar oxygen abundance for planet b (S. Barat et al. 2024). However, oxygen-rich solids can be locked into the central core and/or deep interior without being dredged up all the way to high atmospheric altitudes where we can probe the atmospheric composition. Furthermore, rapid inward migration would have prevented the planet from attaining a water-enriched envelope (E. Vlahos et al. 2024). Alternatively, in situ gas accretion for planet b is possible if the gas metallicity is as low as 0.01× solar value, which is more than an order of magnitude below the value reported by S. Barat et al. (2024), so we do not favor this scenario. By contrast, planet d is more robustly consistent with practically in situ gas accretion. Similar to Kepler-51 planets, these results—in situ gas accretion of planet d and large-scale migration of planets b and e—are robust to the uncertainties in the estimated age within 1σ error (see the Appendix).

At the measured radius and the upper limit on mass from atmospheric constraint ≲20 M_⊕ (P. C. Thao et al. 2024), HIP 67522 b has a bulk density that is below 0.1 g cm⁻³, the lower limit below which very few planets are seen in nature (see the bottom panel of Figure 1 of D. P. Thorngren et al. 2023). At its high irradiation flux (i.e., short orbital period), such low density is expected to spell runaway mass loss through an adiabatic expansion of the envelope during its photoevaporative loss.

Given the extremely young age of HIP 67522 b (t_* ∼ 17 Myr), we may be observing the planet in a semi-runaway state, having not yet lost its entire envelope. Figure 3 of D. P. Thorngren et al. (2023) illustrates planets at various ages undergoing mass loss and thermal evolution and shows that at 10 Myr, small planets (∼20 M_⊕, ∼11 R_⊕) have begun to lose mass but still have most of their envelopes, consistent with the existence and the measured properties of HIP 67522 b. This planet may be an ideal target to study the ongoing mass-loss process.

Recent radial velocity measurements from the High Accuracy Radial Velocity Planet Searcher (M. Mayor et al. 2003) find that TOI-837 b has a mass of ∼120 M_⊕ (O. Barragán et al. 2024), which is in agreement with our proposed mass range with the 2σ metallicity constraint. This mass estimate would imply that TOI-837 b is a young (∼35 Myr), hot Jupiter that migrated to its close-in orbit (∼8.32 days) over a short timescale, potentially through disk-induced migration. Given the rarity of hot Jupiters around young stars (≲100 Myr), however, we argue that disk-induced migration is likely not the dominant channel of hot Jupiter formation.

In general, we find that planets larger than 10 R_⊕ are much more likely to have masses ≳100 M_⊕, consistent with being gas giants. In terms of system age, among the particularly young planets (≲40 Myr), with the exception of TOI-837 b, all the planets in our sample (HIP 67522 and V1298 Tau) have either measured or upper limit on masses that place them well within the Neptune to sub-Neptune regime.

We additionally analyze K2-33 b (5.04${}_{-0.37}^{+0.34}{R}_{\oplus }$, 5.425 days around a 1.02 M_⊙ star aged 9 Myr; A. W. Mann et al. 2016), DS Tuc A b (5.7 ± 0.17 R_⊕, 8.138 days around a 1.01 M_⊙ star aged 45 Myr; E. R. Newton et al. 2019, metallicity from S. Benatti et al. 2019), and TOI-1227 b (9.572${}_{-0.583}^{+0.751}{R}_{\oplus }$, 27.364 days around a 0.17 M_⊙ star aged 11 Myr; A. W. Mann et al. 2022) that were excluded by our target selection (K2-33 b and DS Tuc A b are smaller than 6R_⊕ while TOI-1227 b is around a very low-mass star so its insolation flux was too low to be considered even a warm Jupiter). We find that the range of masses consistent with the 1σ metallicity constraint are 12–30.60 M_⊕, 14.38–18.40 M_⊕, and 38.26–222.65 M_⊕ for K2-33 b, DS Tuc A b, and TOI-1227 b, respectively. Allowing for a 2σ metallicity constraint increases the upper limit to 148.77 M_⊕, 73.29 M_⊕, and 486.37 M_⊕. These results are largely consistent with our main result showcased in Figure 7. The smaller planets K2-33 b and DS Tuc A b are more consistent with being puffed-up Neptunes, whereas the larger TOI-1227 b has a higher likelihood of being a gas giant. We note, however, that in terms of insolation flux, TOI-1227 b is considered a cold Jupiter and the disk temperature (Y. Chachan & E. J. Lee 2023, their Equation (18)) at its current location of ∼200 K is cold enough to allow for rapid gas accretion in situ (GCR > 1 over 0.1 Myr; see Equation (10)). Given that TOI-1227 b likely does not require either disk or high-eccentricity migration to explain its property, we do not consider it as part of our sample of young hot/warm Jupiter.

We further consider the possibility of young hot and warm Jupiters being inflated to a size large enough to be misclassified as an eclipsing stellar binary, and we consider this unlikely. Following the same procedure as D. P. Thorngren et al. (2023), we find that young giants would cool down to ≲1.5 Jupiter radius within ∼20–40 Myr even with extra heating and so they would have been classified as a giant planet.

In the sample we study (plus K2-33 b and DS Tuc A b), we see no planet larger than 10R_⊕ of age younger than ∼100 Myr. This paucity of very young, large, and massive planets provides tentative evidence that the dominant origin channel for hot Jupiters is likely high-eccentricity migration. If disk gas migration were dominant, we would find Jupiter-mass planets around stars as young as ∼10 Myr, but we instead find that the youngest hot/warm Jupiter candidates are inflated lower mass planets. The one exception of TOI-837 b may be a rare case of a hot Jupiter that arrived at its present location by disk migration. Our result is consistent with obliquity measurements from C. Spalding & J. N. Winn (2022), who find evidence that most hot Jupiters arrive late (≳100 Myr) based on tidal theory and stellar evolution models.

We argue that our finding is statistically significant against the general rarity of hot and warm Jupiters seen in the mature system, which reports a cumulative occurrence rate of giants at ∼4%–5% out to ∼200 days (e.g., F. Fressin et al. 2013; E. A. Petigura et al. 2018). We find ∼66 systems younger than 1 Gyr across K2, Kepler, and TESS detections at the time of writing this paper. If the hot/warm Jupiter occurrence rate is the same as that of main-sequence stars, we expect to find approximately three giants. We find more than three (albeit without any correction for detection bias/sensitivity, transit probability, etc.), with a clear trend toward lack of giants at ages ≲100 Myr.