Hyades Member K2-136c: The Smallest Planet in an Open Cluster with a Precisely Measured Mass

K2-136 is a late-K dwarf (0.742 ± 0.039 M ⊙) in the Hyades open cluster with three known, transiting planets and an age of 650 ± 70 Myr. Analyzing K2 photometry, we found that planets K2-136b, c, and d have periods of 8.0, 17.3, and 25.6 days and radii of 1.014 ± 0.050 R ⊕, 3.00 ± 0.13 R ⊕, and 1.565 ± 0.077 R ⊕, respectively. We collected 93 radial velocity (RV) measurements with the High-Accuracy Radial-velocity Planet Searcher for the Northern hemisphere (HARPS-N) spectrograph (Telescopio Nazionale Galileo) and 22 RVs with the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) spectrograph (Very Large Telescope). Analyzing HARPS-N and ESPRESSO data jointly, we found that K2-136c induced a semi-amplitude of 5.49 ± 0.53 m s−1, corresponding to a mass of 18.1 ± 1.9 M ⊕. We also placed 95% upper mass limits on K2-136b and d of 4.3 and 3.0 M ⊕, respectively. Further, we analyzed Hubble Space Telescope and XMM-Newton observations to establish the planetary high-energy environment and investigate possible atmospheric loss. K2-136c is now the smallest planet to have a measured mass in an open cluster and one of the youngest planets ever with a mass measurement. K2-136c has ∼75% the radius of Neptune but is similar in mass, yielding a density of 3.69−0.56+0.67 g cm−3 (∼2–3 times denser than Neptune). Mass estimates for K2-136b (and possibly d) may be feasible with more RV observations, and insights into all three planets’ atmospheres through transmission spectroscopy would be challenging but potentially fruitful. This research and future mass measurements of young planets are critical for investigating the compositions and characteristics of small exoplanets at very early stages of their lives and providing insights into how exoplanets evolve with time.


INTRODUCTION
The timescales on which planets and planetary systems evolve are far longer than any feasible timescale of scientific observations.The only way to learn about how planets form and evolve is to collect snapshots at different stages of their development and assemble these snapshots into a cohesive framework.This is where open clusters prove particularly useful.Open clusters, close collections of young, recently formed stars, are excellent laboratories for studying the early lives of stars, because all of the stars in an open cluster, regardless of size, temperature, metallicity, or location, have a shared formation history, and therefore the ages of the stars can be very tightly constrained.This logic can also be applied to planets; if they form very quickly after the coalescence of their host star (Raymond & Morbidelli 2020), it is possible to determine the age of a planet orbiting an open cluster star, thereby capturing one of the early snapshots required to assemble the framework of a planet's evolution.
In this paper we characterize K2-136c, a sub-Neptune planet in the Hyades open cluster.Orbiting a late K dwarf, this planet is one of the three known, transiting planets in the system.The system was originally observed in K2 Campaign 13 for 80 days (2017 March 8 -2017 May 27) and was proposed for observation by seven guest observer teams: GO13008, GO13018, GO13023, GO13049, GO13064, GO13077, and GO13090.All three planets were originally discovered by Mann et al. (2018) (hereafter M18) and Ciardi et al. (2018) (a parallel analysis published simultaneously).Shortly thereafter a subsequent analysis was completed by Livingston et al. (2018).All three papers are in broad agreement regarding stellar and planetary parameters, but M18 established the tightest constraints on orbital period for all three planets.
As for the stellar age, there are a number of estimates available.Perryman et al. (1998) found the Hyades open cluster to be 625 ± 50 Myr.Gossage et al. (2018) found an age of ∼ 680 Myr while Brandt & Huang (2015) determined a slightly older age of 750 ± 100 Myr.The age we use throughout this paper comes from Martín et al. (2018), who determined the Hyades to be 650 ± 70 Myr old.We thus assume that K2-136 and the three orbiting planets share that approximate age.We chose this age because it is a relatively recent result, it compares and combines results using both old (Burrows et al. 1997) and new (Baraffe et al. 2015) standard evolutionary models, and it also agrees broadly with other, previous estimates.The young age of the system was our primary reason for pursuing K2-136c as a target: there are very few young, small planets with mass measurements.According to the NASA Exoplanet Archive (accessed 2023 Mar 12; NASA Exoplanet Science Institute 2020), there are only 13 confirmed exoplanets with R p < 4 R ⊕ , a host star age < 1 Gyr, and a mass measurement (not an upper limit): HD 18599b (Desidera et al. 2022), HD 73583b and c (Bar-ragán et al. 2022); K2-25b (Stefansson et al. 2020);L 98-59b, c, and d (Demangeon et al. 2021); Kepler-411b and Kepler-411d (Sun et al. 2019); Kepler-462b (Masuda & Tamayo 2020); Kepler-289b and Kepler-289d (Schmitt et al. 2014); and K2-100b (Barragán et al. 2019).Of these, only the Kepler-411, K2-100, HD 73583, K2-25, and HD 18599 systems have an age constraint tighter than 50% (Sun et al. 2019;Barragán et al. 2019Barragán et al. , 2022;;Stefansson et al. 2020;Desidera et al. 2022).
We analyzed photometry of the K2-136 system in order to measure the radii, ephemerides, and other transit parameters of each planet.We also collected spectra of the K2-136 system and measured radial velocities (RVs) as well as stellar activity indices.Then, by modeling these RVs (following Rajpaul et al. 2015), we determined the mass of K2-136c and placed upper limits on the masses of the other two planets.We used this system to investigate the nature, environment, and evolution of young, small exoplanets.
This paper is organized as follows.In Section 2 we discuss our observations.Then we detail our method of stellar characterization in Section 3. Next, in Section 4 we describe our RV and photometry models, data analysis, model comparison, and parameter estimation.In Section 5 we present and discuss our results.Finally, we summarize and conclude in Section 6.

OBSERVATIONS
2.1.K2 Photometric observations of the K2-136 system were collected with the Kepler spacecraft (Borucki et al. 2008) through the K2 mission during Campaign 13 (2017 Mar 08 to 2017 May 27).K2 collected long-cadence observations of this system every 29.4 minutes.

TESS
Photometric observations of the K2-136 system were also collected with the TESS spacecraft (Ricker et al. 2015) during Sector 43 (2021 Sep 16 to 2021 Oct 10) and Sector 44 (2021 Oct 12 to 2021 Nov 06).TESS collected long-cadence full frame image observations of this system every 10 minutes in Sector 43 and shortcadence observations every 20 seconds in Sector 44.

HARPS-N
We collected 93 RV observations using the HARPS-N spectrograph (Cosentino et al. 2012, Cosentino et al. 2014) on the Telescopio Nazionale Galileo (TNG).The first 88 spectra were collected between 2018 August 11 and 2019 February 7 (programs A37TAC 24 and A38TAC 27, PI: Mayo), and the final 5 spectra were collected between 2020 September 18 and 2020 October 31 by the HARPS-N Guaranteed Time Observation program.RVs and additional stellar activity indices were extracted using a K6 stellar mask and version 2.2.8 of the Data Reduction Software (DRS) adapted from the ESPRESSO pipeline.Spectra had an average exposure time of 1776.5 seconds and the average SNR in the order around 550 nm was 51.1.The RV standard deviation was 6.9 m s −1 and the RV median uncertainty was 1.6 m s −1 .Stellar activity indices also extracted and reported in this paper include the cross-correlation function (CCF) bisector span inverse slope (hereafter BIS), the CCF full width at half maximum (FWHM), and S HK (which measures chromospheric activity via core emission in the Ca II H and K absorption lines).The observation dates, velocities, and activity indices are provided in Table 2.

ESPRESSO
We collected 22 RV observations using the ESPRESSO spectrograph (Pepe et al. 2021) on the Very Large Telescope (VLT) between 2019 November 1 and 2020 February 27 (program 0104.C-0837(A), PI: Malavolta).RVs and additional stellar activity indices were extracted using a K6 stellar mask and the same pipeline as the HARPS-N observations (DRS version 2.2.8).Typical exposure time for spectra was 1800 seconds and the average SNR at Order 111 (central wavelength = 551nm) was 79.9.The RV standard deviation was 7.8 m s −1 and the RV median uncertainty was 0.70 m s −1 .These observations and indices are also provided in Table 2.

