The Radius Distribution of M dwarf-hosted Planets and its Evolution

Stars in binary or multiple systems with separations $\lesssim$50 au are known to host a planet population different from that around single stars (Kraus et al., 2016). Furthermore, binary stars with separations ($<$100 au) have different rotational histories than their single counterparts (e.g., Messina et al., 2017; Stauffer et al., 2018), meaning that a rotation-age relation developed for single stars cannot be applied. 90% of M dwarfs in the Kepler field are more distant than 135 pc, thus these physical separations correspond to angular separations $\lesssim$1". Companions or any nearby sources unresolved or partially resolved by Kepler (which had 4" pixels and a 3.1-7.5" PSF) can also dilute transit signals and make identification of the host star itself – and thus planet properties – ambiguous. The $K_{s}$-band photometry used to estimate radius and mass will be inaccurate for binaries that are unresolved by the 2MASS survey (resolution of $\sim$3"). Thus we identified and excluded KOIs around binary stars with separations $<4$". Binary/multiple systems or candidate systems were identified and excluded using (1) Gaia astrometry of resolved sources; (2) Gaia astrometric error relative to single-star astrometric solutions; (3) adaptive optics imaging; and (4) over-luminosity relative to a single-star locus in a color-absolute magnitude diagram.

Gaia resolved companions: We identified companions resolved by the Gaia satellite ($\gtrsim$1 arcsec, Ziegler et al., 2018a) based on similarity of parallax (distance) and proper motion. Rather than fixed criteria for similarity (e.g., El-Badry et al., 2021), we adopted a Bayesian approach. We identified candidate companions with cone searches of the DR3 catalog with a radius of 74" (corresponding to about $10^{4}$ au at the median distance of Kepler M dwarfs) around each Kepler M dwarf. We computed the probability of companionship as $p_{\rm binary}=P_{c}/(P_{c}+P_{b})$, where $P_{c}$ and $P_{b}$ are the likelihoods that the source is a companion or background star, respectively. We evaluated these assuming normally-distributed errors using:

and

where $\Delta\mu$ and $\sigma_{\Delta\mu}$ are the absolute scalar difference in proper motion between the star and source and its standard error, $\Delta\pi$ and $\sigma_{\Delta\pi}$ are the difference in parallax and its standard error, and $\mu_{\rm orb}$ is a “softening" term that accounts for the angular orbital motion of a binary. This last term we set to 5 mas yr^-1 based on Eqn. 4 in El-Badry & Rix (2018). Constant factors that divide out are omitted from Eqns. 1 and 2. The difference in the terms involving parallax and proper motion reflects the latter’s two-dimensionality. To improve the statistical accuracy by averaging over many stars, the summation for $P_{b}$ is performed over the $S$ stars in all the $N-1$ other cone searches (a few thousand stars) that cannot be companions. We imposed no weighting with angular separation other than the generous search radius. Based on the observed distribution of $P_{\rm binary}$ values relative to those calculated using stars from other fields, we identified a threshold of $P_{\rm binary}>0.98$ to select companions.

Gaia astrometric error: Binaries with separations $<1$" cannot be resolved by Gaia (Ziegler et al., 2018b) but can often be identified based on the poor fit of Gaia astrometry to a single-star solution (Belokurov et al., 2020). This is quantified as the Reduced Unit Weight Error (RUWE, Lindegren et al., 2018) and values significantly greater than unity are indicative of a multiple star system. We selected a threshold of 1.4 based on Lindegren et al. (2021).

AO imaging: AO imaging of KOIs was obtained of all but 8 systems with the NIRC2 infrared imager and AO system on the Keck-2 telescope (McLean & Chaffee, 2000). Observations used the narrow camera (10 mas pixel scale) and either the $K_{cont}$ or $K_{p}$ filters. Images were processed, sources extracted, and astrometry performed using the pipeline described in Kraus et al. (2016). We selected those sources with separations $<$2" and contrast $\Delta K<4$ as highly likely to be bound companions based on Kraus et al. (2016).

