The Near-infrared Transmission Spectra of TRAPPIST-1 Planets b, c, d, e, f, and g and Stellar Contamination in Multi-epoch Transit Spectra

We derived wavelength solutions based on the position of the target in the direct images. We adopted up-to-date wavelength calibration coefficients from Wilkins et al. (2014). The centroids of the target point source in the direct images were determined by fitting two-dimensional Gaussian profiles. In Program 1, the direct-image frame had a different aperture from the spectroscopic frame. For those observations, we adjusted the direct-image coordinates accordingly.

We adopted a third-order polynomial function in wavelength for the flat-field correction. For each visit, we separately calculated a wavelength calibration–dependent flat-field correction, i.e., differences introduced primarily by the variations in the position of the target on the direct images. We next applied the correction to individual nondestructive reads. While applying the flat-field correction, we also identified and corrected for low data quality pixels: pixels in the flat-field frame that deviated more than 20% from unity were flagged (most of these also had nonzero data quality flags, i.e., were also flagged by the CalWFC3 pipeline). We also flagged any additional pixels identified by the CalWFC3 pipeline as bad pixels, hot pixels, pixels with unstable response, or bad or uncertain flat values.⁸ Finally, we replaced the flagged pixels by interpolating over the neighboring unflagged pixels.

Correcting for CRs was a crucial step in our reduction, particularly for the orbits heavily affected by SAA passage. We identified seven orbits that suffered from SAA passage using the fits header flag and the number of CR hits. On average, there were at least 20 visually apparent CR hits per frame in these orbits. We developed a suite of custom algorithms to identify and remove the CRs and evaluate the efficiency of the CR corrections.

First, we applied iterative bidirectional median filtering to each nondestructive read and identified pixels as CR-affected if they exceeded the median-filtered image level by a threshold of 11σ (determined through the CR removal assessment described below). For each iteration, pixels that were previously marked as CRs were excluded from median filter calculations. We typically repeated this iterative filtering procedure at least three times (see discussion below for the connection between the algorithm's performance and the number of iterations used).

For each identified CR, we replaced the CR-affected pixel value with the weighted average (f_replace) of the same pixel in the exposures preceding (f₋₁) and following (f₁) the image, as described by the relation

$$\begin{eqnarray}&&{f}_{\mathrm{replace}}=\displaystyle \frac{{t}_{0}-{t}_{-1}}{{t}_{1}-{t}_{-1}}{f}_{1}+\displaystyle \frac{{t}_{1}-{t}_{0}}{{t}_{1}-{t}_{-1}}{f}_{-1},\end{eqnarray} \tag{ 1 }$$

in which t₀ refers to the time of the CR-affected exposure and the subscripts ±1 denote the two adjacent (in time) exposures. The weights for the averaging are effectively the inverse of the time difference between the exposures. Program 1 adopted a bidirectional scanning mode; i.e., forward and reverse scanning directions were applied alternately. In this case, the preceding and following exposures had slightly different scanned image regions due to the upstream/downstream effect (McCullough & Mackenty 2012). To account for this effect, we corrected the CR hits separately for images taken with different scanning directions. Figure 1 shows example images before and after CR removal.

Zoom In Zoom Out Reset image size

We assessed our CR correction algorithm quantitatively by injecting and removing CR hits with a CR template. We constructed the CR template using one uncorrected frame (observation identifier: iddea1meq) that was taken during an SAA crossing and contaminated by over 1000 CR hits of different sizes and amplitudes. We removed the portion of the image with the stellar spectrum (pixel coordinates [135:170, 50:200]) and replaced it with randomly selected copies of regions out of the spectrum. We set all pixels with a flux below 2000 e⁻ (∼3× sky background) as zero and used the resulting image as the CR template. The CR template had 1150 CR hits, representing the most severe CR-affected case. We added the template to 20 cleaned frames and applied our CR removal algorithm. We evaluated the performance of our CR identification and correction algorithm based on three criteria: CR identification rate, false-positive rate, and correction efficiency (η_c). The latter we defined on a pixel-to-pixel basis as

$$\begin{eqnarray}&&{\eta }_{c}=1-\left|\displaystyle \frac{{f}_{\mathrm{corrected}}-{f}_{\mathrm{clean}}}{{f}_{\mathrm{dirty}}-{f}_{\mathrm{clean}}}\right|,\end{eqnarray} \tag{ 2 }$$

in which f_clean, f_dirty, and f_corrected refer to the pixel fluxes in the input, CR-injected, and output images, respectively.

Three parameters influence the effectiveness of the algorithm: the size of the median filter, the CR identification threshold, and the number of iterations. We optimized the algorithm through a three-dimensional grid search, with the size of the median filter ranging from 3 to 25 pixels, the thresholds ranging from 3σ to 20σ, and the number of iterations ranging from 1 to 10. The most effective combination included an 11 pixel median filter, an 11σ threshold, and a minimum of three iterations. Our algorithm yields identification rates of 98.5% on average, false-positive rates below 1.5%, and correction efficiencies above 90% for over 99.9% of the pixels across the entire image for non-SAA exposures. Within the image region containing the stellar spectrum, the identification rate dropped slightly to approximately 90%. Typically, there were less than 40 pixels identified as CR hits in the spectrum region in an exposure obtained out of SAA passage, indicating that the total numbers of missed CR pixels and false-positive pixels are 5 and 1 pixels, respectively. Since the region we used for the spectral extraction has a total size of 60 × 140 = 8400 pixels, we found that unidentified and false-positive CR hits had a negligible influence on our results.

In addition, SAA passage severely affected the pointing accuracy of HST for a few orbits in Program 2. For example, during Orbit 5 of Visit 4, the pointing shifted by ∼14 pixels (2''). Shifts with similar amplitudes were observed in the first four orbits of Visit 1 and the first two orbits of Visit 3 in Program 2, which all indicate failures in the HST fine guidance. In these orbits, the light curve also drops by ∼1%. We identify two potential causes for these systematics. First, after a pointing shift, the spectrum moved to the part of the detector that was not previously illuminated. This increased the ramp effect induced by charge trapping, because charge traps in these newly illuminated pixels have not yet been filled (Zhou et al. 2017). Second, WFC3's IR flat field has intrinsic uncertainties of ∼1.0% (Deustua et al. 2016), which introduce light-curve systematics for such large pointing shifts. The ramp-effect correction would, in principle, be correctable using the RECTE model (Zhou et al. 2017) with additional free parameters describing the image drifts. However, the existing science data and calibrations did not provide a viable option for alleviating the increased flat-field uncertainty. Therefore, we excluded the first four orbits in Visit 1, the first two orbits in Visit 3, and the last orbit in Visit 5 from the remainder of our analysis.

We consulted transit times listed in the NASA Exoplanet Archive (Akeson et al. 2013) and identified three transits that may lie within the discarded orbits. For completeness, we list these transits in Table 2.

We identified and removed the sky background using a σ-clipping algorithm. Pixels within 5σ of the median image level after 10σ clip iterations were considered as background. The median value of the background pixels was then subtracted from the image.

The remaining observations suffered from HST pointing drifts in both the x and y directions at levels of 0.05–0.1 pixels per orbit. Such drifts, especially in the wavelength dispersion direction, introduced systematic slopes in the spectrally binned light curves when left uncorrected. To correct for the drifts, we measured the shifts between each image and the reference image (first image in each data set) by cross correlation. We then used bicubic interpolation to shift and align the images to the reference image.

We generated white-light and spectrally binned light curves. We summed the CR-cleaned, background-subtracted, and aligned subexposures back to a total exposure image. We created 12 10 pixel wide bands from the scanned area ranging from 1.1 to 1.70 μm. We then obtained the light curves by summing pixels inside a 60 pixel wide window for every band. The uncertainty of each point on the light curve includes photon noise, dark current, and readout noise. As wavelength solutions differ slightly for each visit, the central wavelengths of the 10 pixel wide bins vary at levels of (∼0.01 μm) for the different visits. We list the central wavelengths of the bands in Table 3. In this way, five raw light curves were derived, each with 12 bins.

We note that the bins applied here are slightly different from those used by de Wit et al. (2016, 2018), who used 11 (0.05 μm wide) bins in the 1.15–1.7 μm range for Transits 1 and 2 and 10 (also 0.05 μm wide) bins in the 1.15–1.65 μm range for the subsequent transits. In general, determining bin sizes with an integer number of pixels is more widely adopted in HST/WFC3 transmission spectroscopic studies (Deming et al. 2013; Mandell et al. 2013; Kreidberg et al. 2014). Considering that the 0.05 μm bin size is not an integer multiple of the spectral resolution unit of the G141 grism, without knowing the exact binning and interpolations applied there, we could not use identical bins in the light-curve extraction steps. Therefore, to compare our results with those of de Wit et al. (2016, 2018), we instead interpolated our pixel-binned transmission spectra. The interpolation had a negligible effect due to the coarse wavelength resolution of the spectra. We discuss this point further in Section 4.3.

The raw light curves show prominent ramp-effect systematics (Figure 2), typical of HST/WFC3/IR time-resolved observations (e.g., Berta et al. 2012; Apai et al. 2013). This systematic is caused by two populations of charge carriers that are trapped and then, with some delay, released by impurities in the HgCdTe detectors (Zhou et al. 2017). We corrected these systematics using the RECTE model (Zhou et al. 2017), which models the history of illumination, trapping, and release for each pixel. This model offers a consistent solution to correct the ramp-effect systematics from the perspective of the physical cause, instead of fitting empirically determined exponential/polynomial functions, which are used in most HST/WFC3 transiting exoplanet studies to date. The use of this correction is also a major difference between our data reduction and that of de Wit et al. (2016), who used exponential functions to model and correct for the ramp effect. De Wit et al. (2018), which was published after the submission of this paper, used another model bearing more resemblance to the RECTE model to remove the ramp effects, and we compare them in Section 7.3.

Zoom In Zoom Out Reset image size

Zhou et al. (2017) described the charge-trapping processes with six parameters (E_s,tot, E_f,tot, η_s, η_f, τ_s, and τ_f⁹ ) representing the trap numbers, trapping efficiency, and charge release time for slow and fast charge-trap populations. Zhou et al. (2017) found these parameters to vary little in different observations and considered them to be intrinsic to the WFC3 detector. Therefore, we fixed these six parameters and provided the adopted values in Table 4. The free parameters that determine the systematic profiles are as follows.

• 1.  
  f: the incoming flux on each pixel as a function of time. We consider it to be a constant here.
• 2.  
  E_s,0 and E_f,0: the initial numbers of trapped charges.
• 3.  
  ΔE_s and ΔE_f: the number of additional charges trapped during interorbit gaps due to unintended detector illumination.
• 4.  
  v: the slope of the visit-long trend.

Table 5.  Parameters Used as Input in Transit Profile Fitting from Gillon et al. (2016, 2017)

Planet	Inclinationa	Tb	a
	(deg)	(days)	(R*)
b	89.65+0.22−0.27	1.5109	20.50
c	89.67 ± 0.17	2.4218	28
d	89.75 ± 0.16	4.04961	39.55
e	${89.86}_{-0.12}^{+0.10}$	6.099615	51.97
f	89.680 ± 0.034	9.20669	68.4
g	89.710 ± 0.025	12.35294	83.2

Notes.

^aGaussian-distributed priors. ^bFixed parameters.

Download table as:  ASCIITypeset image

We found the best-fit RECTE profiles for the light curves of each band. While the RECTE algorithm essentially models charge-trapping processes in individual pixels, it was not feasible (or necessary) to fit the light curves at the single-pixel level. First, the ramp effect at the pixel level is overwhelmed by other systematics, particularly telescope jitter. Second, the accuracy of the single-pixel level ramp-effect correction is negatively influenced by photon noise in these data. As Zhou et al. (2017) found no evidence for the charge-trap parameters varying between pixels, we therefore adopted an average band-level charge-trapping correction instead of a single-pixel level correction.

For observations using bidirectional scanning, exposures conducted in the opposite scanning directions bear an intrinsic flux level difference of ∼0.5%. For these cases, we assumed different f and v values for the light curves observed in different scanning directions but calculated the charge-trapping/release processes for the two scanning directions together.

We found the best-fit parameters using a Markov chain Monte Carlo (MCMC) with 500 walkers for 600 steps, with the first 300 steps discarded as burn-in. The MCMC runs were performed using the emcee package (Foreman-Mackey et al. 2013). Examples of the best-fit RECTE profiles are shown in Figure 2. In most bands, the ratio of the average value of the standard deviation to the average of the photon noise is within the range of 0.8–1.2, i.e., our complete procedure (including CR and charge-trapping corrections) robustly reaches the photon noise level or very near to it.

For each transit, we also derived a broadband light curve by computing the weighted average of all bands, adopting the inverse variance as the weight.

Our final reduction steps were fitting the transit profiles and extracting the transmission spectra. We first fitted the broadband light curves by generating model transit light curves using the Python package batman (Kreidberg 2015), in which the light-curve shape model is based on Mandel & Agol (2002). The fitting procedure was performed with an MCMC algorithm using the emcee package (Foreman-Mackey et al. 2013). The transit profile model contained nine parameters, namely transit mid-time t₀, orbital period T, relative planet size r_p/r_*, semimajor axis a, eccentricity e, inclination i, argument of periapsis ω, and quadratic limb-darkening coefficients (LDCs) u₁, u₂.

We found that due to the lack of ingress or egress data in some transits, the LDCs could not be constrained well by the light curves. Consequently, the LDCs obtained in different transits are not always consistent. It is important to note that models predict that LDCs of late M stars will vary significantly with wavelength in the 1.1–1.7 μm range and that the transit depths derived are anticorrelated with LDCs (Figure 3). Therefore, errors in LDCs may introduce apparent spectral features in the transmission spectra. To carefully examine the effect of LDCs, we experimented with three different LDC treatments: (i) interpolating LDC values and uncertainties provided in de Wit et al. (2016)—derived using PHOENIX stellar models (Husser et al. 2013)—to our bandpasses and using these as Gaussian-distributed priors; (ii) fixing the LDC to the best-fit values reported in de Wit et al. (2016); and (iii) independently deriving the LDCs by fitting PHOENIX specific intensity model stellar spectra¹⁰ (disk-integrated, multiplied by the HST/G141 bandpass, and normalized) to the HST/G141 out-of-transit spectrum and fitting a quadratic limb-darkening law to the limb-darkening profile of the best-fit model. We then fixed the LDCs in the transit fit using the derived values.

Zoom In Zoom Out Reset image size

The comparison of the results based on the different LDCs showed that the final spectra are only weakly affected by the adopted LDCs: in every band, the transit depths derived from the different methods agreed with each other to levels better than 1σ. We adopt the first limb-darkening treatment described above as our nominal procedure and present the results from this approach here.

The limited phase coverage of ingress and egress in some visits that complicated the LDC studies also precluded precise measurements of the transit durations. This uncertainty likewise hampered the precise determination of the orbital inclinations. Therefore, we adopted the distributions obtained by Gillon et al. (2017) as priors in our fitting procedure.

We fixed all eccentricities to be zero, as justified by the small values found by Gillon et al. (2017), which rendered ω irrelevant. We also fixed the orbital periods and semimajor axes to the values given in Gillon et al. (2017). We did not directly adopt external constraints on the mid-transit times, as transit timing variations (TTVs) are large in the system and not yet well understood. The fit parameters are summarized in Table 5. With r_p/r_*, t₀, i, and the LDCs as the only free parameters, we performed an MCMC search. We adopted 500 walkers and ran them for 2000 steps, the first 1000 of which were treated as burn-in. In the evaluation of the fit quality, we did not use the uncertainty estimates directly derived from the pipeline (which are dominated by photon noise), as these did not include the assessment of the residual systematic noise (even though these are found to be very small). Instead, for each light curve in each band of each visit, we opted to calculate the standard deviation of the data in the baseline (pre- and post-transit) and adopted this value as a uniform relative uncertainty applicable to all data points in the light curves.

We present a corner plot of the MCMC posterior distributions in Figure 3 and an example of the transit profile fit in Figure 4.

Zoom In Zoom Out Reset image size

After fitting the broadband transit light curves, we fitted the transits in the individual spectral bins in each transit, keeping the mid-transit times fixed to the values found in the broadband transit fits. To test the reliability of these fits, we carried out a Shapiro–Wilk test on the residuals of the fittings (Shapiro & Wilk 1965). Most of our fits passed the test with p-values exceeding 0.1. We inspected each of the few exceptions (p-values less than 0.1) visually and found that a few outlying data points, probably caused by stellar flares or other activity, were the reason that the transit models did not provide complete fits. We further investigate possible stellar activity in the light curves in Section 5.4.

Finally, we subtracted the best-fit broadband transit depth value (r_p/r_*) from each spectral bin's transit fit (to determine relative, spectrally dependent differences in the transit depths) to derive the transmission spectra of the seven transits.

We measured the WFC3 broadband transit depths of TRAPPIST-1 b, c, d, e, f, and g with an average precision of 123 ppm. We note that overlapping Transits 1 and 2 have larger uncertainties in transit depths than the other transits, which are all due to single planets. Most TRAPPIST-1 planets' transit depths we measured in the WFC3 G141 broadband are consistent with those measured in Spitzer Channel 2 (central wavelength 4.5 $\mu {\rm{m}}$) light curves (Gillon et al. 2017). However, our transit model's r_p/r_* for Transits 1, 2, and 3 (0.0849 ± 0.0012, 0.0879 ± 0.0012, and 0.0622 ± 0.0005) are deeper than the corresponding Spitzer transit depth measurement for the same planets (0.0828 ± 0.0006, 0.0852 ± 0.0005, and 0.0605 ± 0.0015) by over 1σ. Differences between the HST and Spitzer bands should come as no surprise, as these may be introduced either by planetary absorption features or—more likely—by stellar activity and heterogeneity, which we will explore in greater detail in Section 5.4.

Through our transit light-curve modeling, we obtained high-precision mid-transit time measurements, which were converted from Modified Julian Dates (MJD) to Barycentric Julian Dates (BJD_TDB) using the algorithms described in Eastman et al. (2010). The results are summarized in Table 6.

The mid-transit times have a typical uncertainty of 0.0002 days (or 17 s). The uncertainties are typically dominated by the limited phase coverage (due to HST's visibility windows), i.e., light curves missing either the ingress or egress. Different phase coverage of the individual transits causes the quality of the constraints on mid-transit times to vary. Better constraints were achieved for the transits of planets TRAPPIST-1 b and c because they have better-than-typical transit phase coverages, owing to their shorter transit durations; for these planets, the transit mid-times are constrained with uncertainties of only 0.0001 days (or 8 s).

In contrast, Transit 7—corresponding to a transit of TRAPPIST-1 e—lacks both the ingress and ingress, resulting in a greater uncertainty (0.002 day) on the mid-transit time.

We found substantial TTV signals (observed–predicted time differences) in the observed transits when comparing our results with those from the Online Exoplanet Archive, which are calculated assuming strict periodicity (Akeson et al. 2013). In particular, the best-fit mid-transit time for Transit 4 (TRAPPIST-1 g), occurring on 2016 December 9, deviates from that predicted from the best-fit Spitzer transit mid-time (Gillon et al. 2017) and the planet's orbital period by ∼30 minutes, as was noticed and discussed in Wang et al. (2017).

We obtained WFC3/IR G141 transmission spectra for TRAPPIST-1 b, c, d, e, f, and g. The spectra are shown in Figure 5, with the individual spectra offset by arbitrary levels for clarity. Each spectra has 12 bins covering wavelengths from 1.1 to 1.17 μm, corresponding to a spectral resolution of Δλ = 50 Å (or spectral resolving power of R = Δλ/λ = 22–34; Figure 4). Table 13 lists original spectra with no interpolation, measured by transit depths of different bands, and Table 14 lists the spectra with wavelength interpolation.

Slight differences in the wavelength calibrations of the individual visits resulted in small differences of wavelength solutions among them (Table 3). To directly compare our results with those by de Wit et al. (2016) and to combine spectra from multiple visits, we interpolated our single-transit transmission spectra to align our wavelength bins with those of de Wit et al. (2016). We adopted the third-order spline interpolation for both transit depths and uncertainties. We justify the interpolation of transit depths in two respects. First, the adopted bin size is coarse and significantly larger than the average difference of the interpolated wavelengths and original wavelengths. Second, the transmission spectra and their uncertainties have small spectral variations. The interpolation introduced negligible modifications to the transit depths and uncertainties. In order to allow direct comparison with the spectra published in de Wit et al. (2016), we interpolated our transmission spectra to match the wavelength bins used in that study (Figures 5–7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following, we also explore the combined spectra of the TRAPPIST-1 planets, as well as the combined spectra of different subsets. We combined the spectra by summing the transit depths from the six planets in each of the bands (already aligned to those in de Wit et al. 2016), giving the combined transit depth

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{{\rm{p}}}^{2}}{{R}_{{\rm{s}}}^{2}}\right)}_{{\rm{c}}{\rm{o}}{\rm{m}}{\rm{b}}{\rm{i}}{\rm{n}}{\rm{e}}{\rm{d}}}=\displaystyle \sum _{i}\displaystyle \frac{{R}_{{\rm{p}},i}^{2}}{{R}_{{\rm{s}}}^{2}}\end{eqnarray} \tag{ 3 }$$