Hubble Space Telescope
Near-ultraviolet (NUV) observations of K2-136 were taken as part of a broader Hubble Space Telescope (HST) program observing the Hyades (GO-15091, PI: Agüeros).The target was exposed for 1166.88 seconds on 2019 September 13 using the photon-counting Cosmic Origins Spectrograph (COS; Green et al. 2012) in the G230L filter and had no data quality flags.
After initial data reduction through the CALCOS pipeline version 3.3.10,we additionally confirmed that the star was not flaring during observations by integrating the background-subtracted flux by wavelength over 1 and 10 second time intervals in the time-tagged data.No flares above 3σ were identified.
2.6.XMM-Newton K2-136 was the target of an XMM-Newton (XMM) 43 ksec observation on 2018 September 11 (Obs.ID: 0824850201, PI: Wheatley).The observation was processed using the standard Pipeline Processing System (PPS version 17.56 20190403 1200; Pipeline sequence ID: 147121).The source detection corresponding to K2-136 was detected by both the pn and MOS cameras, for a total of 800 source counts in the 0.2-12.0keV energy band.The X-ray source has data quality flag SUM FLAG=0 (i.e., good quality).No variability or pileup were detected for this X-ray source.

STELLAR CHARACTERIZATION
In order to characterize the star, we started by combining all of our collected HARPS-N spectra (from 2018-2019) into a single, stacked spectrum with S/N ∼ 300 (based on signal divided by scatter on continuum segments near 6000 Å; see Section 3.1 of Mortier et al. 2013 for more details).Then we ran the ARESv2 package (Sousa et al. 2015) to obtain equivalent widths for a standard set of neutral and ionised iron lines (Sousa et al. 2011).We refer to Mortier et al. (2013), Sousa (2014), andSousa et al. (2015) for our choice of typical model parameters.Afterward, we calculated stellar parameters using MOOG 1 (Sneden 1973) with ATLAS plane-parallel model atmospheres (Kurucz 1993) assuming local thermodynamic equilibrium.A downhill simplex minimization procedure (Press et al. 1992) was used to determine the stellar photospheric parameters (see e.g.Mortier et al. 2013, and references therein).We determined that the stellar temperature was less than 5200 K, so we reran the minimization procedure with a sublist of lines designed for cooler stars (Tsantaki et al. 2013); we also constrained our line list to those with equivalent widths between 5 and 150 milliAngstroms (m Å), removing 5 lines above 150 m Å and 1 line below 5 m Å (lines within this range tend to be sufficiently strong and well-described by a Gaussian).Finally, we corrected for log g and re-scaled errors following Torres et al. (2012), Mortier et al. (2014), andSousa et al. (2011).The resulting effective temperature, surface gravity, microturbulence, and metallicity are reported in Table 1.
Then we determined the same stellar parameters from the same spectra with a different, independent tool: the Stellar Parameter Classification tool (SPC; Buchhave et al. 2012).SPC interpolates across a synthetic spectrum library from Kurucz (1992) to find the best fit and uncertainties on an input spectrum.In addition to the stellar parameters calculated from ARES+MOOG, this tool also estimated rotational velocity.All atmospheric stellar parameters from ARES+MOOG and SPC were in good agreement (within 1σ).Like ARES+MOOG, all SPC parameter estimates can be found in Table 1.
We then took our estimated effective temperature and metallicity from ARES+MOOG and SPC, the Gaia Data Release 3 (DR3) parallax (Gaia Collaboration et al. 2016, 2020;Gaia Collaboration & Vallenari 2022), and numerous photometric magnitudes (B, V, J, H, K, W1, W2, and W3) and input them into the isochrones Python package (Morton 2015).This package used two different sets of isochrones: Dartmouth (Dotter et al. 2008) and Modules for Experiments in Stellar Astrophysics (MESA) Isochrones and Stellar Tracks (MIST; Choi et al. 2016;Dotter 2016).Comparing two standard, independent models is useful for mitigating systematic errors and revealing discrepancies or issues in the resulting parameter estimates.We used MultiNest (Feroz et al. 2009(Feroz et al. , 2013) ) for parameter estimation, assuming 600 live points and otherwise standard MultiNest settings: importance nested sampling mode, multimodal mode, constant efficiency mode disabled, evidence tolerance = 0.5, and sampling efficiency = 0.8.As stated earlier, K2-136 is a member of the Hyades and therefore has a very tight age constraint of 650 ± 70 Myr (Martín et al. 2018).We applied a much broader age prior of 475 Myr -775 Myr (a 3σ range on the 625 ± 50 Hyades age estimate from Perryman et al. 1998), which was more than sufficient to achieve convergence.This yielded posterior distributions from both input atmospheric parameter sets (ARES+MOOG and SPC) as well as both isochrone sets (Dartmouth and MIST), for a total of four sets of posterior distributions (based on all combinations of input parameters and isochrones).
The posteriors were then combined together (i.e. the posterior samples were appended together) to yield a single posterior distribution for each parameter.Lastly, systematic uncertainties determined by Tayar et al. (2020) were added in quadrature to the combined posteriors to   3.1.Stellar Rotation Period One parameter of special interest is the stellar rotation period, which we include as a parameter in our RV model (see Section 4.5).M18 conducted a Lomb-Scargle periodogram on the K2 light curve and reported a rotation period of 15.04 ± 1.01 days.Ciardi et al. (2018) analyzed the same light curve and found a rotation period of 15.2±0.2days through a Lomb-Scargle periodogram and 13.8±1.0days through an autocorrelation function.Livingston et al. (2018) conducted a Gaussian process (GP) regression, a Lomb-Scargle periodogram, and an autocorrelation function on the light curve and found a corresponding rotation period of 13.5 +0.7 −0.4 d, 15.1 +1.3 −1.2 d, and 13.6 +2.2 −1.5 d, respectively.Note: Given an offset and uncertainties in the photometric data set, a generalized Lomb-Scargle periodogram would be preferred (Zechmeister & Kürster 2009); however, it is not clear from the referenced papers whether this generalized method was used or just a basic Lomb-Scargle periodogram (Scargle 1982).
Notably, the estimates via a Lomb-Scargle periodogram are longer than estimates with other methods.All results are broadly consistent with our findings from our full model results (13.37 +0.13  −0.17 days; see Table 4), except the 15.2±0.2day result from the Lomb-Scargle analysis by Ciardi et al. (2018).They also have the smallest uncertainties of any rotation period estimate, so it is possible that their value is reasonable but the uncertainties are overly optimistic.
A possible explanation of this discrepancy could be differential rotation.Regardless of activity level, starspots, plage, and other activity may be more prominent at different stellar latitudes when the K2 photometry and our HARPS-N spectroscopy were conducted.This hypothesis is also mentioned by Ciardi et al. (2018) to explain a larger than expected vsin i.Following Barnes et al. (2005) and Kitchatinov & Olemskoy (2012), they estimate that the equatorial rotation period of K2-136 could be faster than higher latitudes by ∼ 1 day.Then again, Aigrain et al. (2015) found that claims of differential rotation should be treated with caution even for long baselines of photometry.We may simply be seeing different starspots at different longitudes creating phase modulation, combined with greater or fewer numbers of starspots leading to better or worse constraints on rotation period.