Photometry: Unresolved late K and M dwarf binaries can be identified by their over-luminosity relative to single stars on the main sequence since these stars evolve imperceptibly on the main sequence. Metallicity is the sole confounding factor, and differences in [Fe/H] produce a dispersion in color-absolute magnitude relations.³³3Gravity does not vary significantly for a given $T_{\rm eff}$ among late K- and M dwarf stars and extinction is negligible at the distance of most M dwarfs observed by Kepler. By experimenting with different color-magnitude combinations constructed from ground-based photometry, we found that a color-magnitude locus of Kepler M dwarfs constructed using PanSTARRS $g$-$y$ vs. $M_{y}$ (Tonry et al., 2012) has the smallest dispersion and is thus the least sensitive to metallicity. We identified stars that are more luminous than the best fit by more than two standard deviations as probable binaries.

The posterior distribution in radius of each planet was determined by multiplying each value from the Monte Carlo Markov Chain (MCMC) posterior chains for planet-to-star radius ratio from DR25 Q1-17 (Hoffman & Rowe, 2017) with a value from the stellar posterior distribution in radius calculated above, assuming Gaussian errors in magnitude and parallax. To correct for detection efficiency, we normalized each distribution by the summed probability for detection of an equivalent planet around every star observed by Kepler that is a late K- or early M-type dwarf $M_{K}>4.5$, 3200 $<$ $T_{\rm eff}$ $<$ 4200K) and that appears in the Santos et al. (2019) catalog of stars evaluated for rotation periods. We used the KeplerPORTs code (Burke & Catanzarite, 2017b) to calculate the probability of detection (including the geometric transit probability). The required window functions, one-sigma detection functions (Burke & Catanzarite, 2017a), “dataspans" (timespans of data), and duty cycles were retrieved from the NASA Exoplanet Archive. Non-linear (four-coefficient) limb-darkening coefficients were obtained from the DR25 Stellar Properties Catalog (Mathur et al., 2017). To mimic the removal of binaries from the KOI list, we excluded all target stars with Gaia RUWE values $>$1.4. Finally, KOIs were then matched with entries in the TIME Table of Gaidos et al. (2023) to retrieve gyrochronologic age estimates.

The summed probabilities over all stars represent the expected number of detections if each star in the observed sample had such a planet. Each posterior distribution in radius was divided by this quantity and then summed to give the number of objects per unit radius per star. As a metric of uncertainty due to radius error and finite counting statistics, we calculated 68% confidence intervals by constructing 1000 Monte Carlo realizations using bootstrap sampling of the planet population with replacement.⁴⁴4While multi-planet systems will affect the counting statistics of total planet occurrence (which is not considered here), correlation between the radii of planets in the same system (e.g., Weiss et al., 2018) would impact the significance of individual features.

To validate our methodology for reconstructing the planet radius distribution we constructed a sample of synthetic KOIs the same size as the actual one (115 planets, see Sec. 3), drawing radii from the sum of two Gaussian distributions representing super-Earths and sub-Neptunes. The fractional populations in the two distributions were set to 0.4 and 0.6, the centroids were set to 1.3 and 2.1 $R_{\oplus}$, and the standard deviations were set to 0.15 and 0.25$R_{\oplus}$, respectively. These values produces a reasonable radius distribution that roughly approximates the observed distribution but because this a test of the fidelity of recovery, the precise values do not matter. Orbital periods were selected from a power-law distribution between 1 and 40 days that reproduces the observed distribution (Gaidos et al., 2016). Synthetic planets were randomly assigned to the Kepler late K/early M population with measured rotation periods as defined above. To calculate the planet-star radius ratio, the stellar radii were determined using the metallicity-dependent relation of Mann et al. (2015), with metallicities drawn from the Gaussian distribution found among Kepler M dwarfs Gaidos et al. (2016). However, for added realism these artificial metallicities were not propagated to the analysis of synthetic KOIs since the metallicities of many KOIs are not well determined; instead, the metallicity-independent relation of Mann et al. (2015) was used. The total probability of detection (including the geometric factor) was calculated for each star-planet pair as described above and the pair accepted conditional to a random uniform draw below this probability. In lieu of actual posterior distributions for the observed radius ratio we adopted Gaussians with fractional standard deviation that are 1/SNR, where SNR is the predicted transit signal-to-noise (Howard et al., 2012):

