K2-291b: A Rocky Super-Earth in a 2.2 day Orbit^{* †}

We collected 75 radial velocity measurements of K2-291 (Table 2) with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the Keck I Telescope on Maunakea and the High Accuracy Radial velocity Planet Searcher in the Northern hemisphere (HARPS-N; Cosentino et al. 2012) on the Telescopio Nazionale Galileo in La Palma (Table 2). HARPS-N is an updated version of HARPS at the European Southern Observatory (ESO) 3.6 m (Mayor et al. 2003).

Table 3.  Stellar Parameters

Parameter	Name (units)	Value
Name and Magnitudea
EPIC		247418783
UCAC ID		558-013367
2MASS ID		05054699 + 2132552
Gaia DR2		3409148746676599168
HD		285181
Kp	mag	9.89
R	mag	9.84 ± 0.14
J	mag	8.765 ± 0.032
K	mag	8.35 ± 0.02
V	mag	10.01 ± 0.03
K2		291
Locationb
R.A.	Right ascention (deg)	05 05 46.991
Decl.	Declination (deg)	+21 32 55.021
π	Parallax (arcsec)	0.011076 ± 6.03e–05
d	Distance (pc)	${90.23}_{-0.46}^{+0.51}$
Stellar Properties
Av	Extinction (mag)	0.11740 ± 0.00061
R*	Radius (R⊙)	${0.899}_{-0.033}^{+0.035}$
M*	Mass (M⊙)	0.934 ± 0.038
L*	Luminosity (L⊕)	${0.682}_{-0.016}^{+0.014}$
Teff	Effective temperature (K)	5520 ± 60
log(g)	Surface gravity (cgs)	4.50 ± 0.05
[Fe/H]	Metallicity (dex)	0.08 ± 0.04
vsini	Rotation (km s−1)	<2.0
log(age)	Age (yr)	${9.57}_{-0.49}^{+0.30}$
log $({R}_{\mathrm{HK}}^{{\prime} })$	Chromospheric activity	−4.726

Notes.

^aMAST. ^bGaia DR2 (Gaia Collaboration et al. 2018).

Download table as:  ASCIITypeset image

We obtained 50 measurements with HIRES between 2017 August and 2018 February. These data were collected with the C2 decker with a typical S/N of 150/pixel (125k on the exposure meter, ∼10 minute exposures). An iodine cell was used for wavelength calibration (Butler et al. 1996). We also collected a higher resolution template observation with the B3 decker on 2017 September 6 with 08 seeing. The template was a triple exposure with a total S/N of 346/pixel (250k each on the exposure meter) without the iodine cell. See Howard et al. (2010) for more details on this data collection method.

We obtained 25 measurements with HARPS-N between 2017 November and 2018 March as part of the HARPS-N Collaboration's Guaranteed Time Observations (GTO) program. The observations follow a standard observing approach of one or two observations per GTO night, separated by 2–3 hr. The spectra have S/Ns in the range of S/N = 35–99 (average S/N = 66), seeing and sky transparency dependent, at 550 nm in 30 minute exposures. This separation was designed to well sample the planet's orbital period and to minimize the stellar granulation signal (Dumusque et al. 2011).

The HIRES data reduction and analysis followed the California Planet Search method described in Howard et al. (2010). The HARPS-N spectra were reduced with version 3.7 of the HARPS-N Data Reduction Software, which includes corrections for color systematics introduced by variations in seeing (Cosentino et al. 2014). The HARPS-N radial velocities were computed with a numerical weighted mask following the methodology outlined by Baranne et al. (1996) and Pepe et al. (2002). The resultant radial velocities are presented in Table 2 and in Figure 6.

The HIRES data were collected with three consecutive exposures of 10 minutes each to well sample the stellar p-mode (acoustic) oscillations which occur on a timescale of a few minutes. The HARPS-N data were collected in single observations. Multiple exposures per night were frequently taken, separated by a few hours, to better sample the planet orbital period.

We derived the stellar parameters by combining constraints from spectroscopy, astrometry, and photometry (Table 3). The methodology is described in detail in Fulton & Petigura (2018) and summarized in the following paragraphs. We used the HIRES template spectrum to determine the parameters described below. A comparison analysis performed on the HARPS-N data resulted in 3σ consistent parameters.