Binarity of K2-136
One of the planet discovery papers, Ciardi et al. (2018), reported a binary companion to K2-136.In addition to their K2 photometric analysis, they collected spectra from the SpeX spectrograph (Rayner et al. 2003(Rayner et al. , 2004) ) at the 3-m NASA Infrared Telescope Facility and the HIRES spectrograph (Vogt et al. 1994) at the Keck I telescope, as well as AO observations with the NIRC2 instrument at the Keck II telescope and the P3K AO system and PHARO camera (Hayward et al. 2001) on the 200" Hale Telescope at Palomar Observatory.The AO observations at both facilities detected an M7/8V star separated from the primary star by ∼ 0.7", corresponding to a projected separation of ∼ 40 AU; the spectroscopic observations did not detect this companion, and no further companions were found by any of the above observations.(Notably, this angular separation is more than enough for HST to resolve; see Sections 4.9 and 2.5.)Gaia DR3 did not detect the binary companion, leaving the issue of boundedness unresolved.However, Ciardi et al. (2018) compared the current position of K2-136 against observations from the 1950 Palomar Observatory Sky Survey (POSS I) and noted that the star had moved 6" in the intervening time with no evidence of background stars.These POSS I observations show that the stellar companion is likely bound.
Further, Gaia DR3 reported K2-136 to have an astrometric excess noise of 96 µas and a Renormalised Unit Weight Error (RUWE) of 1.23, a mild departure from a good single-star model.At the separation and brightness of the companion, this excess variability in the astrometry is unlikely to be due to pollution from its light contribution: at ∼ 0.7" separation a companion can be detected only with a G magnitude difference of 2 (Gaia Collaboration et al. 2021).There is therefore a mild indication of astrometric variability due to unmodeled orbital motion.Using the formalism of Torres (1999), at the distance of the system, given its angular separation, and for the mass range of an M7/8 star, the median astrometric acceleration is expected to be ∼ 25 µas yr −2 , with maximum value close to 40 µas yr −2 , indicating that in addition to simple astrometric noise, the bulk of the astrometric variability could be caused by the detection of the acceleration due to the companion.
It is worth considering whether flux from the companion could bias the measured RVs of the primary K dwarf.Both the HARPS-N and ESPRESSO band passes are approximately 380nm -690nm and centered on the V band.According to Ciardi et al. (2018), the M dwarf companion is at least 10 magnitudes fainter than the primary in the V-band.We can use ∆m = 10 as a worst-case scenario and similarly assume the companion star was well-centered on the fiber for all observations (because the companion and primary are separated by 0.7", the companion would not be well-centered and the actual flux contamination from the companion would be less).Cunha et al. (2013) explored the RV impact of flux contamination from a stellar companion: for a K5 dwarf and an M dwarf (M3 or later) with ∆m = 10, the maximum impact on RVs is < 10 cm s −1 , and therefore negligible for our level of RV precision.Ciardi et al. (2018) explored whether the transit signals may originate from the M dwarf.They found that in order to match the observed transit depth of K2-136c, the M dwarf would have to be a binary system itself that exhibits significant and detectable secondary eclipses, which have not been observed.Further, the transit duration of K2-136c is inconsistent with a transit of an M dwarf.Finally, K2-136c has already been validated by Ciardi et al. (2018) and all three planets have been independently validated by M18 and Livingston et al. (2018).Therefore, it is very unlikely that the planets are false positive signals or planetary signals from the M dwarf.
However, in the Kepler band pass, the companion M dwarf is 6.5 magnitudes fainter, which leads to a very small dilution effect on the planet transit depths.Following Ciardi et al. (2015), we include this effect for the sake of robustness in our final planet radius estimates (which enlarges each planet by ∼ 0.13%).
This binary companion may also cause a long-term RV trend, which we discuss further in Section 4.5 and test in Section 4.10.

DATA ANALYSIS
We analyzed the RV data in conjunction with a number of other common stellar activity indices that are calculated with the ESPRESSO DRS 2.2.8 pipeline, specifically the CCF BIS, CCF FWHM, and S HK .In this section, we first explore the data set by generating a periodogram and conducting a correlation analysis between the RVs and other data types.Then we discuss the transit and RV components of our model, the param-eter estimation process, and how we compare our models.Finally, we conduct tests on our results and discuss the implications of a binary companion in the system.

Periodogram Analysis
In order to investigate periodic signals in our data (planetary or otherwise), we created Generalized Lomb-Scargle periodograms (Scargle 1982;Zechmeister & Kürster 2009) of our RVs and our stellar activity indices.As a point of reference, we also included the window function of our data (built from constant, non-zero values at each of the timestamps of our observations).It is used to determine the regular patterns in the periodograms due to the sampling and gaps in the time series.These periodograms are presented in Fig. 1.
To ascertain the robustness of any apparent signals, we also estimated each periodogram's False Alarm Probability (FAP), the likelihood that an apparent signal of a given strength will be detected when no underlying signal is actually present.The FAP was estimated with the bootstrap method: sampling the observations randomly with replacement while maintaining the same timestamps.We repeat this process 100000 times, each time constructing a periodogram and determining the maximum peak.This reveals how often a given signal strength will appear due only to noise, from which the FAP is calculated.
The strongest RV signals in our combined periodogram (both HARPS-N and ESPRESSO) are at the orbital period of K2-136c, the rotation period of the star, and near 0.017 d −1 (i.e.half the length of the ESPRESSO data baseline of 117.7 days), although none are significant at the 1% level.Among the stellar activity periodograms, the strongest signals are the rotation period signal in the ESPRESSO FWHM power and a long-period signal in the HARPS-N FWHM power which can be attributed to the window function.Also, the strongest peak in the (combined data) window function periodogram above 0.01 days −1 (near 0.03 days −1 ) does not correspond to any significant signals or aliases in any of the four data types.As for low-frequency signals < 0.01 days −1 , we discuss possible long-term trends in Sections 3.2 and 4.5.All of this indicates that K2-136c has a more detectable RV signal than the other two planets, as expected given its size.

Correlation Analysis
We also examined the relationship between RVs and our stellar activity indices.Because RVs are measured from small shifts in spectral absorption lines, and stellar activity can change the shape of absorption lines, stellar activity can significantly affect RV observations (e.g.Queloz et al. 2001;Haywood 2015;Rajpaul et al. 2015).Scatter plots between RVs and all three activity indices are presented in Fig. 2. At least in the case of S HK and FWHM, there are notable correlations with RVs according to the p-values for the Spearman correlation coefficient (which captures nonlinear, monotonic correlations and uses the same -1 to 1 range as the Pearson correlation coefficient).According to the Spearman coefficient, there is a correlation of 0.26 between RV and S HK , a negative correlation of −0.18 between RV and BIS, and the strongest correlation of 0.39 between RV and FWHM.4. Blue data points correspond to HARPS-N observations, orange data points to ESPRESSO.In the top-left corner of each subplot is the Spearman correlation coefficient, which can capture nonlinear, monotonic correlations.(Coefficients were calculated with data values but not uncertainties.) For a data set of this size, these coefficients correspond to p-values of 0.004 for RV and S HK , 0.052 for RV and BIS, and 0.001 for RV and FWHM.Therefore, for RV and S HK and especially RV and FWHM, there appears to be a statistically significant correlation.In other words, there is good reason to believe that this data set includes correlated and structured stellar activity.In fact, the Spearman correlation coefficient is likely an underestimate of the correlation between RVs and stellar activity indices, since there can be a phase shift in the amplitude variation from one data type to another (Santos et al. 2014;Collier Cameron et al. 2019).

K2 Transit Photometry
We cleaned and flattened the photometric data from K2 using the exact same procedure as was used originally in M18.Their procedure follows the self-flatfielding (SFF) method developed in Vanderburg & Johnson (2014) to perform a rough removal of instrumental variability followed by a simultaneous fit to a model consisting of Mandel & Agol (2002) transit shapes for the three planets, a basis spline in time to describe the stellar variability, and splines in Kepler's roll angle to describe the systematic photometric errors introduced by the spacecraft's unstable pointing (Vanderburg et al. 2016).After performing the fit, they removed the best-fit systematics and stellar variability components, isolating the transits for further analysis.The interested reader should refer to M18 for additional detail of the full procedure.
We modeled the flattened and cleaned M18 light curve with the BATMAN Python package (Kreidberg 2015), based on the Mandel & Agol (2002) transit model.Our model included a baseline offset parameter and white noise parameter for our K2 Campaign 13 photometry as well as two quadratic limb-darkening parameters (parameterized using Kipping 2013).Each planet was modeled with five parameters: the transit time, orbital period, planet radius relative to stellar radius, transit duration, and impact parameter.All parameters were modeled with either uniform, Gaussian, Jeffreys, or modified Jeffreys priors (Gregory 2007).Only the photometric white noise parameters used modified Jeffreys priors, with a knee located at the mean of the photometric flux uncertainty for that particular campaign or sector.All priors are listed in Table 3.The raw and flattened data can be seen in Fig. 3.
We also applied a Gaussian prior on stellar density by comparing the spectroscopically derived stellar density to the stellar density found via the following equation (Seager & Mallén-Ornelas 2003;Sozzetti et al. 2007): where orbital period (P ) and the semi-major axis (a/R * ) are derived directly from the light curve model.
All together, our full transit model includes 15 planetary parameters (5 per planet: time of transit, orbital period, ratio of planet radius to stellar radius, transit duration, and impact parameter), 2 quadratic limb darkening parameters, 1 photometric noise parameter, and 1 photometric baseline offset parameter for a total of 19 parameters.
The only parameters in common between our transit model and our RV model (as explained in further detail below) are transit times and orbital periods.

TESS Transit Photometry
No pre-processed light curves were available for the TESS Sector 43 observation of K2-136, so we extracted the photometry from the full frame image (FFI) pixel level.Following Vanderburg et al. (2019), we constructed 20 different apertures (10 circular, 10 shaped like the TESS point spread function) and selected the one that best minimized photometric scatter.As for TESS Sector 44 observations, we used the simple aperture photometry (SAP) light curve produced by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016).Light curves from both sectors were then flattened in the same way: a basis spline fit was performed iteratively on the photometry (with breakpoints every 0.3 days in order to adequately model stellar variability) and 3σ outliers were removed until convergence (this too, aside from the breakpoint length, follows Vanderburg et al. 2019).Finally, we conducted a simultaneous fit of the low-frequency variability and the transits in order to determine the best-fit low-frequency variability.
TESS photometry is not incorporated into our final photometric model, although we did run exploratory joint transit models on K2 and TESS photometry simultaneously.The transit signals of the two smaller planets, K2-136b and d were too small to reliably detect in the TESS photometry: individual transits were indistinguishable in depth and quality from temporally adjacent stellar activity.However, the transits of K2-136c were easily identifiable individually in TESS photometry, so we ran a joint transit model on all K2 photometry and TESS photometry to explore the resulting improvements in the parameters of K2-136c.This joint transit model included all parameters listed in Section 4.3 as well as additional baseline offset and white noise parameters for the two TESS Sectors (43 and 44) and two additional quadratic limb-darkening parameters for TESS photometry for six additional parameters total.The fit resulted in consistent values for all planet and system parameters as well as a dramatically more precise ephemeris for K2-136c: P c = 17.307081 +0.000014 −0.000013 days and t 0,c = 8678.07179+0.00067 −0.00063 (BJD-2450000).For comparison, this period and transit time have uncertainties that are both approximately 15x tighter than those resulting from K2 transit modeling alone (see Table 3).We report these values here to minimize ephemeris drift and facilitate planning of future transit observations of K2-136c.

