K2-79b and K2-222b: Mass measurements of two small exoplanets with periods beyond 10 days that overlap with periodic magnetic activity signals

We present mass and radius measurements of K2-79b and K2-222b, two transiting exoplanets orbiting active G-type stars. Their respective 10.99d and 15.39d orbital periods fall near periods of signals induced by stellar magnetic activity. The two signals might therefore interfere and lead to an inaccurate estimate of exoplanet mass. We present a method to mitigate these effects when radial velocity and activity indicator observations are available over multiple observing seasons and the orbital period of the exoplanet is known. We perform correlation and periodogram analyses on sub-sets composed of each target's two observing seasons, in addition to the full data sets. For both targets, these analyses reveal an optimal season with little to no interference at the orbital period of the known exoplanet. We make a confident mass detection of each exoplanet by confirming agreement between fits to the full radial velocity set and the optimal season. For K2-79b, we measure a mass of 11.8 $\pm$ 3.6 $M_{Earth}$ and a radius of 4.09 $\pm$ 0.17 $R_{Earth}$. For K2-222b, we measure a mass of 8.0 $\pm$ 1.8 $M_{Earth}$ and a radius of 2.35 $\pm$ 0.08 $R_{Earth}$. According to model predictions, K2-79b is a highly irradiated Uranus-analog and K2-222b hosts significant amounts of water ice. We also present an RV solution for a candidate second companion orbiting K2-222 at 147.5d.


