Scaling K2. VII. Evidence For a High Occurrence Rate of Hot Sub-Neptunes at Intermediate Ages

Praesepe (NGC 2632, M44) is a 600–800 Myr old open cluster located about 180 pc away. The age, distance, and metallicity of the cluster have been studied extensively (see, e.g., Brandt & Huang 2015; Gossage et al. 2018 and references therein). The cluster has slightly supersolar metallicity, with estimates ranging from [Fe/H] = 0.12 ± 0.04 (Boesgaard et al. 2013) to 0.21 ± 0.01 (D'Orazi et al. 2020), depending on the calibration and methods used. Fortuitously lying in the ecliptic, it was observed by K2 in campaigns 5, 16, and 18 (the latter of which had the same pointing as campaign 5). We use the member list compiled by Rebull et al. (2017) in their analysis of the rotation periods of Praesepe members in K2 data, which originally included the targets observed in campaign 5 (and 18), and was subsequently augmented by additional members observed in campaign 16. Figure 1 shows the spatial distribution of the 1030 Praesepe targets observed by K2 across the sky.

Zoom In Zoom Out Reset image size

As part of our search of the full K2 data set, presented in Zink et al. (2021), the Praesepe targets were processed through our standard K2 data reduction and transit detection pipeline detailed in Zink et al. (2020a); in brief, light curves are extracted using EVEREST (Luger et al. 2016, 2018), detrended using Gaussian process regression and harmonic removal, searched for periodic transit signals using TERRA (Petigura et al. 2013), and the resulting signals uniformly vetted using the EDI-Vetter tool (Zink et al. 2020a). Although there is some overlap between the campaigns, for this analysis we treat each campaign separately, rather than combine light curves from multiple campaigns; this limits our period sensitivity to below 40 days but greatly simplifies the completeness and reliability corrections in the subsequent analysis. We compared the noise properties of the Praesepe targets to the field K2 field star data set (after removing the known young stars in Upper Scorpius, Pleiades, Hyades, and Praesepe)—young stars rotate more rapidly and therefore, although the Praesepe members studied here are largely on the main sequence, their light curves could potentially have increased correlated noise. This could make transit signals more difficult to detect than would be reflected by the detection efficiency previously measured for the pipeline across all targets, and potentially necessitate an updated detection efficiency measurement specifically for the Praesepe targets. We examine the Combined Differential Photometric Precision (CDPP) values (Christiansen et al. 2012) established for Kepler light curves. The left panel of Figure 2 shows the 8 hr rms CDPP values as a function of the Kepler magnitude for the Praesepe targets in blue, overlaid on an orange heatmap showing the 8 hr rms CDPP values for our full K2 data set. The Praesepe targets largely follow the broad features of the overall sample, except at the very faintest magnitudes. Zink et al. (2021) noted that stars with an 8 hr rms CDPP above 1200 ppm do not contribute meaningfully to the derived occurrence rates, as they have no sensitivity to the planet parameter range under study. ¹³ Applying the same CDPP cut as Zink et al. (2021) removes 429 targets, largely very faint M dwarfs smaller than 0.33 R_⊙. The 601 remaining Praesepe targets are shown in the right-hand panel of Figure 1.

Zoom In Zoom Out Reset image size

In addition to the overall scatter in a light curve, we can also measure the extent of the correlated noise. Burke & Catanzarite (2017) showed that calculating the CDPP "slope," by measuring how the rms CDPP values change as a function of the timescale over which they are calculated, could capture this information. For white noise, the rms CDPP values should decrease with increasing duration as $\sqrt{{N}_{\mathrm{short}}/{N}_{\mathrm{long}}}$, where N_short and N_long are the shorter and longer durations over which the CDPP is calculated, respectively. The right panel of Figure 2 shows the normalized distributions of the CDPP slopes for the 2 and 8 hr rms CDPP values of the 601 Praesepe targets in blue and the remainder of the K2 targets in orange. White noise would produce a CDPP slope of $\sqrt{2/8}=0.5$, and both populations peak at just above that value (and drop off rapidly below it), indicating largely well-behaved light curves. The Praesepe stars have a slightly broader distribution up to a CDPP slope of ∼1, indicating a marginally higher proportion of light curves with correlated noise than the full data set, although above ∼1.1 (where binning over a longer timescale actually increases the noise, indicating highly correlated noise) the fractional proportion of the two populations is very similar. These comparisons indicate that the noise properties of the Praesepe targets are similar to the full K2 data set and that the detection pipeline completeness and reliability measured on the full data set in Zink et al. (2020b) are applicable here.

The pipeline detected 10 planet candidates orbiting nine of the 601 Praesepe targets, originally published in the full K2 catalog provided in Zink et al. (2021). The candidates are listed in Table 1—they range from 1.2 to 3.9 R_⊕, with periods from 1.67 to 21.2 days, typical of the large population of warm sub-Neptunes and super-Earths revealed by Kepler. Given the interest in detecting and characterizing planets earlier in their formation and evolution histories, the K2 Praesepe data have been extensively searched for individual planets, and seven of the 10 candidates are already confirmed in the literature, as indicated in Table 1. However, this paper presents the first statistical sample of planet candidates from Praesepe for population analysis. We note that one of the unconfirmed planet candidates (EPIC 211885995.01) has a consistency score falling below our cut-off of 0.75 (see Section 3.13 of Zink et al. 2020a), and indeed a closer examination of the light curve indicates that the transit signal is likely related to stellar variability; we remove this target from further analysis, leaving nine remaining candidates shown in Figure 3.

Zoom In Zoom Out Reset image size

Table 1. Summary Parameters of the 15 Planet Candidates and their Host Stars Identified in Praesepe and Hyades

EPIC ID	C#	CP Name	Rp	Orbital Period	R*	Teff	[Fe/H]	Sp.T
			(R⊕)	(days)	(R⊙)	(K)	(dex)
Praesepe
211822797.01	16	K2-103 b a	${1.819}_{-0.090}^{+0.101}$	${21.171661}_{-0.003290}^{+0.003413}$	0.598 ± 0.018	3806 ± 138	0.110 ± 0.235	M1
211885995.01	18	f	${3.087}_{-0.131}^{+0.142}$	${9.863718}_{-0.006848}^{+0.005667}$	0.605 ± 0.019	3676 ± 138	0.098 ± 0.235	M1
211913977.01	16	K2-101 b a	${2.006}_{-0.091}^{+0.086}$	${14.677689}_{-0.001147}^{+0.001111}$	0.755 ± 0.019	4694 ± 39	0.109 ± 0.036	K3
211916756.01	18	K2-95 b b	${3.402}_{-0.134}^{+0.134}$	${10.135158}_{-0.000635}^{+0.000698}$	0.415 ± 0.014	3576 ± 138	0.014 ± 0.235	M3
211922849.01	5	⋯	${1.440}_{-0.059}^{+0.066}$	${9.443274}_{-0.002157}^{+0.002052}$	0.517 ± 0.016	3652 ± 138	0.175 ± 0.235	M2
211964830.01	16	K2-264 b c	${2.153}_{-0.105}^{+0.101}$	${5.839636}_{-0.000527}^{+0.000477}$	0.473 ± 0.015	3594 ± 138	0.239 ± 0.235	M3
211964830.02	16	K2-264 c c	${2.503}_{-0.118}^{+0.117}$	${19.662285}_{-0.002070}^{+0.002742}$	0.473 ± 0.015	3594 ± 138	0.239 ± 0.235	M3
211969807.01	5	K2-104 b a	${1.773}_{-0.078}^{+0.085}$	${1.974302}_{-0.000192}^{+0.000168}$	0.505 ± 0.017	3693 ± 138	0.125 ± 0.235	M2
211990866.01	18	K2-100 b a	${3.899}_{-0.103}^{+0.101}$	${1.673847}_{-0.000042}^{+0.000040}$	1.208 ± 0.020	6116 ± 16	0.274 ± 0.013	G1
212035441.01	18	⋯	${1.198}_{-0.066}^{+0.063}$	${2.714743}_{-0.000276}^{+0.000277}$	0.415 ± 0.014	3634 ± 138	0.310 ± 0.235	M3
Hyades
210897587.01	13	K2-155 b d	${1.785}_{-0.065}^{+0.066}$	${6.34399}_{-0.000050}^{+0.000044}$	0.584 ± 0.019	3893 ± 56	−0.831 ± 0.051	K7
210897587.02	13	K2-155 c d	${2.088}_{-0.075}^{+0.081}$	${13.852825}_{-0.000280}^{+0.000442}$	0.584 ± 0.019	3893 ± 56	−0.831 ± 0.051	K7
246711015.01	13	⋯	${1.233}_{-0.094}^{+0.116}$	${13.284219}_{-0.001880}^{+0.002846}$	0.395 ± 0.017	3771 ± 138	−0.112 ± 0.235	M1
246711015.02	13	⋯	${0.505}_{-0.084}^{+0.088}$	${24.213118}_{-0.003407}^{+0.003408}$	0.395 ± 0.017	3771 ± 138	−0.112 ± 0.235	M1
247589423.01	13	K2-136 c e	${3.150}_{-0.095}^{+0.099}$	${17.306369}_{-0.000152}^{+0.000171}$	0.725 ± 0.018	4231 ± 41	−0.112 ± 0.235	K5

Notes. C# = campaign in which the signal was detected by our pipeline. CP name = confirmed planet name in the literature.

^a Mann et al. (2017). ^b Obermeier et al. (2016). ^c Rizzuto et al. (2018). ^d Hirano et al. (2018). ^e Mann et al. (2018). ^f Removed from further analysis due to low consistency score.

Download table as:  ASCIITypeset image

In campaigns 4 and 13, K2 also observed the Hyades open cluster, a smaller, closer, near-twin of Praesepe. Age estimates range from 625 ± 25 Myr (Perryman et al. 1998), to 680 Myr (Gossage et al. 2018), to 800 Myr (Brandt & Huang 2015; David & Hillenbrand 2015). It also has a slightly supersolar metallicity of [Fe/H] = 0.14 ± 0.05 dex (Perryman et al. 1998) and has also been meticulously combed for planets. A list of potential Hyades members was assembled from membership lists in the literature (e.g., Röser et al. 2011, 2019; Goldman et al. 2013; Douglas et al. 2016; Lodieu et al. 2019, among others), the list of Hyades members submitted to the original K2 proposal call, and a run of BANYAN Σ (Gagné et al. 2018), then compared to the list of objects with K2 light curves (L. M. Rebull et al. 2023, in preparation), for a total of 300 targets. This is significantly smaller than our Praesepe target list, in part because Hyades is closer and therefore subtends a larger area on the sky, a smaller fraction of which falls in the K2 field of view; many Hyades cluster members were not observed by K2. Applying the same noise and stellar parameter cuts as before leaves 101 Hyades targets to add to our sample. We examined the noise properties of these targets, as above for Praesepe, and found them to be comparable. Similarly to our Praesepe sample, the Hyades targets are on average cooler than the full K2 sample.

Our pipeline recovered five planet candidates from these 101 Hyades targets, listed in Table 1. Three are known planets, K2-155 b and c (Hirano et al. 2018) and K2-136 c (Mann et al. 2018), although in both cases our pipeline missed known smaller or longer-period additional known planets in the system. The final two are a pair of small (0.5 R_⊕ and 1.2 R_⊕) candidates orbiting the faint (m_Kep = 15.635) target EPIC 246711015, which have thus far proven difficult to confirm.

In order to calculate occurrence rates, we use the ExoMult forward modeling software (Zink et al. 2019, 2020b), which, in brief, produces synthetic K2 planet samples drawn from an underlying parametric planet distribution function, and simulates their detection throughput using the completeness measured in Zink et al. (2021). The synthetic population of planets and the empirical sample were compared using the two-sample Andersen–Darling test statistic and then optimized for each relevant observable. As the aim of the analysis is to compare the intermediate-age planet sample to the field K2 sample and to isolate any age correlation, careful attention must be paid to the form of the parametric model to ensure it captures the differences between the samples to the extent possible.

Given its significantly larger data set, we begin by examining the Praesepe sample. The Praesepe candidates orbit dwarf stars from G1 to M3, so we restrict our population analysis to GKM stars spanning the range 0.33 < R_* < 1.5 R_⊙, 3500 < T_eff < 6500 K, −0.5 < [Fe/H] < 0.5 dex, and log g > 4.0. Our stellar parameters were uniformly computed for the K2 sample as described in Hardegree-Ullman et al. (2020), using a random forest regression on photometric colors trained on stars with LAMOST spectra, ultimately providing spectral types, effective temperatures, surface gravities, metallicities, radii, and masses. Where the machine learning parameters were unavailable, they were supplemented with parameters from isochrone grid matching, as per Zink et al. (2023). For the well-behaved Praesepe targets, this reduces the data set to nine planet candidates and 487 stars. For the well-behaved (8 hr CDPP < 1200 ppm) targets in the full K2 field star data set (with the stars in the four identified young clusters removed), this leaves 520 planet candidates and 108,502 stars, referred to as the "field K2 sample." The rest of the analysis will build on these two samples.

The two stellar samples have significantly different distributions in stellar temperature and metallicity, as shown in Figure 4. There are well-established correlations between the frequencies of planets over these stellar parameters. The median iron abundance of our Praesepe sample is 0.146 dex, consistent with the slightly supersolar metallicity of the cluster. Although sub-Neptunes do not have as steep a metallicity correlation as giant planets (Petigura et al. 2018; Zink et al. 2023), there is still a positive correlation that could lead to an enhancement in the frequency of small planets in Praesepe compared to the field (slightly subsolar) K2 sample if not taken into account. Similarly, the field K2 sample includes a much larger fraction of higher-mass stars than the Praesepe sample, and stellar mass has been shown to be anticorrelated with the frequency of small, short-period planets (e.g., Mullally et al. 2015). To account for these differences we include host star metallicity and stellar effective temperature terms in our parametric model.

Zoom In Zoom Out Reset image size

For the parameter space probed by the nine Praesepe candidates (1–4 R_⊕, 1–30 days), previous analyses using the older and much larger Kepler sample have found that the population is well represented by a broken power law in orbital period with a break at ∼10 days (Mullally et al. 2015), and a single power law in planet radius (e.g., Mulders et al. 2018; Petigura et al. 2018). With only nine Praesepe planet candidates, our statistical power to simultaneously fit many parameters is low, and indeed an attempt to fit the full parametric model with a broken power law in orbital period and single-power-law planet radius was unable to meaningfully constrain the overall occurrence rate, let alone the exponents of the various power laws. However, by concentrating on hot sub-Neptunes, we can simplify the parametric fit in order to make more robust measurements. For the following analysis, we examine candidates with orbital periods 1–10 days, inwards of the expected period break and therefore able to be modeled with a single power law, and planet radii 1.7–4 R_⊕, capturing the sub-Neptune population above the planet radius valley (Fulton et al. 2017; Van Eylen et al. 2018; Hardegree-Ullman et al. 2020). These planets are particularly susceptible to mass-loss mechanisms, due to their proximity to their host stars, and their weakly bound gaseous upper atmospheres. In addition, the same parameter range is investigated in Zink et al. (2023), which allows for a straightforward comparison to their results.

In this space, we fit a parametric model n of the form

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{4}n}{d\mathrm{log}P\ d\mathrm{log}R\ {{dT}}_{\mathrm{eff}}\,d[\mathrm{Fe}/{\rm{H}}]\,}\\ \quad =\,f\ \cdot {P}^{\beta }\ \cdot {R}_{p}^{\alpha }\ \cdot {10}^{[\mathrm{Fe}/{\rm{H}}]\cdot \lambda +\tfrac{{T}_{\mathrm{eff}}}{1000\,{\rm{K}}}\cdot \gamma }\end{array}\end{eqnarray} \tag{ 1 }$$