RV Model
We modeled the RV signal of the orbiting planets and the stellar activity simultaneously.We assumed non-interacting planets with Keplerian orbits.We used RadVel (Fulton et al. 2018) to model the RV signal from each planet with 5 parameters: reference epoch, orbital period, RV semi-amplitude, eccentricity, and argument of periastron.The latter two parameters, eccentricity and longitude of periastron, were parameterized as √ e cos w and √ e sin w.As explained in Eastman et al. (2013), this reparameterization avoids a boundary condition at zero eccentricity that may lead to eccentricity estimates that are systematically biased upward.We conducted trial simulations with circular versus eccentric orbits for all three planets and found excellent agreement in all parameters (less than 1σ).We opted to keep eccentricity and argument of the periastron as parameters in order to constrain or place upper limits on each planet's eccentricity.Additionally, we prevented system configurations that would lead to orbit crossings of any two planets, as well as overlaps of planetary Hill spheres.For each planet's reference epoch and orbital period, we applied a Gaussian prior based on the transit parameters deter-mined from M18.We analyzed the K2 transit photometry with and without TESS photometry and verified these M18 values (see Table 3).
We also applied additional prior limits on the eccentricities of K2-136b and K2-136d.Based on preliminary modeling of all three planets with uninformed prior eccentricity constraints (and everything else identical to our final model), we found that we could determine the eccentricity of K2-136c but not its siblings.Thus, we decided to set eccentricity constraints by using the Stability of Planetary Orbital Configurations Klassifier (SPOCK; Tamayo et al. 2020), an N-body simulator that employs machine learning to improve performance.We input into SPOCK our stellar mass posterior and our orbital period posteriors for all three planets (from our preliminary simulations).We also input the mass, eccentricity, and argument of the periastron posteriors for K2-136c (again, from our early simulations).Lastly, we input uniform distributions for mass, eccentricity, and argument of the periastron for K2-136b and K2-136d.Planet mass ranged from 0 up to that of approximately 100% iron planet composition (3 M ⊕ for K2-136b and 30 M ⊕ for K2-136d; Fortney et al. 2007).Eccentricity ranged from 0 to 0.6 and argument of the periastron ranged from 0 to 360 degrees.We took the subset of our sample that had a > 90% chance of surviving for > 10 9 orbits (of the innermost planet) and then determined the 3σ upper limit on the eccentricities of K2-136b and K2-136d for that subset.We then used those values, e b,max = 0.35 and e d,max = 0.37, as eccentricity upper limit priors in all subsequent simulations.
Given the M dwarf companion to our host star (Ciardi et al. 2018), we wanted to include the potential for an RV trend caused by this companion.For further discussion of binarity and linear trends, see Section 3.2.Our planetary RV signals and the RV trend therefore take the following form: where K is the induced RV semi-amplitude, ω is the argument of the periastron, f is the true anomaly, e is the eccentricity, m is the slope of the RV trend, t is the observation time, E is the eccentric anomaly, M is the mean anomaly, n is the mean motion, τ is the time of periastron passage (as calculated from transit time, orbital period, eccentricity, and argument of the periastron), and P is the orbital period.
The stellar activity was handled via a GP on RVs.We first used a simultaneous model of stellar activity on four different data types: RV, FWHM (a measure of the width of absorption lines), BIS (a measure of line asymmetry), and S HK (an estimate of chromospheric magnetic activity via emission in the cores of the Ca II H & K lines).However, we found that this approach forced the model to include unreasonable amounts of white noise into each data type via a white noise jitter parameter included in our model.Further, it did not lead to a notable improvement in our final parameter constraints.We conducted numerous tests exploring the excess white noise preferred by the model, including modeling different instruments, different numbers of planets, altering the data reduction pipeline (e.g.trying the ZLSD pipeline; Lienhard et al. 2022), and creating synthetic data sets to model against.The most reasonable hypothesis we could find is that for the K2-136, our data (especially for ESPRESSO) is of a sufficiently high quality, with small enough uncertainties, that the model we used from Rajpaul et al. (2015) to relate the stellar activity indices to each other and to the RVs was not complex enough to account for the correlated structure of the stellar activity of K2-136.
Our RV data and the final RV model fit can be seen in Fig. 4. Despite only modeling RVs without any stellar activity indices, a GP is still robust and allows us to separate planet-induced RVs from stellar activity to the extent that disentanglement is possible.We followed the method laid out in Rajpaul et al. (2015), but we constrained the model only to the portions relevant to RVs rather than additional stellar activity indices.They dictate that the RVs are related to the stellar activity as follows: In this equation, G(t) corresponds to the underlying stellar activity GP, while V c and V r correspond to the RV amplitudes of the convective blueshift and rotation modulation effect, respectively.It is important to include both rotation modulation and convective blueshift for data types impacted by both, as one phenomenon may play a larger role than the other depending on the data type and star.In fact, our final results (see Table 4) show that for K2-136, rotation modulation has an outsized effect on RVs compared to convective blueshift; however, we retain both terms in our model since V c is still inconsistent with zero.We follow Rajpaul et al. (2015) further by using a quasi-periodic kernel to establish the covariance matrix of our GP.A quasi-periodic kernel is a good choice to capture the stellar variability of a star because quasi-periodicity describes well the variability exhibited by a rotating star with starspots that come and go.Following Giles et al. (2017), we find a predicted starspot lifetime for K2-136 of 38 +20 −13 days (versus a rotation period of 13.37 days); this estimate is consistent with our GP evolution time-scale result (48 +19 −11 days, see Table 4).In other words, stellar activity during the time frame of a single rotation period is likely to look similar (though not identical) to stellar activity during the previous and subsequent stellar rotation periods, since the starspots for K2-136 are likely to be longer lived than the stellar rotation period.Thus, the quasi-periodic kernel is defined as follows: where h is the GP amplitude (which is folded into the amplitude parameters described above), P * is the stellar stellar rotation period, λ e is a decay timescale proportional to the starspot lifetime, and λ p is a smoothness parameter that captures the level of variability within a single rotation period.t i and t j are any two times between which the covariance is being calculated; for a given time series of N observations, all N 2 combinations of time pairs create the N xN covariance matrix.This covariance matrix (plus a mean model) is a normal multivariate distribution; G(t) can be explored by sampling from this distribution.We refer the reader to Rajpaul et al. (2015) for a more detailed description of this method, as well as Mayo et al. (2019) for an application of this method to the sub-Neptune Kepler-538b.
Finally, we include a jitter parameter (added in quadrature to RV uncertainties) and a baseline offset parameter for each instrument (HARPS-N and ESPRESSO).All together, our full model includes 21 planetary parameters (5 per planet: time of transit, orbital period, RV semi-amplitude, √ e sin ω, and √ e cos ω), 5 GP parameters (2 GP amplitudes corresponding to V c and V r from equation 6 as well as P * , λ p , and λ e from equation 7), 1 linear trend parameter, 2 RV noise parameters (1 per instrument) and 2 RV baseline offset parameters (1 per instrument), for a total of 25 parameters.
We used a uniform prior for the RV semi-amplitude for each planet.However, we conducted a trial simulation to compare a uniform versus log uniform prior for RV semi-amplitude and found no discernible difference in our results (all parameters agreed to within 1σ).

