Updated Parameters and a New Transmission Spectrum of HD 97658b

HD 97658b was observed on 2013 December 19 (visit1) and 2014 January 7 (visit2) under the HST program 13501 (PI: Knutson), and then on 2016 April 20 (visit3) and 2017 January 31 (visit4) under the HST program 13665 (PI: Benneke). A spatial scan mode is adopted for all four visits to accumulate abundant photons on the detector. We use the "round trip" scan method, which means the scan is carried out in two opposite directions alternatively (one direction per exposure) during the observations, and during each scan (one exposure), the image on the detector is read out several times. Each visit of program 13501 was observed with a 256 × 256 pixel subarray, consisting of around one hundred 14 s exposures in each scan direction and three read outs in each exposure, while each visit of program 13665 was observed with a 512 × 512 pixels subarray, consisting of around one hundred 23 s exposures in each scan direction and two read outs in each exposure. Two typical raw images from program 13501 and program 13665 are shown in Figure 1, and raw data from all visits are processed with a standard pipeline described as follows.

Zoom In Zoom Out Reset image size

We start the data analysis from the HST/WFC3 IMA files. For each exposure, we mask out bad pixels that were identified with flags from the WFC3/IR bad pixel table (Hilbert 2012). Then we select a "clean" rectangle on the image—where detected flux roughly flattens spatially—as the background and take the median flux in that rectangle as our background flux level f_bkg. Next, we subtract f_bkg from the whole image and correct the image with the data cube from STScI.

Next, mean flux along the scanning direction is calculated for each wavelength pixel in each exposure; thus, a raw spectrum is produced for each exposure. And finally, we correct for the wavelength shift on the detector over time by picking the spectrum from the first exposure as the template and shifting all subsequent exposures along the dispersion direction to match the template.

The wavelength solutions, which translate pixel values on the detector to wavelengths, are calculated using the wavelength calibration coefficients from STScI, and we select the range of pixels in the dispersion direction such that our wavelength range coincides with that chosen by Knutson et al. (2014) for easy comparison. We then obtain the raw white light curve of each visit by summing up all flux in the dispersion range in each exposure.

Lastly, we correct for the difference in flux levels between the two scan directions of each visit, which results from the dependence of the total flux on the vertical position of the spectrum relative to the middle row of the detector (the "up-stream/down-stream" effect; McCullough & MacKenty 2012). Following a similar process to the one described in Tsiaras et al. (2016), we fit the flux profile in the scan direction of each exposure with a box shape to extract its read-out length and mid-position, which should be related by two different linear relations for the two scan directions. Then, we fit the linear relations and choose the first exposure as our reference point to scale the flux of all subsequent exposures to the level where they should be if they all had the same read-out mid-position as the first exposure. Our final white light curve of each visit is shown in Figure 2.

Zoom In Zoom Out Reset image size

Major systematic effects in the HST/WFC3 observation data include: (1) the "ramp" effect in each HST orbit, which is thought to be caused by free charge carriers trapped in the depletion regions of the detector; (2) visit-long trends, which is a quasi-linear trend across the entire observation period; and (3) "HST breathing effect", caused by the spacecraft temperature variation during each orbital period of HST (Wakeford et al. 2016).

2.2.1. Orbital Ramp Effects

2.2.2. Residual Systematics and Noise Modeling

2.2.3. Visit-long Trend

Visit-long trend corrections can be combined into the base-level flux parameter f in the RECTE model, which becomes (f₀ + bt) in a linear-shape visit-long trend case, where f₀ and b are free parameters, and t represents the time stamp of each exposure. Although most previous works have used a simple linear trend (Knutson et al. 2014; Evans et al. 2018), in this data set of HD 97658b, we observe evidence of a long-term trend that deviates from a linear form from visual examination of the white light curves (Figure 2), which may have been induced or magnified by the high brightness of the host star that is close to the WFC3 detector saturation limit.

Therefore, we try four common functional forms to correct for the visit-long trend:

• (1)  
  linear, with the form (f₀ + bt);
• (2)  
  quadratic, with the form (f₀ + bt + ct²);
• (3)  
  exponential, with the form $({f}_{0}-b\exp (-t/c));$
• (4)  
  logarithmic, with the form $({f}_{0}-b\mathrm{log}(t+c))$.

When applying these functional forms, we make use of the first orbit to help anchor the trend by modeling the second half of its data points, which are less affected by the ramp effect because all charge traps should have been filled. To sum up, the orbits after the first orbit are modeled with function F₁ = f_ramp × f_trend × f_transit, where f_ramp, f_trend, and f_transit represent the ramp-effect model, the visit-long trend model, and the transit model, respectively, while the second half of the first orbits are modeled with the function F₀ = f_trend × f_transit. The data from the first half of the first orbits are discarded.

Simultaneously with instrumental systematic models, we fit the transit signal with models generated using the BAsic Transit Model cAlculatioN (BATMAN) package provided by Kreidberg (2015). In our joint fit of the white light curves, the transit depth, inclination, and orbital semimajor axis are tied to be the same free parameters for all four visits, while the mid-transit time of each visit is a separate free-floating parameter. We adopt the fixed stellar parameters reported in Gaia Data Release 2 (Gaia Collaboration et al. 2018), and the two-parameter limb-darkening coefficients are extracted from Claret & Bloemen (2011) assuming these stellar parameters are fixed in our models.

The best-fit transit model and detrended white light curves of all four visits assuming a logarithmic shape visit-long trend are shown in Figure 3, with input parameters and best-fit output parameters shown in Table 2. We can see the model does not fit the second half of the first orbit very well, which is expected because the ramp-effect component was excluded from the modeling of this section of the data, and this section of the data was only in the analysis to provide an anchor for modeling the visit-long trends. Apart from this, there are a few outliers in the detrended visit1 light curve, but they do not bias our overall science results of this planet. The best-fit white light curves assuming other forms of visit-long trends are of similar quality, and we discuss those results as follows.

Zoom In Zoom Out Reset image size

Table 2.  Input and Best-fit Parameters from HST/WFC3 White Light-curve Fitting

Parameter	Symbol	Value	Unit
Input fixed Parameters
Orbital period	P	9.4903	days
Eccentricity	e	0.078
Argument of periapsis	ω	90.0	degree
Quadratic limb-darkening coefficients	u	[0.246, 0.203]
Best-fit of output parameters
Radii ratio	Rp/Rs	0.0293 ± 0.0001
Mid-transit time visit1	T0,visit1	2456646.4829 ± 0.0011	BJD
Mid-transit time visit2	T0,visit2	2456665.4621 ± 0.0012	BJD
Mid-transit time visit3	T0,visit3	2457491.0312 ± 0.0011	BJD
Mid-transit time visit4	T0,visit4	2457785.2021 ± 0.0011	BJD
Semimajor axis ratio	a/Rs	26.7 ± 0.4
Inclination	i	89.6 ± 0.1