and uncertainty

$$\begin{eqnarray}&&\delta {\left(\displaystyle \frac{{R}_{{\rm{p}}}^{2}}{{R}_{{\rm{s}}}^{2}}\right)}_{\mathrm{combined}}=\sqrt{\displaystyle \sum _{i}{[\delta ({R}_{{\rm{p}},i}^{2}/{R}_{{\rm{s}}}^{2})]}^{2}.}\end{eqnarray} \tag{ 4 }$$

For planet e, for which we have two transits, we used the inverse variance–weighted average of the two as the transmission spectrum.

Combined spectra may provide information about shared spectral features: combining seven spectra could increase the signal-to-noise of the spectra by up to a factor of 2.6.

The combined spectra of all seven transits are shown in the top panel of Figure 6.

Given the differing quality of the spectra, we also explored combinations of different subsets of the transits. As noted above, the transits of planets b and c overlapped in time, and their spectra show an apparent anticorrelation. Additionally, the parameters of Transit 7 (planet e) are determined less precisely than those of the other transits due to its lack of ingress or egress coverage. To allow the assessment of the impacts of these data quality differences and possible systematics on the combined spectra, we show two additional spectral combinations in Figure 6. The middle panel shows the combined spectra of planets d, e, f, and g (with spectra for planets b, c, and e from Transits 1, 2, and 7 excluded). The lower panel in Figure 6 shows the combined spectrum of Transits 1 and 2, i.e., planets c and b, respectively.

Water absorption is the most prominent expected spectral feature in the planetary atmospheres in this wavelength range. As such, we show in Figure 6 a model of water absorption for comparison (gray lines). For this comparison, we adopted a water transmission model calculated using ExoTransmit (Kempton et al. 2017), with an amplitude scaled to provide the best fit to the observed spectra. None of the three combined spectra in Figure 6 resemble the water absorption spectrum: while the water absorption would result in larger r_p/r_* values around 1.4 μm, all three spectra suggests a decrease in r_p/r_* values. We conclude that no water absorption or other planetary molecular absorption features are visible in any of the individual or combined spectra. In Sections 5.3 and 6, we explore the upper limits we can place on planetary water absorption, as well as the impact of stellar contamination on these spectra.

