A Sub-Neptune-sized Planet Transiting the M2.5 Dwarf G 9-40: Validation with the Habitable-zone Planet Finder

G 9-40 was observed by the Kepler spacecraft as part of Campaign 16 of the K2 mission. It was proposed as a K2 Campaign 16 target by the following programs: GO16005_LC (PI: Crossfield), GO16009_LC (PI: Charbonneau), GO16052_LC (PI: Stello), and GO16083_LC (PI: Coughlin). The star was monitored in long cadence mode (30 minute cadence) for 80 days from 2017 December 7 to 2018 February 25 and was originally identified as a candidate planet host by Yu et al. (2018). We used the Everest pipeline (Luger et al. 2016, 2018) to detrend and correct for photometric variations seen in the K2 photometry that are due to imperfect pointing of the spacecraft with only two functioning reaction wheels. The unbinned 6.5 hr photometric precision before detrending was 193 ppm, and it decreased to 30 ppm after detrending with Everest.

We obtained four visits of G 9-40 with the HPF Spectrograph with the goal to measure its RV variation as a function of time. HPF is a high-resolution (R ∼ 55,000) NIR spectrograph recently commissioned on the 10 m Hobby–Eberly Telescope (HET) in Texas covering the information-rich z, Y, and J bands from 810 to 1280 $\,\mathrm{nm}$ (Mahadevan et al. 2012, 2014). HPF is actively temperature stabilized, achieving ∼1 $\,\mathrm{mK}$ temperature stability long term (Stefansson et al. 2016). The HET is a fully queue-scheduled telescope with all observations executed in a queue by the HET resident astronomers (Shetrone et al. 2007). In each visit, we obtained two 945 s exposures, yielding a total of eight spectra with a median signal-to-noise ratio (S/N) of 123 per extracted 1D pixel at λ = 1000 $\,\mathrm{nm}$.

HPF has an NIR laser-frequency comb (LFC) calibrator that has been shown to enable $\sim 20\,\mathrm{cm}\ {{\rm{s}}}^{-1}$ calibration precision and $1.53\,{\rm{m}}\ {{\rm{s}}}^{-1}$ RV precision on sky on the bright and stable M dwarf Barnard's Star (Metcalf et al. 2019). We elected to not use the simultaneous LFC reference calibrator for these observations to minimize any possibility of scattered LFC light in the target spectrum. Instead, the drift correction was performed by extrapolating the wavelength solution from other LFC exposures from the night of the observations. As discussed in detail in the Appendix, HPF's nightly drift is a predictable linear sawtooth pattern with a 10–15 $\,{\rm{m}}\,{\rm{s}}$⁻¹ amplitude (see Figure 14 in Appendix), where dedicated LFC calibration exposures are generally taken nightly a few hours apart, enabling precise wavelength calibration even without simultaneous LFC calibration exposures (Metcalf et al. 2019). To further illustrate this linear behavior, in the Appendix, we show the LFC drift exposures for the four nights we obtained G 9-40 observations. Lastly, in the Appendix, we also estimate the error contribution of the drifts to the total RV errors using a leave-one approach. In doing so, we estimate errors at the <30 cm s⁻¹ level for our G 9-40 observations, which we add in quadrature to the errors estimated from our RV extraction (see Table 6).

The HPF 1D spectra were reduced and extracted with the custom HPF data-extraction pipeline following the procedures outlined in Ninan et al. (2018), Kaplan et al. (2018), and Metcalf et al. (2019). After the 1D spectral extraction, we extracted precise NIR RVs using a modified version of the SpEctrum Radial Velocity Analyzer (SERVAL) pipeline (Zechmeister et al. 2018) following the methodology discussed in Metcalf et al. (2019). In short, SERVAL uses a template-matching technique to derive RVs (see Anglada-Escudé & Butler 2012; Zechmeister et al. 2018). This technique relies on minimizing the differences of the observed spectrum against a high-S/N master template constructed from a coaddition of all available observations of the target star. We generated the master template by performing a best-fit spline regression to the eight as-observed spectra of the target star. Following the methodology described in Metcalf et al. (2019), we masked out all telluric and sky-emission line regions, using a thresholded synthetic telluric-line mask and an empirical thresholded sky-emission line mask, respectively. We generated the telluric-line mask using the telfit (Gullikson et al. 2014) Python wrapper to the Line-by-Line Radiative Transfer Model package (LBLRTM; Clough et al. 2005). For the sky-emission line mask, we used the same empirical blank-sky sky-emission line mask as used in Metcalf et al. (2019). We calculated the barycentric corrections using the barycorrpy Python package (Kanodia & Wright 2018), which uses the barycentric-correction algorithms from Wright & Eastman (2014).

To constrain contamination from background sources and study the possibility that astrophysical false positives could be the source of the transits (e.g., eclipsing binaries or background eclipsing binaries), we performed high-contrast adaptive optics (AO) imaging of G 9-40 using the 3 m Shane Telescope at Lick Observatory on 2018 November 25. We performed the AO imaging using the upgraded ShARCS camera (Srinath et al. 2014) in the K_s bandpass. We observed the star using a five-point dither pattern (see, e.g., Furlan et al. 2017), imaging the star at the center of the detector and in each quadrant. Individual images had an exposure time of 30 s, and we obtained a total of 20 minutes of integration time across all exposures.

Standard image processing (e.g., flat-fielding, sky subtraction), as well as subpixel image alignment, was performed using custom Python software. Using the combined image, we computed the variance in flux in a series of concentric annuli centered on the target star. The resulting 5σ contrast curve is shown in Figure 1. From the images we see that no bright (ΔK_s < 4) secondary companions are detectable within 05. These data show that there is no evidence of blending up to a radius of 25. We therefore infer that the transits observed from both K2 and from the ground (see Section 4) are not contaminated by background sources.

Zoom In Zoom Out Reset image size

To constrain blends within the K2 aperture at distances farther than the 25 field of view (FOV) of the diffraction-limited Lick/ShaneAO images shown in Figure 1, we obtained seeing-limited images using the ARCTIC imager (Huehnerhoff et al. 2016) on the 3.5 m Astrophysical Research Consortium Telescope (Figure 2(b)). On the night of 2019 January 14, we obtained 20 seeing-limited images of G 9-40 using an exposure time of 10 $\,{\rm{s}}$ in the Sloan Digital Sky Survey (SDSS) g' filter. The skies were photometric with a seeing of 075. We median combined the full set of 20 images to create a high-S/N final image shown in Figure 2. Given G 9-40's high proper motion, we compare our median-combined seeing-limited image to archival observations from the POSS1 survey obtained using the astroquery skyview service (Ginsburg et al. 2018) in Figure 2. The POSS1 image shows no evidence for a background star at G 9-40's current coordinates marked by the yellow circle (7'' radius) in both panels in Figure 2. In addition to this, we also queried the Gaia archive (Gaia Collaboration 2018) for nearby sources. In doing so, Gaia reveals one companion at a distance of 20'' from G 9-40, with a ΔG = 5, that is both significantly fainter than G 9-40 and outside the K2 aperture. Given G 9-40's proper motion of ${\mu }_{\alpha }=175\,\mathrm{mas}\ {\mathrm{yr}}^{-1}$ and ${\mu }_{\delta }=-318\,\mathrm{mas}\ {\mathrm{yr}}^{-1}$ and the short time span between the Gaia and K2 observations, we conclude that this star is not a significant source of dilution in the K2 and diffuser-assisted transits.

Zoom In Zoom Out Reset image size