where $\delta$ is the transit depth (approximately the radius ratio squared), $\sigma_{\rm CDPP}$ is the measured combined differential photometric precision of the star over 2 hours (Christiansen et al., 2012, , 2 hr is a typical transit duration), $N$ is the number of transits (equal to the dataspan times the duty cycle divided by the orbital period), and $\tau$ the transit duration (after averaging over a uniformly distributed impact parameter). Host star-planet pairs were generated until the accepted number equaled the size of the actual sample We ran the synthetic KOIs through the same analysis pipeline as the actual KOIs to produced a simulated observed radius distribution. This is plotted along with the input distribution in Fig. 1. The bimodal distribution is largely reproduced, but with a noticeable broadening of the super-Earth distribution, probably because of the absence of measured [Fe/H] "measurements" and the effect of noise on the radii of the smallest detectable planets. The input ratio of super-Earths to Neptunes (divided by the vertical grey line) is 0.74 while the retrieved ratio is 0.78, an increase of 5%. The increase is due to the greater scatter of the more poorly determined super-Earths into the sub-Neptune range compared to vice-versa. The location of the apparent valley is shifted to a slightly larger radius (from 1.64 to 1.75 $R_{\oplus}$) for the same reason.

We identified 117 confirmed/validated or candidate planets around 74 single late K and early M dwarfs with matches to Gaia DR3 sources and rotation periods in Santos et al. (2019). Thirteen have NASA Exoplanet Archive dispositions of "candidate" while the rest are "confirmed" (i.e., validated). Orbital periods, radii and ages of KOIs are provided in Table 1. Of these, 114 (in 73 systems) have rotation-based ages: nearly all are from Gaidos et al. (2023) but two missing systems (KOIs 2542 and 2705) were assigned ages using the methods of Gaidos et al. (2023). Four planets in three other systems (KOIs 430, 4957, and 6863) have rotation periods in Santos et al. (2019) that are shorter than the calibrated range in Gaidos et al. (2023) and are not age-dated.

The two largest objects in the sample, excluded from further analysis, are the hot Jupiter KOI-254.01 (Kepler-45b, Johnson et al., 2012), and the $5.45\pm 0.27$ $R_{\oplus}$ KOI-2715.01 (Kepler-1321b, Morton et al., 2016) orbiting a late K dwarf (neither shown in Fig. 5). The next largest planet is the $3.20\pm 0.14$$R_{\oplus}$  sub-Neptune KOI-253.01 (Kepler-2000b, Valizadegan et al., 2023), highlighting how Neptune-size and larger planets are very rare around M dwarfs. Figure 2 plots the radius vs. orbital period of the 115 small planets, along with a Gaussian kernel density estimator (KDE) to highlight structure. Planets with assigned ages are plotted as black points; while those lacking ages are plotted in grey.

We calculated the detection efficiency (i.e. fraction of detected planets with a given radius and period if distributed uniformly among the late K/early M type stars), in among 3200-4200K stars in Santos et al. (2019) using the procedures described in Sec. 2. This is represented as the grey contours in Fig. 2. As expected, detection efficiency is lowest for smaller planets on wider orbits, which explains the absence of detections in the lower right corner of the plot. There is a single object, KOI-571.05, which is in a region of particularly low detection efficiency ($<$0.002) and would normally warrant scrutiny as a potential false alarm, but this was validated as a member of the multi-planet Kepler-186 system by Quintana et al. (2014).

As was shown for solar-type host stars (Fulton et al., 2017) and subsequently for late K- and M-type dwarfs (Hirano et al., 2018; Cloutier & Menou, 2020; Van Eylen et al., 2021), the planet distribution of this sample contains two distinct clusters at about 1.3$R_{\oplus}$ (“super-Earths") and 2.2$R_{\oplus}$ (“sub-Neptunes"), with an intervening gap or “valley" at 1.6-1.7$R_{\oplus}$. Around solar-type hosts there is a deficit of (sub-)Neptunes on short–period (or high–irradiance) orbits, the so-called “Neptune desert" (Szabó & Kiss, 2011; Mazeh et al., 2016). There are no universally accepted boundaries for this desert but it is expected to depend on stellar mass (McDonald et al., 2019). Adopting the conservative criterion of a period of $<$2 days (corresponding to irradiances $S\gtrsim$ 200$S_{\oplus}$), the “desert" has a single unambiguous denizen in our sample: KOI–952.05 (Kepler–32f), a member of a multi-transiting planet system (Swift et al., 2013), has an orbital period of 0.74 days and radius of $2.04\pm 0.11$$R_{\oplus}$, consistent with but more precise than that determined by Morton et al. (2016) and Berger et al. (2018) but significantly larger than previous estimates (e.g., Swift et al., 2013; Mann et al., 2013). There are five other Earth- and super-Earth-size Ultra Short Period (USP, $P<1$ d) planets in the sample. The smallest object in the sample is the $0.52^{0.04}_{0.05}$$R_{\oplus}$ candidate KOI–7793.01.