INTRODUCTION
The HARPS radial velocity (RV) survey first revealed the existence of a sub-population of super-Earth and Neptunelike planets in tight orbits (Mayor & Udry 2008;Lovis et al. 2009). A few years later, NASA's Kepler mission found that these small planets, with sizes between that of Earth and Neptune, are the most abundant type of exoplanets with periods less than ≈100d (e.g. Borucki et al. 2011;Batalha et al. 2013;Fressin et al. 2013). These types of planets are not found in the Solar System, and their existence was not predicted by contemporary planet-formation models (Ida & Lin 2004a,b;Mordasini et al. 2009a,b). Subsequent studies, using a greater number of small exoplanets from Kepler, re-vealed relatively few planets with radii between 1.5R ⊕ and 2R ⊕ , now commonly referred to as the radius valley (Fulton et al. 2017;Zeng et al. 2017a,b;Van Eylen et al. 2018). There also appears to be a scarcity of exoplanets in the sub-Saturnian desert, with radii larger than 4R ⊕ but smaller than a gas giant (Zeng et al. 2018). Exoplanets in the sub-Saturnian desert appear separated from Sub-Neptunes by the radius cliff, a steep drop-off in exoplanet occurrence starting near 3R ⊕ and attributed to the non-linear solubility of hydrogen in magma at relatively high pressures (Kite et al. 2019).
The exoplanet community has increased efforts to measure masses of small exoplanets to develop a more complete picture of the processes involved in their formation. Bulk densities measured for exoplanets on both sides of the radius valley suggest that those above the gap have very different compositions from those below (Rogers 2015;Zeng et al. 2018Zeng et al. , 2019. Exoplanets below the gap have relatively high densities, consistent with rocky compositions. Exoplanets above the gap have lower densities that are degenerate, consistent with rocky cores surrounded by envelopes of liquid and/or gas.
Several mechanisms have been proposed in the literature to explain the observed small planet radius valley: Photoevaporation, internal heat-powered mass loss, and water-rich vs. water-poor formation. Photo-evaporation refers to the shedding of a planet's envelope from prolonged exposure to its host star (Hansen & Murray 2012;Mordasini et al. 2012b,a;Alibert et al. 2013;Chiang & Laughlin 2013;Chatterjee & Tan 2014;Coleman & Nelson 2014;Lee et al. 2014;Raymond & Cossou 2014;Lee & Chiang 2016). Like photo-evaporation, core-powered mass-loss also causes an exoplanet to shed its envelope, but is instead powered by radiation from the hot rocky core as it cools and winds in the hot gaseous envelope (Owen & Wu 2016;Ginzburg et al. 2016). Successful models of photo-evaporation, corepowered mass-loss, and atmospheric "boil off" thus far have assumed that all planets above the gap host Hydrogen/Helium (H or H/He) envelopes, excluding water from their simulations. However, simulations including water and ices have been carried out to form water-rich exoplanets with rocky or rocky/icy cores that accrete significant amounts of water and gaseous envelopes if massive enough (Zeng et al. 2018(Zeng et al. , 2019. Comparing small planets of varying surface gravity across a range of stellar irradiation levels can help identify driving factors behind the observed population of small exoplanets (Lopez & Rice 2018;Cloutier & Menou 2020). So far, the majority of small planets with mass measurements occur at relatively short orbital periods (P orb < 10 days), and therefore are subject to high irradiation levels. Mass and radius measurements of small planets at longer periods are therefore essential for developing a more complete understanding of small exoplanet formation.
Measuring small exoplanet masses at longer orbital periods presents challenges, including hard-to-detect smallamplitude Doppler signals, and potential interference between periodic signals from stellar activity and P orb . The second problem is even more challenging than initially recognized because P orb need not overlap perfectly with P rot or one of its multiples/harmonics for the stellar activity and Doppler signals to interfere (Nava et al. 2020). Mortier et al. (2016) tackled the feat of determining the mass on an exoplanet with P orb near P rot , but generally these cases have been set aside to measure masses of exoplanets with less challenging interference from stellar magnetic active regions.
In this paper we present a method to mitigate these potential interference effects when radial velocity (RV) and activity indicator observations are available for a target over multiple observing seasons. We apply our method to two exoplanets discovered by the K2 mission, K2-79b and K2-222b, whose orbital periods overlap with observed activity signals from their host stars. K2-79b was first detected and validated by Crossfield et al. (2016), and then re-detected by Mayo et al. (2018) and Kruse et al. (2019). K2-222b was detected as a planet candidate by Mayo et al. (2018) and Petigura et al. (2018). It was validated by Mayo et al. (2018) and then later re-detected by Livingston et al. (2018). Given their radii of 4.09R ⊕ and 2.35R ⊕ , respectively, measuring masses for these two planets is particularly interesting, since they lie in the sub-Saturnian desert and near the upper edge of the radius valley where their compositions could be water-rich or waterpoor.
In Section 2, we describe the data sets utilized in our analysis, and in section 3, we detail the characterization of the two host stars. In section 4 we describe our investigations of potential stellar activity signals in the light curve (LC) and RV data sets of each target. Section 5 describes our MCMC fits to determine the exoplanet parameters. We report our results and discuss the scientific implications for K2-79b and K2-222b in Section 7.

K2 Photometry
We analyzed 3161 and 3424 photometric K2 observations of K2-79 and K2-222, respectively, collected in long cadence (29.4 minutes) mode (Figures 1a and 1b). K2-79 was observed between February 10, 2015and April 20, 2015, and K2-222 between January 6, 2016and March 23, 2016. The K2 spacecraft operated with only two of the four original reaction wheels, causing the spacecraft to drift over time, with intermittent thruster firings to keep desired targets in the telescope's field of view. This led to long-term trends in the light curve (LC) of a given star, due to the stars drifting over time across pixels with different quantum efficiencies. We corrected for thruster firings first, followed by systematic instrumental effects and other low frequency signals, according to the technique developed by Vanderburg & Johnson (2014). The top panels in Figures 1a and 1b show the LCs of each target after correcting for those effects. After deriving a firstpass systematics corrected light curve, we refined the correction by simultaneously fitting for spacecraft systematics, low-frequency variability, and exoplanet transits following Vanderburg et al. (2016). The bottom panels in the same figures show the resulting flattened LCs with exoplanet transits isolated. We set the error of each point in a given LC data set equal to the standard deviation of the out-of-transit points in the flattened LC. In the case of the brighter star, K2-222, the second half of the observations have increased scatter due to uncorrected K2 systematics. We therefore calculate errors for that LC in two segments separated at BJD -2400000 = 57438 days, where the start of the increased scatter begins (see Figure 1b, bottom panel).

HARPS-N Spectroscopy
We collected 79 spectra of K2-79 over four seasons between November 4, 2015 and December 29, 2019 (Figure 2), and 63 spectra of K2-222 over three seasons between August 14, 2016 and December 23, 2019 ( Figure 3). All spectra were collected with HARPS-N, the high precision spectrograph mounted on the Telescopio Nationale de Galileo at the Observatorio del Roque de los Muchachos in La Palma, Spain (Cosentino et al. 2012). The spectrograph covers wavelengths in the range 383 -690 nm, with average resolution R = 115,000. Observations for K2-79 were taken with 30 minute exposure times, yielding a mean signal-to-noise ratio (SNR) of 32.67 in order 60 (627.077 -634.166 nm) of the spectrum. For K2-222, observation exposure times ranged between 15 -30 minutes, yielding a mean SNR of 74.11 in order 60. The observations were collected as part of the HARPS-N collaboration's Guaranteed Time Observations (GTO) program.
We reduced the spectra and calculated associated RVs and activity indicators using the HARPS-N data reduction software (DRS; Lovis & Pepe 2007). Activity indicators calculated include the contrast, FWHM, cross-correlation function area (CCF area = contrast * FWHM), bisector inverse span (BIS), and S-index (Figures 2 and 3).

STELLAR CHARACTERIZATION
We derived stellar parameters for each target using the Stellar Parameter Classification Tool (SPC; Buchhave et al. 2012Buchhave et al. , 2014 to analyze the K2-79 and K2-222 HARPS-N spectra and estimate stellar effective temperature (T eff ), surface gravity (log g), metallicity ([Fe/H]), and projected rotational velocity (v sin i) for each target. Spectra with normalized peak heights of less than 0.9 in the cross-correlation function were discarded due to low signal-to-noise ratios. We aligned all individual spectra to a common RV frame and calculated the maximum stellar rotation period associated with the lower limit of v sin i using the equation P rot,max = 2πR * /v sin i. Final values from SPC are summarized in Table  6.
We also applied the ARES+MOOG equivalent width method to the K2-79 and K2-222 HARPS-N spectra to calculate T eff , log g, [Fe/H], and microturbulence (χ turb ) (Sousa 2014) for each target. Spectra associated with 3-sigma RV outliers were discarded. As with SPC, all used spectra were shifted to a common RV frame before combining them for analysis. For accuracy, the value of log g was corrected following Equation 3 in Mortier et al. (2014), and its error bars were increased according to Sousa et al. (2011). Final values measured with ARES+MOOG and the combined results with SPC are summarized in Table 6.  Table 6), we used the Isochrones (Morton 2015) software to run four isochrone analyses for each target: Dartmouth models (Dotter et al. 2008) and MIST models (Choi et al. 2016) each paired with both the SPC and ARES + MOOG obtained T eff and metallicity values. The Gaia parallax and magnitudes in various filter are listed in Table 6. We used 400 live points for the Multinest algorithm and combined posteriors from all four isochrone runs, as first presented in Malavolta et al. (2018), to obtain final values for T eff , log g, [Fe/H], extinction (A(V)), stellar luminosity, radius, mass, density, distance, and age. For both targets, parameter values between our four isochrone runs were in agreement with each other within 1.5-sigma. The final isochrone values listed in Table 6 were obtained by combining results from the four isochrone fits. When applicable, we inflate errors on values obtained with any method to 2% in temperature and radius, 4% in luminosity, 5% in mass, and 20% in age, according to Tayar et al. (2020)'s approximations for the systematic uncertainty floor associated with current measured parallaxes and bolometric fluxes. The final derived isochrone results had 1-sigma agreement with relevant pa-  K2-222 (right). The top panels show the LCs after removing the thruster firings and other high-frequency systematic signals. The teal dashed lines show the best-fit third-degree polynomials, and the second panels show the LCs after subtracting the respective polynomials (discussed in Section 4.1). The third panels show the Everest-reduced LCs (also discussed in Section 4.1). The bottom panels show the flattened LCs after removing all remaining low-frequency signals. In the case of K2-222, both the scatter of flux values and related errors increase in a later sub-set of the data, beginning at approximately BJD-2400000 = 57438d.

STELLAR ACTIVITY ANALYSIS
We performed analyses of the K2-79 and K2-222 LCs and the HARPS-N RVs and activity indicators to search for prominent stellar activity signals. The next two subsections detail these analyses.

Light Curve Stellar Activity Analysis
To probe prominent stellar activity signals, we calculated the Generalized Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009) on the LC of each target with only K2 thruster firings removed (top panels in Figure 1). We also performed auto-correlation function (ACF) analysis on each LC. We utilized the Astropy.timeseries Lomb-Scargle package to calculate all GLS periodograms using the de-     fault auto method which selects the best periodogram algorithm based on the input data. We calculate all periodograms with a sampling of N = 50,000 for periods between 0.1 days and one-half the baseline time span of the LC, or corresponding frequencies (Astropy Collaboration et al. 2013. We calculate a 1% false alarm probability (FAP) level using the Astropy.timeseries Lomb-Scargle package's default false alarm level tool, which calculates the upperlimit to the alias-free probability using the approach outlined in Baluev (2008). The resulting periodograms of the transitmasked K2-79 and K2-222 LCs are shown in the top panels of Figure 4.
To calculate each window function, we produce a signal with the same time stamps as the observations, replacing the LC fluxes with ones. We then calculate a GLS periodogram on the signal with the same period and frequency limits used to calculate the LC periodogram, this time applying the scipy periodogram algorithm which treats individual points without error. The resulting window functions of the K2-79 and K2-222 LCs are shown in the bottom panels of Figure 4.
In the cases of both targets, the top panel periodograms are dominated by long-term instrumental systematics that can be seen by eye in the LCs in the top panels of Figure 1. Therefore, for each target, we perform periodogram analyses on two more LCs: one with the best-fit third-degree polynomial subtracted (second panels in Figures 1 and 4) and the other reduced with Everest (Luger et al. 2016), an opensource pipeline for removing instrumental noise signals from K2 LCs (third panels in Figures 1 and 4). We truncated the Everest LCs to match the baseline of the observations used in our analysis and removed extreme outliers. After these cuts, the Everest-reduced LCs still include some points excluded in our reduction, but the window functions for the two sets are identical.
Two statistically significant (FAP < 1%) signals found in both targets' periodograms are likely due to remaining highfrequency instrumental systematics in the LCs. In the K2-79 and K2-222 polynomial-and Everest-reduced LCs, respectively, the first signal occurs at 12.9d, 12.7d, 12.9d, and 12.6d, and the second at 10.3d, 9.33d, 10.7d and 10.6d. We attribute the scatter in values to the different reduction methods being applied, but the set of peaks can be spotted by-eye in both targets' periodograms.
We calculated the ACF by applying discrete shifts to the polynomial-reduced and Everest-reduced LCs and crosscorrelating the shifted LCs with the original. The ACF shows repeating peaks separated by P quasi , a timescale often related to P rot , with the correlation power of each peak dropping off at a rate related to the evolution timescale of magnetic active regions, τ. We utilized the scipy non-linear least squares function to fit an under-damped simple harmonic oscillator curve to each ACF, described by: where A is the ACF amplitude, and y 0 is the oscillator offset from zero. The fits provide estimates of P quasi and τ for each star (McQuillan et al. 2014;Giles et al. 2017). For each fit we also tried a model accounting for the common scenario of two magnetic active regions evenly spaces around the star.
We did this be including a second cosine term, B cos 2πt P quasi , inside the brackets of Equation 1. However, all of the best fits returned and amplitude consistent with zero for B. We also place an upper bound of 100 days on τ, as it will suffice to identify a stable (τ >> P quasi ) signal with only 35-40 day ACF baselines. Table 4. K2-79 HARPS-N RV and stellar activity indicator (contrast, cross-correlation function (CCF) area, full-width at half maximum (FWHM), bisector inverse span (BIS), and S-index) observations ( Figure 2). Table 4 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.  Table 5. K2-222 HARPS-N RV and stellar activity indicator (contrast, cross-correlation function (CCF) area, full-width at half maximum (FWHM), bisector inverse span (BIS), and S-index) observations ( Figure 3). Table 5 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. The results of the ACF analysis are shown in Figure 6a for the K2-79 LC and in Figure 6b for K2-222. In the case of K2-79, the fit to the ACF of the polynomial-reduced LC reveals the strongest periodicity at P quasi = 26.6 ± 0.4 days, with an evolution timescale of τ = 72 ± 6 days. The ACF of the Everest-reduced LC reveals the strongest periodicity at P quasi = 26.1 ± 0.2 days, with an evolution timescale of τ = 100 ± 10 days. In the case of K2-222, the fit to the ACF of the polynomial-reduced LC reveals the strongest periodicity at P quasi = 28.2 ± 0.2 days, with an evolution timescale of τ = 100 ± 10 days. The ACF of the Everest-reduced LC returns an amplitude consistent with zero, dominated by the leftover long-period trend that can be seen by-eye in the LC ( Figure  1b, third panel). The P quasi values returned by successful fits for both targets are similar to each other and the strongest signals in all LC periodograms ( Figure 4). The τ values are also all quite long relative to P quasi . It is possible but unlikely that both stars coincidentally have similar rotation periods with very slow evolution of magnetic active regions. Therefore, we cannot rule out the potential that the ACF signal is the result of remaining systematics for either target. We will rely on our analysis of the HARPS-N RVs in the next section to further search for periodic stellar activity signals.

RV Stellar Activity Analysis
We first estimated values of P rot for both targets using the Noyes et al. (1984) correlation between P rot and R HK , and applying a reddening correction of E(B-V) = 0.20 in the case of K2-79 (Green et al. 2019). This method yields a value of P rot,Noyes ≈ 21.3d in the case of K2-79 ( S-index = 0.159) and P rot,Noyes ≈ 8.0d in the case of K2-222 ( S-index = 0.160). We note that lower metallicity stars have higher R HK values and therefore shorter estimated values of P rot , and vice versa, when using the Noyes et al. (1984) correlation. Hence, for K2-79 P rot,Noyes may be slightly overestimated ([Fe/H] = Table 6. Stellar parameters for K2-79 and K2-222. For values with multiple estimates, we indicate the adopted value with asterisks. Abbreviations for analyses applied to the HARPS-N spectra are: A+M = ARES + MOOG (Sousa 2014), SPC = Stellar Parameter Classification tool (Buchhave et al. 2012(Buchhave et al. , 2014. The adopted isochrone values were obtained by combining results from four different isochrone fits detailed in Section 3. Tables 9 and 10 list stellar parameters for K2-79 and K2-222 obtained with fits using the EXOFAST v2 software (Eastman et al. 2019 . After applying the Everest reduction (third panel), additional statistically significant peaks occur at 25.9d (f = 0.0386d −1 ), 16.7d (f = 0.0599d −1 ), and 7.0d (f = 0.1427d −1 ). The strongest peak in each is consistent with the estimated P rot,max (23.8d ± 4.4d) within 0.5σ, and may be related to stellar rotation.
(b) K2-222: After removing the best-fit third-order polynomial (second panel), additional statistically significant peaks occur at 28.2d (f = 0.0355d −1 ) and 7.6d (f = 0.1316d −1 ). After applying the Everest reduction (third panel), additional statistically significant peaks occur at 25.7d 1190d −1 ), and 6.8d (f = 0.1471d −1 ). The strongest peaks at 28.2d and 25.7d in the polynomial-and Everest-reduced LCs, respectively, are barely within 4-sigma and 3-sigma of the estimated P rot,max (17.5d ± 4.4d). Other pronounced peaks unrelated to stellar rotation estimates could still be related to stellar activity and are investigated further in RVs.

Figure 5. Auto-correlation function (ACF) curves for the polynomial-reduced (top) and
Everest-reduced (bottom) K2 LCs of K2-79 (left) and K2-222 (right). The dark purple points show the ACF. The solid orange line shows the result of the under damped simple harmonic oscillator function fit to the ACF curves, as described in Section 4.1. The ACF shows repeating peaks separated by P quasi , a timescale often related to P rot , with the correlation power of each peak dropping off at a rate related to the evolution timescale of magnetic active regions, τ. The P quasi values returned by successful fits for both targets are similar to each other and the strongest signals in all LC periodograms ( Figure 4). The τ values are also all quite long relative to P quasi . We therefore cannot rule out, for either target, the potential that the ACF signal is the result of instrumental systematics.
(a) K2-79: The fit to the ACF of the polynomial-reduced LC reveals the strongest periodicity at P quasi = 26.6 ± 0.4 days, with an evolution timescale of τ = 72 ± 6 days. The ACF of the Everest-reduced LC reveals the strongest periodicity at P quasi = 26.1 ± 0.2 days, with an evolution timescale of τ = 100 ± 10 days.
(b) K2-222: The fit to the ACF of the polynomial-reduced LC reveals the strongest periodicity at P quasi = 28.2 ± 0.2 days, with an evolution timescale of τ = 100 ± 10 days. The ACF of the Everest-reduced LC returns an amplitude (A) consistent with zero, dominated by the leftover long-period instrumental trend that can be seen by-eye in the LC ( Figure  1b, third panel).
-0.05), and that of K2-222 may be underestimated ([Fe/H] = -0.25). We therefore consult other methods in further attempts to estimate P rot from the RVs.
To investigate the overall strength of stellar activity signals in the RVs, we calculated Spearman correlation coefficients (r s ) between the RVs and each of the measured activity indicators. In order to identify changes from season to season, we calculated three sets of correlation coefficients for each target: one for the full RV data set, one with just the first season (S1) RVs, and one with just the second season (S2) RVs. Other observing seasons contained too few points for a meaningful calculation. Figures 6 and 7 show correlation plots and associated r s values between the RVs and each of the activity indicators for K2-79 and K2-222, respectively. The p-value, a measure of how probable it is that the observed correlation would be produced by a null dataset, is also shown. To compare correlation strengths, we establish the following terminology corresponding to the absolute value of r s : • < 0.20 = no correlation In the case of K2-79, the full RV set has weak correlations with the contrast, CCF area, and BIS. The S1 RVs have no correlation with any activity indicators. The S2 RVs have weak correlations with the contrast and CCF area and a moderate anti-correlation with the BIS, where we see the strongest peak at 11.0d. In the case of K2-222, the full RV set has a moderate correlation with the contrast, where we see the strongest peak at 39.2d, and a weak correlation with the CCF area. The S1 RVs have weak correlations with the contrast and S-index. The S2 RVs have moderate correlations with the contrast and CCF area, where we see the strongest peak at 41.8d in both. In the cases of both targets, the S2 RVs show greater overall correlation with activity indicators than the S1 RVs, which we also see in the periodogram analysis.
We calculate GLS periodograms for the RVs and activity indicators of K2-79 and K2-222 using three sets of data for each: The full dataset, the first season (S1) dataset, and the second season (S2) dataset. The strongest periods in stellar activity signals can change significantly from season to season due the sizes, locations, and distribution of stellar features evolving over time. Table 3 provides more detailed information about the full and seasonal data sets. We also calculate the window functions using the method described in Section 4.1. The resultant sets of periodograms for each target are shown in Figures 8 and 9. In the periodograms we only show up to f = 0.20d −1 because no statistically significant peaks occur beyond that frequency in the case of K2-79, and, in the case of K2-222, statistically significant peaks beyond that frequency are aliases of the 1d peak in the window function.
We address statistically significant peaks in the periodograms, as well as those that occur at the known orbital periods of the transiting exoplanets. Although none of the peaks at the known orbital periods meet the traditional criteria for statistical significance (FAP < 1%), the probability of the peaks being false decreases beyond traditional FAP levels with the confirmed existence of exoplanets at those periods.
For both K2-79 and K2-222, magnetic activity signals generally show a stronger potential presence in S2 than in S1. This is apparent both in the relative strengths of the S1 and S2 seasonal activity indicator periodograms and the S1 and S2 Spearman correlation coefficients, especially considering the much lower number of data points in S2. For each target, the strongest power in the S2 activity indicator periodograms is greater than or comparable to the strongest in S1, in spite of far fewer points in the S2 data set (Figures 8 and 9). In the case of K2-79, none of the S1 activity indicators correlate with the RVs, while three of the S2 activity indicators show at least weak correlation ( Figure 6). In the case of K2-222, two activity indicators correlate with the RVs in both S1 and S2. There is a correlation between the S1 S-index values and RVs (r s = 0.22) where no correlation exists in S2, but it is much weaker than the correlations between contrast and RVs for both S1 (r s = 0.28) and S2 (r s = 0.57), where the S2 Spearman coefficient is significantly higher (Figure 7).
For both targets, at P b specifically, seasonal periodograms also show more potential for activity interference in S2 than in S1. Even though the K2-79 and K2-222 seasonal RV periodograms reveal stronger signals near P b in S1, the seasonal activity indicator periodograms show stronger signals at P b in S2. In the case of K2-79, the S1 RV peak near P b has a FAP ≈ 10%, while the S2 RV peak has no statistical significance (FAP < 99%). Yet, near P b , the S1 activity indicators contain no peaks of significance, while the S2 BIS contains a peak with FAP < 10%, and the S2 contrast contains a peak with FAP ≈ 50%. In the case of K2-222, the S1 RV peak near P b is far stronger (FAP < 10%) than the S2 peak (FAP > 50%), which is expected given the relative sizes of the seasonal data sets. However, near P b , the S1 activity indicators again contain no peaks of significance, while the S2 contrast contains a peak with FAP ≈ 50%.
In the case of K2-79 (Figure 8), no statistically significant peaks occur in the full or either of the seasonal RV periodograms. The two strongest peaks occur at the known orbital period of the planet (P b = 10.99d; f = 0.09099d −1 ) and at f = 0.0882d −1 , an alias of P b with the strongest peak in the full window function ( f w = 0.0028d −1 ). The signal at P b could potentially be combined with an activity signal, based on the corresponding peaks with FAP ≈ 10% and ≈ 50% in the S2 BIS and contrast periodograms, respectively, and its location less than one day from the first harmonics of the estimated P rot,max and P rot,Noyes .
In the case of K2-222 (Figure 9), the strongest peak in the full RV periodogram occurs at f = 0.007d −1 (P = 146.6d), notably not at the period of the known exoplanet. An alias of this signal with the strongest peak in the full window function ( f = 0.0026d −1 ) also occurs at f = 0.0042d −1 . The peak at 146.6d does not appear to be related to the window function or activity indicator periodograms; peaks at nearby frequencies do occur in the activity indicators, but not at the exact same frequency, and the signal is much more significant in the RVs. The 146.6d signal remains the most significant peak (FAP < 1%) when we remove the last three solitary observations and recalculate the periodogram. We therefore consider this signal as a potential second object (P c? ).
The periodogram of the full K2-222 RV set also shows peaks at the period of the known transiting exoplanet (P b = 15.39d; f = 0.06498d −1 ) and at f = 0.0624, an alias of the exoplanet signal with the strongest peak in the full window function ( f w = 0.0026d −1 ). The signal at P b could potentially be combined with an activity signal, based on the corre- Figure 6. Correlation plots and strengths (Spearman r coefficients, r s ) between the K2-79 RVs and various activity indicators (contrast, crosscorrelation function area, full-width at half max., bisector inverse span, and S-index). All RVs and error-bars are plotted in dark purple, with season 1 RVs over-plotted in turquoise and season 2 RVs in magenta. P-value is a measure of how probable it is that the observed correlation is due to random chance. sponding peak with FAP ≈ 10% in the contrast periodogram. We analyze an additional statistically significant peak at f = 0.1229d −1 (P act = 8.137d) that occurs in both the full and S1 periodograms. We believe the signal is most likely a result of stellar activity because of a nearby peak with FAP ≈ 10% in the FWHM periodograms. The signal at P act,1/2 with FAP ≈ 10% may also be related to activity, as it occurs at approximately one-half of P act and less than 1d away from a peak with FAP < 50% in the contrast periodogram. The four prominent peaks surrounding the peak at P act,1/2 are aliases with the largest peak in the S2 window function.
As discussed above, for both K2-79 and K2-222, the S2 RVs suffer more potential interference from stellar activity at P b than the S1 RVs. Since each target's S1 and S2 observations were not taken in the same seasons as the other's, we can rule out any seasonal instrumental effect being the cause of more interference in S2. For K2-79 and K2-222, the small number of observations in S2 will already make for a less precise mass determination, but constructive or destructive interference with stellar activity in the S2 RVs could also affect the accuracy of the exoplanet mass measured from the full RV set. Therefore, in addition to full RV fits, we perform fits to the S1 and S2 RV sub-sets to confirm agreement between exoplanet masses estimated using the full result and at least the S1 result.

SIMULTANEOUS RADIAL VELOCITY/TRANSIT EXOPLANET FIT
Attempts to fit the K2-79 and K2-222 LCs and RVs separately revealed that for both targets, neither of the data sets is sufficient to constrain the orbital eccentricity of the transiting exoplanet. We therefore perform a simultaneous fit to the RVs and flattened LCs (top panels Figures 2 and 3, bottom panels of Figures 1a and 1b) using the publicly available Differential Evaluation Markov-Chain Monte-Carlo (MCMC) Figure 7. Correlation plots and strengths (Spearman r coefficients, r s ) between the K2-222 RVs and various activity indicators (contrast, crosscorrelation function area, full-width at half max., bisector inverse span, and S-index). All RVs and error-bars are plotted in dark purple, with season 1 RVs covered in turquoise and season 2 RVs covered in magenta. P-value is a measure of how probable it is that the observed correlation is due to random chance. software, EXOFAST v2 (Eastman et al. 2019). We run our EXOFAST v2 fits with maxsteps = 50,000 and nthin = 30, running for a maximum of 50,000 recorded steps while recording the 30th step of each walker. The global model used in EXOFAST v2 includes spectral energy distribution and integrated MIST stellar evolutionary models informed by the magnitudes in Table 6 to constrain stellar parameters. In the case of both targets, we utilize isochrones from Choi et al. (2016) to derive a physical star (mass, radius, luminosity, age) based on T eff , [Fe/H], and log g. Table 8 shows the median and 1σ stellar values returned for K2-79 and K2-222.
For the transit fits, we set the longcadence flag to account for the 29.4 min cadence of the K2 data, and we use the default quadratic limb darkening law, based on tables reported in Claret & Bloemen (2011). The transit model includes parameters corresponding to each transiting exoplanet's orbital period (P orb ), time of central transit (T c ), orbital inclination (i), orbital eccentricity (e), orbital argument of periastron (ω), and radius relative to its star (R P /R * ). The RV model is a Keplerian consisting of the following parameters: RV semiamplitude (K), e, P orb , and ω.
The transit fits place strong constraints on P orb and T c . Fit alongside with stellar density, they also drastically reduce the degenerate e-ω parameter space. Pairing the transit and RV fits therefore allows us to break the degeneracy between these two parameters and achieve better estimates of e and ω than were possible with the RVs or transits alone.
For each target, we perform three fits assuming a single planet: one using all available RVs, another with just the S1 RVs, and a third with just the S2 RVs. In the case of K2-222, we perform a second set of these three fits for a 2planet scenario to investigate the nature of the strongest RV periodogram signal at 146.6d. We fix the eccentricity of the potential second exoplanet to zero, taking into account the The full RV set shows weak correlations with the contrast, CCF area, and BIS. The S1 RVs show no correlation with any activity indicators. The S2 RVs show weak correlations with the contrast and CCF area, as well as a strong anti-correlation with the BIS. P b is the period of the known transiting exoplanet. The signal at P b could potentially be combined with an activity signal, based on the corresponding peaks with FAP ≈ 10% and ≈ 50% in the S2 BIS and contrast periodograms, respectively. At P b , the seasonal RVs show a stronger signal in S1, but the seasonal activity indicators show stronger signals in S2, particularly in contrast and BIS. Figure 9. K2-222 GLS periodograms of HARPS-N RVs and activity indicators. The full RV set shows a moderate correlation with the contrast and a weak correlation with the CCF area. The S1 RVs show weak correlations with the contrast and S-index. The S2 RVs show moderate correlations with the contrast and CCF area. P b is the period of the known transiting exoplanet. The signal at P b could potentially be combined with an activity signal, based on the corresponding peak with FAP ≈ 10% in the contrast periodogram. At P b , the seasonal RVs show a stronger signal in S1, but the seasonal activity indicators show stronger signals in S2, particularly contrast. P c? is the period of a potential second companion object, and P act with FAP ≈ 1% is a signal we attribute to stellar activity because of the nearby peaks with FAP ≈ 10% in the FWHM periodograms. The signal at P act,1/2 with FAP ≈ 10% may also be related to activity, as it occurs at approximately one-half of P act and less than 1d away from a peak with FAP < 50% in the contrast periodogram. The four prominent peaks surrounding the peak at P act,1/2 are aliases with the largest peak in the S2 window function. exoplanet's relatively long period and gaps in its phase coverage. Table 7 summarizes the probability distributions of priors applied to the EXOFAST v2 fits for both targets. Priors on stellar parameters were based on the combined result from the Stellar Parameter Classification tool and ARES+MOOG analyses of the HARPS-N spectra (see Section 3 and Table  6). In the case of the seasonal 2-planet fits to K2-222 data, we fix the parameters associated with the potential second exoplanet to the median values returned by the full 2-planet fit. This allows us to confirm agreement between the full and seasonal RV sets of the parameters associated with K2-222b, despite the fact that the seasonal baselines are too short to reliably fit the 146.6d signal of a potential second companion object. Table 8 shows the MCMC posterior median values and 1σ uncertainties of the stellar and exoplanet transit parameters for each target fitted with EXOFAST v2. We also list several other parameters describing the individual transits (depth, duration, impact parameter) and full LCs (baseline flux, variance of transit model residuals). The estimated equilibrium temperature (T eq ) returned by EXOFAST v2 follows Equation 1 of Hansen & Barman (2007) where a is orbital semi-major axis and no albedo and perfect energy redistribution are assumed. The stellar values returned by the EXOFAST v2 fits agree with those adopted from our stellar analysis within 1σ, with the exception of K2-222's stellar radius, which sees a slightly higher discrepancy but agrees well within 1.5σ.
Tables 9 and 10 show the posterior median values and 1σ uncertainties of all other parameters from the EXOFAST v2 fits to K2-79 and K2-222, respectively. In the cases of both targets, the obtained K values for K2-79b and K2-222b from the full RVs and seasonal fits agree within 1σ. In the case of K2-222, the 1-planet and 2-planet fits also return values consistent within 1σ, but median value for the semi-amplitude of the transiting planet, K b , is notably smaller in the 2-planet model.
We find evidence in support of a second companion candidate by comparing the robustness of the 1-planet and 2-planet models using their Bayesian information criterion (BIC) and Akaike information criterion (AIC) values. The BIC value is lower for the 1-planet model full fit, but the AIC value is lower for the 2-planet model by a much larger margin. The AIC and BIC of the 2-planet model are also penalized more due to the higher number of parameters fit than in the 1-planet model. Although the AIC and BIC values provide tentative support in favor of a 2-planet model, we adopt the 1-planet solution for now, considering the significant phase gaps in the coverage of the 147.5d period returned by our fit. Figures 10 and 11 show the phase-folded LCs of K2-79b and K2-222b with the best-fit transit models over-plotted. Figures 12 and 13 show the phase-folded RV curves of each transiting exoplanet with the best fit 1-planet Keplerian model over-plotted. Figure 14 shows the phase-folded RV curves of K2-222b and the candidate K2-222c with their best fit 2-planet Keplerian model over-plotted.
Finally, in addition to the simultaneous RV and LC fits with EXOFASTv2, we also fitted the RV data with the RadVel software (Fulton et al. 2018) to see if we could improve the RV solutions using Gaussian Processes (GP) regression to model the stellar activity (Haywood et al. 2014). We detail these fits in the next section.

RADIAL VELOCITY FIT WITH GAUSSIAN PROCESSES REGRESSION
The following RV fits including GPs were not included in the calculation of any final mass results. However, in this section we detail our analysis and results for completeness.
In addition to the simultaneous RV and LC fits, we also use the RadVel software (Fulton et al. 2018) to perform MCMC fits on just the HARPS-N RVs, but this time including a quasi-periodic Gaussian Processes (GP) regression to model a stellar activity signal in addition to the exoplanet Keplerian (Haywood et al. 2014). We once again perform fits on the full and seasonal RV sets, and in the case of K2-222, with both 1-planet and 2-planet models.
The RadVel model fits for the same Keplerian and instrumental parameters fit by EXOFAST v2. In the 2-planet models we once again fix the eccentricity of the potential second exoplanet and, in the seasonal fits, fix all of its associated parameters to the values returned by the full fit. We place priors on P b , T C,b , and e b for both targets based on posteriors from our EXOFAST v2 fits. We set Jeffreys priors on K b values with an upper limit of 10m/s to explore parameter space around the ≈ 6m/s standard deviation of both RV sets. We place a Gaussian prior on σ HARPS−N , set equal to the distribution of RV errors ( RVerr , σ RVerr ). Table 11 summarizes the priors applied to exoplanet parameters in our RadVel fits.

Symbol
Parameter and Units All S1 S2 Exoplanet Parameters:         where k(t,t') is the correlation weight between observations taken at times t and t'. A is the mean amplitude of the activity signal, τ is a hyper-parameter related to the evolution timescale of activity features, P quasi is related to the stellar rotation period, and η describes the level of high-frequency variation expected within a single stellar rotation and is related to the average distribution of activity features on the surface of the star.
We apply Jeffreys priors on the GP stellar amplitudes with the same limits and motivation applied to the K b priors, and also place a Jeffreys prior on τ investigating a wide range of values between 2d and 100d. As mentioned above, η is physically related to the average distribution of magnetic active regions on the stellar surface. Models by Jeffers & Keller (2009) demonstrated that even highly complex activity distributions will average to just two to three large active regions in a given rotation. We therefore assign a prior distribution of η = 0.25 ± 0.025, which allows for two to three local minima or maxima per rotation (Haywood et al. 2014;Grunblatt et al. 2015;López-Morales et al. 2016).
For each target, we ran fits with various priors on P quasi , first with a Jeffreys prior of 2-40d for each. Based on the posterior distributions returned by those fits, we ran additional fits with Jeffreys priors of 10-40d and 2-10d for K2-79 and K2-222, respectively. We also used several Gaussian priors centered around various P rot estimates with a 20% standard deviation. For K2-79 we ran fits with Gaussian priors on P quasi centered around 26.35d, 23.8d, and 11.0d based on P rot,Noyes and the LC ACFs, the v sin i-calculated P rot,max , and the BIS periodogram, respectively. For K2-222 we ran fits with Gaussian priors on P quasi centered around 28.2d, 17.5d, and 8.0d, based on the polynomial-reduced LC ACF, the v sin i-calculated P rot,max , and P rot,Noyes plus the RV and FWHM periodogram peaks, respectively. Final exoplanet parameters returned by successful fits agreed within 1-σ no matter which priors were used on the GP hyper-parameters. All fits including a Gaussian prior on P quasi just returned the input prior in their posterior, so we report values from the fits including 10-40d and 2-10d Jeffreys priors for K2-79 and K2-222, respectively. Table 11 summarizes the priors placed on GP hyper-parameters in our reported RadVel fits.
The majority of the RadVel fits failed to converge on a solution for the stellar activity signal. RVs for both K2-79 and K2-222 were not optimally scheduled to characterize stellar activity signals, with never more than a single observation taken per night and several multiple-day gaps in each of the analyzed seasons. It is therefore unsurprising that our models have a difficult time characterizing the stellar activity signal. We define successful fits as those that converge for values of A and P quasi , at least. Even in successful fits, posterior distributions show τ pushing toward the smallest possible value, consistent with a highly unstable P quasi signal. This was the case for all fits including a GP-generated model with a quasiperiodic kernel. We include a corner plot for each target in the Appendix for those wishing to view the posteriors from the reported Radvel fits in further detail (Figures 16 and 17). We also tested the performance of GP fits with a simpler squared-exponential kernel, but found they show the same behavior in the τ posteriors as the fits with quasi-periodic kernels and are even less likely to converge on a value for stellar amplitude, A. Table 12 shows final parameter values returned by successful RadVel GP fits. The ≈24d P quasi signal returned for K2-79 is consistent with that returned by the LC ACFs and close in value to the largest peaks in the LC periodogram, P rot,Noyes , and the P rot,max calculated from v sin i. A second smaller peak that shows up in the posterior at 12-13d and the 11d peak in the BIS periodogram both sit near the first harmonic of these signals, as well. We therefore estimate that K2-79 has a P rot in the range 21-26 days. The ≈4d P quasi signal returned for K2-222 is consistent with the first harmonics of the estimated P rot,Noyes and the ≈8d peaks in the RV and FWHM periodograms. While we believe the most likely rotation period of this star is ≈8d, estimates of P rot for K2-222 have varied greatly from method to method, and we are less confident in the estimate than that for K2-79. The reported jitter values are low relative to the average RV error for both targets, giving us concern that the GP models could be overfitting noise that is unrelated to stellar activity. However, the exoplanet semi-amplitudes returned by the successful GP fits all have 1-sigma agreement with those in their corresponding EXOFASTv2 fits. We therefore list the RadVel results as further evidence in support of the measured K2-79b and K2-222b semi-amplitudes, but adopt final results for the exoplanet characteristics from our EXOFASTv2 fits to the full RV data sets. 7. RESULTS AND DISCUSSION Table 13 summarizes the fitted and derived parameters of K2-79b, as well as K2-222b in both the 1-planet and 2-planet cases. Figure 15 shows their location in mass-radius space along with the population of confirmed exoplanets with T eq between 300K and 3000K, 1-sigma mass errors less than 50%, 1-sigma radius errors less than 20%, and orbiting stars with radii between 0.6R and 2.0R . With radii of 4.09 and 2.35 R ⊕ , K2-79b and K2-222b exist at two important locations: The sub-Saturnian desert and at the upper edge of the radius valley, respectively. Zeng et al. (2021, accepted) provides evidence for two new dimensionless parameters that can be used to describe exoplanet compositions. The first, ζ ≡ (R p /R ⊕ )(M p /M ⊕ ) −1/4 , can distinguish between three small exoplanet populations: rocky planets composed of varying ratios of silicates and metals (ζ ≈ 1), water worlds dominated by significant amounts of water and ice (ζ ≈ 1.4), and Neptune/Uranuslike planets with ice-dominated cores hosting small gaseous envelopes (ζ ≈ 2.2) (see histogram in Figure 15). ζ is physically related to the chemical composition of the exoplanet core as illustrated by the ternary plot in Figure 15). The second parameter, z ≡ ∈envelope dP ρ G·M ⊕ R ⊕ , characterizes both the amount and the temperature of the added envelope. We will refer to ζ and z in the following two sub-sections when discussing potential compositions of K2-79b and K2-222b.

K2-79b
For K2-79b, we measure a mass and radius of M b = 11.8 ± 3.6 M ⊕ and R b = 4.09 ± 0.17 R ⊕ , which when combined yield a bulk density of ρ b = 0.17 ± 0.06 ρ ⊕ . K2-79b falls in the sub-Saturnian desert, where potential compositions are degenerate. If the planet formed water-poor, its predicted composition is a rocky core with an approximately 10% H/He envelope, by mass (Lopez & Rice 2018). On the other hand, if the planet formed water-rich, the core is likely a rocky and icy combination surrounded by an envelope of mostly H/He gas, and potentially H 2 O vapor (Zeng et al. 2019). K2-79b receives a relatively high average irradiation level of 180 ± 8 F ⊕ (approximated for circular orbit) and has a low surface gravity of 0.71 ± 0.22 g ⊕ . Due to the low molecular weight and escape velocity of H 2 and He, it would be difficult for this exoplanet to maintain a significant water-poor envelope. According to simulations by Lopez & Rice (2018), K2-79b would not retain a purely H/He envelope over its greaterthan-5Gyr lifetime, supporting the presence of heavier gases .
With a ζ value of 2.21, K2-79b is likely a Uranus-analog. A slightly H 2 O-enriched atmosphere could explain K2-79b's ability to maintain its envelope while receiving such a high level of irradiation. Assuming the exoplanet hosts a noncloudy envelope, a pure hydrogen composition would produce a calculated change in transit depth of ∆d = 236.8 parts-per-million (ppm) over five atmospheric scale-heights in transit spectroscopy signals. Assuming a water-enriched Table 11. Prior probability distributions applied to Radvel MCMC fits to K2-79 and K2-222 RVs. The two numbers in each prior are the mean and standard deviation, respectively, for Gaussian prior distributions, or the lower and upper bounds, respectively, for uniform Jeffreys prior distributions. Upper bounds on the eccentricity priors were taken from the posterior distributions returned by EXOFAST v2 for each respective fit. The Gaussian prior on σ HARPS−N is set equal to the distribution of RV errors ( RVerr , σ RVerr ).  Table 13. Final Derived Exoplanet Characteristics. ζ is physically related to the chemical composition of an exoplanet core (see Section 7 and Figure 15).

K2-222b
As described in Section 5, there appears to be a second signal in the RVs of K2-222 with a period of 147.5d. AIC and BIC values provide tentative support in favor of that signal being Doppler-induced (Tables 9 and 10), but we adopt the 1-planet solution for now, considering the significant phase gaps in the coverage of the 147.5d period returned by our fit. Assuming that the 147.5 signal is Doppler-induced, it would correspond to an object with minimum mass Msin i = 19.3 ± 4.2 M ⊕ . Including that second signal in our RV fits yields a mass for K2-222b of M b = 6.0 ± 1.9 M ⊕ . In the case of the adopted single-planet model, we measure a final mass and radius of M b = 8.0 ± 1.8 M ⊕ and R b = 2.35 ± 0.08 R ⊕ , respectively, for K2-222b. Combined, this mass and radius yield a final bulk density of ρ b = 0.62 ± 0.15 ρ ⊕ . K2-222b sits just above the upper edge of the radius valley where potential compositions are again degenerate. If it formed water poor, its predicted composition is a rocky core and an approximately 1% H/He envelope, by mass (Lopez & Rice 2018). If the exoplanet instead formed water-rich, the core is likely equal parts rock and ice hosting an envelope of H 2 O equal to half the core mass (Zeng et al. 2019). K2-222b receives an irradiation of 95 ± 5 F ⊕ and has a surface gravity of 1.45 ± 0.34 g ⊕ . According to 5 Gyr simulations and its estimated minimum age, K2-222b could potentially retain a pure H/He envelope (Lopez & Rice 2018). However, its ζ value of 1.40 suggests that K2-222b is most likely a water world ( Figure 15). Assuming an envelope composed of pure hydrogen, we estimate a calculated change in transit depth for K2-222b of ∆d = 158.4 ppm over five scale-heights in transit spectroscopy signals.

Utilizing seasonal RV analyses
With our seasonal RV analyses we were able to measure the masses of K2-79b and K2-222b, in spite of the potentially challenging interference from stellar activity near P orb in both cases. By analysing RV and activity indicator periodograms separately for each available season of data, we determined the stellar activity signal to be evolving from season to season, and identified S1 RVs as suffering from less interference at P rot than S2 RVs. Fitting the full and seasonal RVs separately, we find that all obtained planet parameters agree within errors, providing confirmation that our final mass estimates are not being hindered by stellar activity interference in S2 RVs.
Many RV surveys are now optimizing observations and data reduction approaches to mitigate the effect of stellar activity in the measurement of planet-induced RVs (e.g. López-Morales et al. 2016;Haywood et al. 2020;Miklos et al. 2020;de Beurs et al. 2020;Collier Cameron et al. 2020). However, the problem of overlapping exoplanet and stellar activity signals remains difficult to solve in systems where the orbital period of the exoplanet overlaps with a periodic signal from stellar rotation. We suspect the key observing strategy with such targets will be to observe with high-cadence (at least once per night) over at least two separate seasons, knowing that the planet signal does not change, but the average spot distribution, and therefore stellar activity signal, does. In cases of highly stable activity signals, several seasons of separation between observations may be favorable in order to observe different average spot distributions. We hope that future simulations will confirm or rule-out these potential strategies. However, much of the archival RV data was collected without the inclusion of optimal strategies for mitigating stellar activity signals. In those cases, an analysis approach like the one described in this paper might help to extract exoplanet signals. Figure 15. The scatter plot (left) shows the placement of K2-79b and K2-222b in mass-radius space among other confirmed small exoplanets (colored points) with T eq between 300K and 3000K, 1-sigma mass errors less than 50%, 1-sigma radius errors less than 20%, and orbiting stars with radii between 0.6R and 2.0R (NASA Exoplanet Archive). Uranus (U) and Neptune (N) are indicated by grey dots. The red curve represents Fe-metal compositions of 100% iron, the orange curve represents pure silicates (MgSiO 3 ) composition, and the blue curve represents pure H 2 O compositions. The contours above the blue curve depict different values of the z parameter discussed in Section 7, which characterizes both the amount and the temperature of the added envelope (Section 7; Zeng et al., 2021 (accepted)). The two sets of contours between the three colored curves represent equally spaced mass-proportions of the two compositions they lie between. The histogram (top right) shows the distribution of ζ values, also discussed in Section 7, for the same population of planets in the scatter plot, with K2-79b and K2-222b shown in blue. The histogram shows clearly a bi-modal distribution, with the first peak (ζ ∼ 1) corresponding to purely rocky compositions, while the second peak (ζ ∼ 1.4) is consistent with ice-dominated compositions. The histogram also shows a potential tri-modal distribution, with the third peak (ζ ∼ 2.1 − 2.4) consistent with Neptune-like planets. K2-222b is consistent with the second peak, and K2-79b is consistent with the third peak. The ternary plot (bottom left) shows how ζ relates to various potential exoplanet compositions (Zeng et al., 2021 (accepted)). The solid diagonal lines are contours of ζ, which range from slightly less than 0.9 in one extreme (pure Fe-metals) to slightly more than 1.4 in the other extreme (pure-Ices). The dashed straight line is a fixed ratio of Fe-metals/Silicates, which is a theoretical prediction from maximum collision stripping models. Figure 16. Corner plot showing posterior distributions of hyper-parameters from the successful K2-79 Radvel MCMC fits including a stellar activity signal modelled with quasi-periodic Gaussian Processes regression. The Radvel hyper-parameters 1-4 represent A, τ, P quasi , and η. The ≈24d P quasi signal is consistent with that returned by the LC ACFs, as well as close in value to the largest peaks in the LC periodogram, P rot,Noyes , and the P rot,max calculated from v sin i. The second smaller peak that shows up in the posterior at 12-13d and the 11d peak in the BIS periodogram both sit near the first harmonic of the 24d signal. We therefore estimate that K2-79 has a P rot in the range 21-26 days. Radvel MCMC fits including a stellar activity signal modelled with quasi-periodic Gaussian Processes regression. The Radvel hyper-parameters 1-4 represent A, τ, P quasi , and η. We only show the corner plot for the fit utilizing all RVs because that from the fit to only S1 RVs contains no notable differences. The ≈4d P quasi signal returned for K2-222 is consistent with the first harmonics of the estimated P rot,Noyes and the ≈8d peaks in the RV and FWHM periodograms. While we believe ≈8d is the most likely rotation period of this star, estimates of P rot for K2-222 have varied greatly for certain methods, and we are less confident in this estimate than that for K2-79.