One of the key differences between our study and that of de Wit et al. (2016) is that we used the RECTE model to correct for the HST/WFC3 IR ramp effect, while the de Wit et al. (2016) study discarded the first orbits of each visit and relied on an empirical fit to correct the systematics in the subsequent orbits, implicitly assuming those to be identical to the ramp seen in Orbit 2. Here we compare the results of the two data reduction approaches.

First, by using RECTE, we successfully corrected the ramp effect in each visit's first orbit, which was discarded in de Wit et al. (2016), as well as in almost all published HST/WFC3 transit spectroscopic observations. As a result, we increased the available useful data and effectively improved the efficiency of the HST observations by about 25%. This increase translated to a better orbital phase coverage and an improved accuracy of the transit baseline levels. Furthermore, the additional baseline observations also enabled a more thorough exploration of the stellar activity and spectral changes (see Section 5.4) than would have been possible had those orbits been discarded.

Second, in addition to the increase in the data quantity and efficiency, we seek to compare the resulting data quality. At this point, any such comparison must be limited to Visit 1 (from Program 1), the only visit for which reduced data had been published at the time of submission for the current work.¹¹ With an uncertain astrophysical signal underlying possible residual systematics, such comparisons are not trivial. We proceed here by assuming that, to first order, the astrophysical signal is well understood and purely consists of the planetary transit that follows the analytic models of Mandel & Agol (2002); we will then discuss the limitations of this approach.

Under the assumptions laid out, the residuals of the observed and modeled light curves contain no systematics and should be photon noise–limited. Therefore, we use the standard deviation of the light-curve fit residuals as a metric to compare the data quality between our reduction and that of de Wit et al. (2016). De Wit et al. kindly provided their detailed results, enabling accurate comparison studies. We note that de Wit et al. (2016) reported a larger broadband residual standard deviation (240 ppm) than that calculated from data in Figure 1 of de Wit et al. (2016; 215 ppm). We conservatively adopt the latter for the comparison. The first three rows of Table 7 compare the standard deviations of residuals between the two studies for Orbits 1, 2, and 4, as well as several orbit combinations.

Table 8.  Planetary Parameters for Water Absorption Calculations

Planet	Mass	Teff	Δd	Radius	AH	${ \mathcal U }({A}_{H})$
	M⊕	K	ppm	R⊕
b	0.79 ± 0.27	400.1 ± 7.7	−413 ± 266	1.086 ± 0.035	−9.2 ± 6.8	11.2
c	1.63 ± 0.63	341.9 ± 6.6	−212 ± 245	1.056 ± 0.035	−11.1 ± 13.6	29.7
d	0.33 ± 0.15	288.0 ± 5.6	17 ± 113	0.772 ± 0.030	0.2 ± 1.0	3.2
e	0.24 + 0.56−0.24	251.3 ± 4.9	−209 ± 129	0.918 ± 0.039	−1.9 ± 4.6	11.9
e	${0.24}_{-0.24}^{+0.56}$	251.3 ± 4.9	−303 ± 126	0.918 ± 0.039	−2.8 ± 6.6	17
f	0.36 ± 0.12	219.0 ± 4.2	−130 ± 190	1.045 ± 0.038	−2.3 ± 3.5	8.2
g	0.566 ± 0.0038	198.6 ± 3.8	77 ± 197	1.127 ± 0.041	2.5 ± 6.6	22.3

Download table as:  ASCIITypeset image

Several points are notable about this comparison. If our underlying assumptions were correct, the residuals in all orbits should be the same. However, differences in Orbits 1, 2, and 4 are visible even within the same reductions. We attribute these differences in part to the scatter of the standard deviations themselves and in part to the fact that physical processes other than the transit are present in the data, as discussed below.

The immediate comparison of the standard deviations of the residuals between the two studies (rows 1–3 of Table 7) have overall very similar levels, considering the limited number of data points from which the standard deviations are calculated. Data from our reduction have lower standard deviations in Orbit 4 than those from de Wit et al. (2016; 136 ppm versus 159 ppm) but larger residuals in Orbit 2 (342 ppm versus 257 ppm). For the combined residuals of Orbits 2 and 4 (row 5), our results still have a slightly larger standard deviation than that of de Wit et al. (2016; 266 ± 25 ppm versus 215 ppm), corresponding to about a 1σ level difference (when assuming similar uncertainties for the standard deviations themselves). We note that Orbit 2's correction quality is worse in both reductions than that of Orbit 4, suggesting that the inherent data quality or other processes also play an important role.

Thus, based on a superficial comparison, one may conclude that (1) both reductions reach similar precision, (2) underlying differences in the data quality are more important than the type of ramp-effect correction applied, and (3) the empirical ramp-effect correction performs sometimes slightly better than the physical model for the charge trapping (not considering the benefits of the nondiscarded first orbit).

Closer inspection of the residuals in Figure 9 reveals isolated groups of outlying data points. These may mark potential small stellar flare events (discussed further in Section 5.4), such as an event with 4σ outliers above the baseline in Orbit 2 of Visit 0. It is instructional to inspect this event: in fact, the outliers contribute the most to the increased standard deviation in the residuals in this orbit. By excluding the four data points around suspected flare events, the standard deviations of the combined Orbit 2 and 4 data decreased by 40% to 195 ppm (last row of Table 7). If the same data points are excluded from the de Wit et al. (2016) data, the resulting standard deviation is, within the uncertainties, the same as that resulting from the RECTE fit. We conclude that the empirical ramp correction may occasionally provide a lower standard deviation because it can potentially fit the ramp-effect correction and stellar flares together.

Zoom In Zoom Out Reset image size

Thus, in our RECTE-based reduction, the higher standard deviations are an indication that the model applied (charge-trapping model + transit) is incomplete; in contrast, the empirical systematics fit and transit have the capability to absorb different sources of astrophysical signal and instrument systematics without distinguishing these, resulting in a slightly lower standard deviation of the residuals. Compared to empirical corrections, the physically motivated RECTE-based correction is less likely to be skewed by astrophysical events. This comparison highlights another advantage of the RECTE model over traditional empirical fits: given the well-determined detector response, it will be less likely to overcorrect and remove astrophysical processes (or other types of instrumental systematics). In fact, in Orbit 4, which has no signs of stellar flares, our RECTE model indeed resulted in a residual standard deviation that is the lowest of all measured here (136 ppm).

The transmission spectra from our reduction and that of de Wit et al. (2016) agree within the uncertainties (Figure 7). For planets b and c, we performed Kolmogorov–Smirnov (K-S) tests on the spectra from this study and de Wit et al. (2016). The tests resulted in p-values of 0.9094 and 0.8286; i.e., they did not show any evidence for the data points being drawn from statistically different parent samples.

The fact that our data sets contain two transit events (Transits 5 and 7) for planet TRAPPIST-1 e offers an opportunity to examine the similarities between the transmission spectra from the two epochs. Figure 8 shows the two spectra of TRAPPIST-1 e interpolated to the same wavelength bins. While the two spectra share some features, such as average transit depth and nondetection of planetary water absorption, they differ in the blue parts (<1.35 μm) of the spectra. There appears to be a difference in the overall slope between the two transit spectra.

To quantify the agreement between the two spectra, we calculated the bin-by-bin difference between them and report an average difference of 2.4σ. A K-S test yields a p-value of 0.01, i.e., supporting the conclusion that the two data sets are drawn from different parent populations. Manual inspection confirms that the blue parts of the spectra contribute most of the difference (Figure 8). This difference may be an indication of a time-evolving stellar contamination signal, which we discuss further in Sections 6.3 and 7.2.

In our data, all six planets (TRAPPIST-1 b–g) have spectra that show no obvious absorption features in the 1.1–1.65 μm wavelength range. We examine the possibility of planets sharing similar spectral features by combining the spectra from the six planets.

We also compared our observed transmission combined spectra with a scaled water transmission model (Kempton et al. 2017), as shown in Figure 6. The deviations, as expressed in χ², are 2.2, 1.3, and 1.0, which are all worse than those of the flat model (0.85, 0.56, and 0.45, respectively). We attribute this to the contamination from stellar heterogeneity, which we will discuss in more detail in Section 6.

Water is of the most interest among molecules with absorption features in the WFC3 G141 bandpass. Since no spectral features were observed, we estimate the upper limits of water absorption for each of the planets instead. Following Fu et al. (2017), we quantified the water absorption amplitude A_H, i.e., the transit depth difference in and out of the water absorption band in terms of scale heights H, as

$$\begin{eqnarray}&&{A}_{H}={\rm{\Delta }}d\,{R}_{* }^{2}/2{R}_{{\rm{p}}}H.\end{eqnarray} \tag{ 5 }$$

We estimated the scale height of each planet based on TTV masses from Wang et al. (2017). We set the host star radius to R_* = 0.114 ± 0.006 R_☉ (Gillon et al. 2016) to calculate the absolute radii of the planets and the host star effective temperature to 2559 K (Gillon et al. 2017). For the mean molecular weight, we adopted 18 amu, corresponding to that of water. We calculated ${ \mathcal U }({A}_{H})$, the 3σ upper bounds of A_H, for planets TRAPPIST-1 b–g.

The results, along with essential information, are presented in Table 8. We do not find evidence for a planetary water absorption feature in any of the transmission spectra; i.e., A_H is consistent with or less than 0 at the 1σ level in each case. However, the wide uncertainties on A_H allow for the possibility of significant water absorption, as illustrated by the upper limits. Improved precisions on the transmission spectra may place tighter constraints on planetary absorption features, though additional observations will have to contend with the time-resolved activity and photospheric heterogeneity of the host star.

We examined TRAPPIST-1's out-of-transit broadband light curves and searched for possible flares and dimming events. We identified the baseline using the scipy UnivariateSpline routine, which smooths the light curves with third-order splines. For the spline fit, the numbers of knots were optimized with inverse standard deviation weighting. We then computed the weighted average of deviations from unity within each window. Windows that had deviations above 3σ were selected as candidates for flare or dimming events. In total, two potential flares were found. We summarize the key properties of these potential events in Table 9 and show their broadband light curves in Figure 9. This figure also shows the spectrum of each event.

The possible events were sampled in only a few readouts, and their signal-to-noise ratio remains low, inhibiting their robust classification. We only note here a naïve expectation for these spectra, which is that (micro-)flares, representing hotter-than-photospheric plasma, would show an overall "bluer" continuum. This expectation may be met by these events, but the generally low quality of the spectra does not allow for meaningful characterizations.

In total, 20 orbits are used in our analysis, amounting to approximately 850 minutes of observations. We estimate the stellar occurrence rate of marginally detectable (micro)flares to be on the order of 1/425 minutes, i.e., an event every seven hours, on average. We note here that the identification of the events leading to this statistic benefited from the use of the RECTE correction, which allowed us to include data from four additional orbits in our analysis and provided a ramp-effect correction that was not affected by the flare events.

With this framework, we investigated the possible contribution of photospheric heterogeneities to the observed transmission spectra of the TRAPPIST-1 planets using the composite photosphere and atmospheric transmission (CPAT) model (Rackham et al. 2017). We used an MCMC approach developed with the PyMC (Patil et al. 2010) Python package to fit the CPAT model to our observations. The free parameters of the model and their priors are given in Table 10. For the prior on the photosphere temperature T_phot, we adopted the stellar effective temperature T_eff from Gillon et al. (2017), with a width of 5× of the reported uncertainty to allow the algorithm to thoroughly explore the parameter space. We adopted the same uncertainty for priors on the spot and facula temperatures, T_spot and T_fac, with means given by T_spot = 0.86 × T_phot and T_fac = T_phot + 100 K, following Rackham et al. (2018). For priors on the spot and faculae covering fractions, both within the transit chord and for the whole disk, we adopted normalized estimates of the covering fractions found for TRAPPIST-1 by Rackham et al. (2018).