Cloutier & Menou (2020) found that the planet radius distribution hosted by stars with masses 0.08–0.42M_⊙ (corresponding to spectral types M3-M8) differs from that of more massive stars in lacking a significant population of (sub)Neptunes. Ment & Charbonneau (2023) also found a very pronounced dearth of sub-Neptunes around mid-to-late type M dwarfs. Our sample does not substantially overlap in stellar mass with these analyses and falls off below 0.5$M_{\oplus}$, although there is at least one KOI hosted by a mid-M dwarf just above the radius valley (Fig. 3).

The location of the radius valley and its dependence on orbital period (a proxy for irradiance) is used to test models of escape of planetary atmospheres/envelopes (see Sec. 4). Following past work (e.g., Ho & Van Eylen, 2023, and references therein), we performed a regression on the full sample with support vector machines (SVM) as implemented in the Python sci-kit learn package, with a choice of regularization parameter of $C=10$, following Van Eylen et al. (2019). Super-Earth and sub-Neptune classification was initially assigned initially based on radii smaller or larger than 1.68$R_{\oplus}$ (the value predicted for this sample from the relation in Ho et al. (2024)), then classification was iteratively redone (10 times) using the SVM. Confidence intervals were calculated with $10^{4}$ Monte Carlo simulations where planet radii were randomized according to their uncertainties. We plot the final assignments and converged regression ($\log R_{p}=-0.03^{+0.01}_{-0.03}\log P+0.26^{+0.04}_{-0.01}$) as the black line in Fig. 4, along with previously published regressions.

Our regression is only weakly period-independent and contrasts with the negative slopes for M dwarf host stars derived by Van Eylen et al. (2021, blue dashed line in Fig. 4), from combined solar- and late-type Kepler host stars derived by Ho et al. (2024) and evaluated at the 0.59M_⊙ median sample mass of our sample (dashed magenta line), and for a sample of TESS-detected planets around M dwarfs with mass measurements (Bonfanti et al., 2024, dashed orange line). It is also inconsistent with the positive slope found by Cloutier & Menou (2020) for planets detected around late-type hosts by Kepler/K2 (dashed cyan line). However, the slope is consistent with that derived from a small sample of M dwarf planets with RV-determined masses by Luque & Pallé (2022). The cause of the disagreements between these estimates is not known, but could arise from sensitivity of the slope to objects with the shortest or longest orbital periods combined with refinement of radii by exclusion of binaries and improved stellar properties.

The detection efficiency-corrected radius distribution constructed for the entire small-planet sample is shown in the top panel of Fig. 5. The grey region spans the 68% confidence interval, and the individual planet radii are marked by the ticks at the bottom. The total occurrence is the integral over the distribution and is 1.84 planets per star. This presentation also shows the super-Earths and sub-Neptunes sa peaks separated by an intervening gap or “valley" around 1.6-1.7$R_{\oplus}$. We computed a valley depth, i.e. the ratio of the mean of the adjacent maximum to the minima, of 3.85. This is similar to that found for solar-type stars by Ho & Van Eylen (2023), and inconsistent with the trend of shallower valley with lower stellar mass found by Ho et al. (2024).

While Fig. 5 appears to resolve the super-Earth population into two overlapping peaks, this structure is not significant compared to the uncertainties in the reconstruction (represented by the shaded area). To verify this, for each of the 1000 bootstrap reconstructions we assessed whether the radius interval 1-1.7$R_{\oplus}$ contained one or two maxima and, in the latter case, calculated the ratio of the intervening minimum to the lowest of the two maxima. 75% of the reconstructions contain two maxima (at $\approx$1.2 and 1.4 $R_{\oplus}$), but among these, the ratio of the minimum to lowest maximum is continuously distributed between about 0.6 to unity, with most values near unity. This indicates that the two observed peaks are a random occurrence from a continuous distribution of possible outcomes.