Parameter estimation
We conducted parameter estimation of our model with the observed data using the Bayesian inference tool MultiNest (Feroz et al. 2009(Feroz et al. , 2013)).We set MultiNest to constant efficiency mode, importance nested sampling mode, and multimodal mode.We used a sampling efficiency of 0.01, 1000 live points, and an evidence tolerance of 0.1.Constant efficiency is typically off, but we turned it on since it allows for better exploration of parameter space in higher-dimensional models such as our own.Further, sampling efficiency is usually set to 0.8 and the number of live points is usually set to 400.Decreasing sampling efficiency and increasing the number of live points leads to more complete coverage of parameter space, at the cost of a typically longer simulation convergence time.Finally, the evidence tolerance is usually set to 0.8; reducing the evidence tolerance causes the simulation to run longer but increases confidence that the simulation has fully converged.In other words, the standard MultiNest settings would likely lead to reliable results, but our choice of settings increases the trustworthiness of our parameter estimation and model evidence results.Unif(0, 1) Kepler/K2 quadratic limb-darkening parameter q 2,Kepler -0.56 +0.24 −0.18
One of the strengths of MultiNest is that it automatically calculates the Bayesian evidence of the selected model, making model comparison very easy.We compared the model evidences of eight different models (RV only, no photometry), based on all possible combinations of planets b, c, and d.The results of our model comparisons are listed in Table 5.We find that the most preferred model is the one that contains only planet c, and not planets b or d.In fact, using the Bayes factor interpretation of Kass & Raftery (1995), we find that almost every other combinations of planets can be decisively ruled out (i.e. the Bayes factor of the planet c model to any other model is > 100).The only excep-tions are the model with planets b and c and the model with planets c and d, which are only strongly disfavored.This tells us that neither K2-136b nor K2-136d are unambiguously detected, so including either in the model quickly worsens the model evidence.However, it is notable that the model with planets b and c is better than the model with planets c and d, and is nearly in the more likely "Disfavored" category rather than "Strongly disfavored".This makes sense, as K2-136b (unlike K2-136d) has a non-zero peak in its semi-amplitude posterior distribution (see Fig. 5).In fact, although only upper limits are reported for the mass of K2-136b in Table 4, the median mass is actually non-zero at the 2.0σ level (M b = Unif(0, 100) GP RV rotation modulation amplitude Vr m s −1 22.0 +12.9 −8.0 Unif(0, 100) GP stellar rotation period P * day 13.37 +0.13   −0.17  2013).e K2-136b and K2-136d had additional eccentricity prior upper limits of 0.35 and 0.37, respectively; see Section 4.5.
Our results tell us that we can be confident that K2-136c has been detected in our observations.In contrast, the RV signals of K2-136b and K2-136d fall below the threshold of detection, at least given the quantity and quality of our specific data set.Continued radial velocity monitoring, particularly with high precision instruments and facilities may be able to measure their masses, especially K2-136b (which appears to already be near the threshold of detection with our current data set).
Although the model with only K2-136c is the most favored, we still use a model with all three planets as our canonical model for parameter estimation for three reasons.First, we already know from transit photometry that planets b and d exist.Accordingly, the goal of the model comparison exercise described above is not to question the existence of these planets but to examine whether the RV signals from each planet can be detected  in our data set.By adopting the three-planet model, we incorporate the uncertainties introduced by the unknown masses and eccentricities of planets b and d.Second, using the three-planet model allows us to determine upper limits for the masses of planets b and d, which is useful for constraining planet compositions and providing guidance for any future attempts to constrain the mass of either planet.Third, both models agree very closely: all parameters are consistent at 1σ or less, and all uncertainties on parameters are similar in scale.For the RV semi-amplitude of K2-136c, our key parameter of interest, our canonical model returned K c = 5.49 +0.54  −0.52 m s −1 while the one-planet model returned K c = 5.17 +0.56  −0.51 m s −1 .

Model Reliability Tests
In order to rigorously assess the accuracy of our results, we conducted tests to analyze different components of our model.Specifically, we removed the GP portion of the model, and we also injected and recovered synthetic planets into the system to compare input and output RV semi-amplitudes.
For our test models we chose not to include planets b and d, as well as photometry, and we then compared against the model with only K2-136c (hereafter referred to as the "one-planet reference model") rather than the three-planet model; this is despite already selecting the three-planet model as the canonical model to report our results (see Section 4.7).This was done for a few reasons.First, the one-planet model is the preferred model according to the Bayesian evidences, so it is a reasonable point of comparison.Second, as stated earlier, the one-planet reference model and the three-planet model agree very closely.Therefore, any test model parameters found to be highly consistent with the one-planet reference model results will also be highly consistent with the three-planet model results.Third, as a practical matter, including only K2-136c in our test models significantly reduced computational complexity, allowing us to test a wider variety of models.

No GP
GPs are very versatile and can fit highly variable and correlated signals.Therefore, it is reasonable to ask whether a GP may, in the process of modeling stellar activity, "steal" part of the RV signal from a planet due to overfitting of the data.In order to address these concerns, we ran a model without a GP, and no alternative method to handle stellar activity.We found the resulting parameters were broadly consistent.The noise parameters for each data type in the no-GP model were notably larger, but that is to be expected given no mitigation of the stellar activity.All other parameters agreed with the one-planet reference model parameters to within 1σ; for the RV semi-amplitude of K2-136c, we determined a value of K c = 5.92 +0.89  −0.91 m s −1 (compared to K c = 5.17 +0.56  −0.51 m s −1 for the one-planet reference model).Finally, it is worth noting we also found that ∆ log 10 (evidence) = −18.1 compared to the one-planet reference model, decisively ruling out the no-GP model (Kass & Raftery 1995).In other words, a GP accounts for the stellar activity satisfactorily, whereas ignoring stellar activity is clearly insufficient.

Synthetic Planet Injections
We also conducted planet injection tests to determine how robustly we could recover the injected signals.Accurate recovery of such signals builds confidence in the accuracy of the RV signal recovered for K2-136c as well as the upper limits placed on K2-136b and K2-136d.
We ran four separate tests in which we injected a 5.5 m s −1 RV signal of a planet on a circular orbit with a period of 4, 12, 20, and 28 days.5.5 m s −1 was chosen because it is approximately the same semi-amplitude as the signal induced by K2-136c, allowing us to directly test our confidence in the recovered RV signal of K2-Figure 5. Left: Phased RV plots for all three K2-136 planets.For each subplot, we used our best fit model parameters to remove stellar activity and the presence of the other two planets.In each subplot, blue data points (HARPS-N) and orange data points (ESPRESSO) are unbinned RV observations.The black line is the median fit and the gray region around that line is the 1σ confidence interval.Right: Posterior mass distribution plots for all three K2-136 planets.It is visually apparent that while there is a strong mass detection for K2-136c (Rp = 3.00 ± 0.13 R ⊕ ), there is at best only a marginal detection for K2-136b (Rp = 1.014 ± 0.050 R ⊕ ) and no evidence of a detection for K2-136d (Rp = 1.565 ± 0.077 R ⊕ ).Therefore, therefore we report only place upper limits on the masses of K2-136b and K2-136d (see Table 4).

Mayo et al.
136c specifically.Our set of orbital periods was selected in order to 1) span the range of known periodic signals in the system (the orbital period of the three known planets and the stellar rotation period), 2) avoid close proximity to those signals (none are within 1.5 days of the injected signals), and 3) be equally spaced in order to uniformly test the encompassed period range.
All four of the recovered signals agree with the injected signal of 5.5 m s −1 to within 1σ.In all four tests, the recovered signal of K2-136c also agrees with our oneplanet reference model (K c = 5.28 ± 0.56 m s −1 ) within 1σ, lending further confidence to our results.