The CPAT model describes the emergent disk-integrated spectrum of the photosphere as the sum of three distinct components—the immaculate photosphere, spots, and faculae—each covering some fraction of the projected stellar disk. We utilized the grid of DRIFT-PHOENIX model stellar spectra (Hauschildt & Baron 1999; Woitke & Helling 2003, 2004; Helling & Woitke 2006; Helling et al. 2008a, 2008b; Witte et al. 2009, 2011) with solar metallicity ([Fe/H] = 0.0) to generate spectra for each component. We parameterized the three components by their temperatures (T_phot, T_spot, and T_fac for the photosphere, spots, and faculae, respectively), which we allowed to vary, and linearly interpolated between models with different temperatures to produce the component spectra. For all components, we linearly interpolated between models with log g = 5.0 and $\mathrm{log}g=5.5$ to produce spectra matching the surface gravity of TRAPPIST-1 ($\mathrm{log}g=5.21$), which we calculated from the star's mass and radius (Gillon et al. 2017).

While fitting the transmission spectra, we required that the stellar parameters also produced a disk-integrated stellar spectrum matching the median observed out-of-transit stellar spectrum (λ = 1.15–1.70 μm) of TRAPPIST-1. We computed the disk-integrated spectrum as

$$\begin{eqnarray}&&{S}_{\lambda ,\mathrm{disk}}=(1-{F}_{\mathrm{spot}}-{F}_{\mathrm{fac}}){S}_{\lambda ,\mathrm{phot}}+{F}_{\mathrm{spot}}{S}_{\lambda ,\mathrm{spot}}+{F}_{\mathrm{fac}}{S}_{\lambda ,\mathrm{fac}},\end{eqnarray} \tag{ 8 }$$

ignoring projection effects owing to the positions of photospheric heterogeneities, which are not constrained by our model. Both the observed and model stellar spectra were normalized to the median flux between 1.27 and 1.31 μm for comparison. To account for discrepancies between the high-precision HST observations and stellar models, we multiplied the observational uncertainties by an error inflation factor η, which was allowed to vary between 1 and 100 (see Table 10).

We utilized two variations of this modeling framework. In the first, which we call the flat model here, we set the active-region covering fractions within the transit chord to the whole-disk values (i.e., f_spot = F_spot and f_fac = F_fac). Thus, _λ = 1 (see Equation (7)), and there was no stellar contribution to the observed transmission spectra. The achromatic transit depth D solely determined the transmission spectra model, while the stellar parameters in the model only affected the fit to the out-of-transit stellar spectrum. In the second framework, which we call the contamination model, we allowed the active-region covering fractions within the transit chord to differ from the whole-disk values. In this case, both planetary and stellar parameters affected the model fit to the observed transmission spectra (Equation (6)), and the stellar parameters (save the transit-chord covering fractions) determined the fit to the stellar spectrum. Flat models have seven free parameters, contamination models have nine, and each fit includes 144 data points—14 for the transmission spectrum and 130 for the stellar spectrum (see Section 6.2).

For each transmission spectrum that we considered, we performed an MCMC optimization procedure for both of these modeling frameworks. In each procedure, we marginalized over the log likelihood of a multivariate Gaussian. We ran three chains of 5 × 10⁵ steps with an additional 5 × 10⁴ steps discarded as the burn-in. We checked for convergence using the Gelman–Rubin statistic $\hat{R}$ (Gelman & Rubin 1992) and considered chains to be well mixed if $\hat{R}\lt 1.03$.

We initially performed this analysis using only the HST transit depths presented here and the Spitzer 4.5 μm transit depths provided by Delrez et al. (2018b). After the submission of this manuscript, Ducrot et al. (2018) presented K2 (0.42–0.9 μm) transit depths and I + z (0.8–1.1 μm) transit depths from the SPECULOOS-South Observatory (SSO; Delrez et al. 2018a) for each of the planets discussed here. In the following analysis, we consider the full K2+SSO+HST +Spitzer transmission spectra. Thus, each transmission spectrum includes 14 transit depths: 11 HST/WFC3 depths from this analysis, along with a K2, SSO, and Spitzer 4.5 μm depth.

We note that Ducrot et al. (2018) examined the impact of stellar contamination from TRAPPIST-1 by comparing the K2+SOO+HST+Spitzer spectra of TRAPPIST-1 b and c to a model from a pre-peer-reviewed version of this work (see Figure 4 in Ducrot et al. 2018). They found a 20σ discrepancy between the observed K2 transit depth and the model prediction and claimed that the stellar contamination model can be "firmly discarded." For the model that was used in the study of Ducrot et al., we assumed that spots and faculae were present only in the nonocculted stellar disk and not in the transit chords; i.e., the spectra of the transit chords were the same as the immaculate photosphere spectrum. During the reviewing process and prior to the publication of Ducrot et al. (2018), we updated the model to a more general form in which spots and faculae affect both the nonocculted disk and the transit chords (see Equation (7)). The analysis of the K2+SOO+HST+Spitzer data set using the more general stellar contamination model is presented here. Appendix C provides the results of this same analysis considering only the HST+Spitzer transmission spectra, along with the accompanying predictions for the K2 and I+z transit depths. The results of the analyses for the two data sets are fully consistent. Most notably, the best-fit contamination models from the HST+Spitzer analysis offer accurate predictions of the K2 and I+z transit depths for the combined transmission spectra.

We also note that transit depth determinations depend on i and a/R_s, which may vary between analyses of transit data. In our analysis, we fixed these values, which are hard to constrain with the HST observations, to those reported by Gillon et al. (2017). As the K2, SSO, and Spitzer observations all covered a longer baseline, however, Delrez et al. (2018b) and Ducrot et al. (2018) both benefited from larger data sets including many repeated transits, which allowed them to fit for more parameters. To be more specific, both these works adopted impact parameter as a free parameter in MCMC analyses. They also included the stellar mass and radius as parameters, either fixed or with a prior distribution, which combine with the orbital periods reported by Gillon et al. (2017) to yield the system scale a/R_s for each planet. For a comparison of the orbit parameters, we refer the reader to Table 5, which contains the parameters that we adopted from Gillon et al. (2017), and Table 1 of Delrez et al. (2018b), which details the fitted parameter values from that analysis.¹² In short, we find that the inclinations and system scales that we adopt differ, on average, from those of Delrez et al. (2018b) by 0.9σ and 0.2σ, respectively. We note that these subtle differences might lead to minor differences in the following analysis.

We first fit CPAT models to each of the seven single-transit transmission spectra presented in this work, as well as the weighted mean spectrum of TRAPPIST-1 e from Transits 5 and 7. We included the corresponding K2, SSO (Ducrot et al. 2018), and Spitzer 4.5 μm transit depths (Delrez et al. 2018b) in the fitting procedure. Table 11 summarizes the posterior distributions of the fitted parameters and provides the Akaike information criterion (Akaike 1974) corrected for small sample sizes (AIC_c; Sugiura 1978), Bayesian information criterion (BIC; Schwarz 1978), χ², and its corresponding p-values for each fit. We use the information criteria (ICs) to evaluate the efficacy of the increased model complexity of the contamination model compared to the flat model. Following convention (e.g., Liddle 2007), for both ICs—i.e., AIC_c and BIC—we interpret ΔIC > + 5 and ΔIC > + 10 relative to the best model as "strong" and "decisive" evidence against the current model, respectively.

Table 11.  Results of Stellar Contamination Model Fits to Single-planet Spectra

Data Set	Model	Fitted Parameter									AICc	ΔAICc	BIC	ΔBIC	χ2	p
		D (%)	Tphot (K)	Tspot (K)	Tfac (K)	Fspot	Ffac	fspot	ffac	η
b	flat	${0.744}_{-0.006}^{+0.005}$	${2118}_{-127}^{+87}$	${1962}_{-131}^{+111}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−736.0	6.2	−716.0	0.9	150.0	0.21
b	cont.	${0.677}_{-0.023}^{+0.022}$	${2433}_{-245}^{+240}$	${2009}_{-99}^{+150}$	${2950}_{-25}^{+50}$	${39}_{-10}^{+10}$	${49}_{-10}^{+7}$	${9}_{-9}^{+5}$	${48}_{-7}^{+8}$	${23}_{-2}^{+1}$	−742.2	⋯	−716.9	⋯	136.4	0.45
c	flat	${0.705}_{-0.005}^{+0.005}$	${2117}_{-125}^{+87}$	${1961}_{-136}^{+106}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-1}^{+1}$	−746.0	1.9	−726.1	−3.5	141.8	0.37
c	cont.	${0.663}_{-0.020}^{+0.021}$	${2271}_{-205}^{+147}$	${1960}_{-165}^{+91}$	${2964}_{-20}^{+36}$	${32}_{-10}^{+10}$	${47}_{-7}^{+6}$	${8}_{-8}^{+4}$	${47}_{-7}^{+6}$	${23}_{-2}^{+1}$	−747.9	⋯	−722.6	⋯	133.4	0.52
d	flat	${0.388}_{-0.003}^{+0.003}$	${2118}_{-127}^{+86}$	${1962}_{-133}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−723.2	24.6	−703.2	19.3	180.9	0.01
d	cont.	${0.309}_{-0.016}^{+0.015}$	${2551}_{-157}^{+253}$	${2000}_{-78}^{+107}$	${2937}_{-27}^{+63}$	${48}_{-9}^{+11}$	${49}_{-12}^{+9}$	${8}_{-8}^{+4}$	${52}_{-10}^{+12}$	${23}_{-2}^{+1}$	−747.8	⋯	−722.5	⋯	151.4	0.16
e (T5)a	flat	${0.477}_{-0.004}^{+0.004}$	${2117}_{-126}^{+86}$	${1961}_{-131}^{+110}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−759.8	−3.2	−739.9	−8.7	141.4	0.38
e (T5)a	cont.	${0.498}_{-0.020}^{+0.021}$	${2125}_{-120}^{+84}$	${1982}_{-124}^{+118}$	${2972}_{-15}^{+28}$	${16}_{-12}^{+8}$	${47}_{-6}^{+5}$	${12}_{-12}^{+5}$	${43}_{-6}^{+5}$	${23}_{-2}^{+1}$	−756.6	⋯	−731.2	⋯	139.5	0.38
e (T7)b	flat	${0.504}_{-0.004}^{+0.004}$	${2116}_{-125}^{+87}$	${1960}_{-135}^{+107}$	${2975}_{-14}^{+25}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−743.3	12.6	−723.4	7.1	158.9	0.10
e (T7)b	cont.	${0.449}_{-0.021}^{+0.021}$	${2559}_{-118}^{+163}$	${2032}_{-75}^{+100}$	${2937}_{-27}^{+63}$	${46}_{-9}^{+9}$	${51}_{-12}^{+8}$	${10}_{-10}^{+4}$	${38}_{-8}^{+11}$	${23}_{-2}^{+1}$	−755.9	⋯	−730.5	⋯	139.7	0.37
ec	flat	${0.493}_{-0.003}^{+0.003}$	${2116}_{-127}^{+85}$	${1960}_{-133}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−767.5	−5.6	−747.5	−10.9	142.0	0.37
ec	cont.	${0.480}_{-0.023}^{+0.022}$	${2267}_{-209}^{+149}$	${1981}_{-172}^{+102}$	${2965}_{-20}^{+35}$	${30}_{-11}^{+10}$	${48}_{-8}^{+6}$	${9}_{-9}^{+4}$	${43}_{-6}^{+6}$	${23}_{-2}^{+1}$	−761.9	⋯	−736.6	⋯	140.6	0.35
f	flat	${0.641}_{-0.005}^{+0.005}$	${2117}_{-125}^{+87}$	${1961}_{-130}^{+111}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−746.7	0	−726.7	−5.3	144.0	0.33
f	cont.	${0.627}_{-0.027}^{+0.026}$	${2395}_{-254}^{+199}$	${2022}_{-106}^{+159}$	${2959}_{-23}^{+41}$	${35}_{-11}^{+11}$	${50}_{-9}^{+6}$	${10}_{-10}^{+5}$	${40}_{-7}^{+8}$	${23}_{-2}^{+1}$	−746.7	⋯	−721.4	⋯	135.8	0.46
g	flat	${0.774}_{-0.006}^{+0.006}$	${2116}_{-126}^{+87}$	${1960}_{-131}^{+109}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−750.4	-7.8	−730.4	−13.1	137.3	0.48
g	cont.	${0.755}_{-0.027}^{+0.027}$	${2176}_{-173}^{+109}$	${1987}_{-133}^{+122}$	${2969}_{-17}^{+31}$	${23}_{-13}^{+12}$	${48}_{-6}^{+5}$	${10}_{-10}^{+4}$	${47}_{-7}^{+6}$	${23}_{-2}^{+1}$	−742.6	⋯	−717.3	⋯	139.1	0.39