There are other, more ambiguous features in the radius distribution, including a “shoulder" at 3$R_{\oplus}$, representing four planets of similar size (Fig. 2). The radius distribution also hints at the existence of a population of Earth-size objects distinct from super-Earths (second peak from left in Fig. 5). This is not an artefact of detection bias, since the efficiency only decreases smoothly with radius in this regime (Fig. 2). The prominent peak at 0.64$R_{\oplus}$ in Fig. 5 is produced by KOI–314.03 (Kepler–138b), which has a low detection probability due to its comparatively wide (10.3 day) orbit. It is a member of a three-planet system and was discovered only after analysis with optimized, automated routines following the discovery of the other, larger members (Batalha et al., 2013). Its presence in the sample could be just a statistical fluke, or it could represent a large population of “sub-Earths" (Sinukoff et al., 2017).

We divided the planets with assigned ages into halves that are younger and older than the median gyrochronologic age of 3.86 Gyr (mean ages of 2.6 and 5.4 Gyr, respectively). These are plotted red and blue curves in the bottom panel of Fig. 5. The sub-Neptune peak is much less prominent in the older sample while the super-Earth is larger, with the ratio between the two (integrated over 1.68–4$R_{\oplus}$ and 0.7–1.68$R_{\oplus}$, respectively) declining from 2.38 to 0.58 from the younger to older moieties. This is consistent with the evolution of planets by the loss and/or contraction of low molecular-weight envelopes over Gyr. Indeed, the integrated probabilities over mini-Neptunes (1.7–4$R_{\oplus}$) fall from 0.59 to 0.31 and is roughly equivalent to the increase in the Earth- and super-Earth (0.7–1.7$R_{\oplus}$) from 0.25 to 0.53. There is no evidence for a decrease in the mean radius of the sub-Neptunes with age nor a deepening of the radius valley as has been seen in other analyses of solar-type stars (Petigura et al., 2022; Chen et al., 2022) but the valley is indistinct in the older sample due to the small number of sub-Neptunes. A two-sided Kolmogorov-Smirnov test provide weak statistical support for the young and old samples to be drawn from different populations ($p=0.08$).

We investigated four systematic effects that could potentially mimic radius evolution: (1) an artefact due to a sensitivity of the radius distribution to the choice of divisor between younger and older sub-samples; (2) a variation in stellar properties with age that influences detection efficiency; (3) a difference in host star properties brought about by selecting for stars with detectable rotation period that in turn affects the detected planets; and (4) a variation in stellar properties (i.e. $T_{\rm eff}$ or mass) with age combined with a correlation between the planet population and stellar properties. We find no evidence that any of these is responsible for the observed difference in radius distribution between the younger and old sub-samples, and details are provided in Appendix A.

Figure 6 compares the period-radius distribution of the younger (red points) and older (blue points) sub-samples and hints at non-uniformity in this evolution. While there are both younger and older sub-Neptunes ($>1.7$$R_{\oplus}$, above the horizontal line) at periods $<25$ days, nearly all those at longer periods ($>$25 days) are younger. This is discussed further in Sec. 4.3.

Models that explain the dichotomy between super-Earths and sub-Neptunes predict a dependence of the location of the radius valley on orbital period or stellar irradiation similar that can be compared to observations (e.g., see Table 6 in Ho & Van Eylen, 2023). Both photoevaporation and core-powered mass loss models predict negative slopes with period, in the former due to stellar X-ray and/or EUV radiation as the driver (e.g., Owen & Jackson, 2012), and in the latter a consequence of the role of stellar bolometric irradiation (Gupta & Schlichting, 2021). For example, Gupta & Schlichting (2019) derived an index of -0.11 for core-powered loss, independent of stellar mass, while Owen & Kollmeier (2017) derive -0.16 for photoevaporative loss. Such models are consistent with the slopes reported by Van Eylen et al. (2021), Ho et al. (2024), and Bonfanti et al. (2024) (Fig. 4). In contrast, a model of planet formation in a gas-depleted inner disk predicts a positive slope of +0.11 (Lopez & Rice, 2018), consistent with the findings of Cloutier & Menou (2020). None of these are consistent with the near-zero slopes found by Luque & Pallé (2022) and our analysis ($-0.03^{+0.01}_{-0.03}$; Fig. 4). Weak dependence on period as well as on stellar luminosity (Fig. 3), suggests a mechanism largely independent of the host star, e.g. differences in water abundance (see Burn et al., 2024, for a discussion of some possible scenarios).