HST NUV observations
To compare the UV quiescent activity of K2-136 with other K stars Hyades members, we measured the surface flux of the Mg II h (2796.35Å) and k (2808.53Å) lines.The Mg II lines are the strongest emission lines in the NUV and correlate strongly with the chromospheric activity of the star.For accurate emission measurements, we subtracted the NUV continuum by fitting the data outside of the Mg II integration region using the astropy module specutils.We then integrated over 2792.0 − 2807.0Å to measure the Mg II emission flux.To convert from observed flux to surface flux, we estimated the radii of the stars using the relationships between age, effective temperature, and radius given by Baraffe et al. (2015), except for K2-136, where we use the radius reported in this work.
Figure 7 shows the Mg II surface flux as a function of the rotation period for K2-136 and 13 other observed K star members of the Hyades from GO-15091.K2-136 has a surface Mg II flux of (8.32 ± 0.17) × 10 5 erg s −1 cm −2 , whereas the median of the sample is (9.66 ± 0.24) × 10 5 erg s −1 cm −2 .
Additionally, Richey-Yowell et al. ( 2019) measured the NUV flux densities of 97 K stars in the Hyades using archival data from the Galaxy Evolution Explorer (GALEX; Martin et al. 2005).The median NUV flux density of their sample of Hyades stars at a normalized distance of 10 pc was 1.89 × 10 3 µJy.We measure a GALEX NUV magnitude for K2-136 of 19.51 mag, which corresponds to a flux density at 10 pc of 1.52 × 10 3 µJy, well within the interquartiles of the total sample from Richey-Yowell et al. (2019).

X-ray observations
We detected a total of 800 EPIC counts for K2-136 in the XMM observation and can therefore extract an Xray spectrum for the star.We used a one-temperature APEC model to fit the spectrum, which is appropriate for representing the hot plasma in stellar coronae.We combined this model with the ISM absorption model tbabs using photoelectric cross-section from Balucinska-Church & McCammon (1992) to account for the neutral hydrogen column density N H .We set N H to 5.5 × 10 18 cm −2 , derived using E(B−V ) = 0.001 for Hyades (Taylor 2006), R V = 3.1, and the relation N H [cm −2 /A v ] = 1.79 × 10 21 (Predehl & Schmitt 1995); allowing N H to float did not improve the fit.The spectra for each of the three XMM cameras are shown in Figure 8 and the best-fit parame- ters, which are obtained from a simultaneous fit to the three, are provided in Table 6.
Using the total EPIC energy flux from the spectral best fit, we obtained an X-ray luminosity L X = (1.26± 0.19)×10 28 erg s −1 and L X /L * = (1.97 ± 0.30)×10 −5 (0.1−2.4 keV energy range).These values are within 1σ of those found by Fernandez Fernandez & Wheatley (2021) using the same XMM observation4 .For the sample of 89 K dwarfs with X-ray detections in the Hyades, the median values for L X and L X /L * are 4.5 +7.9 −3.3 × 10 28 erg s −1 and 4.0 +20.8  −1.6 × 10 −5 , respectively (Núñez et al. in prep.).K2-136, therefore, appears somewhat less luminous in X-rays than most of its coeval K dwarf brethren.A narrower comparison, against late-K (K5 and later) dwarf Hyads, shows that the L X and L X /L * values for K2-136 are within one standard deviation of the median for that cohort.
In addition to X-ray luminosity, we also estimated extreme ultraviolet (EUV) luminosity using stellar age and Equation 4 from Sanz-Forcada et al. (2011) and found L EUV = (22.6 +7.8 −5.7 )×10 28 erg s −1 .Combining L X and L EUV and using the semi-major axis of K2-136c, we are able to estimate the X-ray and UV flux incident on K2-136c to be F XUV = (6.0+2.0 −1.4 )×10 3 erg s −1 cm −2 .Then, using this incident flux value (along with M c and R c ) we estimate an atmospheric mass loss rate with Equation 1 from Foster et al. (2022).This yields a current atmospheric mass loss rate for K2-136c of Ṁc = (3.4+1.3 −0.9 )×10 9 g s −1 = (17.8+6.7  −5.0 )×10 −3 M ⊕ Gyr −1 .This rate is based on current values, so the mass loss rate in the past or future may differ.At this current rate, with a H 2 -He envelope mass fraction of ∼ 5%, it would take 51 +23 −16 Gyr Figure 8. X-ray spectra of K2-136 from the XMM pn (left panel), MOS1 (center), and MOS2 (right) EPIC cameras.X-ray counts are binned by 20 in the pn camera, and 15 in the MOS cameras.The orange lines are the best fits using a one-temperature APEC model and assuming a fixed neutral hydrogen column density of 5.5 × 10 18 cm −2 , typical for Hyads (see Sec. 4.9.2).The residuals of each fit are shown in the bottom panels.The best-fit parameters are presented in Table 6.
to fully evaporate the atmosphere.Even the 95% lower limit evaporation time is still 28 Gyr, longer than the age of the universe.In other words, in ∼ 4 Gyr, when the K2-136 system is as old as the Solar System currently is, we expect K2-136c will have likely only lost 5 − 10% of its current atmosphere.
We also calculated the Rossby number R o of K2-136, which is defined as the star's rotation period P * divided by the convective turnover time τ .We used the (V −K s )log τ empirical relation in Wright et al. (2018) (their equation 5) to obtain τ = 22.3 d for K2-136.Using our measured P * value (see Sec. 3.1) gives R o = 0.6.This R o puts K2-136 well within the X-ray unsaturated regime, in which the level of magnetic activity decays follows a power slope as a function of R o (see figure 3 in Wright et al. 2018).For the sample of 51 K dwarf rotators with X-ray detection in the Hyades, the median R o = 0.46 +0.06 −0.08 (Núñez et al, in prep.), which suggests that the lower levels of X-ray emission from K2-136 relative to its fellow Hyades K dwarfs can be explained by its slower rotation rate.
In conclusion, K2-136 does not appear unusually active in either the NUV or the X-ray relative to its fellow Hyades K dwarfs.

Considering the Nearby Stellar Companion
As discussed in Section 3.2, prior observations of this system revealed a likely bound M7/8V companion with a projected separation of approximately 40 AU.The presence of a nearby stellar companion can easily cause a trend in RV observations.We wanted to estimate the range of possible trend amplitudes for our system using some reasonable assumptions.Using the distance of 58.752 +0.061  −0.072 pc from Bailer-Jones et al. ( 2021) and the projected angular separation of 0.730 ± 0.030" from adding in quadrature the R.A. and Dec separation components in Ciardi et al. (2018), we found a projected separation of 42.9 ± 1.7 AU.We considered the possibil- ity of additional radial separation by folding in a uniform distribution on radial separation between 0 and twice the median projected separation to estimate an overall separation (this broad range was chosen to include radial separations of approximately the same scale as the projected separation).As an approximation, we treat this overall separation as the semi-major axis.
Next, the stellar companion was reported in Ciardi et al. (2018) to have a spectral type consistent with M7/8 (we were unable to find uncertainties associated with this result, but nearby spectral types were never mentioned).Taking a conservative approach, we assumed a stellar companion mass between 1 M Jup (for a low-mass brown dwarf) and 0.2 M (for a mid to late M dwarf).
Then, we combined our separation distribution and stellar mass distributions (primary and companion) via Kepler's Third Law to get a broad orbital period estimate of 520 +320 −180 years.Including a wide range of eccentricities (Unif(0,0.9)),we estimated an RV semi-amplitude of 490 +340 −320 m s −1 .On the timescale of our observations, this centuries-long sinusoidal signal would manifest as a linear trend, with a maximum (absolute) slope of 5.3 +7.5 −3.6 m s −1 yr −1 .This is a very rough estimate with many assumptions, but it serves to demonstrate that a drift of only a few m s −1 each year or less is very reasonable.Indeed, from our model of the RV data we found an RV trend of 0.1±2.5 m s −1 yr −1 , highly consistent with both a zero trend as well as our estimate calculated here.