Notes. Posterior distributions are provided for 16 MCMC optimization procedures resulting from two model frameworks and eight data sets. The model frameworks are the flat model, in which the transmission spectrum is unaffected by photospheric features, and the contamination model (identified as cont. here), in which the covering fractions of spots and faculae are allowed to differ from the whole-disk covering fractions. The eight data sets are the seven individual transits and the combined TRAPPIST-1 e data set. Medians and 68% confidence intervals of the fitted parameters are quoted. The AIC_c corrected for small sample sizes, BIC, χ², and corresponding p-value for each model are provided. Each model has 144 data points (14 for the transmission spectrum and 130 for the stellar spectrum). Flat models have 137 degrees of freedom, and contamination models have 135.

^aTransit 5. ^bTransit 7. ^cWeighted mean of Transits 5 and 7.

Download table as:  ASCIITypeset image

Of the single-planet spectra, both ICs only indicate decisive evidence against the flat model for one case: the TRAPPIST-1 d data set. Considering its χ² value, the flat model for this data set is ruled out at 99% confidence (p = 0.01). For TRAPPIST-1 e, of the two transits that were observed in two visits separated by 12 days, one shows evidence to support the contamination model (Transit 7), but the other (Transit 5) does not. This difference result for two transits of the same planet could suggest temporal variability of stellar contamination. For the remaining data sets, both ICs generally support the same model—the exception being the TRAPPIST-1 c data set—though they do not both rise to the level of decisive evidence. In general, the additional complexity of the contamination models results in lower χ² values, as expected, though the ICs show that the additional complexity is not decisively warranted by the data for any of the data sets besides that of TRAPPIST-1 d.

The observed effect of a stellar contamination signal scales with the transit depth (Equation (6)). The impact of stellar contamination is therefore more readily observable in the spectra of exoplanets with deeper transit depths. In the case of TRAPPIST-1, assuming that a steady-state stellar contamination signal similarly affects all the individual transmission spectra, we can coadd the individual transmission spectra to increase the signal-to-noise ratio. An examination of the combined transit spectrum can reveal whether the regions probed by the transit chords have different spectra from the average spectrum of the stellar disk.

Thus, in addition to the single-planet spectra, we fit both model frameworks to seven combinations of TRAPPIST-1 transmission spectra. The first combination is the sum of transit depths for all TRAPPIST-1 planets observed with HST, b–g, using the weighted mean spectrum of TRAPPIST-1 e from Transits 5 and 7. The resulting spectrum utilizes all of the available data and is, in effect, what one would observe if TRAPPIST-1 b–g transited simultaneously. This approach probes for shared spectral features, similar to the analysis of the double transit of TRAPPIST-1 b and c by de Wit et al. (2016). In this case, we are primarily interested in a stellar contribution that affects all transmission spectra similarly, such as surface active regions that are outside of all the planetary transit chords. This can be more easily studied in the combined spectra because the stellar contamination signal combines multiplicatively with any planetary transmission spectrum. The remaining combinations exclude the contribution from one of the six planets in turn, allowing us to examine the effect of the individual planets on the combined result. For each combination, we also included the sum of the corresponding K2 and SSO (Ducrot et al. 2018) and Spitzer 4.5 μm transit depths (Delrez et al. 2018b) in the fitting procedure and added the uncertainties of the individual transit depths in quadrature.

Figure 10 shows the TRAPPIST-1 b–g combined transmission spectrum, the out-of-transit stellar spectrum, and the best-fit flat and contamination models. The K2, SSO, and Spitzer transit depths are 3.6σ, 1.0σ, and 7.1σ below the mean of the combined HST transit depth, respectively. The combined HST transmission spectrum displays a notable decrease in transit depth around 1.4 μm, which coincides with a strong water absorption band. This "inverted" water feature is the opposite of the water absorption signature commonly observed in transiting exoplanet atmospheres (e.g., Sing et al. 2016). The offsets between instruments and the apparent 1.4 μm decrease are also evident in all five-planet combined transmission spectra (Figure 11), which illustrates that it is not due solely to the spectrum of an individual planet. We find that the offsets between instruments and the apparent 1.4 μm decrease are both reproduced well by the contamination model for each combined spectrum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 12 provides the complete results of the model fits to the combined spectra. For each combination, the AIC_c and BIC both prefer the contamination model. According to the ICs, each data set provides decisive evidence against the flat model (ΔAIC_c > + 10, ΔBIC > + 10). In each case, the ICs show that the data warrant the inclusion of the additional parameters in the contamination models.

Table 12.  Results of Stellar Contamination Model Fits to Combined Spectra

Combination	Model	Fitted Parameter									AICc	ΔAICc	BIC	ΔBIC	χ2	p
		D (%)	Tphot (K)	Tspot (K)	Tfac (K)	Fspot	Ffac	fspot	ffac	η
b–g	flat	${3.748}_{-0.013}^{+0.013}$	${2117}_{-128}^{+84}$	${1961}_{-132}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−687.7	39.1	−667.8	23.6	177.8	0.01
b–g	cont.	${3.467}_{-0.059}^{+0.058}$	${2425}_{-178}^{+168}$	${2006}_{-93}^{+127}$	${2957}_{-25}^{+43}$	${38}_{-8}^{+8}$	${48}_{-8}^{+6}$	${10}_{-10}^{+4}$	${45}_{-6}^{+6}$	${23}_{-2}^{+1}$	−726.8	⋯	−701.4	⋯	130.7	0.59
c–g	flat	${3.000}_{-0.011}^{+0.012}$	${2116}_{-123}^{+89}$	${1960}_{-128}^{+113}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−705.0	26.2	−685.0	20.8	164.8	0.05
c–g	cont.	${2.794}_{-0.051}^{+0.054}$	${2419}_{-192}^{+179}$	${2013}_{-98}^{+133}$	${2954}_{-25}^{+46}$	${38}_{-9}^{+8}$	${49}_{-9}^{+6}$	${10}_{-10}^{+5}$	${45}_{-6}^{+6}$	${23}_{-2}^{+1}$	−731.2	⋯	−705.9	⋯	130.4	0.59
b, d–g	flat	${3.043}_{-0.012}^{+0.012}$	${2117}_{-125}^{+86}$	${1961}_{-134}^{+107}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-1}^{+1}$	−697.0	29.1	−677.1	23.6	172.0	0.02
b, d–g	cont.	${2.808}_{-0.054}^{+0.053}$	${2480}_{-162}^{+206}$	${2033}_{-85}^{+119}$	${2949}_{-27}^{+51}$	${39}_{-8}^{+8}$	${50}_{-10}^{+7}$	${11}_{-11}^{+5}$	${45}_{-6}^{+7}$	${23}_{-2}^{+1}$	−726.1	⋯	−700.7	⋯	134.5	0.50
b–c, e–g	flat	${3.366}_{-0.013}^{+0.012}$	${2118}_{-129}^{+84}$	${1961}_{-136}^{+106}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−700.0	26.8	−680.1	21.3	166.4	0.04
b–c, e–g	cont.	${3.165}_{-0.055}^{+0.055}$	${2425}_{-211}^{+173}$	${2014}_{-101}^{+146}$	${2954}_{-25}^{+46}$	${37}_{-9}^{+8}$	${49}_{-9}^{+7}$	${10}_{-10}^{+5}$	${44}_{-6}^{+6}$	${23}_{-2}^{+1}$	−726.8	⋯	−701.4	⋯	131.2	0.58
b–d, f–g	flat	${3.254}_{-0.012}^{+0.012}$	${2117}_{-124}^{+87}$	${1961}_{-133}^{+107}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−693.6	34.5	−673.6	29.1	172.7	0.02
b–d, f–g	cont.	${2.999}_{-0.052}^{+0.052}$	${2444}_{-183}^{+181}$	${2015}_{-92}^{+125}$	${2952}_{-26}^{+48}$	${39}_{-9}^{+8}$	${49}_{-9}^{+7}$	${10}_{-10}^{+5}$	${46}_{-6}^{+6}$	${23}_{-2}^{+1}$	−728.1	⋯	−702.7	⋯	130.1	0.60
b–e, g	flat	${3.105}_{-0.012}^{+0.012}$	${2117}_{-125}^{+87}$	${1960}_{-133}^{+106}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−693.6	36.8	−673.6	31.4	174.6	0.02
b–e, g	cont.	${2.849}_{-0.051}^{+0.051}$	${2407}_{-191}^{+147}$	${2001}_{-95}^{+125}$	${2956}_{-24}^{+44}$	${38}_{-9}^{+8}$	${48}_{-8}^{+6}$	${9}_{-9}^{+4}$	${46}_{-6}^{+6}$	${23}_{-2}^{+1}$	−730.4	⋯	−705.0	⋯	130.3	0.60
b–f	flat	${2.970}_{-0.012}^{+0.011}$	${2117}_{-124}^{+90}$	${1960}_{-133}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−690.6	40.3	−670.6	35.9	178.5	0.01
b–f	cont.	${2.721}_{-0.049}^{+0.051}$	${2390}_{-164}^{+141}$	${1987}_{-96}^{+119}$	${2959}_{-23}^{+41}$	${38}_{-8}^{+8}$	${47}_{-7}^{+6}$	${9}_{-9}^{+4}$	${46}_{-6}^{+5}$	${23}_{-2}^{+1}$	−731.9	⋯	−706.5	⋯	129.8	0.61

Note. The combined spectra are listed in the first column. The combination b–g utilizes all of the current data, including the weighted mean of the two transmission spectra of TRAPPIST-1 e. The others are five-planet combinations, which in turn exclude the spectrum of a single planet. The remaining table elements are the same as those in Table 11.

Download table as:  ASCIITypeset image

Table 15.  Results of Stellar Contamination Model Fits to Single-Planet Spectra Using Only HST and Spitzer Data

Data Set	Model	Fitted Parameter									AICc	ΔAICc	BIC	ΔBIC	χ2	p
		D (%)	Tphot (K)	Tspot (K)	Tfac (K)	Fspot	Ffac	fspot	ffac	η