Stellar radius is derived from the Stefan Boltzman Law given an absolute bolometric magnitude M_bol and an effective temperature. We derived stellar effective temperature T_eff, surface gravity log(g), and metallicity [Fe/H] by fitting our iodine-free template spectrum using the Spectroscopy Made Easy³⁶ (SME) spectral synthesis code (Valenti & Piskunov 2012) following the prescriptions of Brewer et al. (2016). Stellar mass is then calculated using the package isoclassify³⁷ (Huber et al. 2017). We then derived bolometric magnitudes according to

$$\begin{eqnarray}&&{M}_{\mathrm{bol}}={m}_{K}-{A}_{k}-\mu +{BC},\end{eqnarray} \tag{ 1 }$$

where m_K is the apparent K-band magnitude, A_k is the line-of-sight K-band extinction, μ is the distance modulus, and BC is the K-band bolometric correction. In our modeling, constraints on m_K come from 2MASS (Skrutskie et al. 2006) and constraints on μ come from the Gaia DR2 parallax measurement (Gaia Collaboration et al. 2018). We derived BC by interpolating along a grid of T_eff, log g, [Fe/H], and A_V in the MIST/C3K grid³⁸ (Choi et al. 2016; Dotter 2016; C. Conroy et al. 2019, in preparation). To find A_k, we first estimate A_v from a 3D interstellar dust reddening map by Green et al. (2018), then convert to A_k using the extinction vector from Schlafly et al. (2018).

The stellar rotation velocity vsini, is computed using the SpecMatch-Syn code (Petigura 2015). Due to the resolution of the instrument the code has been calibrated down to 2 km s⁻¹, and values smaller than this should be considered as an upper limit. Although we measured a value of 0.2 km s⁻¹, we adopt vsini < 2 km s⁻¹. To determine the chromospheric activity measurement log$({R}_{\mathrm{HK}}^{{\prime} })$, we measured the flux in the calcium H and K lines relative to the continuum as described in Isaacson & Fischer (2010). We measured the flux in the Calcium H and K lines relative to the continuum. Small differences are noted as S_HK and are tracked to determine if the stellar activity is influencing the radial velocity data.

We searched for stellar companions and blended background stars to K2-291 since these stars could contaminate the stellar flux in the K2 aperture, resulting in an inaccurate planet radius and affecting our radial velocity data if bound.

We searched for secondary spectral lines with the ReaMatch algorithm (Kolbl et al. 2015). This algorithm searches for faint orbiting companion stars or background stars that are contaminating the spectrum of the target star. There are no companions detected down to 1% of the brightness of K2-291 with a radial velocity offset of less than 10 km s⁻¹.

We further looked for stellar companions to K2-291 with AO. We observed K2-291 on 2017 August 3 UT with NIRC2 on the Keck II AO system (Wizinowich et al. 2000). We obtained images with a three-point dither pattern in the Br-γ and J_cont filters at an air mass of 1.71. We do not detect any companions down to ΔBr-γ = 6.41 at 103 as shown in Figure 2.

Zoom In Zoom Out Reset image size

Complementary follow-up observations were taken on 2017 September 7 UT with PHARO-AO on the Hale telescope (Hayward et al. 2001). We obtained images with a five-point dither pattern in the Br-γ filter at an air mass of 1.04. The conditions of our observations allowed us to be sensitive down to ΔBr-γ = 8.05 at 105 as shown in Figure 2 and confirm that we detect no companions to K2-291 above our limits; this also suggests that the transit signal detected is not by a background eclipsing binary.

Stars produce intrinsic radial velocity variations due to their internal and surface processes that can be mistaken as planetary signals. The timescales of these radial velocity variations range from a few minutes or hours (p-modes and granulation) to days or years (stellar rotation and large-scale magnetic cycle variations; Schrijver & Zwaan 2000).

We examine the K2 light curve periodicity (Figure 3) with a Lomb–Scargle periodogram from scipy (Jones et al. 2001) and attribute the clear signal at 18.1 days to rotational modulation of stellar surface features (e.g., spots). There is a secondary peak at half of the strongest peak, and no other significant peaks.

Zoom In Zoom Out Reset image size