We observed the transit of G 9-40b using the ARCTIC imager (Huehnerhoff et al. 2016) on the ARC 3.5 m Telescope at Apache Point Observatory on the night of 2019 April 13 local time (04:50 UT–07:30 UT April 14). The night was photometric with a seeing of ∼075 at the beginning of the night. The target rose from air mass 1.55 to air mass 1.05 during the observations.

We observed the transit using the Engineered Diffuser on the ARCTIC imager, which has been described in detail in Stefansson et al. (2017, 2018a). In short, the Engineered Diffuser molds the focal-plane image of the star into a broad and stable top-hat shape, allowing us to increase our exposure times to gather more photons per exposure while minimizing correlated errors due to point-spread function variations and guiding errors. The increased exposure time allows us to further average over scintillation errors by increasing the duty cycle of the observations (Stefansson et al. 2017). We used a newly commissioned, custom-fabricated narrowband filter from AVR Optics²⁵ and Semrock Optics,²⁶ centered on 857 nm with a 37 nm width in a region with minimal telluric absorption. The filter, along with its as-measured throughput curve, has been further discussed in Stefansson et al. (2018b). The filter is 5 inches in diameter, covering the full ARCTIC beam footprint in the optical path, resulting in minimal vignetting and allowing the observer to make use of the full field of view (FOV) of ARCTIC. The combination of the Engineered Diffuser and the narrowband filter allowed us to make use of a bright reference star HD 76780 (V = 7.63, I = 6.7; G 9-40 has an SDSS i' magnitude of 11.994) within the FOV at a high observing efficiency of 76.9% with a cadence of 18.2 s between successive exposures (exposure time of 14 s).²⁷

We reduced the on-sky transit observations along with the 25 bias and 25 dome flat-field calibration frames using AstroImageJ (Collins et al. 2017) following the methodology in Stefansson et al. (2017). We calculate the Barycentric Julian Date (BJD_TDB) time using the Python package barycorrpy. After experimenting with different aperture radii and annuli, we found that the aperture setting that yielded the smallest residuals was a value of 30 $\,\mathrm{pixels}$ (69) for the target star aperture, 40 $\,\mathrm{pixels}$ (92) for the inner annulus, and 80 $\,\mathrm{pixels}$ (184) for the outer annulus for background estimation. This resulted in peak ADU/pixel counts of ∼1000 and ∼42,900 for the target and reference star, respectively. After removing the best-fit transit model and six additional significant (>3σ) outliers present in the data, we obtained an unbinned precision of 1225 ppm, which is further discussed in Section 4. The transit is shown in Figure 7, where the data are shown with error bars estimated following the methodology in Stefansson et al. (2017) and accounting for errors from photon, dark, read, background, digitization, and scintillation noise.

In this subsection, we discuss our implementation of the SpecMatch-Emp algorithm described in Yee et al. (2017) to derive precise estimates of the spectroscopic effective temperature (T_eff), metallicity ([Fe/H]), and surface gravity (log g) of G 9-40 through comparing our high-resolution NIR HPF spectra of G 9-40 to a library of high-S/N as-observed HPF spectra with well-characterized stellar parameters. As this is the first description of this implementation for this algorithm for HPF spectra, we provide an initial description here on its performance. As we assemble a larger HPF spectral library, we plan to further detail the performance and applicability of this algorithm on the larger high-resolution NIR M dwarf library. In the next subsections we further discuss our spectral library, the algorithm, and its performance in recovering the stellar parameters of the stars internal to the library using a cross-calibration scheme. We then discuss our best-estimate stellar parameter values of G 9-40, and finally we discuss future work we plan to perform to further improve the performance of the code framework described here.

3.1.1. Spectral Library

3.1.2. Empirical SpecMatch Algorithm

Following Yee et al. (2017), we use a χ² metric to compare the goodness of fit of a target spectrum S_target to a given S_reference star spectrum in the spectral library,

$$\begin{eqnarray}&&{\chi }_{\mathrm{initial}}^{2}=\sum _{i=1}^{N}{\left({S}_{\mathrm{target},i}-{S}_{\mathrm{reference},i}(p,v\sin i\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where S_target,i is the deblazed and normalized flux of the target star, and ${S}_{\mathrm{reference},i}(p,v\sin i)$ is the deblazed and normalized flux of the library reference star at the ith element of the resampled wavelength array (N in total). We rotationally broadened the reference star flux, S_reference,i(p, v sin i), by a projected rotational velocity v sin i and multiplied by a fifth-order Chebyshev polynomial denoted by the six-element vector p to remove any low-order residual variations in the deblazed spectra. For the v sin(i) rotational broadening, we adopt the broadening kernel from Gray (1992), following the implementation in Yee et al. (2017). Following Yee et al. (2017), we specifically do not scale the χ² value in Equation (1) by the estimated photon-noise error bars, as the residuals of the high-S/N spectra are completely dominated by systematic astrophysical or instrumental differences rather than our estimate of the photon noise. Using the χ² metric from Equation (1), we loop through all of the spectra in the library to optimize for a minimum χ² value as a function of the Chebyshev polynomial vector p and v sin i. We used the Nelder-Mead simplex ("Amoeba") algorithm (Nelder & Mead 1965) as implemented in the scipy.optimize package. Following Yee et al. (2017), we enforce a limited range of allowed v sin i values between $0\,\mathrm{km}\ {{\rm{s}}}^{-1}$ (no broadening) and 15 $\,\mathrm{km}\,{\rm{s}}$⁻¹ to minimize excursions to artificially high v sin i values.

Following the ${\chi }_{\mathrm{initial}}^{2}$ loop, we then select the top five lowest-${\chi }_{\mathrm{initial}}^{2}$ spectra to form a new composite spectrum, ${S}_{\mathrm{composite}}\,={\sum }_{j=1}^{5}{c}_{j}{S}_{j}({p}_{j},{[v\sin i]}_{j})$, where the five coefficients $c=\{{c}_{1},{c}_{2},{c}_{3},{c}_{4},{c}_{5}\}$ are optimized to further minimize the χ² of the target star and the composite spectrum according to the following χ² metric:

$$\begin{eqnarray}&&{\chi }_{\mathrm{composite}}^{2}=\sum _{i=1}^{N}{\left({S}_{\mathrm{target},i}-\sum _{j=1}^{5}{c}_{j}{S}_{j,i}({p}_{j},{[v\sin i]}_{j}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where p_j and [v sin i]_j are the best-fit values determined from the optimization step in Equation (1) for the jth best reference spectrum ($j\in \{1,2,\,\ldots 5\}$). We fit for the first four c coefficients in the optimization step and then set

$$\begin{eqnarray}&&{c}_{5}=1-\sum _{j=1}^{4}{c}_{j}\end{eqnarray} \tag{ 3 }$$

to ensure that all five coefficients sum up to unity. All of the c coefficients are further constrained to have a value between 0 and 1.

3.1.3. Cross-validation

To check the performance of the algorithm described above on the NIR HPF spectra, we performed a cross-validation procedure consisting of removing a given spectrum from the library and comparing the recovered best-fit stellar parameter to its known library value. We then repeated this comparison for all of the stars in the library, and we computed the standard deviation (σ) of the residuals between the recovered best-fit stellar parameters and the known library value for the three stellar parameters considered (i.e., computing ${\sigma }_{{T}_{\mathrm{eff}}},{\sigma }_{\mathrm{Fe}/{\rm{H}}}$, and σ_logg). To compare the performance of different HPF orders in recovering the known parameter values, we ran this cross-validation procedure on six different HPF orders that are relatively clean of telluric absorption features. Table 1 summarizes the comparison between the different orders, showing the resulting ${\sigma }_{{T}_{\mathrm{eff}}}$, ${\sigma }_{\mathrm{Fe}/{\rm{H}}}$, and σ_logg values. From Table 1, we see that the different orders overall show similar scatter, with order 5 (wavelengths: [8670–8750 Å]) having the lowest ${\sigma }_{{T}_{\mathrm{eff}}}$ value and order 17 (wavelengths: [10460–10570 Å]) having the lowest ${\sigma }_{\mathrm{Fe}/{\rm{H}}}$ value. In addition, Table 1 shows the individual best-fit stellar parameter point estimates for G 9-40 for the same orders. The derivation of the point estimates is further discussed in Section 3.1.4.

Table 1.  Results of Stellar Parameter Estimation for G 9-40 for the HPF Orders Considered (see Section 3.1.4), with the Library Cross-validation Residual Standard Deviations (see Section 3.1.3) for the Same Orders

0	1	2	G 9-40 Best-fit Values			Library Cross-validation
Order	Wavelength Region	S/N	Teff	Fe/H	log g	${\sigma }_{{T}_{\mathrm{eff}}}$	${\sigma }_{\mathrm{Fe}/{\rm{H}}}$	σlogg
	(Å)		(K)	(dex)	(dex)	(K)	(dex)	(dex)
5	[8670, 8750]	85	3426	−0.084	4.896	64	0.15	0.05
6	[8790, 8885]	90	3420	−0.031	4.860	83	0.18	0.06
14	[9940, 10055]	126	3366	−0.048	4.878	100	0.17	0.06
15	[10105, 10220]	131	3378	−0.063	4.873	86	0.17	0.05
16	[10280, 10395]	134	3381	−0.061	4.898	69	0.15	0.05
17	[10460, 10570]	138	3405	−0.078	4.909	73	0.13	0.05

Notes. The cross-validation standard deviations are the standard deviation of the residuals between the recovered and known library values of all of the stars in the stellar library. We assign the cross-validation standard deviations as our error on the corresponding stellar-parameter estimate. For the spectroscopic stellar-parameter point estimates, we see that all of the orders yield consistent parameters within the error bars determined from the library cross-validation. We adapt the boldface values as the stellar parameters of G 9-40. Figure 3 graphically depicts the individual library residuals for order 17, our highest S/N order, and we assign the stellar parameters derived from order 17 as our final parameters for G 9-40, as is further discussed in Section 3.1.3.

Download table as:  ASCIITypeset image

Although all of the orders considered perform similarly in the cross-validation, we elect to use order 17 for our stellar parameter point estimate, as that order has the highest S/N and this order has some of the fewest telluric absorption lines of the orders considered, minimizing systematic errors from telluric absorption. Further, this order empirically shows a slightly better performance in recovering the known metallicity. Although we do not observe G 9-40 to exhibit significant chromospheric Ca infrared (IRT) activity from inspection of the HPF spectra, we elect not to use order 5 for our final spectral parameter point estimate (although it formally shows the lowest ${\sigma }_{{T}_{\mathrm{eff}}}$ value in our cross-validation), as that order contains one of the Ca IRT lines that can be in emission for active stars and could adversely affect the χ² comparison. As such, using the cross-validation results for order 17, we assign ${\sigma }_{{T}_{\mathrm{eff}}}=73\,{\rm{K}}$, σ_Fe/H = 0.13 $\,\mathrm{dex}$, and ${\sigma }_{\mathrm{log}g}=0.05\,\mathrm{dex}$ as our best estimates of the 1σ errors for our point estimates of T_eff, [Fe/H], and log g, respectively. Figure 3 graphically shows the cross-validation residuals for order 17.

Zoom In Zoom Out Reset image size

3.1.4. G 9-40 Parameter Estimation

We ran the spectral-matching algorithm described above on the highest-S/N spectrum of G 9-40, which had an S/N of 148 at ∼1.1 μm and an S/N of 137 in order 17. The top panels in Figure 4 visualize the ${\chi }_{\mathrm{initial}}^{2}$ values as calculated using Equation (1) for all of the library stars in the T_eff and [Fe/H] plane (left) and the T_eff and log g plane (right). The five best-matching stars (GJ 251, GJ 581, GJ 109, GJ 1148, and GJ 105 B) are highlighted in red in the upper panels in Figure 4. We then used these five best-matching stars for the second linear combination step optimizing the ${\chi }_{\mathrm{composite}}^{2}$ value in Equation (2). The lower panel in Figure 4 compares the spectra of these five stars, along with the best-fit linearly combined composite spectrum. The final composite spectrum and the corresponding residuals from the target star are shown at the bottom. Using the weights, we determine the following stellar parameters for G 9-40 using order 17: ${T}_{\mathrm{eff}}=3405\pm 73\,{\rm{K}}$, Fe/H = −0.078 ± 0.13, and log (g) = 4.909 ± 0.05, where our best estimate of the error bars is derived from our cross-validation step for order 17 in Table 1. We note that running the algorithm independently on the other orders in Table 1 results in consistent stellar parameter estimates within the 1σ uncertainties.

Zoom In Zoom Out Reset image size

3.1.5. Future Work

To estimate the mass, radius, and age of G 9-40, we use EXOFASTv2 to perform a spectral energy distribution (SED) and isochrone fit to the available literature photometry (see Table 2), using the Gaia parallax and the spectroscopically determined stellar parameters discussed in the previous subsection as Gaussian priors. EXOFASTv2 uses the BT-NextGen Model grid of theoretical spectra (Allard et al. 2012) and the Mesa Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016) to fit the SED and derive model-dependent stellar parameters. To test the consistency of our spectroscopic log g measurement from our spectral-matching algorithm, we use the same log g = 4.909 value as a starting point for our SED analysis but do not impose a strict prior, to minimize biases on the inferred evolutionary state, mass, and radius of the star (see, e.g., Torres et al. 2012). Figure 5 shows the resulting SED fit. The derived model-dependent stellar parameters show good agreement with the spectroscopic values derived from the HPF spectra from the previous subsection.

Zoom In Zoom Out Reset image size

Table 2.  Summary of Stellar Parameters

Parameter	Description	Value	References
Main Identifiers:
EPIC	⋯	212048748	Huber
LSPM	⋯	J0858+2104	Lepine
2MASS	⋯	J08585232+2104344	Huber
Gaia DR2	⋯	684992690384102528	Gaia
Equatorial Coordinates, Proper Motion, and Spectral Type:
αJ2000	Right ascension (R.A.)	08:58:52.32	Gaia
δJ2000	Declination (decl.)	+21:04:34.20	Gaia
μα	Proper motion (R.A., $\,\mathrm{mas}\ {\mathrm{yr}}^{-1}$)	175.512 ± 0.103	Gaia
μδ	Proper motion (decl., $\,\mathrm{mas}\ {\mathrm{yr}}^{-1}$)	-318.469 ± 0.067	Gaia
Spectral type	⋯	M2.5	Reid
Optical and Near-infrared Magnitudes:
B	APASS Johnson B mag	15.462 ± 0.074	APASS
V	APASS Johnson V mag	13.823 ± 0.040	APASS
g'	APASS Sloan g' mag	14.626 ± 0.047	APASS
r'	APASS Sloan r' mag	13.217 ± 0.061	APASS
i'	APASS Sloan i' mag	11.994 ± 0.061	APASS
Kepler-mag	Kepler magnitude	12.771	Huber
J	2MASS J mag	10.058 ± 0.022	2MASS
H	2MASS H mag	9.433 ± 0.023	2MASS
KS	2MASS KS mag	9.190 ± 0.018	2MASS
WISE1	WISE1 mag	9.032 ± 0.023	WISE
WISE2	WISE2 mag	8.885 ± 0.020	WISE
WISE3	WISE3 mag	8.774 ± 0.029	WISE
WISE4	WISE4 mag	8.597 ± 0.411	WISE
Spectroscopic Parametersa:
${T}_{\mathrm{eff}}$	Effective temperature in $\,{\rm{K}}$	3405 ± 73	This work
[Fe/H]	Metallicity in dex	−0.078 ± 0.13	This work
log(g)	Surface gravity in cgs units	4.909 ± 0.05	This work
Model-dependent Stellar SED and Isochrone Fit Parametersb:
Teff	Effective temperature in $\,{\rm{K}}$	3348 ± 32	This work
[Fe/H]	Metallicity in dex	${0.04}_{-0.11}^{+0.10}$	This work
log(g)	Surface gravity in cgs units	4.926 ± 0.027	This work
M*	Mass in M⊙	0.290 ± 0.020	This work
R*	Radius in R⊙	${0.3073}_{-0.0061}^{+0.0059}$	This work
ρ*	Density in $\,{\rm{g}}\ {\mathrm{cm}}^{-3}$	${14.11}_{-0.92}^{+0.99}$	This work
Age	Age in Gyr	${9.9}_{-4.1}^{+2.6}$	This work
L*	Luminosity in L⊙	${0.01069}_{-0.00029}^{+0.00028}$	This work
Av	Visual extinction in mag	${0.077}_{-0.046}^{+0.034}$	This work
d	Distance in pc	27.928 ± 0.045	Gaia
π	Parallax in mas	35.807 ± 0.059	Gaia
Other Stellar Parameters:
Prot	Rotational period in days	${29.85}_{-0.94}^{+1.01}$	This work
vsin i*	Stellar rotational velocity in km s−1	<2	This work
RV	Absolute radial velocity $\,\mathrm{km}\ {\rm{s}}$−1 (γ)	14.65 ± 0.26	This work

Notes.

^aDerived using our empirical spectral-matching algorithm. ^bEXOFASTv2-derived values using MIST isochrones with the Gaia parallax and spectroscopic parameters in a) as priors.

References: Huber (Huber et al. 2016), Lepine (Lépine & Shara 2005), Reid (Reid et al. 2004), Gaia (Gaia Collaboration 2018), APASS (Henden et al. 2015), UCAC2 (Zacharias et al. 2004), 2MASS (Cutri et al. 2003), WISE (Cutri et al. 2014).

Download table as:  ASCIITypeset image

Dressing et al. (2019) provide spectroscopic and photometric estimates of the stellar parameters of G 9-40 that are independent of the parameters presented here. Their spectroscopic estimates are T_eff = 3264₋₁₂₁⁺¹³⁷$\,{\rm{K}}$, ${R}_{* }={0.328}_{-0.054}^{+0.064}\,{R}_{\odot }$, ${M}_{* }={0.169}_{-0.153}^{+0.142}\,{M}_{\odot }$, and $[\mathrm{Fe}/{\rm{H}}]=-0.202\pm 0.084$. Their photometric parameters are T_eff = 3310 ± 57 $\,{\rm{K}}$, ${R}_{* }\,=0.313\pm 0.009\,{R}_{\odot }$, and ${M}_{* }=0.295\pm 0.007\,{M}_{\odot }$, where they adapt the photometric parameters as the best estimate of the stellar parameters. Their adapted parameters are in good agreement with both our spectroscopic stellar parameters and our EXOFASTv2 model-dependent parameters as summarized in Table 2.

We measure the projected rotation velocity v sin i of the star using the high-resolution spectra from HPF, following the method in Reiners et al. (2012). This method compares the CCF of the target spectrum to the resulting CCFs of artificially broadened spectra of a slowly rotating template star. To calculate the CCF, we use the CCF routines from the publicly available Collection of Elemental Routines for Echelle Spectra code base (Brahm et al. 2017), which calculate the CCF by cross-correlating a given spectrum to a fixed binary mask. We elected to use Barnard's Star as the template star as it is known to be a slow rotator (see, e.g., Ribas et al. 2018) and has a similar spectral type of M4 to G 9-40 of M2.5. To generate the binary mask for the CCF, we used a peak-finding algorithm to find single unblended lines in the Barnard's Star spectra. Blended lines have been shown to introduce substantial systematic errors in the determination of the v sin i for relatively small differences in the spectral characteristics of input stars (Reiners et al. 2018). To minimize the impact of tellurics, we only run the peak-finding algorithm on the eight HPF orders cleanest of tellurics (the same eight orders as in the RV analysis in this work). Although relatively clean of tellurics, we further filtered the resulting list to produce a list of lines with minimal telluric overlap using the same telluric mask used for the RV analysis in this work (see Section 2.2). To estimate the scatter in our determination of the v sin i value, we independently estimated the v sin i for each of the eight orders separately, resulting in eight independent point estimates of the v sin i. In doing so, all independent orders resulted in a v sin i at the lower limit of our v sin i measurement precision (i.e., the width of the target star CCF was fully consistent with the width of the slow-rotator calibrator CCF), suggesting that the v sin i is below the resolution of HPF. At the resolution of HPF (R = 55,000, FWHM = 6 $\,\mathrm{km}\ {{\rm{s}}}^{-1}$), we estimate that we can detect rotation velocities of about $2\,\mathrm{km}\ {{\rm{s}}}^{-1}$ or more. We thus conclude that G 9-40 has a $v\sin i\lt 2\,\mathrm{km}\ {{\rm{s}}}^{-1}$ from the resolution of HPF. This agrees with the slow rotational velocity observed in the K2 data, as is further discussed in the following subsection.

To measure the absolute RV of G 9-40, we fit a Gaussian profile to the CCF discussed above. Analogous to the v sin i determination, we fitted each order independently, yielding eight point estimates of the absolute RV. All eight point estimates yielded a consistent value, and we used the standard deviation of the eight different point estimates as our estimate of the error. This resulted in an absolute RV of ${RV}=14.65\,\pm 0.26\,\mathrm{km}\ {{\rm{s}}}^{-1}$, which is also summarized in Table 2.

We use the K2 photometry to estimate the rotation period of G 9-40. Looking at the photometry from K2 (Figure 7(a)), we see that there are small-scale (<1%) modulations with a period between 25 and 35 days seen on top of a long-term trend after a sharp flux decrease after the first three days. We argue that both the initial sharp flux decrease and the long-term trend are likely not astrophysical (e.g., due to potential systematics including thermal settling of the spacecraft). Therefore, to estimate the rotation period, we remove the first three days and then remove a simple linear trend from the resulting photometry, yielding the light curve shown in Figure 7(b). We then fit the rotation period using two methods as described below.

First, we use the first peak of the autocorrelation function (ACF) as an estimate of the rotation period (see, e.g., McQuillan et al. 2013). To generate the ACF, we used the acf function in the statsmodels Python package. Before calculating the ACF, we masked the data points during transit, removed any >3σ outliers, and then interpolated the masked values to give a uniformly sampled light curve at the ∼30 minute cadence of K2. Figure 6 shows the resulting ACF for the linear-trend-removed photometry in Figure 7(b) along with the location of the highest peak, yielding a rotation period of 27.5 $\,\mathrm{days}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Second, we derived an additional estimate for the rotation period by using a Gaussian process. Angus et al. (2018) show that a quasi-periodic covariance function is effective for making probabilistic measurements of the rotation period even if the data is sparsely sampled. Exact Gaussian process modeling has the disadvantage of the run time scaling as ${ \mathcal O }\left({N}^{3}\right)$. For this purpose, we use the juliet analysis package (Espinoza et al. 2019), which models the photometry using the formalism presented in Foreman-Mackey et al. (2017), where a simple function is constructed that mimics the properties of the quasi-periodic covariance function and implemented with the celerite Python package. The period is the only hyperparameter, and there are three nuisance parameters. Using the same data that generated the ACF, we calculate a rotation period of ${29.85}_{-0.94}^{+1.01}$ days, which is in broad agreement with our rotational period estimate from the ACF method (within 3σ). This results in a maximum stellar equatorial rotational velocity of ${v}_{\mathrm{eq}}\sim 500\,{\rm{m}}\ {{\rm{s}}}^{-1}$, in good agreement with the v sin i constraint in the previous subsection. From the K2 photometry, the amplitude of the modulation is ∼0.5%, which could be measured through extensive long-baseline ground-based photometric observations and used to confirm the P_rot value presented here. Although they do not present a rotation period estimate, Pepper et al. (2008) classify G 9-40 as a long-period variable using photometric data from the Kilodegree ELT, which is consistent with the rotation period estimate presented here.

Although Yu et al. (2018) provide an ephemeris for the transits in the K2 data for EPIC 212048748, we describe our independent transit search here. To search for transits, we flattened the K2 light curve of G 9-40 using the best-fit GP model from Everest. We then ran a Fortran and Python implementation²⁸ of the box least squares (BLS) algorithm (Kovács et al. 2002) on the flattened light curve. After finding a significant periodic signal using the BLS algorithm, we masked 2 × T₁₄ transit-time windows surrounding the calculated transit midpoints and then reran the BLS algorithm to look for potential further transits. We found no strong evidence for other transits in the system. The final K2 photometry from Everest is shown in Figure 7 along with the identified transits.

After identifying the transit locations using the BLS algorithm, we ran three different fits: a fit of the K2 transits only, a fit of the ground-based transit only, and finally a joint fit with the K2 and ground-based transits simultaneously to put further constraints on the orbital period. We performed all transit fits in a Markov Chain Monte Carlo (MCMC) framework to obtain parameter posteriors following the methodology described in Stefansson et al. (2017, 2018a). In short, for all three fits, we used the batman Python package for the transit model, which uses the transit model formalism from Mandel & Agol (2002). Before starting the MCMC runs, we found the global best-fit solution using a differential-evolution global optimization package called PyDE,²⁹ maximizing the log-posterior probability (sum of the log-likelihood function and log prior probabilities). We initialized 100 walkers distributed in a uniform N-dimensional sphere (where N is the number of jump parameters used) close to the maximum-likelihood solution using the emcee affine-invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013). We computed 50,000 steps for each Markov Chain walker, resulting in a total of 5,000,000 computed samples. We removed the first 5000 steps in each chain as burn-in, resulting in 4,500,000 samples used for final parameter inference. The Gelman–Rubin statistics for the resulting chains were all <3% from unity, which we considered well mixed (see, e.g., Ford 2006).

We used the following five MCMC jump parameters describing the planet transit: transit center (T_C), period (log(P)), inclination (cos(i)), radius ratio (R_p/R_*), and semimajor axis (log(a/R_*)), along with one parameter for the out-of-transit baseline flux. We explicitly fix the eccentricity and the argument of periastron to be equal to zero. We impose broad uniform priors on all five jump parameters that are summarized in Table 3. For T_C, log(P), and R_p/R_*, we center the priors on the values determined by the BLS algorithm, but as we do not get constraints on the values of log(a/R_*) and cos(i) from the BLS algorithm, we imposed wide uniform priors on these parameters. We impose a Gaussian prior centered on unity for the baseline flux parameter. We used a quadratic limb-darkening law to describe the limb-darkening using the triangular-sampling parameterization described in Kipping (2013). Following Kipping (2013), we fully marginalize over the complete parameter space of the q₁ and q₂ limb-darkening coefficients (from 0 to 1) to minimize biases on our planet parameter values due to inaccuracies in the stellar parameters.

Table 3.  Summary of Priors Used for the Three MCMC Transit Fits Performed

Parameter	Description	K2-only	Ground-only	Joint
log(P) (days)	Orbital period	${ \mathcal U }(0.7592,0.7594)$	${ \mathcal N }(0.75936,0.00002)$	${ \mathcal U }(0.7592,0.7594)$
TC	Transit midpoint $({\mathrm{BJD}}_{\mathrm{TDB}})$	${ \mathcal U }(2458095.54,2458095.56)$	${ \mathcal U }(2458497.76,2458497.80)$	${ \mathcal U }(2458497.76,2458497.80)$
${({R}_{p}/{R}_{* })}_{K2}$	Radius ratio	${ \mathcal U }(0,0.1)$	⋯	${ \mathcal U }(0,0.1)$
${({R}_{p}/{R}_{* })}_{\mathrm{ground}}$	Radius ratio	⋯	${ \mathcal U }(0,0.1)$	${ \mathcal U }(0,0.1)$
cos(i)	Transit inclination	${ \mathcal U }(0,0.2)$	${ \mathcal U }(0,0.2)$	${ \mathcal U }(0,0.2)$
log(a/R*)	Normalized orbital radius	${ \mathcal U }(0.9,2.0)$	${ \mathcal U }(0.9,2.0)$	${ \mathcal U }(0.9,2.0)$
frawK2	Transit baseline for K2 data	${ \mathcal U }(0.9,1.1)$	⋯	${ \mathcal N }(1.00002,0.0002)$
frawGround	Transit baseline for ground-based data	⋯	${ \mathcal U }(0.9,1.1)$	${ \mathcal N }(1.001,0.001)$
DLine	Ground detrend parameter: line	⋯	${ \mathcal U }(-0.1,0.1)$	${ \mathcal N }(0.00104,0.00001)$

Notes. ${ \mathcal N }(m,\sigma )$ denotes a normal prior with mean m and standard deviation σ; ${ \mathcal U }(a,b)$ denotes a uniform prior with a start value a and end value b. The eccentricity was assumed to be zero for all fits. For all fits, we uniformly sampled quadratic limb-darkening parameters q₁ and q₂ from 0 to 1 using the formalism from Kipping (2013). Priors on T_eff and R_* are adapted from Table 2, but no prior is placed on the stellar density ρ.

Download table as:  ASCIITypeset image

Table 3 summarizes the priors used for all of the different fits. We report the median values along with the 16th and 84th percentile error bars derived from the posteriors for all three fits in Table 4. Our joint-fit median values collectively provide a good description of the transit (i.e., no obvious bimodality is seen in the resulting posteriors). We further note that our resulting best-fit transit parameters are consistent within the 1σ error bars with the parameters obtained by Yu et al. (2018). We further compare our resulting transit ephemeris to the ephemeris obtained by Yu et al. (2018) in Section 6.

Table 4.  Median Values and 68% Confidence Intervals for the Transit Fit Parameters for Our K2-only, Ground-based-only, and Joint K2 and Ground-based MCMC Analyses

Parameter	Description	K2	Ground	Joint Fit (Adopted)
TC(BJDTDB)	Transit midpoint	${2458095.55737}_{-0.00030}^{+0.00030}$	${2458497.77751}_{-0.00036}^{+0.00036}$	${2458497.77747}_{-0.00033}^{+0.00032}$
P (days)	Orbital period	${5.745951}_{-0.00004}^{+0.00004}$	${5.74596}_{-0.00026}^{+0.00026}$	${5.746007}_{-0.000006}^{+0.000006}$
${({R}_{p}/{R}_{* })}_{K2}$	Radius ratio (K2)	${0.059}_{-0.0026}^{+0.0035}$	⋯	${0.0605}_{-0.0028}^{+0.0026}$
${({R}_{p}/{R}_{* })}_{\mathrm{ground}}$	Radius ratio (ground)	⋯	${0.0635}_{-0.0047}^{+0.0031}$	${0.0605}_{-0.0033}^{+0.0032}$
${R}_{p,K2}({R}_{\oplus })$	Planet radius (K2)	${1.981}_{-0.095}^{+0.12}$	⋯	${2.025}_{-0.097}^{+0.096}$
${R}_{p,\mathrm{ground}}({R}_{\oplus })$	Planet radius (ground)	⋯	${2.13}_{-0.16}^{+0.12}$	${2.03}_{-0.11}^{+0.11}$
${\delta }_{p,K2}$	Transit depth (K2)	${0.00348}_{-0.00030}^{+0.00043}$	⋯	${0.00365}_{-0.00033}^{+0.00032}$
${\delta }_{p,\mathrm{ground}}$	Transit depth (ground)	⋯	${0.00403}_{-0.00058}^{+0.00041}$	${0.00366}_{-0.00039}^{+0.00039}$
$a/{R}_{* }$	Normalized orbital radius	${29.8}_{-7.6}^{+5.8}$	${22.4}_{-2.9}^{+6.4}$	${27.0}_{-3.7}^{+5.4}$
a (au)	Semimajor axis (from $a/{R}_{* }$ and R*)	${0.0425}_{-0.011}^{+0.0082}$	${0.032}_{-0.0042}^{+0.0092}$	${0.0385}_{-0.0053}^{+0.0078}$
${\rho }_{* ,\mathrm{transit}}$ (g cm−3)	Density of star	${15.1}_{-8.9}^{+11.0}$	${6.4}_{-2.2}^{+7.2}$	${11.2}_{-4.0}^{+8.3}$
i(◦)	Transit inclination	${88.88}_{-0.97}^{+0.79}$	${87.98}_{-0.49}^{+0.85}$	${88.57}_{-0.47}^{+0.63}$
b	Impact parameter	${0.58}_{-0.38}^{+0.23}$	${0.788}_{-0.20}^{+0.064}$	${0.672}_{-0.22}^{+0.100}$
e	Eccentricity	0.0 (adopted)	0.0 (adopted)	0.0 (adopted)
ω (◦)	Argument of periastron	0.0 (adopted)	0.0 (adopted)	0.0 (adopted)
${T}_{\mathrm{eq}}$ (K)	Equilibrium temp. (assuming a = 0.3)	${304.0}_{-26.0}^{+48.0}$	${350.0}_{-42.0}^{+26.0}$	${319.0}_{-28.0}^{+25.0}$
Teq (K)	Equilibrium temp. (assuming a = 0.0)	${434.0}_{-37.0}^{+69.0}$	${500.0}_{-59.0}^{+37.0}$	${456.0}_{-40.0}^{+35.0}$
S (S⊕)	Insolation flux	${5.9}_{-1.8}^{+4.7}$	${10.4}_{-4.1}^{+3.4}$	${7.2}_{-2.2}^{+2.5}$
T14 (days)	Transit duration	${0.0546}_{-0.0018}^{+0.0026}$	${0.0583}_{-0.0025}^{+0.0026}$	${0.0557}_{-0.0017}^{+0.0019}$
τ (days)	Ingress/egress duration	${0.0044}_{-0.0015}^{+0.0044}$	${0.0086}_{-0.0040}^{+0.0035}$	${0.0056}_{-0.0019}^{+0.0023}$
${T}_{S}({\mathrm{BJD}}_{\mathrm{TDB}})$	Time of secondary eclipse	${2458098.43034}_{-0.00028}^{+0.00028}$	${2458500.65049}_{-0.00038}^{+0.00039}$	${2458500.65047}_{-0.00033}^{+0.00032}$

Notes. For our joint fit, we fit for the transit depth (δ) and the planet radius ratio (R_p/R_*) separately for K2 and the diffuser-assisted observations, resulting in slightly different values derived for the planet radii and transit depths in the two different bands. These values are denoted by (K2) and (ground), respectively. We fixed the eccentricity and argument of periastron to zero for all fits.

Download table as:  ASCIITypeset image

4.2.1. K2 Light Curve Analysis

4.2.2. Ground-based Light Curve Analysis

For the ground-based analysis, we estimate the photometric uncertainty, including scintillation, following the discussion surrounding Equation (3) in Stefansson et al. (2018a), yielding error bars that account for photon, dark, readout, background, digitization, and scintillation noise. For the scintillation error calculation, we assumed that the single reference star used was fully uncorrelated to the scintillation error from the target star. The jump parameters, along with the associated priors for this fit, are summarized in Table 3. In addition to the transit jump parameters, in this MCMC fit we also simultaneously detrend with a line to account for a low-amplitude linear slope observed in the photometry.

In Figure 8, we show the photometric precision (or the rms of the best-fit residuals) as a function of bin size using the MC3 package from Cubillos et al. (2017), which estimates error bars on the rms values assuming they follow an inverse gamma distribution. Additionally, shown in red in Figure 8 is the expected bin-down behavior assuming Gaussian white noise. From Figure 8, we see that in our ground-based observations, we achieve a photometric precision of ${\sigma }_{\mathrm{unbin}}=1225\,\mathrm{ppm}$ unbinned (18.2 s cadence), which bins down to ${\sigma }_{1\min }\,=689\pm 37\,\mathrm{ppm}$ and ${\sigma }_{30\min }={88}_{-27}^{+52}\,\mathrm{ppm}$, in 1 minute and 30 minute bins, respectively. For comparison, we note that the Gaussian expected precision in 30 minutes is 138 ppm. We observe that the rms values largely follow the expected white-noise behavior, although we note that slight excursions below the Gaussian expected values are observed at the largest bin sizes, which still are largely consistent with the Gaussian expected precision within the reported error bars. Similar and larger departures below the Gaussian expected precision have been reported by a number of groups in the literature (e.g., Blecic et al. 2013; Stefansson et al. 2017), and Cubillos et al. (2017) show that these excursions are not statistically significant after taking into account the increasingly skewed inverse gamma distribution of the rms values at the largest bin sizes. Following our previous work (Stefansson et al. 2017, 2018a), we argue that excursions much below the Gaussian expected precision are likely an overestimate of the actual precision achieved, and we conservatively say that we achieve 138 ppm precision in 30 minute bins for these transit observations.

Zoom In Zoom Out Reset image size

4.2.3. Joint K2 and Ground-based Light Curve Analysis

To enable better future scheduling of transit observations, we fit our ground-based diffuser transit and the K2 transits together to provide an updated ephemeris of G 9-40b. Figure 9 compares errors on the transit center as a function of time into the beginning of the JWST era. From our K2-only fit (blue curve), we see that the expected errors are ∼13 minutes at the start of the JWST era in 2021. However, from our updated joint-fit ephemeris (green curve), we reduce this error by a factor of 8 down to ∼1.7 minutes. We further see that our ground-based transit was observed at the upper edge of our 1σ error bar estimate from our K2-only fit.

Zoom In Zoom Out Reset image size

In Figure 9 we additionally compare our ephemerides to the ephemerides presented in Yu et al. (2018) calculated from K2 Campaign 16 data only (purple curve). We see that there is a discrepancy between the ephemerides derived in this work and the ephemerides presented in Yu et al. (2018). As mentioned above, our ground-based transit was observed at the upper edge of our 1σ error bar estimate from our K2-only fit, whereas our ground-based transit was observed ∼23 minutes later than expected from the Yu et al. (2018) ephemeris (no uncertainty estimate is reported in Yu et al. 2018 on the ephemeris). We attribute the discrepancy to the different photometric reduction and detrending methods, and we prefer the ephemeris reported here for the planning of future transit observations, given the consistency of the K2-only fits and our ground-based transit observations.

We used our joint transit model to get an independent estimate of the stellar density, using the equation from Seager & Mallén-Ornelas (2003),

$$\begin{eqnarray}&&{\rho }_{* ,\mathrm{transit}}=\displaystyle \frac{3\pi }{{{GP}}^{2}}{\left(\displaystyle \frac{a}{{R}_{* }}\right)}^{3},\end{eqnarray} \tag{ 4 }$$

where G is the gravitational constant, and the eccentricity is assumed to be zero. From Table 2, we see that our stellar density derived from our stellar radius and mass (${\rho }_{* }\,={14.11}_{-0.92}^{+0.99}\,{\rm{g}}\ {\mathrm{cm}}^{-3}$) is consistent with the stellar density derived from our transit observations (${\rho }_{* ,\mathrm{transit}}={11.2}_{-4.0}^{+8.3}\,{\rm{g}}\ {\mathrm{cm}}^{-3}$).

We predict the most-likely mass of G 9-40b using two different methods. First, we use the Forecaster Python package (Chen & Kipping 2017), which adopts a broken-power-law model to model the exoplanet MR relation and is capable of predicting the masses of exoplanets given their radii. Using Forecaster, we obtain an expected mass of $M={5.04}_{-1.91}^{+3.79}{M}_{\oplus }$, which yields an expected RV semiamplitude of $K={4.11}_{-1.56}^{+3.08}\,{\rm{m}}\ {{\rm{s}}}^{-1}$ assuming a circular orbit. Figure 10 shows the expected mass and RV-semiamplitude posteriors as calculated using Forecaster. Further, using these posteriors, Forecaster is capable of classifying the planet as Terran, Neptunian, Jovian, or stellar. From the posteriors in Figure 10, Forecaster classifies G 9-40b as Neptunian with 99.7% confidence.

Zoom In Zoom Out Reset image size

Second, we compare our predicted mass from Forecaster to the value from the MRExo toolkit. As is discussed in Kanodia et al. (2019), MRExo offers two different MR relations for forecasting: a Kepler planet–sample MR relation, and an M dwarf–planet MR relation. For this fit, we used the M dwarf–planet MR relation given the M2.5 spectral type of the host star, yielding a predicted mass of $M={3.94}_{-3.3}^{+11.5}{M}_{\oplus }$. Although both values are within the formal error bars of each other, we note that the median value from MRExo is smaller than the $M={5.04}_{-1.91}^{+3.79}{M}_{\oplus }$ value we obtain from Forecaster. Further, the MRExo value has a larger spread than the value from Forecaster. This larger spread is noted in Kanodia et al. (2019) as being due to the nonparametric fitting technique and the low number of M dwarf-only planets in the M dwarf–MR relation used by MRExo, resulting in a less-precise estimate.

Assuming a median semiamplitude of $4\,{\rm{m}}\ {{\rm{s}}}^{-1}$, this shows that G 9-40b can be followed up with a number of high-precision RV instruments already online or under construction (Fischer et al. 2016; Wright & Robertson 2017). As noted by Kanodia et al. (2019), there are rising statistical trends that suggest a difference between the exoplanet MR relation between M dwarfs and exoplanets orbiting earlier-type stars. However, the comparison of the two populations is still limited by the low number of M dwarf planets with both precise radii and masses. With a measurement of its mass, G 9-40b will yield further direct insights into those statistical trends.

To further constrain possible false-positive binary scenarios and to provide an upper bound of the mass of G 9-40b, Figure 11(a) shows the RVs of G 9-40 as a function of time as observed by HPF (Table 6 lists the RVs in tabular format). The blue points are individual 945 s exposures (${\sigma }_{\mathrm{unbinned}}\,=6.49\,{\rm{m}}\ {{\rm{s}}}^{-1}$), whereas the red points show weighted-average RVs per individual HET track (${\sigma }_{\mathrm{binned}}=5.32\,{\rm{m}}\ {{\rm{s}}}^{-1};$ effective exposure time of 2 × 945 s = 31.5 minutes). The associated error bars as derived from the SERVAL pipeline are 9.13 $\,{\rm{m}}\,{\rm{s}}$⁻¹ and 6.23 $\,{\rm{m}}\,{\rm{s}}$⁻¹ for the unbinned and binned points, respectively. Figure 11(b) shows the phase-folded RVs using our best-fit ephemeris from our joint fit in Table 4, showing that our RVs are predominantly located in phases between 0 and 0.4. Further, for comparison, also shown in Figures 11(a) and (b) is the predicted RV curve using our best prediction of the mass of G 9-40b given its radius using the Forecaster package ($M={5.04}_{-1.91}^{+3.79}{M}_{\oplus }$ and $K={4.11}_{-1.56}^{+3.08}\,{\rm{m}}\ {{\rm{s}}}^{-1}$) assuming a circular orbit.

Zoom In Zoom Out Reset image size

To estimate the upper bound on the mass of G 9-40b, we modeled the RVs in Figure 11 using the radvel Python package. For our RV model, we fixed the orbital period and transit center to the joint-fit transit ephemeris in Table 6, and given the small number of RV points, we additionally fixed the orbital eccentricity and argument of periastron to zero. This resulted in a two-parameter RV model using only the orbital semiamplitude (K) and a zero-point offset (γ) as the MCMC jump parameters. For the modeling, we placed minimal uninformative priors on the two parameters, only constraining K to be positive. To sample the parameter space, we initialized 100 MCMC walkers in the vicinity of the expected best-fit solution using the radvel MCMC interface. We removed the first marginally well-mixed iterations as burn-in, and we ran the MCMC chains until the Gelman–Rubin statistic was within 0.1% of unity, which we consider indicative of well-mixed chains. Figure 11(c) shows the resulting posteriors of both K and γ. From our fit, we see that the best-fit semiamplitude for the orbital semiamplitude is $K={0.8}_{-0.7}^{+2.1}\,{\rm{m}}\ {{\rm{s}}}^{-1}$ (black solid curve in Figures 11(a) and (b)), which is consistent with zero within the ∼1σ lower error estimate. Although the median estimated semiamplitude from this analysis is lower than our predicted semiamplitude from both Forecaster and MRExo as discussed above, it is within the 2σ error estimates. We attribute the low semiamplitude we estimate here to a combination of the small number of RV points and their sparse phase coverage of the full RV phase curve, and we argue that more RVs are needed to further test the robustness of this measurement. The semiamplitude posteriors reported here could be biased by uncorrected astrophysical systematics (e.g., stellar activity), instrumental systematics (e.g., persistence in the HPF H2RG detector; see Metcalf et al. 2019 or Ninan et al. 2019 for a discussion of systematics in NIR RVs from H2RG detectors), or model-mismatch systematics (e.g., an eccentric system, or another unknown planet could be present in the system). With these caveats in mind, we use the resulting posteriors on K to say that, given the current data set and assuming a circular orbit, the semiamplitude of G 9-40b is $K\lt 9.0\,{\rm{m}}\ {{\rm{s}}}^{-1}$ at 99.7% (3σ) confidence. Using our M = 0.290 ± 0.020M_⊙ mass estimate of the host star from Table 2, this translates to an upper limit mass constraint of M < 11.7M_⊕ at 99.7% confidence assuming a circular orbit. If instead we adopt the same priors as discussed above, but let the eccentricity and argument of periastron float (sampled as $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$ in radvel), our best-fit semiamplitude increases to $K={21}_{-12}^{+19}\,{\rm{m}}\ {{\rm{s}}}^{-1}$, due to the fit favoring a high-eccentricity solution of e ∼ 0.67 due to the sparse phase sampling of the RVs. Even in the unlikely scenario of such a high-eccentricity orbit, this translates to an upper mass limit of ∼0.5M_J at 99.7% confidence, showing that G 9-40 is a planetary-mass object. We prefer the e = 0 mass constraint, as we expect G 9-40 to be fully circularized. We estimated a circularization timescale of ∼150 Myr for G 9-40b, following the methodology in Bodenheimer et al. (2001), assuming a tidal quality factor of Q ∼ 10³ for a mini-Neptune mass planet (extrapolating between Q ∼ 10² and 10⁶ for Earth and Jupiter in the solar system; Goldreich & Soter 1966), which is substantially smaller than our SED estimated age of ${9.9}_{-4.1}^{+2.6}\,\,\mathrm{Gyr}$, as summarized in Table 2.

Table 6.  Radial Velocities of G 9-40 from HPF

$\,\mathrm{BJD}$TDB	RV	σRV	RV Drift Error	S/N
	(${\rm{m}}\,{\rm{s}}$−1)	$({\rm{m}}\ {{\rm{s}}}^{-1})$	$({\rm{m}}\ {{\rm{s}}}^{-1})$	@1100 $\,\mathrm{nm}$
2458543.6173586124	−7.78	6.23	0.02	140.03
2458543.6284200638	−10.90	5.96	0.02	148.49
2458544.622642529	0.94	8.99	0.31	107.54
2458544.6348675573	−7.47	8.94	0.18	106.49
2458562.782749085	3.98	7.65	0.20	131.85
2458562.794474252	11.27	9.04	0.22	113.95
2458590.6907325243	11.00	10.43	0.10	103.96
2458590.7019076305	−1.04	9.85	0.11	108.93

Notes. All exposures are 945 s in length. The derivation of the RV drift error—the contribution of the RV drift of HPF to the total estimated RV error—is described in the Appendix. The total RV error (σ_RV) is composed of the error bar given by the SERVAL RV pipeline added in quadrature to our RV drift error estimate.

A machine-readable version of the table is available.

Download table as:  DataTypeset image

Figure 12 compares G 9-40b to other known planets from the NASA Exoplanet Archive (Akeson et al. 2013) in the exoplanet radius, and exoplanet distance from the solar system plane. At a distance of $d=27.928\pm 0.045\,\mathrm{pc}$, G 9-40b is the second-closest transiting planet discovered by K2 to date, behind the K2-129 system from Dressing et al. (2017), which has a Gaia DR2 distance of ${27.796}_{-0.074}^{+0.075}\,\mathrm{pc}$.

Zoom In Zoom Out Reset image size

Being a nearby system, G 9-40 is bright at both visible and NIR wavelengths; for example, G 9-40 is brighter in the I band (I = 10.9) than GJ 1214 (I = 11.1), and only slightly fainter than GJ 1214 (J = 9.8) in the J band with a J magnitude of J = 10.1. Due to its relatively large transit depth and the brightness of its host star at NIR wavelengths, G 9-40b is amenable to transmission spectroscopic observations that could provide insight into the planet's bulk composition and formation history through measuring the elemental composition of its atmosphere and overall metal enrichment. To compare the applicability of performing precision transmission spectroscopic observations of G 9-40 during transit, we calculate the transmission spectroscopy metric (TSM) from Kempton et al. (2018) for G 9-40b and compare it to the corresponding TSM metric for other known exoplanets. Specifically, the TSM metric from Kempton et al. (2018) is defined as

$$\begin{eqnarray}&&\mathrm{TSM}=(\mathrm{Scale}\ \mathrm{Factor})\times \displaystyle \frac{{R}_{p}^{3}{T}_{\mathrm{eq}}}{{M}_{p}{R}_{* }^{2}}\times {10}^{-{m}_{J}/5},\end{eqnarray} \tag{ 5 }$$

where R_p is the radius of the planet in Earth radii, M_p is the mass of the planet in Earth masses, R_* is the radius of the host star in solar radii, T_eq is the equilibrium temperature in Kelvin of the planet calculated for zero albedo and full heat redistribution, and m_J is the magnitude of the host star in the J band. Kempton et al. (2018) define the scale factor for different radius bins where

$$\begin{eqnarray*}&&\begin{array}{l}(\mathrm{Scale}\ \mathrm{Factor})=0.19\qquad \mathrm{for}\qquad R\lt 1.5{R}_{\oplus },\\ (\mathrm{Scale}\ \mathrm{Factor})=1.26\qquad \mathrm{for}\qquad 1.5{R}_{\oplus }\lt R\lt 2.75{R}_{\oplus },\\ (\mathrm{Scale}\ \mathrm{Factor})=1.28\qquad \mathrm{for}\qquad 2.75{R}_{\oplus }\lt R\lt 4.0{R}_{\oplus }.\end{array}\end{eqnarray*}$$

Further, Kempton et al. (2018) define the TSM metric specifically for JWST/NIRISS as NIRSS gives more transmission spectroscopy information per unit observing time compared to the other JWST instruments for a wide range of planets and host stars (Batalha & Line 2017; Howe et al. 2017). For the planets with an unknown mass, we forecast the mass using the prescription suggested in Kempton et al. (2018), using the MR relationship of Forecaster as implemented by Louie et al. (2018).

Using this metric—with G 9-40b falling in the "small sub-Neptune" (1.5R_⊕ < R_p < 2.75R_⊕) planet size bin—G 9-40b has a high TSM metric of TSM = 96₋₄₁⁺⁶⁴, assuming a predicted mass of M = 5.04M_⊕ using Forecaster, and exceeds the TSM > 90 threshold for high-quality priority atmospheric targets suggested by Kempton et al. (2018). Figure 13 compares the TSM of G 9-40b to the TSM for other currently known planets orbiting M dwarfs (${T}_{\mathrm{eff}}\lt 4000\,{\rm{K}}$ ³¹ ) with radii less than 4.0R_⊕ as a function of the host star J magnitude. This plot shows that G 9-40b is currently among the top best small (<4.0R_⊕) M dwarf planets for transmission spectroscopic observations, after GJ 1214b (Charbonneau et al. 2009; TSM ∼ 615), the newly discovered mini-Neptune (R = 1.57R_⊕) L98-59d orbiting a nearby M3 dwarf (Cloutier et al. 2019; TSM ∼ 230), the Neptune-sized (R = 3.5R_⊕) M dwarf planet K2-25b in the Hyades (Mann et al. 2016; TSM ∼ 140), and the newly discovered sub-Neptune (R = 2.42R_⊕) TOI 270c (Günther et al. 2019; TSM ∼ 130).

Zoom In Zoom Out Reset image size

We note that even though G 9-40b has a high TSM metric, the atmosphere of G 9-40b has the possibility of having damped atmospheric features due to its low equilibrium temperature of T_eq = 319 $\,{\rm{K}}$ (a = 0.3). Past work has shown that hotter planets tend to have larger spectral features, while cooler planets tend to have less prominent features due to hazes (Crossfield & Kreidberg 2017). This has been attributed to the fact that that at temperatures below ∼1000 K, methane is abundant and can easily photolyze to produce carbon hazes that mute atmospheric features. Even if the atmosphere of G 9-40 is observed to be featureless, as noted by Kempton et al. (2018), the study of a large sample of sub-Neptune exoplanet atmospheres is especially important because these planets do not exist in our solar system, and therefore no well-studied benchmark objects exist. Further, Kempton et al. (2018) argue that a large sample of sub-Neptune planets is needed for precise atmospheric characterization, to adequately study the high degree of diversity expected for their atmospheric compositions.