b	flat	${0.745}_{-0.006}^{+0.006}$	${2117}_{-128}^{+85}$	${1961}_{-134}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−707.8	−8.6	−687.9	−14.0	148.3	0.21
b	cont.	${0.675}_{-0.023}^{+0.021}$	${2210}_{-241}^{+190}$	${1980}_{-136}^{+119}$	${2964}_{-20}^{+36}$	${29}_{-16}^{+16}$	${47}_{-8}^{+5}$	${10}_{-10}^{+4}$	${52}_{-9}^{+11}$	${23}_{-2}^{+1}$	−699.2	⋯	−673.9	⋯	151.5	0.13
c	flat	${0.705}_{-0.005}^{+0.005}$	${2118}_{-126}^{+87}$	${1961}_{-134}^{+107}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−718.6	−6.8	−698.8	−12.2	139.5	0.38
c	cont.	${0.663}_{-0.020}^{+0.020}$	${2227}_{-236}^{+154}$	${1973}_{-160}^{+98}$	${2965}_{-19}^{+35}$	${29}_{-13}^{+15}$	${47}_{-7}^{+5}$	${9}_{-9}^{+4}$	${48}_{-9}^{+9}$	${23}_{-2}^{+1}$	−711.8	⋯	−686.6	⋯	140.4	0.31
d	flat	${0.388}_{-0.003}^{+0.003}$	${2117}_{-127}^{+86}$	${1961}_{-132}^{+109}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−696.4	25.7	−676.5	20.4	178.5	0.01
d	cont.	${0.309}_{-0.016}^{+0.015}$	${2653}_{-98}^{+151}$	${2030}_{-62}^{+91}$	${2931}_{-31}^{+62}$	${53}_{-8}^{+11}$	${54}_{-16}^{+12}$	${9}_{-9}^{+4}$	${42}_{-18}^{+19}$	${23}_{-2}^{+1}$	−722.1	⋯	−696.9	⋯	149.5	0.16
e (T5)a	flat	${0.478}_{-0.004}^{+0.004}$	${2117}_{-125}^{+87}$	${1961}_{-133}^{+109}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−732.6	−4.2	−712.8	−9.6	139.6	0.38
e (T5)a	cont.	${0.493}_{-0.021}^{+0.021}$	${2123}_{-120}^{+82}$	${1982}_{-120}^{+119}$	${2973}_{-15}^{+27}$	${15}_{-12}^{+8}$	${47}_{-6}^{+5}$	${12}_{-12}^{+5}$	${44}_{-7}^{+6}$	${23}_{-2}^{+1}$	−728.4	⋯	−703.2	⋯	139.0	0.34
e (T7)b	flat	${0.506}_{-0.004}^{+0.004}$	${2117}_{-125}^{+86}$	${1960}_{-136}^{+105}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-5}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−721.4	6.6	−701.5	1.3	152.0	0.15
e (T7)b	cont.	${0.447}_{-0.022}^{+0.021}$	${2559}_{-135}^{+187}$	${2030}_{-81}^{+113}$	${2941}_{-27}^{+59}$	${46}_{-10}^{+11}$	${51}_{-12}^{+9}$	${10}_{-10}^{+4}$	${39}_{-12}^{+14}$	${23}_{-2}^{+1}$	−728.0	⋯	−702.8	⋯	139.2	0.34
ec	flat	${0.494}_{-0.003}^{+0.003}$	${2116}_{-124}^{+88}$	${1960}_{-133}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−742.8	−7.5	−723.0	13.0	137.7	0.42
ec	cont.	${0.479}_{-0.023}^{+0.021}$	${2214}_{-178}^{+127}$	${1970}_{-161}^{+98}$	${2968}_{-18}^{+32}$	${27}_{-11}^{+11}$	${47}_{-7}^{+5}$	${9}_{-9}^{+4}$	${45}_{-7}^{+6}$	${23}_{-2}^{+1}$	−735.3	⋯	−710.0	⋯	139.1	0.34
f	flat	${0.643}_{-0.005}^{+0.006}$	${2116}_{-126}^{+85}$	${1961}_{-132}^{+109}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−728.5	−8.3	−708.6	−13.7	134.8	0.49
f	cont.	${0.625}_{-0.025}^{+0.025}$	${2196}_{-193}^{+116}$	${1984}_{-141}^{+112}$	${2969}_{-17}^{+31}$	${24}_{-13}^{+12}$	${48}_{-6}^{+5}$	${10}_{-10}^{+4}$	${46}_{-7}^{+7}$	${23}_{-2}^{+1}$	−720.2	⋯	−694.9	⋯	137.1	0.39
g	flat	${0.776}_{-0.006}^{+0.006}$	${2117}_{-128}^{+86}$	${1961}_{-131}^{+111}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−726.3	-4.7	−706.4	-10.3	134.7	0.49
g	cont.	${0.747}_{-0.028}^{+0.027}$	${2138}_{-138}^{+90}$	${1981}_{-124}^{+118}$	${2971}_{-15}^{+29}$	${18}_{-13}^{+9}$	${47}_{-6}^{+5}$	${11}_{-11}^{+5}$	${49}_{-7}^{+7}$	${23}_{-2}^{+1}$	−721.6	⋯	−696.3	⋯	134.2	0.45

Notes. The table elements are the same as those in Table 11. Each model has 142 data points (12 for the transmission spectrum and 130 for the stellar spectrum). Flat models have 135 degrees of freedom, and contamination models have 133.

^aTransit 5. ^bTransit 7. ^cWeighted mean of Transits 5 and 7.

Download table as:  ASCIITypeset image

Table 16.  Results of Stellar Contamination Model Fits to Combined Spectra Using HST+Spitzer Data Only

Combination	Model	Fitted Parameter									AICc	ΔAICc	BIC	ΔBIC	χ2	p
		D (%)	Tphot (K)	Tspot (K)	Tfac (K)	Fspot	Ffac	fspot	ffac	η
b–g	flat	${3.754}_{-0.013}^{+0.013}$	${2119}_{-125}^{+88}$	${1961}_{-133}^{+110}$	${2974}_{-14}^{+26}$	${16}_{-14}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−671.9	30.5	−652.1	25.0	169.0	0.03
b–g	cont.	${3.474}_{-0.059}^{+0.059}$	${2408}_{-252}^{+235}$	${2002}_{-134}^{+142}$	${2955}_{-24}^{+45}$	${38}_{-10}^{+9}$	${48}_{-10}^{+7}$	${10}_{-10}^{+5}$	${44}_{-7}^{+9}$	${23}_{-2}^{+1}$	−702.4	⋯	−677.1	⋯	131.3	0.52
c–g	flat	${3.005}_{-0.011}^{+0.012}$	${2116}_{-126}^{+86}$	${1960}_{-135}^{+106}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−687.9	17.2	−668.0	11.9	157.2	0.09
c–g	cont.	${2.796}_{-0.054}^{+0.054}$	${2351}_{-250}^{+175}$	${1987}_{-157}^{+120}$	${2958}_{-22}^{+42}$	${36}_{-10}^{+9}$	${48}_{-9}^{+6}$	${9}_{-9}^{+4}$	${45}_{-6}^{+8}$	${23}_{-2}^{+1}$	−705.1	⋯	−679.9	⋯	132.8	0.49
b, d–g	flat	${3.048}_{-0.012}^{+0.012}$	${2118}_{-126}^{+88}$	${1962}_{-130}^{+112}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−680.4	20.2	−660.6	14.7	163.9	0.05
b, d–g	cont.	${2.810}_{-0.055}^{+0.054}$	${2377}_{-275}^{+198}$	${1997}_{-163}^{+117}$	${2957}_{-23}^{+43}$	${37}_{-10}^{+10}$	${48}_{-9}^{+7}$	${9}_{-9}^{+4}$	${45}_{-7}^{+9}$	${23}_{-2}^{+1}$	−700.6	⋯	−675.3	⋯	136.7	0.40
b–c, e–g	flat	${3.372}_{-0.013}^{+0.013}$	${2116}_{-126}^{+88}$	${1959}_{-132}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−687.1	14.4	−667.2	9.0	154.7	0.12
b–c, e–g	cont.	${3.168}_{-0.054}^{+0.057}$	${2294}_{-223}^{+152}$	${1956}_{-165}^{+97}$	${2962}_{-21}^{+38}$	${33}_{-10}^{+8}$	${47}_{-8}^{+5}$	${8}_{-8}^{+4}$	${45}_{-6}^{+7}$	${23}_{-2}^{+1}$	−701.5	⋯	−676.2	⋯	133.5	0.47
b–d, f–g	flat	${3.259}_{-0.013}^{+0.013}$	${2116}_{-126}^{+85}$	${1961}_{-131}^{+109}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−675.9	25.8	−656.0	20.4	165.6	0.04
b–d, f–g	cont.	${3.007}_{-0.054}^{+0.054}$	${2392}_{-269}^{+223}$	${1997}_{-151}^{+128}$	${2956}_{-24}^{+44}$	${37}_{-10}^{+10}$	${48}_{-10}^{+7}$	${9}_{-9}^{+5}$	${45}_{-7}^{+9}$	${23}_{-2}^{+1}$	−701.7	⋯	−676.4	⋯	133.0	0.48
b–e, g	flat	${3.110}_{-0.012}^{+0.012}$	${2115}_{-122}^{+88}$	${1959}_{-132}^{+108}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−673.1	30.9	−653.2	25.5	170.1	0.02
b–e, g	cont.	${2.850}_{-0.053}^{+0.051}$	${2406}_{-238}^{+230}$	${2001}_{-127}^{+140}$	${2954}_{-25}^{+46}$	${38}_{-10}^{+9}$	${48}_{-9}^{+7}$	${9}_{-9}^{+4}$	${45}_{-7}^{+9}$	${23}_{-2}^{+1}$	−704.0	⋯	−678.7	⋯	131.9	0.51
b–f	flat	${2.976}_{-0.012}^{+0.012}$	${2118}_{-126}^{+87}$	${1961}_{-128}^{+112}$	${2974}_{-14}^{+26}$	${16}_{-15}^{+7}$	${46}_{-6}^{+5}$	⋯	⋯	${23}_{-2}^{+1}$	−672.2	31.3	−652.3	26.0	171.4	0.02
b–f	cont.	${2.733}_{-0.052}^{+0.050}$	${2502}_{-153}^{+259}$	${2028}_{-91}^{+151}$	${2950}_{-25}^{+50}$	${41}_{-10}^{+10}$	${50}_{-10}^{+8}$	${10}_{-10}^{+5}$	${43}_{-8}^{+11}$	${23}_{-2}^{+1}$	−703.5	⋯	−678.3	⋯	132.7	0.49

Note. The table elements are the same as those in Table 12. Each model has 142 data points (12 for the transmission spectrum and 130 for the stellar spectrum). Flat models have 135 degrees of freedom, and contamination models have 133.

Download table as:  ASCIITypeset image

6.4.1. Impact of the Stellar Spectrum

For all model fits, both to the individual and combined transmission spectra, the fitted value of η is ${23}_{-2}^{+1}$ (Tables 11 and 12). Notably, this is the case for all flat models, in which the stellar parameters are determined by the out-of-transit stellar spectrum alone and have no effect on the fit to the transmission spectra. This implies that the observational errors of the HST/WFC3 G141 spectrum of TRAPPIST-1 must be largely inflated to match the DRIFT-PHOENIX model spectra, even when using one composed of multiple component spectra, which should provide more flexibility to the model fits. Inspection of the right panel of Figure 10 shows that the data and model agree on the general shape of the spectrum but disagree on many details and the continuum level for the bluest and reddest wavelengths. For further clarity, Figure 12 illustrates these same spectra but with the true uncertainty on the observations shown. Also shown is a DRIFT-PHOENIX model spectrum for a star with TRAPPIST-1's effective temperature and surface gravity (T = 2559 K, $\mathrm{log}g=5.21;$ Gillon et al. 2017). This single-component model disagrees with the data significantly between roughly 1.34 and 1.6 μm, which demonstrates why even the flat models prefer multicomponent stellar spectra, as shown by the fitted values of the F_spot and F_fac parameters in Tables 11 and 12.