One must consider these timescales when planning radial velocity data collection and analysis to adequately average out or monitor these signals (Dumusque et al. 2011). As described in Section 3.1, we chose the exposure time, spacing, and number of exposures to reduce the effects of p-modes and granulation. We investigated the potential radial velocity signal from the stellar rotation by examining the calcium II H and K lines (S_HK, Table 2) in the HIRES and HARPS-N data (Isaacson & Fischer 2010).

We found a clear signal in both the S_HK and radial velocity data that matches the timescale of the rotation period of K2-291 (Figure 4), as determined from the K2 light curve; therefore we need to account for this signal in our radial velocity analysis.

Zoom In Zoom Out Reset image size

We then estimated the correlation coefficient between the measured radial velocity and activity indexes. Due to different zero-points in both radial velocity and S_HK, we performed the analysis for the two instruments, HARPS-N and HIRES, independently. From the calculation of the correlation coefficient value and the knowledge of the sample size, p-value analysis is often used to reject the null hypothesis of non-correlation at a given significance level. We calculated the p-value for both data sets using scipy.stats.pearsonr (Jones et al. 2001). The HARPS-N radial velocity and S_HK data have a p-value of 0.01 allowing us to reject the null hypothesis, therefore suggesting a correlation. The HIRES data, however, have a p-value of 0.45 which does not support a correlation.

To check any potential flaws in the p-value test we also used the Bayesian framework described in Figueira et al. (2016) that allows us to estimate the probability distribution of the coefficient, providing important insight on the correlation presence. This framework calculates the Pearson's correlation coefficient to test for the presence of a linear correlation, and the Spearman's rank to test for the presence of a monotonic correlation.

On HARPS-N data we obtain a Pearson's correlation coefficient of 0.56 with a 95% highest probability density (HPD) between the values [0.29, 0.79], and a Spearman's rank of 0.63 with 95% HPD of [0.39, 0.83]. This shows that not only the correlation coefficient is large but that its distribution populates essentially positive correlation values. As such, the correlation is strong and significant, both in linear and monotonic terms. On the other hand, for HIRES we obtain an average value of 0.10 with 95% HPD of [−0.17, 0.35] and 0.13 with 95% HPD of [−0.12, 0.39] for Pearson's correlation coefficient and Spearman's rank, respectively. The correlation coefficients are low in absolute value and distributed from negative to positive values; its distribution does not support the presence of a correlation. Different instrument properties, such as wavelength ranges and resolution, may explain the differences in the S_HK values and correlation strengths.

We first searched for K2-291 b in the combined HIRES and HARPS-N data sets without any priors from our transit analysis to provide an independent planet detection. The radial velocity data sets from HIRES and HARPS-N are merged using the γ values reported in Table 4 to adjust for their different zero-points in this search. The 75 data points thus obtained are then analyzed in frequency (Figure 5) using the Iterative Sine-Wave fitting (Vaníček 1971), by computing the fractional reduction in the residual variance after each step (reduction factor). This is an iterative process, and peaks should be directly compared within an iteration but not between them. The power spectrum immediately supplies the rotational period at f = 0.055 day⁻¹ (top panel), corresponding to P_rot = 18.1 days. The light curve is very asymmetrical (Figure 3) and therefore signals are visible at the harmonics values, f, 2f, 3f, and 4f. We were successful in detecting the expected frequency of the planet signal at f = 0.45 day⁻¹ after including the stellar rotational frequencies in a simultaneous fit (middle panel). We also searched for any other additional signals, but we did not detect any clear peaks (bottom panel). Indeed, the interaction of the noise with the spectral window (insert in the top panel) prevents any reliable further identification.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Radial Velocity Fit Parameters

Parameter	Name (Units)	Value
Planet Parameters
Pb	Period (days)	${2.225172}_{-7{\rm{e}}\mbox{--}05}^{+6.9{\rm{e}}\mbox{--}05}$
Tconjb	Time of conjunction	${2457830.0616}_{-0.0010}^{+0.0011}$
	(BJDTDB)
eb	Eccentricity	≡0.0
ωb	Argument of periapse	≡0.0
	(radians)