RESULTS AND DISCUSSION
The results of our stellar and planet analyses are listed in Tables 1, 4, and 3. Phase plots of all three planets can be seen in Fig. 5.After conducting our analysis and tests, we find that K2-136c has a mass of 18.1 +1.9 −1.8 M ⊕ and a radius of 3.00 ± 0.13 R ⊕ .This radius is consistent with and slightly larger than the value estimated in M18 (2.91 +0.11  −0.10 R ⊕ ).This is because we find a stellar radius value slightly larger than M18 (by about 3%).
Using planet mass and radius we find K2-136c has a density of 3.69 +0.67  −0.56 g cm −3 (or 0.67 +0.12 −0.10 ρ ⊕ ).For comparison, Neptune 2 is roughly similar in mass (17.15 M ⊕ ) but larger in radius (3.883 R ⊕ ); as a result, K2-136c is more than twice as dense as Neptune (2.25 +0.41  −0.34 ρ ).Similarly, Uranus 3 is slightly less massive (14.54 M ⊕ ) but still larger in radius (4.007 R ⊕ ); thus, K2-136c is nearly three times as dense as Uranus (2.90 +0.52  −0.44 ρ ).This is visually apparent in Fig. 6, which shows K2-136c almost perfectly on the 100% H 2 O composition line.It is important to remember that mass and radius alone do not fully constrain a planet's composition.Although K2-136c may have a density similar to that of a large ball of water (an unrealistic reference composition), it is also consistent with a gaseous sub-Neptune with a massive core or metal-rich atmosphere.
Unfortunately, it is very difficult to determine compositional properties of a planet without atmospheric characterization, especially sub-Neptunes (since there are no analogs in our own Solar System).On the one hand, sub-Neptunes may include ocean worlds with H 2 O abundance fractions not seen in our Solar System (Mousis et al. 2020).And indeed, water vapor has already likely been detected in the atmosphere of the sub-Neptune exoplanet K2-18b (Benneke et al. 2019b).On the other hand, sub-Neptunes like K2-136c may instead be composed of a rocky, Earth-like core composition, very little water, and an atmosphere close to solar metallicity and thus primarily hydrogen and helium (Van Eylen et al. 2018;Benneke et al. 2019a).
As a valuable point of comparison, there are three confirmed planets that share a similar mass and radius to K2-136c to within 10%: Kepler-276c, Kepler-276d (Xie 2014), and TOI-824b (Burt et al. 2020).The masses of the planets in the Kepler-276 system were measured via transit timing variations (TTVs) in a TTV catalog paper.Unfortunately, because they were characterized alongside so many other systems, there is no discussion regarding the formation or composition of those two specific planets.TOI-824b, however, was characterized in a standalone paper that investigated the nature of the planet thoroughly.Unlike K2-136c, TOI-824b is near 2 https://nssdc.gsfc.nasa.gov/planetary/factsheet/neptunefact.html 3 https://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html the hot-Neptune desert (Mazeh et al. 2016) with a very short orbital period (1.393 d).Despite its proximity to its host star, TOI-824b still retains a H 2 -He atmosphere, which the authors estimate has a mass fraction of ≥ 2.8%.They hypothesize that the larger than average mass of the planet (compared to planets of a similar radius) helps the planet retain its atmosphere.K2-136c, with a similar mass and radius and a lower insolation flux, would therefore be able to retain a H 2 -He atmosphere even more easily.
We can go further and form a picture of a reasonable composition for K2-136c.We may assume the planet has a rocky core surrounded by a gaseous H 2 -He envelope.Going a step further, we may also assume the rocky core is similar to that of Earth, namely a core-mass fraction (CMF) of 0.325, i.e. a rock-iron composition of 32.5% Fe and 67.5% MgSiO 3 (Seager et al. 2007).Following the theoretical models of Howe et al. (2014), we find that with an Earth-like rocky core and a H 2 -He envelope, the measured mass and radius of K2-136c are most consistent with a H 2 -He mass fraction of ∼ 5%.
With this rough envelope mass fraction estimate and the measured mass and radius of K2-136c, we wanted to investigate the potential for past or ongoing atmospheric mass loss.After consulting the theoretical models presented in Lopez et al. (2012), Lopez & Fortney (2013), Lopez &Fortney (2014), andJin et al. (2014), we concluded that if there is any historical or contemporary mass loss for K2-136c, it is minimal: likely somewhere between 0 − 10% loss of the H 2 -He envelope across the entire lifespan of the planet.
As suggested in Mann et al. (2016), young planets may be puffier than older planets due to an early-age atmospheric mass loss phase.And yet, K2-136c is not particularly puffy, in fact being notably denser than Uranus or Neptune.However, this planet may indeed have an extended atmosphere but also a lower atmospheric mass fraction than Uranus, Neptune, and other lower-density planets.In other words, as this system ages, the atmosphere of K2-136c may settle to some extent, reducing the planet radius and increasing planet density.Without atmospheric characterization, further insights into the planet's composition are very limited.
Bayesian model comparison proved that we could not conclusively detect K2-136b or K2-136d in our data set (see Fig. 5).Even so, we also conduct similar analyses for K2-136b and K2-136d in order to report upper mass limits and corresponding upper density limits.With 95% confidence, K2-136b is no denser than 24 g cm −3 (a largely unhelpful limit, since that would be much more dense than pure iron) and K2-136d is no denser than 4.3 g cm −3 , corresponding to semi-amplitudes of 1.7 m s −1 and 0.78 m s −1 , respectively.Referring to Fig. 6, we can see that unlike the middle planet K2-136c, the other two planets K2-136b and K2-136d could have a wide variety of densities and compositions.K2-136d could range from a low-density gas planet to Earth composition.K2-136b has an even wider array of possible compositions and theoretically could range from very gaseous to pure iron.
The RV signals of K2-136b and K2-136d may be detectable with more data from next generation spectrographs.K2-136d would be particularly interesting, since its radius places it near the planet radius gap (Fulton et al. 2017).As for K2-136b, the peak of the planet mass posterior distribution is already non-zero, and a marginal detection may already be noted at the 2.0σ level (2.4±1.2M ⊕ ), suggesting a firm planet mass measurement may be within reach with further observations.
To check this, we followed the mass-radius relationship laid out in Wolfgang et al. (2016) and found predicted masses for K2-136b and K2-136d to be 1.17 +0.79  −0.72 M ⊕ and 4.1 +1.9 −1.8 M ⊕ , respectively.These correspond to densities of 7.7 +3.7 −4.7 g cm −3 (i.e.1.40 +0.66 −0.85 ρ ⊕ ) and 8.4 +4.1 −3.7 g cm −3 (i.e.1.52 +0.74  −0.66 ρ ⊕ ), respectively.We note that these mass and density estimates do not make use of the upper limits determined in this paper.By folding in our stellar mass, orbital period, and eccentricity posteriors as well as the orbital inclinations determined in M18, we found the estimated masses of K2-136b and K2-136d correspond to semi-amplitudes of 0.5±0.3m s −1 and 1.1±0.5 m s −1 , respectively.The current RV upper limit on K2-136b (1.7 m s −1 ) is much larger than the estimated semi-amplitude and therefore fully consistent.As for K2-136d, we ac-knowledge that the estimated semi-amplitude is smaller than the upper limit (0.80 m s −1 ), which perhaps suggests K2-136d has a density on the lower end of the range predicted from the Wolfgang et al. (2016) relationship.
There are very few young and small exoplanets that also have measured masses.As can be seen in Fig. 9, the vast majority of young exoplanets do not have a firm mass measurement.There are some young planets with both robust radius measurements and notable upper mass limits, such as the low-density planet TS Duc A b (Benatti et al. 2021), which can be of interest for follow-up study and comparison.However, according to the NASA Exoplanet Archive (accessed 2023 Mar 12; NASA Exoplanet Science Institute 2020), there are only 13 known planets, excluding K2-136c, with R p < 4 R ⊕ , a host star age < 1 Gyr, and a mass measurement (not an upper limit): HD 18599b (Desidera et al. 2022), HD 73583b and c (Barragán et al. 2022); K2-25b (Stefansson et al. 2020);L 98-59b, c, and d (Demangeon et al. 2021); Kepler-411b and Kepler-411d (Sun et al. 2019), Kepler-462b (Masuda & Tamayo 2020), Kepler-289b and Kepler- The gray line in both panels is the modelled atmospheric spectrum assuming low metallicity (5x solar) and no clouds, while the blue data points are the simulated instrument spectra for one observed transit with NIRISS SOSS, and two observed transits with NIR-Spec G395M, including the effects of photon noise, and instrument systematics.
K2-136c is now the smallest exoplanet in an open cluster to have a mass measurement.It is also one of the youngest exoplanets to ever have a mass measurement.The only planets with firm age, radius, and mass measurements (< 25% uncertainties) known to be younger are AU Mic b and c (Klein et al. 2021) as well as Kepler-411b and d (Sun et al. 2019), as can be seen in Fig. 9.In general, measuring the masses of young planets like K2-136c provides an interesting window into the early childhood of planetary systems, allowing us to probe how planet masses and compositions evolve over time.