Zoom In Zoom Out Reset image size

The stellar spectrum contains many more data points than the transmission spectrum. Therefore, disagreement between the observed and model stellar spectra could strongly influence our results. In principle, the η parameter guards against this by effectively deweighting the stellar spectrum in the MCMC optimization procedure. Nonetheless, we investigated how the inclusion of the stellar spectrum in the model framework affects the results using the TRAPPIST-1 b–g data set. Considering only the transmission spectrum, the flat model has 14 data points and one fitted parameter, giving 13 degrees of freedom. The ${\chi }^{2}$ value for this model is 51.8, indicating that it is conclusively ruled out (p < 10⁻⁵). The AIC_c and BIC for this model are −129.3 and −129.0, respectively. By contrast, the contamination model has eight fitted parameters and thus gives 6 degrees of freedom when considering only the 14 data points of the transmission spectrum. The χ² value for the contamination model is 4.2, indicating that it is not ruled out by the data (p = 0.64). Additionally, the AIC_c and BIC for this model are −135.2 and −158.1, respectively. A comparison of the ICs for the two models yields ΔAIC_c = +5.9 and ΔBIC = +29.1 in favor of the contamination model. Thus, relative to the contamination model, the flat model is ruled out decisively by the BIC and strongly by the AIC_c, which more stringently penalizes model complexity. This exercise shows that the additional complexity of the contamination model is warranted by the transmission spectrum alone. For the remainder of this analysis, however, we opt to utilize the additional information in the stellar spectrum and discuss the results taking into account both the transmission and stellar spectra.

6.4.2. Impact of Instrumental Offsets

6.4.3. Physical Interpretation

Returning to the interpretation of the fits to the full K2+SSO+HST+Spitzer transmission spectra, we note that the posterior distributions for the stellar parameters (T_phot, T_spot, T_fac, F_spot, F_fac, f_spot, and f_fac) in the contamination models listed in Table 12 are consistent for all combined transmission spectra. They broadly agree on the same general picture for TRAPPIST-1, namely a heterogeneous photosphere comprised of three components: a hot component or faculae (T ∼ 3000 K) covering ∼50% of the projected stellar disk, a cool component or spots (T ∼ 2000 K) covering ∼40%, and an intermediate component or immaculate photosphere (T ∼ 2400 K) covering the remaining ∼10% of the disk. Within the region transited by the TRAPPIST-1 planets, spots are less prevalent, covering only ∼10% of the transit chord, which results in a spectral mismatch between the disk-integrated stellar spectrum and the light source for the transit spectroscopy observations.

Figure 13 illustrates examples of the joint posterior distributions for both modeling frameworks. For the contamination model, the distributions of two parameters, f_spot and T_fac, deserve special notice. Here f_spot piles on the lower bound of the allowed parameter space in the contamination models and is consistent at 1σ confidence with no spots being present within the transit chord in each case. The facula temperature is notable because it piles on the upper bound of the allowed parameter space, even in the flat models. This suggests that the HST spectrum of TRAPPIST-1 shows evidence for a high-temperature component beyond the limit of the DRIFT-PHOENIX model grid, a result in broad agreement with the bright spots proposed by Morris et al. (2018a) from an analysis of K2 and Spitzer photometry. Given the limitations of the model grid to sample high-enough values for T_fac, the fitted values of F_fac and f_fac are likely overestimated.

Zoom In Zoom Out Reset image size

The discrepancy between f_spot and F_spot suggests that the 28% of the projected stellar disk (or 56% of a half hemisphere) covered by the transit chords (Delrez et al. 2018b) is less spotted than the whole disk and seems to point to the presence of an active region not represented within the transit chord. As we discuss in Section 6.5, this arrangement bears some resemblance to the active high latitudes and circumpolar spot structures observed on fully convective M dwarfs (Phan-Bao et al. 2009; Barnes et al. 2015, 2017) and polar spots observed on active earlier-type stars (see Strassmeier 2009 and references therein). Still, while the ICs show that models that allow for stellar contamination in the transmission spectra are strongly preferred over those that assume no stellar contamination, the tendency of T_fac and f_spot to pile on the boundaries of their allowed ranges highlights the limitations of this analysis (discussed in more detail in Section 6.7) and cautions against overinterpreting the specific component temperatures or covering fractions determined by this analysis.

Nonetheless, these results demonstrate that, in principle, the features in the TRAPPIST-1 transmission spectra can result from a heterogeneous stellar photosphere and not from transmission through the planetary atmospheres. Of course, both stellar and planetary signals can contribute to the observed spectra, and future efforts to jointly constrain the contributions of each source may bear fruit. However, these results urge caution for planetary interpretations of observed features in near-infrared transmission spectra from low-mass host stars. We consider this point more broadly in Section 7.

We present here a relatively simple model to explain a complicated physical phenomenon, namely the heterogeneity of the surface of TRAPPIST-1 and the its effect on the transmission spectra of the TRAPPIST-1 planets. Despite the limitations of the approach, we show that models that allow for imprints of stellar features in the combined TRAPPIST-1 transmission spectra are strongly preferred to those that do not. The model fits point to large spot and faculae covering fractions, suggesting a highly heterogeneous photosphere for TRAPPIST-1. In the following paragraphs, we provide further context for interpreting this result.

Photometric variability shows that M dwarfs have highly heterogeneous surfaces. McQuillan et al. (2013) measured rotation periods and variability amplitudes for 1570 M dwarfs in the Kepler sample with masses between 0.3 and 0.5 M_☉. They found amplitudes ranging from 1.0 to 140.8 mmag in the Kepler bandpass, with a median amplitude of 7.6 mmag or 0.70% (McQuillan et al. 2013, Table 2). Newton et al. (2016) measured rotation periods for 387 field mid-to-late M dwarfs (M < 0.35 M_☉) in the MEarth bandpass (roughly i+z) and found that they typically vary in brightness by 1%–2% as they rotate. Using an ensemble of model M-dwarf photospheres with randomly distributed active regions, Rackham et al. (2018) found that a 1% I-band variability amplitude corresponds to spot covering fractions ${F}_{\mathrm{spot}}={14}_{-7 \% }^{+16 \% }$ and faculae covering fractions ${F}_{\mathrm{fac}}={63}_{-25 \% }^{+5 \% }$ across all M-dwarf spectral types, assuming a typical spot radius of 2°. Importantly, they found that the relation between active-region coverages and variability amplitudes is not linear, as variability monitoring only traces the nonaxisymmetric component of stellar surface features (see also Jackson & Jeffries 2012). Therefore, the typical amplitudes of M dwarfs in the Kepler and MEarth samples indicate that variable M dwarfs can be covered in large part by active regions.

As noted above, TRAPPIST-1 demonstrates a photometric variability of roughly 1% in the I+z bandpass. This has generally been interpreted as active regions rotating into and out of view (Gillon et al. 2016; Vida et al. 2017), though Morris et al. (2018a) suggested that the variability of TRAPPIST-1 may be driven by a magnetic activity timescale rather than a rotation period. Assuming the variability owes to rotation, Rackham et al. (2018) found that it could be caused by spot and faculae with covering fractions of ${F}_{\mathrm{spot}}={8}_{-7 \% }^{+18 \% }$ and ${F}_{\mathrm{fac}}={54}_{-46 \% }^{+16 \% }$, respectively. We used this information as priors in the current analysis, so the results presented here are not completely independent. However, we point out these previous results to show that the active-region coverages found in this analysis are neither unexpected for TRAPPIST-1 given its variability level nor unique among the population of M dwarfs.

The difference between the best-fit values that we find for F_spot and f_spot for the combined transmission spectra argues for variation in the latitudinal distribution of spots for TRAPPIST-1. This should not be surprising: Doppler imaging (Vogt & Penrod 1983) of early M dwarfs provides examples of both stars with spots emerging at preferential latitudes and those with spots at all latitudes (Barnes & Collier Cameron 2001; Barnes et al. 2004). However, the few Doppler images that exist for fully convective M dwarfs suggest that spots may emerge preferentially at high latitudes. To date, six mid-to-late M dwarfs have been studied with this technique: V374 Peg (M4V; Morin et al. 2008), G 164-31 (M4V; Phan-Bao et al. 2009), GJ 791.2A (M4.5V; Barnes et al. 2015, 2017), the binary GJ 65A and GJ 65B (M5.5V and M6V; Barnes et al. 2017), and LP 944-20 (M9V; Barnes et al. 2015). Three of these six targets (G 164-31, GJ 65A, and LP 944-20) only display spots at high-latitude or circumpolar regions. Additionally, GJ 791.2A shows a high-latitude circumpolar spot structure, as well as low-latitude spots. For earlier spectral types, which can be studied more readily with Doppler imaging techniques, polar spots are commonly observed on active stars (see Strassmeier 2009 and references therein) and may be stable and long-lived (Jeffers et al. 2007).

In the case of TRAPPIST-1, the combined transit chords of the seven TRAPPIST-1 planets probe a substantial portion (28%) of the projected stellar disk (Delrez et al. 2018b). However, given the frequent and prominent latitudinal variations in active-region occurrence observed in M dwarfs, it is unsurprising that our analysis suggests significant differences between the 28% equatorial zone and the polar regions of TRAPPIST-1. Therefore, while we cannot precisely constrain the latitudes of active regions in the current analysis, we suggest that features like the active high latitudes observed on fully convective M dwarfs or polar spots observed on active earlier-type stars may account for the different spot covering fractions that we fit to the transit chord and whole disk of TRAPPIST-1.

Spectral template fitting provides another approach for studying the heterogeneity of stellar surfaces. Observations of TiO molecular bands in dwarf stars with photospheres warmer than T ∼ 3500 K have been used to constrain spot sizes and covering fractions (Vogt 1979; Ramsey & Nations 1980; Vogt 1981). Using spectra of inactive G, K, and M dwarfs as templates, Neff et al. (1995) and later O'Neal et al. (1996, 1998) showed that TiO bands in medium-resolution (R ∼ 10,000) spectra of active G and K stars point to covering fractions of cool spots as large as 64%. Gully-Santiago et al. (2017) studied the T Tauri star LkCa 4 and found that spectra features in high-resolution near-infrared spectra were produced by hot (∼4100 K) and cool (∼2700–3000 K) photospheric components, with the cool component covering ∼80% of the stellar surface. The TiO bands in R ∼ 1000 spectra of a large sample of stars in the Pleiades show that spot covering fractions of ∼50% are common among K and M stars at ∼125 Myr (Fang et al. 2016).

In this work, we find that the low-resolution HST spectrum of TRAPPIST-1 is best fit by a template combination of three DRIFT-PHOENIX model stellar spectra. This is the case even for our flat models, in which the spot and faculae parameters are determined by fitting the stellar spectrum alone. While TRAPPIST-1 is both later-type and older than the stars in the examples referenced above, these previous efforts show that the inferred heterogeneity for TRAPPIST-1 is not extreme but rather typical of active stars.