Kb	Semi-amplitude (m s−1)	3.33 ± 0.59
Mb	Mass (M⊕)	6.49 ± 1.16
ρb	Density (g cm−3)	${8.84}_{-2.03}^{+2.50}$
Other Parameters
γHIRES	Mean center-of-mass	−3.5 ± 3.2
	velocity (m s−1)
γHARPS-N	Mean center-of-mass	${25126.2}_{-3.5}^{+3.4}$
	velocity (m s−1)
$\dot{\gamma }$	Linear acceleration	≡0.0
	(m s−1 day−1)
$\ddot{\gamma }$	Quadratic acceleration	≡0.0
	(m s−1 day−2)
σHIRES	Jitter (m s−1)	${1.85}_{-0.37}^{+0.43}$
σHARPS-N	Jitter (m s−1)	${1.43}_{-0.67}^{+0.85}$
η1, HIRES	Amplitude of covariance	${8.45}_{-1.65}^{+2.21}$
	(m s−1)
η1, HARPS-N	Amplitude of covariance	${8.59}_{-1.77}^{+2.23}$
	(m s−1)
η2	Evolution timescale	${26.09}_{-3.62}^{+3.50}$
	(days)
η3	Recurrence timescale	${18.66}_{-0.79}^{+0.95}$
	(days)
η4	Structure parameter	0.41 ± 0.04

Download table as:  ASCIITypeset image

After this initial, transit-blind radial velocity analysis, we analyzed the radial velocity data using RadVel³⁹ (Fulton et al. 2018). RadVel is an open source Python package that models Keplerian orbits to fit radial velocity data by first performing a maximum-likelihood fit to the data and then determining errors through a MCMC analysis. We use the default number of walkers, number of steps, and criteria for burn-in and convergence as described in Fulton et al. (2018).

A single planet at an orbital period of P = ${2.225177}_{-6.8{\rm{e}}\mbox{--}5}^{+6.6{\rm{e}}\mbox{--}5}$ days was found in the K2 photometry (Section 2); we include a Gaussian prior on the orbital period P and time of transit T_conj from the K2 data (Table 1). We first modeled this system using a one-planet fit including a constant offset for each data set γ. This fit results in a semi-amplitude for the 2.2 day planetary signal of K_p = 3.1 ± 1.7 m s⁻¹.

Next, we tested models including an additional trend ($\dot{\gamma }$), curvature ($\ddot{\gamma }$), and eccentricity (e, ω). We used the Bayesian information criteria (BIC) to evaluate if the fit improved sufficiently to justify the additional free parameters; a positive ΔBIC indicates an improved fit. The trend is the only additional parameter which has a noticeable ΔBIC (ΔBIC = 8.29); the trend is $\dot{\gamma }$ = 0.07 ± 0.02 m s⁻². There is nearly no change for the curvature (ΔBIC = 0.84) or eccentric (ΔBIC = −1.90) cases. All three additional parameters result in semi-amplitudes within 1σ of the circular fit.

Stellar activity of K2-291 has an appreciable effect on our measured radial velocities. As discussed in Section 3.4, there is a periodic signal in the radial velocity data that matches both the stellar rotation period determined from the K2 data and the periodicity in the calcium H and K lines (S_HK). We modeled this stellar signal simultaneously with our planet fit using a GP with the default GP model available in RadVel (S. Blunt et al. 2019, in preparation). Gaussian process (GP) regression is a nonparametric statistical technique for modeling correlated noise in data. GP regression enables the determination of physical parameter posterior distributions with uncertainties that reflect the confounding effects of stellar activity noise (e.g., Haywood et al. 2014; Grunblatt et al. 2015; López-Morales et al. 2016).

Stellar noise characteristics in GP models are controlled by a kernel function with one or more hyperparameters, but radial velocity data are often too sparse to confidently determine the values of these hyperparameters (see Faria et al. 2016 for a counterexample). To address this problem, authors in the literature use other data sources to constrain the values of the hyperparameters, then incorporate this information into the radial velocity fit as priors on the hyperparameters (e.g., Haywood et al. 2014; Rajpaul et al. 2015). In this paper, we constrain the values of the hyperparameters in our GP model using K2 photometry.

We modeled the correlated noise introduced from the stellar activity using a quasi-periodic GP with a covariance kernel of the form