The depth of the radius valley relative to the adjacent super-Earth and sub-Neptune peaks has also been compared to model predictions. In particular, Ho et al. (2024) found that the valley became shallower (filled in) at lower host star masses. A trend could appear due to error in the radius ratio $R_{p}/R_{*}$ from fitting the transit lightcurve, which, all else being equal, is expected to be higher for shorter duration transits around smaller M dwarf host stars (Petigura, 2020). Astrophysical interpretations of a shallower radius valley include: (a) scatter in the XUV-emission luminosity of young stars at a given mass and age; (b) heterogeneity in the composition of planetary cores; and (c) collisions (Ho et al., 2024). Our failure to reproduce the trend motivates a deeper dive into the systematics that could affect reconstructions of the planet radius distribution.

The decline of the sub-Neptune population and commensurate rise in super-Earths between our younger and older sub-samples (Fig. 5) indicates that some sub-Neptunes have evolved into super-Earths over a timescale of $\sim$3 Gyr. On orbits with $P>20$ day this radius evolution has presumably caused some planets to leave the sample due to the low detection efficiency of super-Earths on wider orbits (grey contours in Fig. 2). The observed change is analogous to that seen around solar-type Kepler host stars (Berger et al., 2020; Sandoval et al., 2021; David et al., 2021; Chen et al., 2022), which has been explained by a scenario where sub-Neptunes are super-Earths with thick hydrogen-dominated envelopes which can escape by XUV-driven photoevaporation or core-powered mass-loss from rocky/icy cores (Owen, 2019; Gupta & Schlichting, 2021; Owen & Schlichting, 2024).

High-energy emission from solar-type host stars is most intense during the first few hundred Myr (Johnstone et al., 2021), so a Gyr timescale for planet evolution has been previously interpreted as supporting core-powered (Berger et al., 2020). However, more quantitative estimates of photoevaporative timescales are longer (King & Wheatley, 2021) and both mechanisms may operate at different times in different regimes (Owen & Schlichting, 2024). Moreover, our sample cannot be used to investigate shorter timescales because it does not include systems with ages much younger than 1 Gyr. This is because low mass stars take longer to converge to an initial condition-independent rotation sequence, where gyrochronology can be applied.

In a scenario where sub-Neptunes are rocky cores with thick envelopes of high molecular weight volatiles, specifically H₂O (Turbet et al., 2020), radius evolution could be at least partly explained by contraction due to cooling (Howe & Burrows, 2015; Vazan et al., 2018), and/or chemical reaction between the envelope and the rocky core (Kite et al., 2020; Kim et al., 2023). These processes have timescales different from that of atmospheric escape mechanisms.

The near-complete absence of sub-Neptunes in the “desert" very close to the host stars is consistent with the canonical picture of envelope loss driven by stellar irradiation. However, the disappearance of sub-Neptunes elsewhere appears to be non-uniform (Fig. 6). There is a striking paucity of older sub-Neptunes on the widest orbits ($>$20 days), where irradiation is lowest and escape should be slower than at intermediate orbits (2–20 days) where older sub-Neptunes are found (Fig. 6).

The top panel of Figure 7 shows planet radius versus stellar irradiance, comparing the younger and older sub-samples. Irradiance is based on bolometric luminosities calculated using interpolations of the corrections to $K$-band luminosities tabulated in Pecaut & Mamajek (2013), and assuming near-circular orbits. Constructed irradiance distributions for younger and older sub-Neptunes ($R_{p}>1.7$$R_{\oplus}$, dotted grey line) are shown in the bottom panel. Both Gaussian kernel density estimators (solid curves) and histograms with Poisson statistics errors (shaded bars) suggest a deficit of older sub-Neptunes at irradiances below several times the present terrestrial value ($S_{\oplus}$). The fraction of sub-Neptunes in the older-sub-sample falls from 53% at irradiances $>$5$S_{\oplus}$ to 17% at lower irradiance (see red and blue histograms). This cannot be due to a differences in detection efficiency since planets on wider orbits should be more readily detected around older, photometrically quieter stars. Nor is there evidence for a combination of correlations between age with host star mass (or luminosity), and stellar mass with planet period distribution (Figs. 14 and 15 in Appendix A) that could mimic this effect. There are no sub-Neptunes in the sample below an irradiance of 1$S_{\oplus}$, although this could be due in part to low detection efficiency at long orbital periods. The two planets that are marginally that are least irradiated are KOIs 812.03 (Kepler–235e) and 952.03 (Kepler–32d), with ages of $7.8\pm 1.43$ and $5.60\pm 0.96$ Gyr, respectively (Gaidos et al., 2023).