where f captures the number of planets per star within the range of our sample, and α, β, γ, and λ are all free parameters, as in Zink et al. (2023). We start by fitting the field K2 sample, allowing all tunable parameters to float. In our targeted stellar and planet parameter space, this leaves 149 planet candidates from Zink et al. (2021). Examining the corner plot shows the convergence of the Markov Chain Monte Carlo parameter chains and highlights their independence. The fitted values are given in Table 2. For the field K2 sample, the occurrence rate is ${16.54}_{-0.98}^{+1.00}$%, and we recover the strong negative radius scaling (α = − 2.6 ± 0.2) and positive period power-law scaling (β = 2.1 ± 0.1) seen in previous studies (e.g., α = − 2.7 and β = 1.97 from Zink et al. (2023). We also see the expected weak but positive dependence on stellar metallicity (0.19 ± 0.09), and the negative dependence on stellar effective temperature ($-{0.072}_{-0.017}^{+0.018}$).

Table 2. The Measured Parameters of the Parametric Model for the Field K2 Sample, the Praesepe Sample, and the Combined Praesepe+Hyades Sample, Over the Two Planet Parameter Ranges Defined in the Analysis

Sample	Rp (R⊕)	Period (days)	f	α (Radius)	β (Period)	λ	γ
K2	1.7–4	1–10	${0.1654}_{-0.0098}^{+0.0100}$	$-{2.565}_{-0.209}^{+0.179}$	${2.103}_{-0.111}^{+0.109}$	${0.191}_{-0.085}^{+0.085}$	$-{0.0718}_{-0.0171}^{+0.0179}$
	1.8–6	1–12.5	${0.2076}_{-0.0108}^{+0.0110}$	$-{3.697}_{-0.153}^{+0.144}$	${1.546}_{-0.082}^{+0.082}$	${0.200}_{-0.079}^{+0.084}$	$-{0.0725}_{-0.0167}^{+0.0163}$
Praesepe	1.7–4	1–10	${0.983}_{-0.253}^{+0.272}$	−2.565	2.103	${1.06}_{-0.64}^{+0.71}$	$-{0.246}_{-0.270}^{+0.176}$
	1.7–4	1–10	${0.7927}_{-0.2482}^{+0.3081}$	−2.565	${2.301}_{-1.116}^{+1.430}$	${1.082}_{-0.708}^{+0.776}$	$-{0.262}_{-0.0285}^{+0.0193}$
	1.7–4	1–10	${1.066}_{-0.260}^{+0.319}$	$-{5.04}_{-1.91}^{+1.63}$	2.103	${1.16}_{-0.66}^{+0.68}$	$-{0.277}_{-0.282}^{+0.185}$
	1.7–4	1–10	${0.9448}_{-0.3182}^{+0.4112}$	$-{5.470}_{-2.008}^{+1.752}$	${3.168}_{-1.727}^{+4.588}$ a	${1.333}_{-0.695}^{+0.729}$	$-{0.301}_{-0.289}^{+0.193}$
	1.7–4	1–10	${0.9999}_{-0.2455}^{+0.2911}$	−2.565	2.103	0.191	0.0718
	1.7–4	1–10	${0.9474}_{-0.2346}^{+0.2695}$ a	−2.565	2.103	0.191	0.0718
	1.8–6	1–12.5	${1.0451}_{-0.2557}^{+0.2544}$	−3.697	1.546	${1.076}_{-0.668}^{+0.658}$	$-{0.233}_{-0.276}^{+0.163}$
	1.8–6	1–12.5	${0.8736}_{-0.2497}^{+0.3304}$	$-{6.242}_{-1.969}^{+1.585}$	${1.289}_{-0.630}^{+0.742}$	${1.159}_{-0.642}^{+0.707}$	$-{0.265}_{-0.279}^{+0.175}$
Prae+Hya	1.7–4	1–10	${1.1072}_{-0.2771}^{+0.2601}$	−2.565	2.103	${0.348}_{-0.535}^{+0.602}$	$-{0.268}_{-0.222}^{+0.156}$
	1.7–4	1–10	${1.2072}_{-0.2678}^{+0.2977}$	$-{6.017}_{-1.772}^{+1.528}$	2.103	${0.455}_{-0.539}^{+0.583}$	$-{0.295}_{-0.236}^{+0.164}$
	1.7–4	1–10	${0.9094}_{-0.2620}^{+0.3052}$	−2.565	${5.226}_{-2.717}^{+5.174}$ a	${0.415}_{-0.600}^{+0.650}$	$-{0.320}_{-0.236}^{+0.192}$
	1.7–4	1–10	${0.9454}_{-0.2934}^{+0.3645}$	$-{6.334}_{-1.962}^{+1.625}$	${4.582}_{-2.218}^{+3.877}$ a	${0.681}_{-0.584}^{+0.685}$	$-{0.239}_{-0.239}^{+0.183}$
	1.7–4	1–10	${1.1631}_{-0.2475}^{+0.2578}$	−2.565	2.103	0.191	0.0718
	1.8–6	1–12.5	${1.1640}_{-0.2432}^{+0.2693}$	−3.697	1.546	${0.310}_{-0.514}^{+0.590}$	$-{0.276}_{-0.226}^{+0.155}$
	1.8–6	1–12.5	${0.8925}_{-0.2547}^{+0.3054}$	$-{6.889}_{-1.953}^{+1.522}$	${2.070}_{-0.761}^{+1.094}$	${0.575}_{-0.599}^{+0.626}$	$-{0.303}_{-0.238}^{+0.172}$

Notes. Bold values indicate parameters that were fixed in each analysis.

^a A repeat of the analysis on the line immediately above without the CDPP < 1200 ppm cut, to confirm that it does not significantly affect the result.

Download table as:  ASCIITypeset image

In the Praesepe sample, we have three candidates in this parameter space (K2-100 b, K2-104 b, and K2-264 b), shown in the red dashed box in Figure 3. Even our simplified model is poorly constrained with only three candidates. We attempt to fit all five parameters (f, α, β, γ, and λ), which gives an occurrence rate of 95%, however, not unexpectedly, the fit does not converge, specifically for β (the period power-law coefficient). This is perhaps driven by the fact that the three candidates occupy only the middle ∼half of the 1–10 day range. We proceed by first fixing the radius power-law coefficient, α, while fitting the remaining parameters, then fixing β while fitting the remaining parameters, and finally fixing both α and β, using the values for α and β from the full K2 field sample. Under the assumption that the radius and period power-law coefficients are independent, which is likely appropriate over this limited region of parameter space, this allows us to explore the extent to which we can test for differences between these values for Praesepe compared to the full K2 data set, and in the overall occurrence rates.

The results for all fits are shown in Table 2. The resulting occurrence rates vary from ${79}_{-25}^{+31}$% to ${107}_{-26}^{+32}$% of GKM stars hosting a hot sub-Neptune. ¹⁴ The uncertainties are large, due to the small size of the sample, but the rates are still 2.5–3.5σ higher than the 17% occurrence rate recovered from the full K2 sample. In each case, the parameters describing the correlations with period, planet radius, stellar effective temperature, and metallicity (α, β, λ, and γ) lie within 1.3–1.6σ of the values from the full K2 sample, indicating that our (small) Praesepe sample cannot demonstrate any significant differences in the strength of these correlations at younger ages. The only significant difference between the two samples is in the overall occurrence rate, which is significantly higher for the young small (1.7–4 R_⊕), short-period (1–10 days) hot sub-Neptunes in Praesepe than for field stars. To test this, we perform a final fit where we fix all but the occurrence rate, f, assuming the null hypothesis that there is no difference in the correlations between intermediate-age and field stars, and find an occurrence rate of 100%.

Given the high occurrence rate of hot sub-Neptunes found in Praesepe, we combine the Praesepe and Hyades targets to create a larger target list and test the result. The stellar effective temperature and metallicity distributions of the Hyades targets are shown in Figure 5. The median metallicity of our 101 Hyades targets is [Fe/H] = −0.002, which is consistent with the cluster being slightly supersolar on average, given the typical uncertainty (0.235 dex) on our metallicity values and the small size of the sample. In the planet parameter space defined in Section 3.1, we add one additional Hyades planet (K2-155 b) to our sample. We repeat the same sequence of fits as above, using the parametric model described in Section 3.1. We first allow all five tunable parameters to float, which again does not converge for β (the period power-law coefficient), the additional 6.3 day Hyades planet being close in period to the 5.8 day Praesepe planet and likely not providing additional constraints over the 1–10 day period range. We then fix the α parameter, the β parameter, both parameters, and finally all four non-f parameters. The fitted parameters α, λ and γ lie within 0.4–1.8σ of the values from the field K2 sample; the β value is poorly constrained. As before, we cannot exclude the null hypothesis that the shape of the underlying planet population is the same for intermediate-age stars as for field stars. The occurrence rates range from ${90}_{-26}^{+31}$% to ${121}_{-27}^{+30}$%, again indicating a significantly (2.8–3.8σ) higher rate of hot sub-Neptunes around these intermediate-age stars.

Zoom In Zoom Out Reset image size

One intriguing property of small, short-period planets is their bimodal distribution in planet radius (Fulton et al. 2017; Van Eylen et al. 2018; Hardegree-Ullman et al. 2020), or perhaps more accurately trimodal in planet density (Luque & Pallé 2022). Many of the ideas proposed to explain this distribution involve mass loss from planets early in their histories, including through photoevaporation (Lopez et al. 2012; Lopez & Fortney 2013; Owen & Wu 2013, 2017; Rogers & Owen 2021) or core-powered mass loss (Ginzburg et al. 2018; Gupta & Schlichting 2019, 2020, 2021). Perhaps the very high rate of hot sub-Neptunes found in Praesepe and Hyades is in the process of reducing to the smaller rate we see in the K2 field stars.

In the photoevaporation scenario, the gaseous atmosphere is stripped away by high energy (XUV) radiation from the host star. The extent of the mass loss, and the timescale over which it occurs, are driven by the initial planetary core mass, the initial envelope mass fraction, the orbital separation of the planet, and the stellar spectral energy distribution. Lopez et al. (2012) investigated the role of photoevaporation mass loss for a specific system, finding that the 1.8 R_⊕ planet Kepler-11 b could have lost a factor of 2–3 in radius and 1.5–5 in mass between ∼100 Myr and 10 Gyr. Lopez & Fortney (2013) extended this study, predicting that planets in the range 1.8–4.0 R_⊕ should become significantly less common on orbits <10 days. Subsequently, Rogers & Owen (2021) concluded that although some super-Earths are born rocky, approximately four times as many are born with large H/He atmospheres that are subsequently stripped by photoevaporation, a number that is tantalizingly in agreement with our finding of 4–5 times as many young sub-Neptunes as old sub-Neptunes if a significant fraction of the deficit evolve to become super-Earths.

In the core-powered mass loss scenario, the gaseous atmosphere can be blown away by both the planet's primordial energy from formation and bolometric heating from the star. If the formation energy is similar to or greater than the gravitational binding energy of the atmosphere, the cooling luminosity of the planet may blow off a H/He atmosphere entirely. For planets with masses between Earth and Neptune, the core temperatures can be 10,000–100,000 K, which can take up to several Gyr to cool (Ginzburg et al. 2016). Gupta & Schlichting (2020) found that the average size of sub-Neptunes decreases significantly with age, and commensurately that the relative occurrence rate of sub-Neptunes decreased by a factor of 2–2.5 over the timescale of the mass loss (see their Figure 10), which is consistent at the 1–2σ level with our results.

One obvious question is whether our finding, that the occurrence rate of hot sub-Neptunes is still very high at ∼700 Myr, before dropping by a factor of up to 5 at several Gyr, agrees with the putative timescales of these mass-loss mechanisms. Owen & Wu (2013, 2017) investigated the timescales of photoevaporation, finding that most atmospheric erosion happens within the first 100 Myr (see, e.g., Figure 2 of Owen & Wu 2017). This would imply that planets shrink in size somewhat sooner than our Praesepe+Hyades result but is not at odds with the younger TESS result. Core-powered mass loss is expected to occur over much longer timescales than photoevaporation, typically 0.5–2 Gyr (Ginzburg et al. 2016; Gupta & Schlichting 2019, 2020), which is much more consistent with our finding that at 0.7 Gyr these planets are still large. Gupta & Schlichting (2020) state: "We expect a drastic change in the planet size distribution between stars that are younger and older than the typical core-powered mass-loss timescale, which is of the order of a Gyr. Due to these long mass-loss timescales, we predict the transformation of sub-Neptunes into super-Earths to continue over Gyr timescales." Indeed, early examination of individual planets being detected in the Praesepe and Hyades K2 observations noted that they appeared to be larger, on average, than planets orbiting older stars of the same mass observed with Kepler (see, e.g., Figure 12 of Rizzuto et al. 2018), which could have been evidence for ongoing thermal contraction of their atmospheres, and/or could have hinted at a preference for core-powered mass loss (Gupta & Schlichting 2020). However, without understanding the completeness of the survey detecting those planets, this could just as easily have been a selection bias, where the noisier K2 light curves limited the detection to larger planets when compared with the more well-behaved Kepler light curves. Berger et al. (2020) and Sandoval et al. (2021) used isochronal ages to divide the Kepler planet sample into coarse age bins, and found that the relative frequency of sub-Neptunes to super-Earths decreases with increasing age, i.e., that either the number of sub-Neptunes was decreasing, the number of super-Earths was increasing, or both. Here we extend the Kepler ages to the young clusters observed by K2 and find that there is indeed a "drastic change" (a factor of ∼4–5) in the absolute frequency of sub-Neptunes between stars that are younger and older than ∼1 Gyr. This result, specific to the Praesepe and Hyades clusters studied here, is more consistent with the published predictions of the core-powered mass loss origin for the planet radius valley rather than a photoevaporation origin, or perhaps points to some missing process in the photoevaporation model that would act to slow down the expected mass loss rate. Indeed, recent work by Owen & Schlichting (2023) conversely theorizes that planets that have undergone core-powered mass loss can then transition to photoevaporation, depending on the initial envelope mass fraction; whichever the order, the bulk of the mass loss is indicated to happen after the first 0.5 Gyr.

Another factor that must be considered is the natal disk opacity, which plays an important role in both photoevaporation and core-powered mass loss processes. Planets that have formed from material that is richer in metals have the ability to maintain their envelopes for a longer period of time due to their increased heat capacity, enabling efficient cooling of their outer envelopes. Given that our cluster planet sample is more metal-rich than the [Fe/H] ∼ −0.2 of Kepler (Dong et al. 2014) and the [Fe/H] ∼ 0.0 of the field K2 sample, it is important to consider how this may affect our conclusions. According to Owen & Murray-Clay (2018), the envelope mass loss rate ($\dot{M}$) dictated by the photoevaporation model scales with the disk metallicity (Z) as $\dot{M}\propto {Z}^{-0.77}$, while the mass loss originating from the residual heat in the planet's core scales as $\dot{M}\propto {Z}^{-1}$ (Gupta & Schlichting 2020). This implies that the removal of sub-Neptune atmospheres will be delayed in a metal-rich environment, increasing the expected timescales of photoevaporation and core-power mass loss processes by 1.4× and 1.6×, respectively. ¹⁶ However, the surplus of hot sub-Neptunes at 600–800 Myr is more than 6× the fiducial timescale for photoevaporative mass loss (∼100 Myr). Therefore this potential delay, arising from the cluster's metallicity enhancement, cannot account for the abundance of hot sub-Neptunes identified in this study.

Finally, our sample covers a wide range of stellar masses, from 0.3 to 1.3 solar masses. The evolutionary timescales discussed above have been derived for an FGK sample of stars. The inclusion of low-mass stars in our sample may modify the time frame over which hot sub-Neptunes are converted to super-Earths. Under the action of the photoevaporation framework, small M dwarfs experience a heightened integrated XUV flux, accelerating and/or increasing the magnitude of the mass loss process (Rogers et al. 2021). Qualitatively, this would reduce the expected timescale needed to dissociate the planetary envelopes of sub-Neptunes, culminating in an overall occurrence reduction in the 600–800 Myr clusters. However, our results show that the opposite is true. In contrast, the core-powered mass loss evolutionary process relies on the star's total luminosity and therefore expects decelerated and/or reduced mass loss around less massive stars (Gupta & Schlichting 2020). The inclusion of M dwarfs should increase the expected mass loss timescale and leave a larger reservoir of short-period sub-Neptunes at 600–800 Myr, matching the observation of this survey. Overall, the complexities of stellar mass and composition unique to our sample are directionally consistent with a core-powered mass loss mechanism, providing further evidence in favor of this formation pathway.

An alternative explanation for the overabundance of hot sub-Neptunes in Praesepe and Hyades compared to the field K2 sample is that there could be systematic evolution of planetary system formation and architectures with different environmental/Galactic conditions over cosmic time. There are a number of known correlations between small, short-period planet occurrence rates and their host star properties—including properties that have evolved over the lifetime of the Galaxy, as the natal conditions for each generation of stars are impacted by previous generations. One such property is the host star metallicity—stars that are born today are typically formed in giant molecular clouds that have been enriched by the supernova remnants of earlier stars, in a more heavy-element-rich environment than in early Galactic times. This is reflected in the slightly supersolar average metallicities of Praesepe and Hyades, but was also explicitly fit as part of our parametric model, so presumably what we are seeing here is an additional effect. Both observations and models of galaxy formation in Milky Way-like galaxies suggest that stars born at earlier times are typically born in more strongly dynamically clustered environments (Lee & Hopkins 2020; Grudić et al. 2023), are much more strongly irradiated externally by elevated galactic star formation rates (Guszejnov et al. 2020; Tacconi et al. 2020), with likely larger ionizing cosmic ray fluxes (Hopkins et al. 2021) and systematically warmer natal gas clouds (Guszejnov et al. 2020). There is a long history of work on how this may or may not influence the low-mass end of the stellar initial mass function, but until recently little work on its possible effects on planet formation. Possibly, early and frequent interactions with nearby stars destabilize planetary systems to the extent that the overall occurrence rate, even for planets on the close-in orbits probed by our analysis, is reduced. Zink et al. (2023) recently discovered a strong negative correlation between the occurrence rate of 1.7–4 R_⊕, 1–10 day planets (the same planet parameter space probed here), and the maximum height above the Galactic plane that the host stars reach as they oscillate through the plane. This height can be used as a proxy for thick disk (older, more metal-poor, more likely to have formed in more clustered stellar environments) versus thin disk (younger, more metal-rich, less likely to have formed in more clustered environment) membership, and as Zink et al. (2023) also accounted for metallicity, points to an additional line of evidence that older stars have fewer small, short-period planets when compared to younger stars, although as here there was no systematic difference identified in the shape of the planet population. Although there has been less work on quantitative predictions for these primordial planet formation models compared to the atmosphere evolution models, nothing in this work disfavors such explanations as a category.