$$\begin{eqnarray}&&k(t,t^{\prime} )={\eta }_{1}^{2}\ \exp \left[-\displaystyle \frac{{\left(t-t^{\prime} \right)}^{2}}{{\eta }_{2}^{2}}-\displaystyle \frac{{\sin }^{2}\left(\tfrac{\pi (t-t^{\prime} )}{{\eta }_{3}}\right)}{{\eta }_{4}^{2}})\right],\end{eqnarray} \tag{ 2 }$$

where the hyper-parameter η₁ is the amplitude of the covariance function, η₂ is the active-region evolutionary timescale, η₃ is the period of the correlated signal, η₄ is the length scale of the periodic component (Haywood et al. 2014; López-Morales et al. 2016).

We explore these hyperparameters for this system by performing a maximum-likelihood fit to the K2 light curve with the quasi-periodic kernel (Equation (2)) then determine the errors through an MCMC analysis. We find γ_K2 = ${1567969.00}_{-1830.87}^{+1766.12}$, σ = 54.60 ± 9.57, η₁ = ${4429.95}_{-673.95}^{+897.65}$, η₂ = ${25.18}_{-3.59}^{3.50}$, η₃ = ${19.41}_{-1.14}^{+0.68}$, and η₄ = ${0.42}_{-0.03}^{+0.04}$. This stellar rotation period (η₃) is consistent with the results of our periodogram analysis in Section 3.4.

We then perform a radial velocity fit including a GP to account for the affects of stellar activity on our measurements. We model our GP as a sum of two quasi-periodic kernels, one for each instrument as HIRES and HARPS-N have different properties, such as wavelength ranges, that could alter the way that stellar activity affects the data. Each kernel includes identical η₂, η₃, and η₄ parameters but allows for different η₁ values.

We inform the priors on these hyperparameters from the GP light curve fit (Section 4.3). η₁ is left as a free parameter as light curve amplitude cannot be directly translated to radial velocity amplitude. η₂ has a Gaussian prior describing the exponential decay of the spot features (25.18 ± 3.59). η₃ has a Gaussian prior constrained from the stellar rotation period (19.14 ± 1.14). η₄ constrains the number of maxima and minima per rotation period with a Gaussian prior (0.42 ± 0.04), as described in López-Morales et al. (2016). We do not include a prior on the phase of the periodic component of the stellar rotation because spot modulation tends to manifest in radial velocity data with a relative phase shift.

The planet parameters derived from our GP analysis are consistent with our original, non-GP fit within 1σ. The uncertainty on the semi-amplitude of the planet signal has decreased by a factor of three to K_p = 3.33 ± 0.59 m s⁻¹. We then investigate the inclusion of additional parameters with our GP fit. All of the tested models increased the BIC value, therefore none of them justified the additional parameters. We adopt the model including the GP with no additional parameters as our best fit, all other models have results within 1σ; our best-fit parameters are listed in Table 4 and the fit is shown in Figure 6.

We choose to include a GP in our analysis to improve the accuracy of our results by including the affects of stellar activity. The GP was able to also improve the precision of the mass measurement by a factor of three since the planet orbital period is far from the stellar rotation period, both periods were well sampled with the data, and the stellar activity is dominated by the rotation signal.

We perform an independent radial velocity analysis using the PyORBIT code⁴⁰ (Malavolta et al. 2016, 2018) with results well within 1σ with respect to those reported in Table 4.

Planet compositional models and radial velocity observations of small Kepler planets have shown a dividing line between super-Earth and sub-Neptune planets at 1.5–2 R_⊕ (Lopez & Fortney 2014; Marcy et al. 2014; Weiss & Marcy 2014; Dressing & Charbonneau 2015; Rogers 2015). Kepler planet radii also display a bimodality in sub-Neptune-sized planets that matches the location of this divide (Fulton et al. 2017). K2-291 b is near the inner edge of the divide (${1.589}_{-0.072}^{+0.095}$ R_⊕), which makes its composition particularly interesting.