The persistence of “warm" sub-Neptunes and disappearance of cool sub-Neptunes at older ages could be due to an inflation mechanism that depends on stellar irradiation or planet equilibrium temperature and hence proximity to the host star. Turbet et al. (2020) showed that primarily rocky planets with water-dominated atmospheres could have sub-Neptune-like radii at irradiances exceeding a “runaway" threshold where H₂O cannot condense into a liquid ocean.

Kopparapu et al. (2013) and Turbet et al. (2019) found the irradiance threshold for a runaway greenhouse on Earth-like planets to be $\approx$1.05$S_{\oplus}$. Around M dwarfs stars the threshold is slightly lower because of reduced Rayleigh scattering and increased absorption by atmospheric H₂O at the longer wavelengths where cooler stars emit more radiation; the value is predicted to be 0.9$S_{\oplus}$ for $T_{\rm eff}$=3900K or M0 spectral type (dashed vertical green line in Fig. 7. Three-dimensional general circulation models find a somewhat higher threshold that depends on planet rotation rate (Leconte et al., 2013; Kodama et al., 2018, ; vertical orange dashed line). The higher Bond albedo of a cloudy steam atmosphere (that of Venus is around 0.76 (Moroz et al., 1985) but measurements of the phase function suggest a value as high as 0.9 (Mallama et al., 2006)⁵⁵5The albedo depends would depend strongly on cloud droplet size and the high reflectivity of Venus is largely due to sulfuric acid aerosols and a similar scenario might not apply to steam atmospheres.), might move the runaway threshold to $\sim 5$$S_{\oplus}$. This is the approximate boundary where sub-Neptunes become persistent (vertical dashed magenta line). Finally, atmospheres containing significant H/He will have a lower threshold (Innes et al., 2023).

While a decline in sub-Neptune occurrence across the water condensation boundary could, in principle, be explained by the collapse of steam atmospheres, it is unclear if H₂O alone can explain the observed distributions with radius and age. In Fig. 7 we plot predictions of Turbet et al. (2020) for H₂O mass fractions between 0.1 and 25% (cyan lines), assuming a 3 $M_{\oplus}$ rocky core mass chosen to reproduce super-Earth radii. This comparison suggests that H₂O mass fractions exceeding 25% are necessary to explain some sub-Neptunes. Models of H₂O-rich planets produce radii up to 3$R_{\oplus}$ (Burn et al., 2024), but Vivien et al. (2022) found that pure water/steam atmospheres on sub-Earths become unstable if the atmospheric scale height exceeds 0.1 times the planet radius. Moreover, some internal governing mechanism (e.g., reaction with a magma ocean, Schaefer et al., 2016) is needed to explain the radius evolution over Gyr since the change in luminosity of these very low-mass host stars is small and the response of an atmosphere nearly instantaneous.

Alternatively, the inflated radii of sub-Neptunes and the differential evolution of cool vs. warm objects could be explained by atmospheres of mixed H/He and H₂O. H/He lowers the mean molecular weight and decreases the irradiance threshold for a runaway water-vapor greenhouse (Innes et al., 2023). Slow escape of H/He would eventually collapse sub-Neptunes to smaller objects. At an elevated irradiances above the H₂O condensation threshold, H₂O could be persist in the upper atmosphere where its radiative cooling the lowers temperature and inhibits H/He escape (Yoshida et al., 2022). In a cooler sub-Neptune, water condenses into clouds deeper in the atmosphere and cannot play that role. The upper atmosphere is therefore hotter, H/He escape faster, and the sub-Neptunes eventually collapses to a rocky core with a geometrically thin atmosphere.

Finally, we searched for irradiance-dependent evolution of sub-Neptunes among Kepler transiting planets around solar-type host stars with isochrone ages determined by Berger et al. (2018), restricting the sample to those with ages more precise than 30%. Although the median age of the sample is 5 Gyr, we use the same division of 3.86 Gyr between young and old populations for comparison. No such effect is obvious in the irradiance distributions (Fig. 8). The ratio of young to old planets is statistically indistinguishable ($\approx 0.44$) for planets experiencing $S>5$$S_{\oplus}$ and $S<5$$S_{\oplus}$. Three explanations for the difference between the two samples are (1) the late K/early M dwarf sample is affected by small-number statistics and is a statistical fluke; (2) the ages of the Berger et al. (2018) sub-sample with relatively cool planets are more uncertain because the host stars are likely to be late G and early K dwarfs which evolve very slowly in a color-magnitude diagram; and (3) there is a fundamental difference in the evolution of the two planet populations. At a given bolometric irradiance, planets around M dwarfs will experience far more X-ray and UV irradiation than their solar-type counterparts at (McDonald et al., 2019) and it is this radiation that drives photoevaporation of H/He atmospheres (Owen, 2019). Thus temperate sub-Neptunes around solar-type stars will have retained H/He atmospheres while their counterparts around M dwarfs have lost most or all their H/He and have radii controlled by the behavior of water-dominated atmospheres.

These different pathways are illustrated in Fig. 9, where the Kepler M dwarf sample spans Zones I, II, and III (but not Zone V due to small sample size), while the Berger et al. (2018) sample spans Zones I, II, and IV. In Zone I, rapid stellar-driven escape of H/He is not suppressed by atmospheric H₂O and leads to planets with extended steam-rich atmospheres or rocky super-Earths (Fig. 9). In Zone II, water-rich (but not water-poor) sub-Neptunes retain thick H/He envelopes. In Zone III (only around very low-mass stars), H/He envelopes are lost, leaving “ocean worlds" and rocky super-Earths, and in contrast to Zones 4 and 5, where these envelopes are retained. Measurements of planet mass by RVs and observations of atmospheres by JWST for examples around much closer, brighter host stars could distinguish between these scenarios.

Biases and selection effects that could possibly mimic the observed radius evolution are summarized in Sec. 3.3 and evaluated in Appendix A. There are other effects that could quantitatively impact our results. We restricted our sample to planet host stars with detected rotational signals in Kepler lightcurves and these are a biased subset of the overall population. Younger stars will be rapidly rotating, magnetically active, and highly spotted and their rotational periods can be readily measured. However, due that rotational variability as well as intense flaring can inhibit the detection of transiting planets (Gaidos et al., 2017). On the other hand, the oldest stars will be slowly rotating, less spotted, and exhibit less rotational variability, and the rotation periods can be similar to the Kepler quarterly roll ($90$ days) and thus more difficult to retrieve (Santos et al., 2019). This is particularly true for the coolest stars in the sample (mid M-dwarfs) because these stars have the longest rotation periods for a given age. These selection effects can manifest themselves in the overall distribution of M dwarf planet host ages (e.g., Gaidos et al., 2023). However, we see no trend of planet properties (i.e., radius) with host star properties (i.e., luminosity, Fig. 3).

The overall occurrence rate estimated from KOIs with rotation periods could be biased upwards to the extent that stellar rotation and planetary orbits are aligned (low stellar obliquity Albrecht et al., 2022). Stars viewed at a high inclination (equator on) will have the strongest rotational signal and their planets will be more likely to transit compared to the overall population. We quantified this effect by determining the fractional numbers of KOIs around K/M dwarfs observed by Kepler with rotation periods in Santos et al. (2019) vs. K/M dwarfs lacking rotation periods. We applied the same cuts ($M_{K}>4.5$, RUWE $<1.4$, and not otherwise identified as binaries) to both samples, but we did not apply a $T_{\rm eff}$ cut since temperatures are only available for the Santos et al. (2019) sample. The KOI fractions are $2.70\pm 0.24$% and $2.50\pm 0.39$%, respectively, with uncertainties from counting statistics. These are thus indistinguishable and we conclude the effect is negligible among these types of stars, perhaps because of the high (%73) success rate of detecting rotation among these stars.

Transit parameters and stellar parameters were not computed self-consistently: the former adopted DR25 values and the latter were based on revised luminosities and new empirical relations. Stellar density is covariant with transit impact parameter which, along with $T_{\rm eff}$ (which controls limb darkening coefficients) is covariant with the planet-to-star radius ratio $R_{p}/R_{*}$. Lastly, our filtering of binaries from the Kepler late K- and early M-dwarf target sample was less thorough because AO imaging is not available for the overall Kepler late K/early M dwarf population as it is for the KOIs.