Download table as:  ASCIITypeset image

Using the raw white light curve from visit4 as an example, we show the four best-fit visit-long trend models in Figure 4. The first thing we examine is whether those models produce consistent transit depths for different visits. Figure 5 shows the best-fit R_p/R_s of each HST visit assuming four different visit-long systematic trends. The logarithmic trend produces the most internally consistent white light transit depths—all well within 1σ uncertainty with respect to each other. And while the other three cases all have maximum transit depth differences around 1σ uncertainty, the quadratic trend model produces the largest transit depth uncertainties and the most discrepant results among visits. This result is expected since the curvature of a quadratic model is much more sensitive to its parameters than that of the three other models, and a transit signal itself could approximately be fitted with a quadratic shape if it is blended at a high systematic level. Similar to this work, Agol et al. (2010) compared long-term systematic trend functions to be used to model Spitzer light curves and found that a quadratic function could bias transit depth measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The white light-curve fitting result comparison does not strictly rule out any of the four visit-long trend models, even though the quadratic model is disfavored. Hence, we proceed to calculate the transmission spectra assuming these four different models and make further comparison of their corresponding spectra in the following sections.

We group pixels along the dispersion direction into 28 channels with the same wavelength boundaries as those in Knutson et al. (2014) for easy comparison and extract raw light curves in each channel. One thing to notice is that the long wavelength end of the dispersed flux ran over the edge of the detector in visit3 and visit4 because of an observation error, as is shown in Figure 1. Therefore, the last two channels of these two visits were not observed. The channel boundaries compared against the sensitivity curve are shown in Figure 6.

Zoom In Zoom Out Reset image size

We use the divide-by-white technique to analyze the spectral light curves. This technique was introduced by Kreidberg et al. (2014a) and has been widely used in HST/WFC3 transmission spectrum measurements (Knutson et al. 2014; Damiano et al. 2017; Tsiaras et al. 2018). First, a common-mode signal is generated by dividing the white light curve by the best-fit transit signal, and then this common mode is injected as a systematic component of spectral light curves in each channel. This systematic component is multiplied by a new transit model with the mid-transit time fixed to the best-fit transit time of the white light curve. And since we notice that the light curves in different channels have different visit-long trends, we also multiply the model of spectral light curves by a parameterized visit-long trend with the same form (linear, exponential, logarithmic, or quadratic) as that of the white light curve.

After extracting the best-fit R_p/R_s from each channel in each visit, we calculate the R_p/R_s spectra over four visits using the weighted mean on each spectral channel. The uncertainty of R_p/R_s in each channel is calculated by taking the larger value of either σ₁(k) or σ₂(k), where k represents the channel number, σ₁(k) represents the standard deviation of the transit depth values from four visits weighted by their uncertainties, and σ₂(k) represents the standard deviation of the mean transit depth of four visits. Formulae to calculate σ₁(k) and σ₂(k) are shown as follows:

$$\begin{eqnarray}&&{\sigma }_{1}(k)=\sqrt{\displaystyle \frac{{\sum }_{i=1}^{N}{w}_{i}(k){\left({D}_{i}(k)-\bar{D}(k\right)}^{2}}{\tfrac{\left(N-1){\sum }_{i=1}^{N}{w}_{i}(k\right)}{N}}},\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{\sigma }_{2}(k)=\sqrt{\displaystyle \frac{{\sum }_{i=1}^{N}{{\sigma }_{i}(k)}^{2}}{N(N-1)}},\end{eqnarray} \tag{ 2 }$$

where N is the total number of visits, D_i(k) is the transit depth in channel k from visit i, σ_i(k) is the uncertainty in D_i(k), $\bar{D}(k)$ is the weighted average of D_i(k) defined as $\bar{D}(k)=\tfrac{{\sum }_{i=1}^{N}{w}_{i}{D}_{i}(k)}{{\sum }_{i=1}^{N}{w}_{i}}$, and w_i(k) is the weight on D_i(k) defined as ${w}_{i}(k)=\tfrac{1}{{{\sigma }_{i}(k)}^{2}}$.

The final averaged R_p/R_s results in each channel, assuming linear, quadratic, logarithmic, and exponential visit-long trend models, are presented in Table 3.

Table 3.  R_p/R_s Averaged over Four Visits in Each WFC3 Bandpass Assuming Four Different Visit-long Trends

Bandpass Center	Linear		Logarithmic		Exponential		Quadratic
(μm)	Rp/Rs	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$	Rp/Rs	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$	Rp/Rs	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$	Rp/Rs	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$
1.132	0.0308	0.0003	0.0292	0.0003	0.0285	0.0003	0.0286	0.0004
1.151	0.0306	0.0004	0.0290	0.0004	0.0280	0.0009	0.0283	0.0006
1.170	0.0305	0.0005	0.0293	0.0004	0.0287	0.0003	0.0282	0.0005
1.188	0.0308	0.0005	0.0294	0.0004	0.0290	0.0003	0.0274	0.0003
1.207	0.0300	0.0004	0.0286	0.0002	0.0278	0.0003	0.0273	0.0003
1.226	0.0303	0.0004	0.0287	0.0003	0.0283	0.0003	0.0277	0.0004
1.245	0.0302	0.0003	0.0288	0.0002	0.0283	0.0003	0.0276	0.0003
1.264	0.0299	0.0004	0.0286	0.0003	0.0282	0.0005	0.0268	0.0006
1.283	0.0301	0.0005	0.0289	0.0003	0.0281	0.0003	0.0278	0.0007
1.301	0.0307	0.0002	0.0293	0.0002	0.0286	0.0002	0.0276	0.0004
1.320	0.0303	0.0002	0.0290	0.0002	0.0282	0.0002	0.0272	0.0005
1.339	0.0300	0.0003	0.0287	0.0002	0.0280	0.0003	0.0268	0.0004
1.358	0.0302	0.0002	0.0291	0.0002	0.0284	0.0003	0.0268	0.0007
1.377	0.0302	0.0003	0.0288	0.0002	0.0280	0.0003	0.0273	0.0005
1.396	0.0308	0.0004	0.0294	0.0005	0.0289	0.0005	0.0275	0.0006
1.415	0.0314	0.0004	0.0302	0.0002	0.0298	0.0004	0.0285	0.0006
1.433	0.0308	0.0002	0.0294	0.0002	0.0288	0.0002	0.0276	0.0003
1.452	0.0306	0.0003	0.0292	0.0002	0.0283	0.0010	0.0281	0.0003
1.471	0.0303	0.0004	0.0291	0.0004	0.0285	0.0002	0.0273	0.0003
1.490	0.0302	0.0002	0.0290	0.0003	0.0280	0.0003	0.0273	0.0007
1.509	0.0305	0.0003	0.0291	0.0003	0.0282	0.0006	0.0267	0.0006
1.528	0.0302	0.0002	0.0288	0.0004	0.0284	0.0004	0.0270	0.0008
1.546	0.0302	0.0002	0.0288	0.0002	0.0283	0.0004	0.0262	0.0003
1.565	0.0304	0.0002	0.0290	0.0003	0.0280	0.0003	0.0258	0.0006
1.584	0.0302	0.0005	0.0295	0.0002	0.0290	0.0003	0.0279	0.0006
1.603	0.0314	0.0002	0.0298	0.0002	0.0297	0.0010	0.0288	0.0006
1.622	0.0311	0.0009	0.0299	0.0007	0.0294	0.0007	0.0277	0.0005
1.641	0.0312	0.0008	0.0299	0.0005	0.0291	0.0005	0.0294	0.0005

Note. We list the results using all four different visit-long trend formulae, but only the transit depth results using the logarithmic visit-long trend are adopted in our atmospheric retrieval and modeling analysis. The reason for this choice and comparisons between results using different visit-long trends are presented in Section 2.2.3.

Download table as:  ASCIITypeset image

We plot the transmission spectra for all four visit-long trends in Figure 7, where all spectra are shifted to have zero mean for easy comparison. It is apparent that the spectra assuming three different visit-long trends have highly consistent shapes, with the only difference between them being their mean (white light) transit depth, whereas the spectral shape resulting from a quadratic visit-long trend is distinctive from the other three. Taking into consideration the white light-curve transit depth comparisons shown in Figure 5, we decide that using a quadratic shape to model the HST/WFC3 visit-long trend systematics of this data set could have deformed the resulting transmission spectrum shape, so we discard this model in the rest of this work. On the other hand, transmission spectral features on the 1.1–1.7 μm wavelength range extracted by assuming the other three visit-long trend models are consistent with each other. Therefore, we adopt the spectrum assuming a logarithmic visit-long trend, which shows the most consistent white light-curve transit depths among different visits (Figure 5), but we set the mean depth of the transmission spectrum on this wavelength range as a free parameter when performing atmosphere retrieval, so that uncertainties from the visit-long trend model selection are included in the final error budget.

Zoom In Zoom Out Reset image size

Knutson et al. (2014) fitted the light curves from visit1 and visit2 with an exponential orbital ramp model and a linear visit-long model. In Figure 8, we compare our transmission spectra (assuming a linear trend) averaged over visit1 and visit2 (pink shaded region) and averaged over all four visits (blue shaded region) with the final spectrum from Knutson et al. (2014, gray shaded region). We can see that our two-visit-averaged spectrum shows a similar rising trend to the spectrum from Knutson et al. (2014). Although there are some shifts in transit depths between our spectra and the previous one, their differences are mostly within 1σ, and the reduced χ² between our two-visit-averaged spectrum and the previous work is 0.98, showing a high consistency between the two results. However, for the spectrum averaged over all four visits, the rising trend is mitigated, except for the possible feature redward of 1.6 μm. Therefore, it is reasonable to speculate that the rising trend of the HD 97658b spectrum presented in Knutson et al. (2014) results from an unidentified systematic effect. Nonetheless, there could be astrophysical features remaining in the spectrum, and we discuss atmospheric retrieval results in Section 5.

Zoom In Zoom Out Reset image size

We observed three transits of HD 97658b with HST/STIS using the G750L grism (0.524–1.027 μm) as part of HST program 13665. Like HST/WFC3, an initial orbit is required to settle the telescope. The detector was purposefully saturated by about a factor of three in the brightest part of the spectrum to increase the signal-to-noise ratio (S/N; Gilliland et al. 1999). Data from each of the three visits was reduced using the same steps as Lothringer et al. (2018). Counts were added along columns to reconstruct the observed flux at each wavelength. The observations were then split into 10 bins of approximately 0.05 μm covering the G750L bandpass, and a transit and noise model were fit to the data assuming orbital parameters from Knutson et al. (2014). Unlike the systematics marginalization procedure that was used in Lothringer et al. (2018), we instead fit transit models to the data using GPs in a similar fashion to Bell et al. (2017) and Evans et al. (2013, 2018), using a squared-exponential kernel. Using GPs instead of a parametric approach allows us to include time as a nonparametric covariate, as it was found that assuming a linear slope in time (as was done in Lothringer et al. 2018) produced worse fits to the data and unrealistically small error bars. By including time as a covariate in the GP, we can better reflect our uncertainty in the baseline flux, which becomes more apparent over four science orbits rather than the three in most other STIS data sets (also see Demory et al. 2015). HST's orbital phase was also used as a covariate to account for orbit-to-orbit systematics.

In Table 4, we report the spectroscopic transit depths and the mid-transit time from each of the three transits observed with STIS, along with their uncertainties for future reference. The STIS data quality limits our ability to achieve consistent transit depths among different observations. As Table 4 shows, the transit depths measured with data from three STIS observations are highly discrepant in the long wavelength part (λ > 0.79 μm) of the STIS bandpass. And with the same uncertainty calculation procedure as described in Section 2.3, we find that the transit depth uncertainties can be as large as 100–200 ppm in the red end channels of the STIS bandpass. This uncertainty level is far from adequate for providing any meaningful constraint to atmospheric properties. Therefore, we do not include the STIS data set in our atmospheric analysis.

Table 4.  Best-fit Spectroscopic Transit Depths for Each Transit Observed with STIS

Bandpass Center	Observation 1		Observation 2		Observation 3
(μm)	${R}_{{\rm{p}}}/{R}_{{\rm{s}}}$	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$	${R}_{{\rm{p}}}/{R}_{{\rm{s}}}$	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$	Rp/Rs	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$
0.553	0.0273	0.0015	0.0255	0.0011	0.0252	0.0028
0.601	0.0282	0.0006	0.0272	0.0007	0.0277	0.0009
0.650	0.0281	0.0004	0.0282	0.0004	0.0290	0.0005
0.699	0.0278	0.0009	0.0281	0.0004	0.0290	0.0020
0.748	0.0282	0.0010	0.0302	0.0011	0.0308	0.0019
0.797	0.0268	0.0015	0.0286	0.0032	0.0315	0.0011
0.845	0.0265	0.0014	0.0294	0.0009	0.0312	0.0008
0.894	0.0246	0.0016	0.0275	0.0014	0.0315	0.0014
0.943	0.0260	0.0026	0.0294	0.0023	0.0326	0.0017
0.992	0.0213	0.0032	0.0346	0.0034	0.0378	0.0049
Mid-transit Times (BJD)	2457149.4191 ± 0.0005		2457196.8642 ± 0.0003		2457206.3537 ± 0.0008

Download table as:  ASCIITypeset image

Six transits were observed with the IRAC 3.6 μm channel, and five transits were observed with the 4.5 μm channel from 2014 July to 2016 April under Spitzer program 11131 (Dragomir et al. 2014). Each transit was observed with 0.08 s exposure per frame and approximately 0.13 s per frame cadence.

We analyze the raw light curves using the "pixel-level decorrelation" (PLD hereafter) technique (Deming et al. 2015), following the same procedure as described in Guo et al. (2019). Before fitting the light curves, we manually remove the first 50–80 minutes of each observation so that the drastic systematic ramp at the beginning of each Spitzer observation does not bias our fitting result. Using the PLD technique, fractional contributions to the total flux at the observational time points from each selected pixel are treated as an eigenvector, and we set the weight of each eigenvector as a free parameter and combine them together to model the total flux variations. Since the PLD technique assumes that the incoming stellar flux is falling on the same set of pixels throughout the entire time series, we include all pixels around the star that contribute more than 1% of the total flux to ensure that the PLD technique is valid to apply. The same pixel selection procedure is also successfully used in Guo et al. (2019). Simultaneously, we fit a transit model generated using BATMAN and defined with a free-floating transit depth and mid-transit time. The rest of the transit parameters are fixed in the same way as described in Section 2.2.2.

The best-fit transit models and detrended light curves of all Spitzer transits are shown in Figure 9. The best-fit parameters of all transits are presented in Table 5.

Zoom In Zoom Out Reset image size

Table 5.  Best-fit R_p/R_s and Mid-transit Times of Each Transit Observed by Spitzer

AOR ID	Channel	T0 (BJD-2400000.5)	${\sigma }_{{T}_{0}}$ (days)	Rp/Rs	${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$
49696512	3.6 μm	56864.23938	0.00027	0.0276	0.0005
49697536	3.6 μm	56883.21892	0.00027	0.0270	0.0004
49698048	3.6 μm	56892.70742	0.00026	0.0280	0.0004
52197888	3.6 μm	57091.98243	0.00034	0.0261	0.0004
52198144	3.6 μm	57082.49324	0.00032	0.0275	0.0006
52198656	3.6 μm	57101.47164	0.00034	0.0277	0.0005
53908736	4.5 μm	57481.04364	0.00031	0.0272	0.0005
53909248	4.5 μm	57471.55445	0.00025	0.0281	0.0004
53909504	4.5 μm	57253.30001	0.00045	0.0284	0.0005
53909760	4.5 μm	57243.81157	0.00022	0.0293	0.0004
53910016	4.5 μm	57234.32179	0.00026	0.0284	0.0004

Download table as:  ASCIITypeset image

Figure 10 shows the R_p/R_s of each transit arranged according to their mid-transit time. Discrepancies among transits of the same channel are around 1σ. To ensure that the scatter in transit depths and the associated error bars of each transit are properly included in our error budget, we again take the larger one between the standard deviation of the weighted average of center values and the weighted reduced error bars as our uncertainties in R_p/R_s. We find that the best-fit values of R_p/R_s in the 3.6 μm and 4.5 μm channels are 0.0273 ± 0.0003 and 0.0284 ± 0.0003, respectively.

Zoom In Zoom Out Reset image size

MOST (Walker et al. 2003; Matthews et al. 2004) is a microsatellite carrying a 15 cm optical telescope that acquires light through a broadband filter spanning the visible wavelengths from 350 to 700 nm. It is in a Sun-synchronous polar orbit with a period of 101.4 minutes, which allows it to monitor stars in a Continuous Viewing Zone (CVZ) without interruption for up to 8 weeks. The CVZ covers a declination range of +36° > δ > −18°.

HD 97658 was observed by MOST in Direct Imaging mode, in which defocused images of the stars were projected directly onto the CCD (Rowe et al. 2006). One, four, and five transits were observed in 2012, 2013, and 2014, respectively. The 2012 and 2013 transits have been previously published in Dragomir et al. (2013), while the 2014 observations are unpublished. For the analysis performed for this paper, we used the three transits observed in 2013 that do not show interruptions (March 10, 19, and 29; see Dragomir et al. 2013), and all five transits that we observed in 2014 (all of which are also continuous).

The exposure times were all 1.5 s, and the observations were stacked on board the satellite in groups of 21 for a total integration time of 31.5 s per data point. Raw light curves were extracted from the images using aperture photometry (Rowe et al. 2008). Outlier clipping and de-trending from the sky background and position on the CCD were performed as described in Dragomir et al. (2013). After these steps, a straylight variation at the orbital period of the satellite remained. This variation was filtered by folding each light curve on this 101.4 minute period, computing a running average from this phased photometry, and removing the resulting waveform from the corresponding light curve.

We fit the eight transits simultaneously using EXOFASTv2 (Eastman 2017), a differential evolution Markov Chain Monte Carlo (MCMC) algorithm that uses error scaling, and we obtained a best-fit R_p/R_s value of ${0.02866}_{-0.00056}^{+0.00054}$. We summarize the ephemeris of these eight transits into two best-fit mid-transit times: 2456361.8050 ± 0.0033 (BJD) in year 2013 and 2456712.9096 ± 0.0024 (BJD) in year 2014. The detrended light curves and the best-fit transit models are plotted in Figure 11.

Zoom In Zoom Out Reset image size

Since 1997 January, we have collected 553 radial velocity measurements with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the Keck I Telescope on Maunakea and 215 measurements with the Levy spectrograph on the Automated Planet Finder at Lick Observatory (APF; Radovan et al. 2014; Vogt et al. 2014). These data were all collected through an iodine cell for wavelength calibration and point-spread function reference (Butler et al. 1996). One set of iodine-free spectra were collected with each instrument to use as a model of the intrinsic stellar spectrum. The HIRES data were often taken in sets of three due to the short ∼2 minute exposures to mitigate the effects of stellar oscillations; this was not necessary for the APF due to the smaller aperture and longer exposure times (∼10–20 minute exposures). The HIRES data from 2005 January to 2010 August were previously analyzed in the discovery paper of HD 97658 b (Howard et al. 2011). The data reduction and analysis followed the California Planet Search method described in Howard et al. (2010). The resultant radial velocities are presented in Table 6 and in Figure 15.

We first investigate the star for signs of stellar activity by examining the Calcium H and K lines (S_HK; Isaacson & Fischer 2010, Figure 13) in the HIRES and APF data. There is a clear periodicity in both the S_HK and the radial velocity data around 3500 days in the HIRES data set (Figure 14). The APF data does not have a long enough baseline to detect such a long signal. In addition, we compare this long-term variation in S_HK and radial velocity signals with the brightness and color variation of HD 97658, which was measured with the Fairborn T8 0.80 m automatic photoelectric telescope (Henry 1999). As is shown in Figure 13, there is a clear correlation in the variations seen in radial velocities, stellar activity data, stellar brightness, and stellar color. This relation implies that the long-term radial velocity variation is actually caused by stellar activity. The length of the signal indicates that it is likely the star's 9.6 yr magnetic activity cycle (slightly shorter than our Sun's 11 yr cycle); we discuss the stellar activity in more detail in Section 4.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We analyze the radial velocity data using RadVel²⁴ (Fulton & Petigura 2017), which models Keplerian orbits to fit radial velocity data by performing a maximum-likelihood fit to the data and then subsequently determining the uncertainties through an MCMC analysis. We use the default MCMC parameters for RadVel of: 50 walkers, 10000 steps, and 1.01 as the Gelman-Rubin criteria for convergence, as described in Fulton & Petigura (2017).

We model this system in RadVel as a two-Keplerian system for planet b and the stellar activity. We include priors on the transit parameters of planet b from Van Grootel et al. (2014). We incorporate this stellar activity signal at around 3500 days into our radial velocity fit as an additional Keplerian signal because it has a sinusoidal shape, and only two cycles of this signal are captured by the data. We use a Gaussian prior on the period (3424 ± 41 days) and reference phase of this signal (2457372 ± 21 BJD) from a RadVel 1-Keplerian fit of the HIRES S_HK data. Our radial velocity fit is shown in Figure 15, and the output parameters are listed in Table 7. Note that the planet mass is calculated assuming our best-fit inclination (i = 89.6) with the HST/WFC3 data set, which shows sin(i) ≈ 1.

Zoom In Zoom Out Reset image size

Table 7.  Radial Velocity Fit Parameters

Parameter	Name (Units)	Value
Planet Parameters
Pb	Period (days)	9.49073 ± 0.00015
Tconjb	Time of conjunction	2456361.805 ± 0.0006
	(BJDTDB)
eb	Eccentricity	≡0.0
ωb	Argument of periapse	≡0.0
	(rad)
Kb	Semi-amplitude (m s−1)	2.81 ± 0.15
Mb	Mass (M⊕)	${7.81}_{-0.44}^{+0.55}$
ρb	Density (g cm−3)	${3.78}_{-0.51}^{+0.61}$
Pactivity	Period (days)	${3652}_{-120}^{+130}$
Tconjactivity	Reference Time	${2457605}_{-89}^{+100}$
	(BJDTDB)
eactivity	Eccentricity	≡0.0
ωactivity	Argument of periapse	≡0.0
	(rad)
Kactivity	Semi-amplitude (m s−1)	1.96 ± 0.21
Other Parameters
γHIRES	Mean center-of-mass	−0.85 ± 0.17
	velocity (m s−1)
γAPF	Mean center-of-mass	$-{0.42}_{-0.34}^{+0.33}$
	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 (ms−1)	${2.93}_{-0.1}^{+0.11}$
σAPF	Jitter (${\rm{m}}\ {{\rm{s}}}^{-1}$)	${1.3}_{-0.35}^{+0.31}$

Download table as:  ASCIITypeset image

We also test a nonzero planet eccentricity for completeness; the resulting eccentricity is small, consistent with zero to 2σ (e_b = ${0.030}_{-0.021}^{+0.034}$), and the results are consistent with the planet parameters of the circular case. Therefore, we adopt the circular fit results. We also test including a GP to model the stellar activity signal with the hyper-parameters constrained from a fit of the HIRES S_HK data. The results are consistent with the Keplerian fit; the baseline covers only two cycles of the activity; therefore, the deviation from a simple sinusoid is small. Since the fit has consistent posteriors, the additional parameters needed for the GP fit do not seem warranted, and we present the Keplerian fit as our final result.

As described above, our Keck-HIRES spectra allow us to identify a 9.6 yr stellar activity cycle for HD 97658. In addition, we also estimate the star's rotation period by calculating a Lomb–Scargle periodogram of the Keck/HIRES and APF activity indices after removing the long-term variations induced by the stellar activity cycle. No single period dominates, but we see an excess of power from 32–36 days in both data sets. We therefore report a stellar rotation period of P_rot = 34 ± 2 days, slightly lower than previously reported (Henry et al. 2011).

HD 97658's rotation period is typical for stars with similar photometric colors and activity levels (as measured by ${R}_{\mathrm{HK}}^{{\prime} }$; Suárez Mascareño et al. 2016). The activity cycle is also typical, with HD 97658 falling nicely on the empirical relation between P_rot and P_activity/P_rot (Baliunas et al. 1996; Suárez Mascareño et al. 2016). Both HD 97658's overall activity level (${R}_{\mathrm{HK}}^{{\prime} }$ and S_HK) and the 2 m s⁻¹ RV variations induced by its activity mark it as a quieter-than-average star for its spectral type (Isaacson & Fischer 2010; Lovis et al. 2011).

Low-mass stars (M or K dwarfs) often have excessive X-ray and UV radiation from their chromosphere and corona, and these energetic emissions can drive photochemistry and ionization processes in atmospheres of the planets orbiting around them. Loyd et al. (2016) put together a catalog of seven M dwarfs and four K dwarfs, including HD 97658, in their MUSCLES Treasury Survey, and they obtained Chandra or XMM-Newton observations of each of them. An interline continuum in the FUV bandpass is detected at 6.3σ significance in HD 97658. No observation of the X-ray bandpass was made on this star; although, integrated X-ray flux was detected higher than 10⁻¹⁴ erg s⁻¹ cm⁻² for all three of the other K dwarfs in the survey. Since UV and X-ray radiation are strongly related to the dissociation of atmosphere molecules including H₂O, CH₄, CO, O₃, etc. (Rugheimer et al. 2013) and production of hazes (Hörst et al. 2018), atmospheric models ought to take the effect of high-energy stellar radiation into account. And to understand the atmospheric compositions and evolutionary history of HD 97658b, more detailed stellar spectroscopy in full wavelength coverage needs to be conducted.

We searched for additional planet candidates orbiting HD 97658 by applying an iterative periodogram algorithm to our radial velocity data. First, we define an orbital frequency/period grid over which to search, with sampling such that the difference in frequency between adjacent grid points is (2πτ)⁻¹, where τ is the observational baseline. Using this grid, we compute a goodness-of-fit periodogram by fitting a sinusoid with a fixed period to the data, for each period in the grid. We choose to measure goodness-of-fit as the change in the Bayesian Information Criterion (BIC) at each grid point between the best-fit one-planet model with the given fixed period and the BIC value of the zero-planet fit to the data. We then fit a power law to the noise histogram (50%–95%) of the data and accordingly extrapolate a BIC detection threshold corresponding to an empirical false-alarm probability of our choice (we choose 0.003).

If one planet is detected, we perform a final fit to the one-planet model with all parameters free, including eccentricity, and record the BIC of that best-fit model. We then add a second planet to our RV model and conduct another grid search, leaving the parameters of the first planet free to converge to a more optimal solution. In this case, we compute the goodness-of-fit as the difference between the BIC of the best-fit one-planet model, and the BIC of the two-planet model at each fixed period in the grid. We set a detection threshold in the manner described above and continue this iterative search until the (n+1)th search rules out additional signals. This search method is known as a recursive periodogram, also described in Anglada-Escudé & Tuomi (2012). Similar to our RV analysis described in the previous section, we use RadVel (Fulton & Petigura 2017) to fit Keplerian orbits, and we use an implementation of this search algorithm known as RVSearch, to be released at a later date (L. Rosenthal et al. 2020, in preparation). For HD 97658, we search from 1.5 days to five times the observational baseline and detect no new planet candidates with significance higher than FAP = 0.01.

PLATON is a forward modeling and atmosphere retrieval tool for exoplanet transmission spectral analysis, developed by Zhang et al. (2019). PLATON uses the same opacity files and part of the same algorithms as the widely used atmosphere forward modeling package Exo-transmit (Kempton et al. 2017), but it is 100–1000 times faster than Exo-transmit, so that an MCMC or a nested sampling retrieval method can be used to extract the posteriors of atmospheric parameters.

We first retrieve the HD 97658b atmospheric parameters with the observed transmission spectrum on the HST/WFC3 bandpass only. The input stellar radius and effective temperature are set to the latest published value in Gaia DR2, with R_* = 0.746 R_⊙ and T_eff = 5192 K, and the input planet mass is set according to our radial velocity results (Section 4) with M_p = 7.81 M_⊕. We let five atmospheric parameters—planet radius R_p, planet surface temperature T_p assuming an isothermal atmosphere, logarithmic metallicity $\mathrm{log}(Z)$, C/O ratio, and logarithmic cloud-top pressure $\mathrm{log}(P)$—and one hyper-parameter "err_multiple", which is applied to the spectroscopic error bars as a scaling factor, float as free parameters during the retrieval process. We apply a uniform prior on T_p and constrain it to between 550 and 950 K, which are estimated assuming a 0–0.67 albedo and zero to full global heat redistribution. The rain-out condensation process is taken into account, and since our wavelength coverage from 1.1 to 1.7 μm is not adequate to identify the Mie scattering shape at the short wavelength limit, we set the base-10 logarithm of the scattering factor to its default value of 0, assuming a uniform opaque cloud deck where the cloud-top pressure is a free parameter in the atmosphere. We test using both MCMC and nested sampling methods to explore the parameter space and find that the nested sampling method is better suited for this task to prevent samples from being trapped in local minima and for faster convergence.

Figure 16 shows the posterior distributions of free parameters in the retrieval, and we present the input fixed parameter, priors, and the output posterior distribution center and 1σ uncertainties in Table 8. Our fits constrain the cloud-top pressure to values greater than 0.01 bars and C/O ratios to above 0.8 (i.e., super-solar), but the posteriors for both parameters are effectively uniform within this preferred range. The planet radius posterior is consistent with our results from the WFC3, MOST, and Spitzer data, the log(Z) posterior peaks around $\mathrm{log}(Z)=2.4$, and the equilibrium temperature posterior favors temperature as high as 900 K. Aside from these, the posterior of err_multiple is constrained at ${1.47}_{-0.13}^{+0.13}$, indicating that the uncertainties in our HST/WFC3 transmission spectrum may be slightly larger than estimated, despite the conservative approach we have adopted when analyzing our uncertainty budget.

Zoom In Zoom Out Reset image size

Table 8.  Parameters for the PLATON Atmospheric Retrieval Using Our WFC3 Spectrum

Parameter Names	Median	Lower Error (1σ)	Upper Error (1σ)
			Fixed Parameters
Rs(R⊙)	0.746	⋯	⋯
Mp(M⊕)	7.81	⋯	⋯
Teff(K)	5192	⋯	⋯
$\mathrm{log}(\mathrm{scattering}\,\mathrm{factor})$	0.0	⋯	⋯
			Fit Prior
Rp(RJup)	${ \mathcal U }(0.19,0.22)$	⋯	⋯
Teq(K)	${ \mathcal U }(550,950)$	⋯	⋯
$\mathrm{log}Z$	${ \mathcal U }(-1.0,3.0)$	⋯	⋯
C/O	${ \mathcal U }(0.05,2.0)$	⋯	⋯
$\mathrm{log}P(\mathrm{bar})$	${ \mathcal U }(-8.99,2.99)$	⋯	⋯
			Fit Posteriors
Rp/Rs	0.0283	0.0004	0.0002
Teq(K)	809	142	103
$\mathrm{log}Z$	2.4	0.4	0.3
C/O	1.3	0.4	0.4
logP(bar)	0.15	1.81	1.89
error multiple	1.47	0.13	0.13

Download table as:  ASCIITypeset image

Additionally, we also combine all transit depth results from WFC3, MOST, and Spitzer to perform a full (0.5–4.5 μm) transmission spectrum analysis using PLATON. During the fitting, the observed mean transit depth of the WFC3 spectrum is adjusted while the relative spectral shape is fixed to account for the uncertainty in visit-long trend model selection when analyzing the WFC3 data (Section 2.3). The retrieved parameter posteriors are consistent with retrieval results using only the WFC3 data set within 1σ uncertainty. We plot the best-fit full-range retrieval model (T_eq = 770 K, $\mathrm{log}(Z)=2.7$, C/O = 1.0, and $\mathrm{log}P(\mathrm{Pa})=3.6$, with a WFC3 spectrum shift of −78 ppm) against our transmission spectrum data sets in Figure 17. And since adding three data points (MOST and Spitzer) in the spectrum does not change our atmospheric analysis result, we leave out the repeated process in the following sections and only present results using the WFC3 transmission spectrum.

Zoom In Zoom Out Reset image size

In addition to PLATON, we use an existing generic forward model generated with ATMO (Tremblin et al. 2015; Goyal et al. 2019) to further explore HD 97658b's atmospheric properties. We adopt this method based on PLATON's assumption that planetary atmospheres can be parameterized using only bulk metallicity enhancement, C/O ratio, and equilibrium chemistry. We first take a subset of the entire grid, which spans across four planet equilibrium temperatures from 600 to 900 K, two surface gravities 10 and 20 m s⁻², five atmospheric metallicities (1×–200× solar), four C/O ratios (0.35–1.0), and four uniform cloud parameters. The equilibrium temperature grid span is constrained by our equilibrium temperature prior, while the other three-dimensional spans are constrained by the original ATMO grid upper and lower limits. At each parameter setting, we also adjust the planet's photospheric radius and find its best-fit value by maximizing the log-likelihood. Then, we calculate the posterior distribution of each parameter by marginalizing maximum likelihoods over all remaining dimensions. In order for the comparison to be effective, we repeat the same evaluation process with PLATON's transmission spectra generation module on the same parameters grid, and the resulting posterior distribution comparisons of three critical atmospheric parameters are shown in Figure 18.

Zoom In Zoom Out Reset image size

A high consistency is achieved for T_p posterior distributions. For the distributions of $\mathrm{log}(Z)$ and C/O, the ATMO and PLATON models give similar shapes. ATMO produces a bigger tail in the low-metallicity end, but ATMO and PLATON both indicate high-metallicity (200× solar) peaks on the $\mathrm{log}(Z)$ plot; on the C/O ratio plot, the ATMO models favor C/O = 0.7 while the PLATON models favor C/O = 1.0 or higher.

In the previous retrieval sections, we showed that no clear peak value of the parameters is found even when a large volume of atmospheric parameter space is explored. We now also pick out a small subset of forward models to compare with our data.

As a baseline, we test similar scenarios that were considered in Knutson et al. (2014), which include a range of models with varying metallicity, cloud-top pressure, and C/O ratios, and a range of models composed of only H₂ and H₂O with varying H₂O number fractions. On top of previously proposed models, we also investigate the effect of adjusting the abundance of other molecular species in the atmosphere. Several typical atmospheric models that we investigated are presented in Figure 19, along with a perfectly flat model (i.e., a wholly featureless transmission spectrum). The equilibrium temperature is set to 900 K in all atmospheric models, and we adjust the mean transit depth of the model spectrum until we find the location with the maximum likelihood. Table 9 summarizes all models that have been compared with our transmission spectrum data, as well as the best model retrieved with PLATON (T_eq = 950 K, $\mathrm{log}(Z)=2.2$, C/O = 0.8, and $\mathrm{log}P(\mathrm{Pa})=7.8$). Their reduced χ² values and the confidence levels with which they can be ruled out are presented in the same table.

Zoom In Zoom Out Reset image size

Table 9.  Forward Modeling of HD 97658b Transmission Spectrum in the HST/WFC3 G141 Bandpass

Model	${\chi }_{\nu }^{2}$	Rule-out Confidence	${\chi }_{\nu }^{2}$ (Error Scaled Up)	Rule-out Confidence (Error Scaled Up)
Flat	2.5	5.4σ	1.2	0.6σ
Best model retrieved with PLATON	2.5	4.9σ	1.2	0.5σ
1× solar, C/O = 0.8, P = 1 bar	4.1	10.2σ	1.9	3.0σ
200× solar, C/O = 0.8, P = 1 bar	2.5	5.0σ	1.2	0.5σ
200× solar, C/O = 0.8, P = 10 mbar	2.6	5.3σ	1.2	0.7σ
200× solar, C/O = 0.8, P = 1 mbar	3.3	7.7σ	1.5	1.7σ
200× solar, C/O = 0.1, P = 1 bar	4.0	9.9σ	1.9	2.8σ
1000× solar, C/O = 0.8, P = 1 bar	2.7	5.8σ	1.2	0.8σ
1%CH4 + 4%CO2 + 10%CO	2.2	4.0σ	1.0	0.1σ
4%CO2 + 10%CO	3.3	8.1σ	1.5	1.8σ
100%H2O	2.4	5.2σ	1.1	0.4σ
50%H2O	2.6	5.6σ	1.2	0.7σ

Download table as:  ASCIITypeset image

The result shows that with this data set, even the best-fit atmospheric model is excluded at a 4σ significance level, indicating that the uncertainties associated with the transmission spectrum are underestimated, as pointed out in Section 5.1, despite the careful error budget treatments we implemented.

As mentioned in Section 5.1, the retrieval with PLATON shows that ${\chi }_{\nu }^{2}\approx 1.0$ can be achieved for best-fit models when the spectrum data uncertainties are scaled up by a factor of around 1.5. We apply this scaling factor to our data and recalculate the forward modeling ${\chi }_{\nu }^{2}$ and the number of σ, and present the results in the last two columns of Table 9. After applying the scale factor, we find that only models with low metallicity (1× solar) or low C/O ratios (C/O = 0.1) can be ruled out with around 3σ significance.

In light of the above analyses, we conclude that the existing atmospheric models can barely be distinguished with the current HST/WFC3 data set; although, the Z = 200× solar, C/O = 0.8, and P = 1 bar or 10 mbar models, which fall into the high posterior regions presented in Sections 5.1 and 5.2, and an H-dominated atmosphere consisting of around 1% of CH₄, a few percent of CO₂, and around 10% of CO are favored over other metallicity, C/O ratio, and cloud-top pressure settings or other atmosphere molecular combinations.

We propose that more observations of HD 97658b are needed in order to more tightly constrain its atmospheric features, and more discussions are provided in Section 6.

Rackham et al. (2018) proposed the transit light source effect, which describes the problem that spots and faculae of M dwarf stars can produce contamination in the transmission spectra of nearby planets more than 10 times larger than the transit depth changes expected from planet atmospheric features. In Rackham et al. (2019), the transit light source effect analysis was extended to F/G/K host stars. They found that while the stellar contamination can bias the transit depth measurement by as much as  1% for late G and K type dwarfs from UV to NIR, the offsets in transmission spectral features, including H₂O, CH₄, O₃, and CO amplitudes, induced by F/G/K stellar contamination assuming both spots only and spots+faculae models on stellar surface are far smaller than the atmospheric features expected in the transmission spectra of the planets around these stars. Therefore, the stellar contamination in the HST/WFC3 wavelength range and the Spitzer IRAC bandpass is not problematic for transmission spectroscopy analysis for typically active F/G/K dwarfs. As is described in Section 4.2, HD 97658b's activity cycle and rotation period are typical of early K type stars, and its overall activity level (${R}_{\mathrm{HK}}^{{\prime} }$ and S_HK) and 2 ms⁻¹ RV variations indicate that HD 97658b is slightly quieter than average K type stars. These facts ensure the reliability of the transmission spectroscopy results presented in this work. Nonetheless, Rackham et al. (2019) show that F/G/K stellar contamination may have larger impacts on transmission spectra at UV and visual wavelengths, in which TiO and VO display significant opacity. Therefore, we note that stellar contamination should be taken care of when analyzing the transmission spectrum of HD 97658b at UV and visual wavelengths in future works.

We assess the transmission spectrum detectability of all currently confirmed small planets with future missions by calculating the transmission spectroscopy metric (TSM) of each planet, as defined in Kempton et al. (2018). The TSM is defined to approximate an S/N of one scale height in transmission spectra when observed with the NIRISS instrument on JWST for 10 hr; therefore, it is proportional to ${R}_{{\rm{p}}}H/{R}_{* }^{2}\times F$, where H = kT_eq/μg is the atmospheric scale height, R_p is the planet radius, R_* is the stellar radius, and F is the stellar flux received on the detectors. For planets with R_p < 2 R_⊕, a mean molecular weight of μ = 18 is chosen, representing pure water atmospheres, and for planets with R_p ≥ 2 R_⊕, a mean molecular weight of μ = 2.3 is chosen, representing H/He-dominated atmospheres. The final formula of TSM is shown as Equation (1) of Kempton et al. (2018) Applying the scale factors specified for different R_p bins in Kempton et al. (2018), the TSM value represents a near-realistic S/N of 10 hr of JWST observations.

Following the above procedure, we calculate the TSM of all confirmed planets with 1 R_⊕ < R_p < 4 R_⊕ and equilibrium temperature T_eq < 1000 K, and we rank them according to their TSM values. The 20 highest-ranked planets are listed in Table 11, and the full table is available online. In Figure 20, we show the TSM value of each planet versus its radius, and each data point is color coded with the planet equilibrium temperature. The names of the top 20 planets with the highest TSM values are labeled. The shift in TSM values that we observe between planets smaller than 2 R_⊕ and larger than 2 R_⊕ is caused by an artificial sudden jump of the mean molecular weight values of these two groups of planets as described above, but it represents the actual TSM trend when we move from small planets to larger planets.

Zoom In Zoom Out Reset image size

Louie et al. (2018) simulates the transmission spectroscopy S/N of 18 known planets with sizes between 0.5 R_⊕ and 4.0 R_⊕ assuming they are observed with JWST/NIRISS. Twelve of the 18 planets are also in our sample. Their treatment of atmospheric mean molecular weight is similar to our method, except that they divide the exoplanets into two categories at 1.5 R_⊕ instead of 2.0 R_⊕, and as a result, there is a similar S/N shift at 1.5 R_⊕ in Figure 5 of Louie et al. (2018). Although we are measuring the detectability with different JWST instruments, and the definition of TSM in this work is slightly different from the simulated S/N in Louie et al. (2018), the resulting TSM and S/N from Louie et al. (2018) of the overlapping planets are highly consistent. Ten of the 12 overlapping planets are ordered the same in this work as in Louie et al. (2018), and the TSM values of all 12 planets are between half and two times of their S/N values in Louie et al. (2018), showing that our analysis is reliable.

Among all planets in our sample, GJ 1214b has the highest TSM value, but the transmission spectrum of GJ 1214b has been revealed to be featureless from 1.1 to 1.6 μm (Kreidberg et al. 2014a), indicating a cloudy atmosphere or no atmosphere. Nevertheless, the opacity of clouds may vary across a broader wavelength range, and emission and reflection spectra could also contain additional features. Therefore, GJ 1214b could still be a valuable target for JWST observations. Ranked second is K2-25b, a Neptune-sized planet orbiting an M4.5 dwarf in the 650 Myr old Hyades cluster (Mann et al. 2016; Thao et al. 2020). The unusually large size (3.43 R_⊕) of K2-25b in comparison to other planets with similar orbital periods (3.485 days) and the fact that it is orbiting a young star suggest that K2-25b could represent an early or intermittent phase of planetary evolution, where young Neptunes are in the process of losing their atmospheres (Mann et al. 2016). These features make K2-25b a target of high scientific value; although, we must consider that a young M dwarf like K2-25 could have large spot and faculae covering fractions, which may generate stellar contamination as large as several times the transit depth changes expected in the transmission spectrum of K2-25b due to its atmospheric features (Rackham et al. 2018). Ranked fifth is HD 97658b, the planet analyzed in this work. Since we are not able to precisely determine the atmospheric composition of HD 97658b with the current HST/WFC3 data as described in previous sections, we simulate JWST observations of HD 97658b using its Near InfraRed Spectrograph (NIRSpec) instrument in the following subsection and analyze quantitatively how well we can characterize the atmosphere of HD 97658b with JWST observations.

Out of all 1404 planets in Figure 20, 515 have TSM > 5.0 and 820 have TSM > 3.0. At least one-third of small planets cooler than 1000 K can be well characterized using JWST, and more valuable targets will be added to the pool with the ongoing TESS mission.

The upcoming JWST mission, with a 6.5 m diameter primary mirror and four near/mid-infrared instruments covering the wavelength range from 0.6 to 28.5 μm, will provide unprecedented opportunities to characterize exoplanets of all sizes and environments. Here, we simulate a transmission spectrum of HD 97658b as observed by JWST's NIRSpec with its G235M filter (1.6–3.2 μm) using the PANDEXO package (Batalha et al. 2017). We assume one transit observation with a 45% in-transit time and 90% saturation level. A reasonably optimistic noise floor of 30 ppm is adopted according to Greene et al. (2016). Figure 21 shows the simulated transmission spectrum observed by JWST assuming a 200× solar metallicity atmosphere with 1 bar cloud-top pressure and C/O = 0.8 and compares it with a flat spectrum and two other atmospheric models with C/O = 0.5 and 0.1, respectively. Our calculation shows that with only one transit observed by JWST, we will be able to distinguish a C/O = 0.8 atmosphere from a C/O = 0.7 atmosphere, C/O = 0.5 atmosphere, and C/O = 0.1 atmosphere with 5σ, 12σ, and 17σ significance, respectively, which will enable further study of the formation environment and formation process of HD 97658b. The same data will also be able to exclude the flat spectrum with 4σ significance.

Zoom In Zoom Out Reset image size

With this result, and considering the fact that HD 97658b is ranked as the fifth-best target for future transmission spectra observations from our previous calculations, we propose that HD 97658b be assigned top priority in JWST exoplanet proposals.