Recently, Morris et al. (2018c) investigated the problem of stellar contamination from TRAPPIST-1 using a novel transit light-curve "self-contamination" technique (Morris et al. 2018b). They determined R_p/R_s for TRAPPIST-1 b–h (individually) using Spitzer transit light curves based on the durations of the ingress/egress and the durations from the mid-ingress to mid-egress, for which stellar heterogeneity has a negligible effect, and compared it to the R_p/R_s values determined by the transit depths. They found consistent R_p/R_s measurements from the two methods and claimed a nondetection of stellar contamination. However, the R_p/R_s values calculated by the light-curve self-contamination method have large uncertainties, e.g., ∼10% versus 1.4% in this work for TRAPPIST-1 b, and therefore do not have the sensitivity to probe the level of stellar contamination in transmission spectra that we study here. In this light, the nondetection of Morris et al. (2018c) should be more properly viewed as a weak upper limit on stellar contamination, rather than evidence for the lack of it. Put simply, none of the individual transit depths (HST, Spitzer, K2, and SOO) that we used in this work would satisfy the criterion used by Morris et al. (2018c) for rejecting the null hypothesis (nondetection of stellar contamination), although we show that the combined transmission spectra strongly favor a stellar contamination model.

While we build upon previous attempts to characterize the contribution of stellar heterogeneity to transmission spectra by imposing the constraint provided by the out-of-transit stellar spectrum, this model has several notable limitations. We utilize disk-integrated stellar spectra for the spectral components, which may neglect unique spectral features that can emerge from magnetically active regions (e.g., Norris et al. 2017). Similarly, we do not consider any chromospheric contribution to the stellar contamination spectrum. We neglect the effect of limb darkening and consider the photosphere to be static during the transit events, though the observed activity and relatively short rotation period of TRAPPIST-1 may allow for photospheric evolution on a timescale important for transit observations. We note that future efforts could improve on this initial analysis in these respects and more.

Nonetheless, despite the model limitations, we find that models including stellar heterogeneity effects provide significantly better fits to the combined TRAPPIST-1 transmission spectra than flat planetary transmission models. Interestingly, the spot and faculae covering fractions that we infer from the transmission spectra have considerably lower uncertainties than those inferred through modeling the star's rotational variability (Rackham et al. 2018). While the rotational variability of TRAPPIST-1 demonstrates the presence of photospheric heterogeneities, its magnitude poorly constrains their covering fractions. We find that the observed transmission spectra place tighter constraints on the stellar heterogeneity.

This interpretation of the spectra also offers testable predictions. We find the best-fit models to display notable decreases in the visual transit depths of the TRAPPIST-1 planets, which represent an observational challenge given the faintness of this ultracool dwarf in the visual, but may be observable with large ground-based telescopes. Broadband photometry over repeated transits can also overcome this challenge, and in Appendix C, we show that stellar contamination models fit to the HST+Spitzer data sets accurately predict the recently measured K2 (0.42–0.9 μm) and I+z (0.8–1.1 μm) combined transit depths (Ducrot et al. 2018). Looking forward, we find "inverse water features" akin to the 1.4 μm feature in the unexplored spectral region between 1.7 and 4 μm, which can be studied with JWST.

The contamination seen in the combined TRAPPIST-1 data emerges due to the transit light source effect: the difference between the baseline stellar spectrum (disk-integrated) and the spectrum of the transit chord (actual light source for the transmission spectroscopy), as described by Rackham et al. (2018) in detail. The presence of ∼200 ppm-level (inverted) water features as stellar contamination in the HST/WFC3 spectra should serve as a red flag for investigators planning major future HST/WFC3 or JWST transit spectroscopy campaigns on TRAPPIST-1 or similar host stars. The possible range of stellar spectral contaminations has been discussed by Rackham et al. (2018), and here we will briefly discuss two aspects not addressed there: connection with atmospheric retrievals and emission spectroscopy.

Atmospheric retrievals use a Bayesian exploration of possible atmospheric models to identify best-fitting models and derive confidence intervals (and posterior probability distributions) for the atmospheric pressure–temperature structures, molecular abundances, and cloud properties. These models have been shown to be powerful in characterizing atmospheres (e.g., Madhusudhan & Seager 2009; Benneke & Seager 2012; Waldmann et al. 2015) and are widely anticipated to be the primary tools for interpreting high-quality transmission spectra from HST and JWST (e.g., Greene et al. 2016; Morley et al. 2017). The stellar contamination seen in the TRAPPIST-1 planets and sub-Neptune GJ 1214b (Rackham et al. 2017), as well as those predicted by Rackham et al. (2018) for most M-dwarf stars, highlight the importance of including this effect in retrievals. Indeed, without modeling and correcting for the relatively strong (∼200 ppm-level) inverted water absorption feature introduced by the stellar contamination, one cannot hope to accurately measure the water abundance (probably <80 ppm levels; Morley et al. 2017) in the transmission spectra of small planets. Recent efforts to include the effects of stellar contamination in retrievals, however, show promise in separating stellar and planetary contributions to transmission spectra (Pinhas et al. 2018, in press; Bixel et al. 2018, submitted; Espinoza et al. 2018, submitted).

We point out, furthermore, that stellar emission (eclipse) spectroscopy will be affected much less by stellar contamination, as the planet is its own light source for emission or eclipse spectroscopy. In fact, for the limiting case of infinitely slowly rotating stars with a constant (nonevolving) starspot/facular pattern, there will be no stellar contamination. Therefore, JWST emission spectroscopy of the TRAPPIST-1 planets is likely to be much easier to interpret than HST and JWST transmission spectroscopy. We caution, however, that for rapidly rotating stars—such as TRAPPIST-1—the evolution of the stellar spectrum due to its heterogeneity (e.g., Apai et al. 2013; Yang et al. 2015) may also pose a nonnegligible astrophysical systematic.

A central question in HST and JWST transiting exoplanet observations is whether data from multiple transits can be combined. If the instrumental systematics can be robustly corrected and the astrophysical systematics are negligible, data can be coadded. This approach has been used in the HST/WFC3 IR transmission spectroscopy study of the sub-Neptune GJ 1214b by Kreidberg et al. (2014), where observations from 12 transits of the planet were coadded. In fact, the overall success of the empirical ramp-effect correction and the physically motivated RECTE model are very encouraging and suggest that photon noise–limited performance could be achievable for large, multi-epoch campaigns.

However, the astrophysical noise—stellar contamination through the transit light source effect (Rackham et al. 2018)—emerges as a major concern. Not only will the stellar contamination features often dominate the planetary absorption features in amplitude (as explored in detail in Rackham et al. 2018), but the contamination itself may change as the star rotates and its starspot/faculae coverage evolves over the course of a transit observation.

With a rotational period of 1.4 ± 0.05 days (Gillon et al. 2016), TRAPPIST-1 is a relatively rapidly rotating host star at the stellar/substellar boundary. It is worthwhile to briefly outline the implications of this rapid rotation by comparing four different timescales: starspot evolution, planetary orbits, stellar rotation, and the lengths of planetary transit observations.

Stellar contamination in transit spectra will depend on the heterogeneity of the projected stellar disk during the observations. A straightforward mechanism that alters the brightness and spectral distribution on the visible hemisphere of a star is the appearance, evolution, and disappearance of starspots and surrounding faculae. Although we have a limited understanding of starspot lifetimes, it seems reasonable to assume that their evolution may occur over timescales of 10–20 days, based on the similar processes observed in the Sun (Bradshaw & Hartigan 2014). This timescale is generally longer than the orbital periods of the TRAPPIST-1 planets. Therefore, we tentatively conclude that starspot evolution may only have a limited impact on stellar contamination variations in the case that consecutive planetary transits can be observed; nonetheless, it will have a major impact in coadding data from transits separated by months or years.

The stellar rotation, however, will introduce a more problematic source of time-varying stellar contamination. The rotational period of TRAPPIST-1 is shorter than the orbital periods of all but the closest planet (b), which means that even during subsequent transits of the same planet, a different (essentially random) stellar longitude will face the observer, i.e., the stellar contamination will be different between even the closest transits of the same planet.

The rapid rotation of TRAPPIST-1 poses another interesting problem: the rotational period (1.4 ± 0.05 days; Gillon et al. 2016) is comparable to the duration of a typical transit observation (five HST orbits, or about 1/3 of a day). This means that even during a single transit, the star will rotate by about 60°. Thus, about 1/3 of the visible stellar hemisphere will be different between the beginning and the end of the transit. This rapid change means that stellar contamination will not only be different between different transits but will change even during a single transit, making potential starspot modeling and correction even more difficult. We note here that the same effect may also complicate planet eclipse spectroscopy, albeit to a lesser extent.

With its rapid rotation, it may be tempting see TRAPPIST-1 as a worst-case scenario in terms of stellar contamination; however, it should be considered as a typical example of very late M dwarf or brown dwarf planet hosts. Photometric monitoring programs have found that many brown dwarfs are much faster rotators, with rotational periods of a few hours being not uncommon (Metchev et al. 2015). Spectroscopic monitoring of brown dwarfs (with the same HST/WFC3 IR grism as used for this study) showed that most, if not all, brown dwarfs have heterogeneous atmospheres (Buenzli et al. 2014), a result that was reinforced by a larger but less sensitive ground-based broadband photometry survey (Radigan et al. 2014) and a large and unbiased Spitzer photometric survey (Metchev et al. 2015). At least some, and perhaps most, of these brightness and spectral modulations are attributable to the rearrangement of condensate clouds driven by global circulation patterns (Apai et al. 2017).

Therefore, rotational modulations and stellar heterogeneity evolution should be expected to be very common for ultracool planet hosts—the coolest host stars and brown dwarfs. The rapidly evolving stellar contamination in the transit spectra of exoplanets orbiting such ultracool hosts will be a major challenge to obtaining the multi-epoch data necessary to build up the data quality required for atmospheric abundance analysis (Morley et al. 2017). In fact, this goal may only be achievable if we can develop a robust understanding of ultracool dwarf heterogeneities and derive a reliable contamination correction method for this astrophysical noise source.

De Wit et al. (2018) was published after the first submission this work. They independently reduced and analyzed data from program GO-14873. They adopted a ramp removal method that is similar to ours and different from that in de Wit et al. (2016). We identify a few minor differences between the reduction processes, none of which led to significant differences in the reduced data, but the differences in the analysis led to different interpretations.

• 1.  
  The most crucial difference between our study and that of de Wit et al. (2018) is the interpretations of the transmission spectra. De Wit et al. (2018) compared observed spectra with planetary models and flat lines, while we compared the observations with flat spectra, planetary atmosphere model spectra with water absorption features, and stellar contamination model spectra. We agree with de Wit et al. (2018) that flat spectra are preferred over spectra with planetary water absorption features, but we also show that stellar contamination models are favored over flat spectra for the combined TRAPPIST-1 transmission spectra.
• 2.  
  Each visit in de Wit et al. (2018) was designed to include two transit events, and the three events identified in Table 2 are genuine transits. While we discarded some data sets heavily affected by severe CR events and the resulting guiding failure due to SAA passages, de Wit et al. (2018) attempted to recover all the data. Although the Visit 1 data—which covered the transit by TRAPPIST-1 f—was successfully recovered, it was not possible to recover the other two events.
• 3.  
  De Wit et al. (2018) aligned every single readout, while we applied the alignment on each ima file frame. As for the CR correction, in de Wit et al. (2018), both spatial and temporal interpolation were used to replace the pixels affected by CRs.
• 4.  
  De Wit et al. (2018) adopted a ramp correction formula that resembles the correction predicted by the RECTE model (Zhou et al. 2017). They reported an average standard deviation of residuals in the broadband light curve of 220 ppm, while we reached 198 ppm—very similar values given the uncertainties. For the individual (narrowband) bins, they reported average uncertainties of 545, 526, 493, and 494 ppm for Visits 1 to 4, respectively; by comparison, we find average uncertainties of 500, 550, 487, and 454 ppm for the same data sets. Thus, the residuals from the two reductions are very similar.