Atmospheric Characterization Prospects
To explore the suitability of the K2-136 system planets for atmospheric characterization, we calculated the transmission spectroscopy metric (TSM) defined in Kempton et al. (2018).K2-136c has a TSM of 32.7 +4.8 −4.1 , which is well below the recommended TSM of 90 for R p > 1.5 R ⊕ .Because K2-136b and K2-136d have unconstrained masses, we followed Zeng et al. (2016) and assumed an Earth-like CMF of 0.325 in order to predict planet masses of 0.85 +0.27  −0.21 M ⊕ and 3.35 +1.08 −0.84 M ⊕ , respectively.For K2-136b, this yields a TSM of 4.20 +0.42  −0.37 , well below the recommended TSM of 10 for R p < 1.5 R ⊕ .As for K2-136d, its radius of R d = 1.565 ± 0.077 R ⊕ is very near 1.5 R ⊕ , where the TSM metric includes a scale factor that jumps dramatically, thus creating a TSM bi-modal distribution.Thus, for R d < 1.5 R ⊕ we find a TSM of 2.36 +0.15  −0.13 and for R d > 1.5 R ⊕ we find a TSM of 13.6 +1.0 −1.1 .In their respective radius ranges, these values are both well below the recommended TSM, so K2-136d is also probably not a good target for atmospheric characterization.
We also calculated the emission spectroscopy metric (ESM) for K2-136b and K2-136d as defined in Kempton et al. (2018) (the metric applies only to "terrestrial" planets with R p < 1.5 R ⊕ , excluding K2-136c).We find K2-136b and K2-136d have ESM metrics of 0.520 +0.049 −0.047 and 0.342 +0.043 −0.041 , respectively, both below the recommended ESM of 7.5 or higher.Therefore, these planets do not appear to be particularly attractive targets for emission spectroscopy or phase curve detection.
The TSM analysis of K2-136c is not very favorable for atmospheric characterization, but we decided to conduct a more thorough transmission spectroscopy analysis.We used the JET tool (Fortenbach & Dressing 2020) to model atmospheric spectra and to simulate the performance of the JWST instruments for certain atmospheric scenarios.We opted for the broad wavelength coverage of combining NIRISS SOSS Order 1 (0.81-2.81 µm) with NIRSpec G395M (2.87-5.18µm), as recommended by Batalha & Line (2017) to maximize the spectral information content.A single instrument, the NIRSpec Prism, can also cover this wavelength range, but brightness limits preclude its use here.We assumed pessimistic pre-launch instrumental noise values (Rigby et al. 2022), so a future JWST program for atmospheric characterization should outperform our conservative expectations.
The JET tool found that for an optimistic, cloudless, low-metallicity atmosphere (5x solar) we can meet a ∆BIC detection threshold of 10 (corresponding to a ∼ 3.6σ detection of the atmosphere compared to a flat line) with 5 free retrieval parameters (i.e., recon level) with only one transit for NIRISS SOSS and two transits for NIRSpec G395M.For a less optimistic, higher metallicity atmosphere (100x solar) with clouds at 100 mbar, we can meet the same detection threshold with two transits for NIRISS SOSS and five transits for NIR-Spec G395M.
With 10 free retrieval parameters (a more typical number), and with the optimistic atmosphere, we can meet a ∆BIC detection threshold of 10 with only one transit for NIRISS SOSS and three transits for NIRSpec G395M.For the less optimistic atmosphere, we can meet the detection threshold with three transits for NIRISS SOSS, but we show no detection for NIRSpec G395M for up to 50 transits considered.This analysis makes the conservative assumption that the instrument noise floor is not reduced by co-adding transits.
The resulting spectra for the low metallicity, cloudless, case are shown in Fig. 10.It seems that atmospheric characterization of K2-136c may be within reach (assuming a relatively low mean molecular weight/low metallicity atmosphere, and low cloud level), but could require a more significant investment of JWST resources if the actual atmospheric properties are less favorable.
It should be noted that given the on-sky position of K2-136, the ability to observe the system with JWST will be limited due to aperture position angle constraints.In addition, the very close (∼ 0.7") stellar companion may

Figure 1 .
Figure1.Periodograms of RV, CCF FWHM, CCF BIS, and S HK for the K2-136 system.In the top panel is the window function (computed from observation times only).Each subplot has the periodogram of HARPS-N and ESPRESSO combined (black), HARPS-N alone (blue), and ESPRESSO alone (orange).The gray region corresponds to the 1σ confidence interval of the stellar rotation period (as determined from our model results); the three vertical, black lines correspond to the orbital periods of K2-136b, c, and d.Finally, the horizontal dashed lines refer to different false alarm probabilities (for HARPS-N and ESPRESSO combined).

Figure 2 .
Figure2.Scatter plots of S HK , BIS, and FWHM against RV for the K2-136 system.All HARPS-N and ESPRESSO data have been separately offset shifted according to the median model offsets listed in Table4.Blue data points correspond to HARPS-N observations, orange data points to ESPRESSO.In the top-left corner of each subplot is the Spearman correlation coefficient, which can capture nonlinear, monotonic correlations.(Coefficients were calculated with data values but not uncertainties.)

Figure 3 .
Figure 3. Transit plot of K2-136.The top and bottom subplots of the top plot are the raw and normalized K2 photometry versus time from Campaign 13, respectively.In the top subplots of the bottom plot are the phase-folded light curves and transit model fits for K2-136b, c, and d.The gray points are the raw data.The best-fit transit model is the orange line and binning is represented by the blue points.The bottom subplots of the bottom plot are the residuals after the best-fit model has been subtracted.

Figure 4 .
Figure 4. K2-136 observations and model fits for RVs (and fit residuals).In each panel, the blue points (HARPS-N) or orange points (ESPRESSO) are the observations while the black line and gray region are the model fit and 1σ confidence interval, respectively.RVs have been mean-subtracted (corresponding to their respective instrument) and planet-induced reflex motion has been subtracted as well.RV errors have been inflated from their original values by adding the model-estimated RV jitter term in quadrature.Note the two time gaps between the first 88 HARPS-N observations, the 22 ESPRESSO observations, and the final 5 HARPS-N observations.

Figure 6 .
Figure 6.Top: mass-radius diagram of transiting planets with fractional mass and radius uncertainties less than 50%.K2-136b, c, and d are plotted in dark orange, with mass uncertainties on K2-136b and d as 68% and 95% upper limits denoted in light orange.Data collected from the NASA Exoplanet Archive (accessed 2023 Mar 22; NASA Exoplanet Science Institute 2020).Venus, Earth, Uranus, and Neptune are also labeled and plotted in blue for reference.Except for the K2-136 system, planets with larger fractional mass and radius uncertainties are fainter.Gray lines correspond to planetary compositions (from top to bottom) of 100% H 2 O, 50% H 2 O, 25% H 2 O, 100% MgSiO 3 , 50% MgSiO 3 + 50% Fe, and 100% Fe, respectively (Zeng & Sasselov 2013; Zeng et al. 2016).Kepler-136c lies closest to the 100% H 2 O composition line, and is similar in mass to Uranus and Neptune although smaller and much more dense.Middle: the same sample plotted in mass versus planet density, with the same solar system references and composition lines (order inverted from top panel).Bottom: the same sample plotted in planet insolation flux versus planet density, with the color of all data points corresponding to planet mass (except K2-136b and d, which only have mass upper limits and thus are plotted in light orange).

Figure 9 .
Figure 9. Age-radius diagram for all planets smaller than Jupiter, younger than 5 Gyr, and with radius and age uncertainties both smaller than 25%.Data point color corresponds to planet density for planets with mass uncertainties smaller than 25%.Data collected from the NASA Exoplanet Archive (accessed 2023 Mar 22; NASA Exoplanet Science Institute 2020).K2-136b, c, and d are labeled in red to the left of the planet symbol.There are only four planets in this figure with a stellar age younger than K2-136 and a mass measurement better than 25%: AU Mic b and c (stellar age = 22 ± 3 Myr; Mamajek & Bell 2014) and Kepler-411 b and c (stellar age = 212 ± 31 Myr; Sun et al. 2019).The only other plotted planet < 1 Gyr with a mass measurement is the open cluster planet K2-25b (Stefansson et al. 2020), which is slightly larger than K2-136c.

Figure 10 .
Figure10.Simulation of JWST transmission spectra for K2-136c using JET(Fortenbach & Dressing 2020).The top and bottom panels correspond to the NIRISS SOSS-Or1 (0.81-2.81 µm) instrument, and the NIRSpec G395M (2.87-5.18µm) instrument, respectively.The gray line in both panels is the modelled atmospheric spectrum assuming low metallicity (5x solar) and no clouds, while the blue data points are the simulated instrument spectra for one observed transit with NIRISS SOSS, and two observed transits with NIR-Spec G395M, including the effects of photon noise, and instrument systematics.

Table 5
Model Evidence Comparisons for K2-136 Planet Configurations Douglas et al. (2019) a function of stellar rotation period for K star Hyades members.The blue star represents K2-136, and the black points are the 13 other K star Hyades observed in a broader HST program (GO-15091, PI: Agüeros).The rotation periods are fromDouglas et al. (2019)and have assumed errors of 10%, except for K2-136, which we determined to be 13.88 +0.17 −0.18 d in this work.K2-136 does not show any distinct chromospheric activity or unique rotation period compared to the rest of the sample.