As shown in the mass–radius diagram (Figure 7), the composition of K2-291 b is consistent with a silicate planet containing an iron core and lacking substantial volatiles (Zeng et al. 2016). We investigated its composition further using Equation (8) from Fortney et al. (2007), which assumes a pure silicate and iron composition, to estimate the mass fraction of each. For our mean mass and radius, the mass fraction of silicates is 0.61 and the mass fraction of iron is 0.39, similar to the 0.35 iron core mass fraction of the Earth. For a high gravity case (1σ low radius, 1σ high mass), the mass fraction of silicates would be 0.39. For a low gravity case (1σ high radius, 1σ low mass), the mass fraction of silicates would be 0.94. In all cases, no volatiles are needed to explain the mass and radius of K2-291 b.

Zoom In Zoom Out Reset image size

We also estimated the maximum envelope mass fraction of K2-291 b through a model grid from Lopez & Fortney (2014). This grid assumes a solar metallicity envelope with a minimum envelope mass fraction of 0.1%. We generated 100,000 random samples of the envelope fraction from our normal distributions on the mass, radius, age, and flux of K2-291 b. From this, we determined that the 3σ upper limit on the envelope fraction is 0.3%.

Similarly, Kepler planets within 0.15 au and smaller than 2 R_⊕ have an envelope fraction less than 1% (Wolfgang & Lopez 2015). Figure 8 shows the relationship between density and stellar insolation for planets smaller than 4 R_⊕. K2-291 b exhibits a density similar to other small, close-in planets.

Zoom In Zoom Out Reset image size

K2-291 b's lack of a substantial volatile envelope could be explained by atmospheric loss. For lower mass planets experiencing a large amount of stellar insolation, photoevaporation (hydrodynamic escape) is the dominant atmospheric loss process. Photoevaporation occurs when high-energy photons from the host star ionize and heat the atmosphere causing it to expand and escape (Owen & Wu 2013).

K2-291 b is potentially the core of a sub-Neptune planet that underwent photoevaporation. We cannot, however, rule out a scenario where K2-291 b formed with a high density from its onset. In that case, perhaps K2-291 b formed after the gas disk had dissipated, or giant impacts by planetesimals stripped the envelope early in its formation. Although these two scenarios cannot yet be distinguished for an individual planet, population studies can be of use. Swain et al. (2018) finds two separate groups of small planets in radius-insolation-density space. One group is consistent with small solar system bodies and likely has an Earth-like formation, the other forms a bulk density continuum with sub-Neptunes and is likely composed of remnant cores produced by photoevaporation. Another large-scale approach is to look for a radius trend among close orbiting planets of different ages; a trend of smaller young planets compared to larger old planets would suggest photoevaporation. David et al. (2018) finds one such planet and Mann et al. (2017) finds seven close orbiting young planets. There is an emerging trend that these young planets are larger but more planets will need to be found to be statistically significant.

We examine here the possibility that K2-291 b formed by photevaporation. Due to the hydrodynamic escape of the envelope for close-in planets, the boundary between complete loss and retention of 1% of the envelope is at 0.1 au for a 6 M_⊕ planet orbiting a solar mass star (Owen & Wu 2013). K2-291 b orbits within this boundary at a = 0.03261 ± 0.00044 au. For the mass (M_p = 6.49 ± 1.16 M_⊕) and stellar insolation (S_inc = ${633}_{-56}^{+59}$ S_⊕) of K2-291 b specifically, all of its hydrogen–helium should have been lost between 100 Myr and 1 Gyr, depending on the original hydrogen–helium mass fraction and mass loss efficiency (Lopez & Fortney 2013). We determined an age from the HIRES spectra of ${3.7}_{-2.5}^{+3.7}$ Gyr, longer than this photevaporation timescale.

We ran additional models using the Lopez & Fortney (2014) model grid to calculate the radius K2-291 b would have with an additional hydrogen–helium envelope. Adding 0.1% H/He by mass would result in a planet radius of R_p = 1.82 R_⊕. Similarly, an additional 1% or 10% would equal a radius of R_p = 2.2 R_⊕ or R_p = 3.7 R_⊕, respectively. Therefore, a small addition of between 1% and 10% H/He would increase the radius of K2-291 b enough to move the planet across the Fulton gap to the sub-Neptune side.

Together, these analyses imply that K2-291 b may have formed as a sub-Neptune with a substantial volatile envelope and transitioned across the Fulton gap to a super-Earth planet through photevaporation.