An HST/STIS Optical Transmission Spectrum of Warm Neptune GJ 436b

Aerosols may form in the atmosphere of sub-Jovian exoplanets as clouds of condensates or as hazes of photochemical products (e.g., hydrocarbons). Aerosols in the atmospheres of sub-Jovian exoplanets have long plagued transit observations due to their ability to mute potential absorption features; the long path length through the atmospheres from the geometry of transit observations amplifies this problem significantly (Fortney 2005). This issue can be especially troublesome for small exoplanets,   where small scale heights and planet-to-star radius ratios can lead to intrinsic variations of transit depth with wavelength that are already small (i.e., on the order of the systematic noise), even in an aerosol-free scenario.   In some exoplanet studies, clouds are defined as a gray (wavelength-independent) opacity source, while hazes are defined as a scattering opacity source that can induce slopes in the spectrum. For this work, we define clouds and hazes based on their formation and physical properties (i.e., clouds are condensed aerosols, while hazes are aerosols formed through photochemical processes), not on their effect on the planet's spectrum. In the models we use to interpret our data, whether an aerosol acts as a gray or scattering opacity depends primarily on the particle size rather than its classification as a cloud or haze.

High-metallicity atmospheres can affect the planet's spectrum in ways similar to aerosols. High metallicity results in a high mean molecular mass atmosphere, causing a reduced atmospheric scale height. This reduced scale height results in a smaller signal in the transmission spectrum, serving to mute spectral features. High metallicity can be a natural outcome of the formation for low-mass planets through core accretion (Fortney et al. 2013; Thorngren et al. 2016; Venturini et al. 2016). Though distinguishing between high-metallicity and aerosol-rich atmospheres is difficult and requires data with high signal-to-noise ratio, Benneke & Seager (2013) provided a framework to do so by measuring line wing steepness and the relative absorption depth of different spectral features. Additionally, the transmission spectra of high-metallicity atmospheres and atmospheres with aerosols begin to diverge at both short and long wavelengths, providing an opportunity to break this degeneracy (see Section 5.2).

Observations at optical wavelengths provide a path forward in studying clouds, hazes, and high-metallicity exoplanet atmospheres in transit. As mentioned above, scattering by small aerosol particles can dominate the optical spectrum by producing a slope toward larger transit depth at shorter wavelengths.  Cloud-free atmospheres, on the other hand, have transmission spectra that reach minimum transit depths near 0.5 μm before Rayleigh scattering dominates shortward of 0.5 μm and molecular opacities dominate longward of 0.9 μm. Meanwhile, atmospheres with large aerosol particles will remain flat at optical wavelengths. Thus the optical spectrum provides a unique way to characterize opacity sources even in an otherwise featureless spectrum. Additionally, absorption from atomic Na and K can shape much of an exoplanet's optical transmission spectrum (e.g., Seager & Sasselov 2000; Nikolov et al. 2014; Sing et al. 2015). Characterization of these features can provide a measurement of the atmosphere at lower pressures than those probed by the infrared spectrum (Sing et al. 2008; Vidal-Madjar et al. 2011a, 2011b; Heng et al. 2015; Wyttenbach et al. 2015). The absence of Na and K in the optical spectrum of an exoplanet may indicate the condensation of these elements into clouds (e.g., KCl and Na₂S; Morley et al. 2013).

Caution is necessary when interpreting slopes at optical wavelengths in transmission spectra because occulted and unocculted star spots on the photosphere of the host star  can produce slopes in the optical transmission spectrum (Berta et al. 2011; Oshagh et al. 2013; McCullough et al. 2014; Oshagh et al. 2014; Rackham et al. 2017, see Section 4.1). Photometric monitoring can reveal changes in the star spot or plage filling factor as well as the overall activity cycle. Brightness modulations from stellar rotation, star spot variability, and activity cycles all need to be accounted for in secondary eclipse and phase curve observations (see Section 2.2). Spectroscopic monitoring provides information on the host star's absolute activity level via activity indicators like the Ca ii H and K lines.

GJ 436b is a 21.4 ${M}_{\oplus }$ (1.25 M_Neptune, 0.0673 M_Jupiter) warm Neptune with a radius of 4.2 ${R}_{\oplus }$ (1.1 R_Neptune, 0.37 R_Jupiter; Turner et al. 2016; Trifonov et al. 2017). Discovered by radial velocity by Butler et al. (2004) as the first Neptune-mass exoplanet, it was subsequently found to transit by Gillon et al. (2007). Transit spectra were first obtained by Pont et al. (2009) using NICMOS on the Hubble Space Telescope (HST), placing upper limits of a few parts per 10,000 on the potential water absorption feature at 1.4 μm. HST/WFC3 transit spectra also revealed a featureless spectrum, ruling out a cloud-free, hydrogen-dominated atmosphere (Knutson et al. 2014a). Using Spitzer transit measurements, Beaulieu et al. (2011) claimed to detect CH₄ and found no evidence of CO or CO₂ due to large absorption measured in the 3.6 and 8 μm IRAC bands. This was refuted by Knutson et al. (2011), who hypothesized that stellar activity caused the spectrum to vary not only in wavelength but also with time. Reanalysis of these data by Lanotte et al. (2014) and Morello et al. (2015) using new detrending techniques found that the transmission spectrum of GJ 436b was constant with wavelength and did not vary between epochs.

While the near-infrared transmission spectrum of GJ 436b is featureless given current precision, large Ly α absorption has been measured. Up to 56% of the flux in the Ly α wings is absorbed beginning two hours before the optical transit and lasting 10–25 hours after the optical transit (Ehrenreich et al. 2011; Kulow et al. 2014; Ehrenreich et al. 2015; Lavie et al. 2017). This has been interpreted to be evidence for an escaping exospheric cloud of atomic hydrogen, shaped by stellar radiation pressure and interactions with the stellar wind (Bourrier et al. 2015, 2016).

GJ 436b's dayside spectrum has been observed through secondary eclipse measurements with Spitzer. GJ 436b's equilibrium temperature of 700–800 K would imply CH₄ is the most abundant carbon-bearing molecule in chemical equilibrium; however, measurements from Stevenson et al. (2010) suggested that CO was in high abundance, rather than CH₄. Additional studies have supported the assertion that GJ 436b is enhanced in CO and CO₂ and deficient in CH₄ (Madhusudhan & Seager 2011; Agúndez et al. 2012; Lanotte et al. 2014). This could potentially be explained by disequilibrium processes like vertical mixing and tidal heating. Photochemistry is likely not the cause of this CH₄ deficiency; Line et al. (2011) estimate that the rate at which CH₄ should be destroyed was much less than could explain the observations.  

A more recent analysis including both self-consistent modeling and retrievals   of the emission and transit spectra placed a 3σ lower limit on the metallicity at 106× solar (Morley et al. 2017). A cloud-free atmosphere is still possible if metallicity is on the order of 1000× solar; for lower metallicities, clouds are needed to help mute spectral features to match observations. Additionally, Morley et al. (2017) found that models with disequilibrium chemistry through quenching of CH₄, CO, and CO₂ with enhanced internal heating, presumably from tides, best matched the data. These results tend to agree with previous modeling of Neptune-sized planets from Moses et al. (2013), who found CO enrichment and CH₄ depletion to be a natural consequence of high metallicity. In this work, we present new observations at optical wavelengths in order to test the conclusions of Morley et al. (2017).

At the present time, GJ 436b is a candidate target for the NIRISS, NIRCam, and MIRI James Webb Space Telescope (JWST) Guaranteed Time Observation (GTO) programs. These programs will observe multiple secondary eclipses of the planet from 0.7 to 11 μm, providing an unprecedented look at the atmospheric composition and structure of GJ 436b's day side (see Greene et al. 2016).   We discuss how complementary transit observations will help distinguish cloudy and high-metallicity scenarios and determine conditions at the planet's terminator.

Here, we present the first space-based measurements of the optical transit spectrum of GJ 436b, interpret the full optical-to-IR spectrum, and constrain the rotation period and activity cycle of the host star. We organize the paper as follows. In Section 2 we explain the observations and data-reduction procedures used. Section 3 describes our light curve fitting techniques. In Section 4, we show our results, including a look at the effects of stellar variability and   the use of different orbital solutions. In Section 5, we compare our results to other sub-Jovian exoplanets and describe a common trend found in their optical transmission spectrum. We close with a discussion of expectations for future JWST transmission spectroscopy of GJ 436b.

Two transits of GJ 436b were observed on 2015 June 10 UT and 2016 June 14 UT using the STIS on the HST with the G750L grism (0.53–1.03 μm) as part of GO 13665 (PI Benneke). The 52 × 2 arsec² slit was used to minimize slit losses, while a 128-pixel subarray mode reduced readout overhead. Each transit consisted of four total HST orbits. The first orbit of a visit often exhibited strong systematic variations, inconsistent with the systematics in subsequent orbits. For this reason, after confirmation of this phenomenon, the first orbit was not included in the analysis. The second and fourth orbits are used to characterize the stellar baseline flux, while the transit occurs during the third orbit. Exposure times were 100 s, resulting in 20 frames per HST orbit. The first frame of each orbit consistently shows anomalous flux, so this frame is not included in the light curve fitting. Both visits were scheduled with the same orbital phasing, with the in-transit HST orbit covering ingress through the transit center to about four or five frames after transit center.

An STIS pipeline was built for the program using existing CALSTIS routines in addition to custom procedures for cosmic-ray identification. This pipeline was tested and validated using previous STIS data sets and is described below.  

2.1.1. Cosmic-ray Identification and Removal

2.1.2. Spectral Extraction

GJ 436 (M2.5V) has been monitored in Strömgren b and y filters for the past 14 years using the Tennessee State University's T12 0.8 m Automatic Photoelectric Telescope (APT) at Fairborn Observatory in southern Arizona (Henry 1999; Eaton et al. 2003; Henry & Winn 2008) in order to better characterize the star and how star spots may affect our transit observations. In these observations, we nod the telescope between GJ 436 and comparison stars of comparable or greater brightness; we use these stars to remove extinction and seeing effects. Of the three comparison stars analyzed, HD 102555 (hereafter C1) and HD 103676 (C2) were the most constant, with very little variation between them over the 14 years. The third star, HD 99518, showed a gradual brightening throughout the 14 years of observations. Altogether we have obtained 1735 measurements, which we present in Table 1 and Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Strömgren Photometric Observations of GJ 436b

Date	(P−C1C2a)b	(P−C1C2)y	(P−C1C2)${}_{(b-y)}$	(P−C1C2)${}_{(b+y)/2}$
(BJDTDB − 2,400,000)	(mag)	(mag)	(mag)	(mag)
52970.9856	4.40370	3.60900	0.79470	4.00635
52971.9748	4.40740	3.61605	0.79135	4.01175
52972.9752	4.40685	3.61350	0.79335	4.01020
52973.9700	4.39960	3.60950	0.79010	4.00455
52977.9596	4.40075	3.60855	0.79220	4.00470
52984.9665	4.40960	3.61610	0.79350	4.01285

Note.

^aC1C2 denotes that the differential magnitudes are computed with respect to the mean brightness of comparison stars C1 and C2.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The photometry of GJ 436 in Figure 1 shows an obvious long-term variation with a period of roughly 7.4 years and a peak-to-peak amplitude of 10 mmag. This coherent signal is consistent in period and amplitude with observations of many low-mass stars that are interpreted as stellar activity cycles (Suárez Mascareño et al. 2016). We therefore conclude that we have measured GJ 436's stellar activity cycle for the first time. This activity cycle is consistent with the finding that early M-type stars have magnetic cycles that are on average 6.0 ± 2.9 years and mid-M-type stars average 7.1 ± 2.7 years (Suárez Mascareño et al. 2016).

After removing this signal and a linear trend, a periodogram analysis reveals a strong signal with a period of 44.1 ± 0.2 days and peak-to-peak amplitude of 3 mmag.¹² To estimate the uncertainty on the rotation period, we split the data into 3-year blocks and computed a Lomb–Scargle periodogram of each data subset. The 44-day signal is not detectable in each individual season's photometry, but a strong peak is visible in the periodogram for each of these five contiguous subblocks. We measured the location of the peak of each of these five periodograms and find the mean and standard deviation of the mean to be 44.1 ± 0.2 days, which we interpret as GJ 436b's rotation period and the uncertainty on that parameter.

  The rotational period we measure for GJ 436 is similar to periods observed in many low-mass stars   (Irwin et al. 2011; Suárez Mascareño et al. 2016). This measurement is also consistent with an early rotation period derived from spectroscopic indicators (Demory et al. 2007), but inconsistent with the 57-day period derived using a 1-year segment of our data (Knutson et al. 2011). When using those same data, we also see a 57-day periodogram peak, but that season is the only one of our 14 whose periodogram shows a peak at that period. We conclude that GJ 436's rotation period is indeed 44.1 ± 0.2 days. This rotation rate is consistent with the interpretation that GJ 436  is not an active star and is at least a few gigayears old (Saffe et al. 2005; Kiraga & Stepień 2007; Maness et al. 2007; Jenkins et al. 2009; Sanz-Forcada et al. 2010). With respect to the relationship between rotational period and activity cycle, GJ 436  fits in with other M stars (Suárez Mascareño et al. 2016).

Figure 2 shows the same photometric observations, but for the individual b and y filters, showing that as GJ 436 first gets dimmer, the star becomes bluer. Most chromospherically active stars get redder as the star gets fainter, implying that dark spots dominate the variability (e.g., Innis et al. 1997); however, another well-known chromospherically active RS CVn variable binary, UX Ari, trends as GJ 436 does, becoming bluer as it becomes dimmer (Padmakar & Pandey 1999; Aarum Ulvås & Henry 2003). For GJ 436, this behavior does not seem repetitive, since during the next dimming phase, the star becomes slightly more red. During the current brightening phase, GJ 436 remains red compared to the 14-year mean color.

UX Ari's behavior of becoming bluer as it gets dimmer is seen both within the binary's orbital period of 6.4 days (Carlos & Popper 1971) and throughout the 25-year activity cycle (Aarum Ulvås & Henry 2003). Two explanations for this behavior include flare and facular activity (Rodono & Cutispoto 1992) or that the relative component of the hotter (and bluer) member of the binary contributes more flux as the cooler one becomes more spotted (Mohin & Raveendran 1989; Raveendran & Mohin 1995). Aarum Ulvås & Engvold (2003) showed that this behavior can be successfully reproduced using a model that includes dark spots surrounded by bright faculae on the active K star. We thus suggest that GJ 436's odd color behavior as activity changes may be due to the interplay between dark star spots and the faculae that surround them. It may be the case that different activity cycles have different proportions of star spots and faculae. GJ 436's stellar activity cycles will be further studied by continuing long-term stellar monitoring of GJ 436 and should be kept in mind when interpreting future high-precision transit spectra of GJ 436b.

Since interepoch stellar variability and inhomogenous stellar surfaces can induce spurious features in transmission spectroscopy (Knutson et al. 2011; Oshagh et al. 2013; Fraine et al. 2014; McCullough et al. 2014), it is imperative to check that the brightness variations of GJ 436 do not bias our atmospheric measurements. In Figure 1, we indicate the time and rotational phase of each STIS/G750L and WFC3/G141 transit observation. The WFC3 visits all occur within about a month and therefore span a range of rotational phases but essentially a single epoch of stellar activity. In contrast, the STIS visits are separated by a year, with the first occurring fairly near stellar minimum, but at otherwise similar magnitude in the rotation period. We discuss the effects of this stellar variability on the transmission spectrum in Section 4.1.

Crucial to the fitting of transit light curves is proper knowledge of the stellar limb darkening since the effects of limb-darkening biases on the calculated transit depth can be on the order of the atmospheric features we expect to measure. Ideally, one could fit for the coefficients that describe the limb darkening (Kreidberg et al. 2014); however, the low signal-to-noise ratio, low temporal sampling, and limited phase coverage of our observations prevent this from being a viable option in our spectral analysis. We thus use stellar models to calculate stellar limb profiles in order to estimate limb-darkening coefficients (LDCs).

We use the "Limb Darkening Toolkit" (LDTk), which uses the PHOENIX stellar models from Husser et al. (2013) to calculate LDCs for a variety of different parameterizations (Parviainen & Aigrain 2015). We use interferometrically determined stellar parameters from von Braun et al. (2012): ${T}_{\mathrm{eff}}=3416$ K, log(g) = 4.843, and Z = 0.02 dex. We also rescaled the μ values such that μ = 0 occurs where $\left|\tfrac{{dI}}{d\mu }\right|$ is at maximum, as suggested in Espinoza & Jordán (2015). We chose the nonlinear limb-darkening law to fit four coefficients to the stellar limb profile. Since we are fixing the LDCs, it is to our advantage to use more coefficients than in the linear or quadratic limb-darkening law, since the resulting fit to the modeled limb intensity profile will be more exact.

We tested LDTk by running a custom PHOENIX model for GJ 436. This has the advantage of using the exact known stellar parameters (Teff, logg, and metallicity) versus interpolating from the grid of models as in LDTk. We found that given the relatively large uncertainties in our data, it made little difference which method was used. We chose to present the LDTk LDCs for ease of repeatability.

We checked the robustness of the model LDCs by testing LDCs from stellar models ±100 K, the current uncertainty for GJ 436's effective temperature. The resulting transit spectra for GJ 436b for the different LDCs are within the error bars of the spectra using LDCs of the best-known stellar parameters. We consider our transit spectra to be robust to uncertainties in the stellar parameters, but caution is warranted as fixed LDCs can lead to biases, especially in the absolute transit depth (Espinoza & Jordán 2015).

We use the analytical transit model from Mandel & Agol (2002) and the transit center, orbital period, inclination, ratio of the orbital distance to the stellar radius ($a/{R}_{* }$), ratio of the radius of the planet to the radius of the star (${R}_{p}/{R}_{* }$), and four parameter nonlinear LDCs to describe the transit. Due to the observed phasing and the absence of data at the transit's egress, the constraints that our STIS data alone can put on the orbital parameters, particularly the inclination and ratio of orbital distance to the stellar radius ($a/{R}_{* }$), are limited.   We therefore chose to fix these values to the values found in Morello et al. (2015) using PHOENIX LDCs. The independent component analysis used in that work has been shown to give repeatable and accurate measurements (Ingalls et al. 2016). The parameters are listed in Table 2. The raw and corrected white light curves are shown in Figure 3, along with the resulting residuals.

Zoom In Zoom Out Reset image size

Table 2.  Transit Properties Used for GJ 436b

Parameter	Value
Transit center (Tc) (BJDTDB)a	2454222.616632 ± 0.00012
Period (days)a	2.6438986 ± 0.0000016
Inclination (°)	86.49 ± 0.12
$a/{R}_{* }$	13.82 ± 0.34
Impact parameter (b)	0.846 ± 0.05
Eccentricityb	0.16 ± 0.02
Longitude of periastron (deg)b	351 ± 1.2

Notes.

^aFrom Cáceres et al. (2009). ^bFrom Maness et al. (2007).

Download table as:  ASCIITypeset image

  The spectra were then split into wavelength bins to see how GJ 436b's transit depth varied as a function of wavelength. As mentioned above, we fixed the orbital parameters, allowing only the transit depth and systematic model parameters to vary. Many different wavelength binning schemes were tested, and we present here an analysis with 10 bins of approximately 500 Å width. Various numbers of smaller bins and locations were analyzed in a search for absorption from Na or K in the atmosphere of GJ 436b at approximately 5895 and 7684 Å, respectively (see Figure 4). To fit a transit model to the data as well as estimate uncertainties, we use both Levenberg–Marquardt and Markov chain Monte Carlo algorithms. Markov chain Monte Carlo, or MCMC, sampling methods are often used in parameter estimation as they provide empirically estimated uncertainties by exploring the posterior space of parameter likelihoods. We use the affine-invariant MCMC sampler implementation "emcee" (Foreman-Mackey et al. 2012). We find that both Levenberg–Marquardt and MCMC give very similar results for the parameters and the uncertainties estimation. This is likely a result of the Gaussianity of the posterior space when orbital parameters and LDCs are not free parameters. Because of the number of fits required, we chose to use the Levenberg–Marquardt for the marginalization analysis described in Section 3.3. The average standard deviation of residuals for our bins is about 1.25× the expectation from photon noise alone, similar to other transit spectroscopy studies with STIS (Sing et al. 2015).

Zoom In Zoom Out Reset image size

The top panel of Figure 3 shows the effect of instrumental systematics throughout the orbit. Wakeford et al. (2016) and Sing et al. (2016) describe a method for accounting for systematic uncertainties in HST data by marginalizing over the results from many different systematic models. Instead of choosing a single systematic model as the "correct" model, the full spectral time-series analysis is repeated for a variety of systematic models. The   measured transit depths are then combined in a way that weights systematic models that describe the data well, while also penalizing models with many free parameters. For this analysis, we choose five different covariates to parameterize the systematics found in STIS data: S1, HST orbital phase; S2, time; S3, position of the spectral trace in the spatial direction; S4, slope of the spectral trace; and S5, position of the spectral trace in the dispersion direction. S1, S3, S4, and S5 are well-known parameters used to describe systematics caused by HST's low Earth orbit; HST is continually cycling through day–night temperature variations throughout its  90-minute orbit. This results in small changes in the telescope's focus as components expand and contract (Hasan & Bely 1994; Sing et al. 2013).   S2 accounts for a visit-long linear change in the measured flux from the star and is similarly thought to be an effect from the telescope itself. S3 and S4 are found through a linear fit to the spectral trace before image rectification. S5 is found through the cross-correlation of the spectrum with a reference spectrum (in this case the fringe flat) in the cross-dispersion direction during 1D spectral extraction.

Zoom In Zoom Out Reset image size

  Each covariate is fit with different orders of polynomial: S1 is fit from first up to fourth order, which previous STIS studies have shown to adequately capture this trend (e.g., Sing et al. 2013; Nikolov et al. 2014). S2 is kept as a first-order linear trend, as we find no evidence of a quadratic variation in flux with time, and such a variation would be degenerate with the transit model. S3, S4, and S5 are all fit up to second order. As we show below, our analysis and previous analyses have not found justification to include any higher orders.   Each systematic model is then used to fit the data as described in Section 3.2, resulting in 108 different fits to the data. Figure 5 shows the difference in the fit using the least versus the most complex systematics model for each visit. The evidence of fit, or the marginal likelihood, of each model is calculated in order to compare each of the models. Models that fit the data well will have a high evidence-of-fit and will therefore be given a greater weight in the final marginalization. However, each model will be penalized according to its complexity. This is commonly done with the Bayesian information criterion (BIC), or as in Wakeford et al. (2016) with the Akaike information criterion (AIC).

Zoom In Zoom Out Reset image size

The AIC evidence function we use is

$$\begin{eqnarray}&&\mathrm{ln}{E}_{q}=-N\mathrm{ln}\sigma -0.5N\mathrm{ln}2\pi -0.5{\chi }^{2}-M,\end{eqnarray} \tag{ 1 }$$

where E_q is the "evidence" for a given systematics model, N is the number of data points being fit, σ is the uncertainty placed on the data, χ is the chi-squared statistic, and M is the number of free parameters being fit (see Equation (13) of Wakeford et al. 2016).

We find that different systematics models are preferred for each visit. Both visits require systematics models with a third-order polynomial of the HST orbit phase (S1). However, the preferred systematics model for visit 1 also includes a second-order polynomial to describe the slope of the spectral trace on the detector. This suggests that different systematics may be affecting the data for different STIS visits, even if the target and phasing are the same. The difference in evidence between systematics models with and without the additional spectral trace covariate was not large (${{\rm{\Delta }}}_{\mathrm{AIC}}=3.5$), however.

The final results for the marginalized transit depths do not depend on a single systematics model. Since the exact origins of many of the systematic trends seen in HST data are not fully understood, it would be impossible to claim a single "correct" systematics model. The advantage that the marginalized results do not depend on a single systematics model is shared by the more complex Gaussian process (GP) technique, but unlike GP, our method is still a parameterization. Figure 6 shows the mean transit depth of all 108 systematics models weighted by their evidence for the individual visits. Figure 7 and Table 3 show the spectra from both visits combined via a weighted average along with several model scenarios and previous observations. We find that the two visits are statistically consistent with each other by calculating the chi-squared between the two visits. We note that the long-wavelength half of the visit 2 transit spectrum appears consistently above the long-wavelength half of the visit 1 spectrum. We show below that this difference is likely not from stellar variability (see Section 4.1).

Zoom In Zoom Out Reset image size

Another complementary approach to removing systematics in-volves correcting for the common-mode wavelength-independent systematics before the removal of wavelength-dependent systematics (Huitson et al. 2012; Nikolov et al. 2014; Sing et al. 2016). The common-mode systematics are found by fitting a full systematics model to the white light curve data and then are "divided out" of the individual spectral bins. The spectral bins are theoretically left with only the wavelength-dependent systematics. This method can allow for the fits to the individual spectral bins to require fewer free parameters, helping to decrease the uncertainty in ${R}_{p}/{R}_{* }$. We implemented the common-mode correction to this analysis, but found it did not significantly improve either the fits or the estimated uncertainties. The common-mode correction did improve the evidence for less complex systematics models, but none increased the evidence enough to justify their use over more complex systematics models.

It was found that some areas of the detector exhibited different systematic trends than the rest of the detector. Figure 8 shows the difference in the measured raw flux from two adjacent 20 Å bins near the stellar Na line during visit 2. While still being described by the same systematics model, the overall trends are different, especially in those systematics associated with HST's orbital phase. Instead of decreasing in flux throughout the orbit, some bins increase throughout the orbit. The same effect is seen in visit 1. Areas like this on the detector will prevent the application of the common-mode systematic technique from being beneficial.

Zoom In Zoom Out Reset image size

An explanation for this effect could be uncorrected systematics in the position of the spectral trace. Since the bins that exhibit different systematics tend to be those in which the stellar spectrum is decreasing with increasing wavelength at the bin edges, shifts in the spectral trace can lead to strong and differing systematics in adjacent bins. However, even after applying the dispersion solution and shifting the spectra to a common rest frame via cross-correlation, this effect is still present.

  Since unocculted star spots and plages can cause slopes at optical wavelengths, we investigate the role stellar activity may play in GJ 436b's transit spectrum (Berta et al. 2011; McCullough et al. 2014; Rackham et al. 2017). The difference between the flux from a stellar photosphere and the flux from a star spot is a wavelength-dependent quantity since star spots are redder due to their lower temperature.   When a planet transits the stellar disk, a slope will be induced in the spectrum because unocculted star spots will contribute a greater proportion of the observed stellar flux at shorter wavelengths, making the planet appear larger at short wavelengths. Stellar plages will have a similar but opposite effect on the transit spectrum.

McCullough et al. (2014) analyzed the effect of unocculted star spots on the spectrum of HD 189733b, concluding that they could explain the observed spectral slope as well as Rayleigh scattering. The apparent transit depth as a function of wavelength may be modeled as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\tilde{{R}_{p}^{2}}}{\tilde{{R}_{s}^{2}}}=\displaystyle \frac{{R}_{p}^{2}}{{R}_{s}^{2}}\displaystyle \frac{1}{1-\delta (1-{F}_{\nu }(\mathrm{spot})/{F}_{\nu }(\mathrm{phot}))},\end{eqnarray} \tag{ 2 }$$

where $\tilde{{R}_{p}^{2}}/\tilde{{R}_{s}^{2}}$ is the observed transit depth, ${R}_{p}^{2}/{R}_{s}^{2}$ is the actual radius ratio squared, δ is the fraction of the star's projected surface area that contains spots, and ${F}_{\nu }(\mathrm{spot})/{F}_{\nu }(\mathrm{phot})$ is the ratio of flux from the star spot to the flux from the photosphere (McCullough et al. 2014).   Since the optical spectrum of M dwarfs can be rich in opacity features, we model ${F}_{\nu }(\mathrm{phot})$ and ${F}_{\nu }(\mathrm{spot})$ using PHOENIX atmosphere models rather than approximating them as blackbodies. Given the results of our photometric stellar monitoring, we estimate below the effect that stellar activity could have on GJ 436b's transmission spectrum.

  The maximum difference in stellar flux from our monitoring of GJ 436 is 13 mmag (∼10 mmag from the stellar activity cycle and 3 mmag from rotational variability).  A difference in stellar magnitude of 13 mmag corresponds to a difference in flux of 1.4%. This flux difference corresponds to variability in the star spot coverage up to   6.3% for photosphere–star spot temperature contrasts of 200 K (${T}_{\mathrm{phot}}=3400$ and ${T}_{\mathrm{spot}}=3200$) and   3.5% for contrasts of 400 K (${T}_{\mathrm{spot}}=3000$).   With these values, we can expect slopes of up to 150 ppm to vary throughout the activity cycle, as seen in Figure 7, if star spots are the main source of stellar variability. If plages are dominating the stellar variability, the slopes induced would be opposite those from star spots (i.e., transit depths would be smaller at shorter wavelengths). Note that, as mentioned above, GJ 436 does not exhibit the expected color variation if star spots dominate the variability, and the interplay between star spots and stellar plages will complicate this picture.

We mention above that the two STIS visits gave statistically consistent spectra, yet appeared to show some correlation in that the long-wavelength half of the visit 2 transit spectrum showed larger transit depths than the long-wavelength half of the visit 1 spectrum (see Section 3.3). If this were due to stellar plages inducing a slope toward smaller transit depths at shorter wavelengths in the spectrum of visit 2, we would expect the visit 1 transit depths to be significantly higher at all wavelengths. Stellar plages only serve to reduce transit depths relative to a homogeneous stellar surface with the same photospheric temperature. Alternatively, if star spots were affecting visit 1, reducing the long-wavelength part of the visit 1 transit spectrum relative to visit 2, we would expect the blue part of the spectrum to diverge even more. In general, if an inhomogenous stellar photosphere were causing any difference between the two STIS visits, we would expect the long-wavelength half of the spectra to agree better relative to the short-wavelength half. We therefore conclude that any difference between visit 1 and visit 2 appears not to be due to stellar variability.

Since we chose to fit the orbital solution to values found in the literature, we compared orbital solutions from Knutson et al. (2014a), Lanotte et al. (2014), and Morello et al. (2015) (see Table 4). Figure 9 shows how the resulting STIS spectra of GJ 436b compare. For each orbital solution, there is a uniform offset. The largest difference is between the Knutson et al. (2014a) and Lanotte et al. (2014) orbital solutions, with an average offset of about 260 ppm. Because the offset is uniform in wavelength, we can be confident that the orbital solution does not affect our nondetection of both a scattering slope and alkali features.

Zoom In Zoom Out Reset image size

Table 4.  Orbital Solutions for GJ 436b

Reference	Inclination (°)	$a/{R}_{* }$	Impact parameter (b)	e	ω
Morello et al. (2015)	86.49 ± 0.12	13.82 ± 0.34	${0.846}_{-0.049}^{+0.050}$	0.16 ± 0.02	351 ± 1.2
Knutson et al. (2014a)	86.774 ± 0.03	14.41 ± 0.10	${0.81}_{-0.012}^{+0.014}$	0.1495 ± 0.016	336 ± 12
Lanotte et al. (2014)	${86.858}_{-0.052}^{+0.049}$	${14.54}_{-0.15}^{+0.14}$	0.7969 ± 0.021	0.1616 ± 0.004	327.2 ± 2.2

Download table as:  ASCIITypeset image

We compare the optical-to-infrared spectrum of GJ 436b to the models described in Morley et al. (2017). These models are self-consistent radiative–convective models in both chemical equilibrium and disequilibrium via quenching. In the models with quenching, the abundances of CH₄, CO, and CO₂ at pressures below 10 bars are fixed to be the equilibrium abundances at 10 bars (i.e., the quench pressure is 10 bars). Interior heating is included from tidal dissipation because of GJ 436b's eccentric orbit ($e\sim 0.15$, Turner et al. 2016). Clouds are modeled through a modified version of the Ackerman & Marley (2001) models (Morley et al. 2012, 2013, 2015). In this cloud model, f_sed represents a sedimentation parameter that determines whether the clouds are thin (${f}_{\mathrm{sed}}\gt \sim 0.3$) or thick (${f}_{\mathrm{sed}}\lt \sim 0.3$). Photochemical hazes are modeled using the results from Line et al. (2011) with the approach from Morley et al. (2013, 2015). This method models the creation of soot precursors (C₂H₂, C₂H₄, C₂H₆, C₄H₂, and HCN) and "converts" them to scattering hazes via a haze efficiency parameter, f_haze (i.e., f_haze is the fraction of soot precursors that actually form hazes; Line et al. 2011).

We allowed for a uniform offset between the STIS spectrum and the WFC3 spectrum in our calculation of the ${\chi }^{2}$ values. The original white light curve depth found with the Knutson et al. (2014a) orbital solution was 7000 ppm. This is above the average transit depth of the STIS spectrum of around 6700 ppm using the Morello et al. (2015) orbit. Since the overall offset between the modeled optical and infrared spectrum is small, we choose to implement a simple offset. Because the WFC3 spectrum was analyzed using a template fitting technique, the differential depths within the spectrum would not change given a uniform offset.

We find that the best-fit model from Morley et al. (2017; 1000× metallicity, T_int = 240 K, f_sed = 0.3 salt/sulfide clouds, and disequilibrium chemistry through quenching of CH₄, CO, and CO₂) is consistent with our new HST/STIS observations. This model has a moderately thin cloud layer of salt/sulfide grains. A lower f_sed (∼0.1) can induce a slope in the transmission spectrum toward lower transit depths at shorter wavelengths, but our observations do not have the sensitivity to constrain this well. As in Morley et al. (2017), clear high-metallicity models (∼1000 × solar) can explain the full transmission spectrum, as can models with lower metallicity but with some clouds (e.g., ∼100× solar metallicity +0.3 f_sed). Since both scenarios are relatively flat in the optical, they are both consistent with our observations. We must look to future observations in the infrared to determine which scenario is taking place on GJ 436b (see Section 5.2).

Our observations help rule out most hazes with low f_haze efficiency.   The observed transit depths around 1.4 μm are inconsistent with models having large mean haze particle radii and low f_hazevalues. Our new STIS observations disfavor smaller mean particle radius hazes as well, indicated by the lack of a strong scattering slope in the STIS bandpass (see Figure 7). Some haze models are still consistent with the data but require a large f_haze value to sufficiently flatten the transit spectrum.

Only about a dozen sub-Jovian exoplanets have been characterized spectroscopically in transmission. While the smallest exoplanets have either flat spectra or very tentative detections in their transmission spectra (Knutson et al. 2014b; Tsiaras et al. 2016), a few exoplanets closer to the size of the solar system's ice giants have detectable water absorption features in their near-IR spectrum with WFC3.  HAT-P-11b and HAT-P-26b have a relatively strong detection of water at 1.4 μm (Fraine et al. 2014; Wakeford et al. 2017), while observations of GJ 436b have not yet revealed a water absorption feature (Knutson et al. 2014a). Stevenson et al. (2016) also interpreted ground-based observations of HAT-P-26b over 0.7 to 1.0 μm to show evidence of water. Differences in the amplitude and existence of the 1.4 μm water absorption feature could point to variations in either the water abundance or aerosol properties. Crossfield & Kreidberg (2017) provide evidence for a correlation between the amplitude of the water absorption feature and either the equilibrium temperature or hydrogen mass fraction. To explain this correlation, they suggest that low-temperature planets may have more optically thick and photochemically produced hazes or that smaller planets have higher metallicity and thus reduced scale heights.

Recently, Wakeford et al. (2017) constrained the metallicity of the Neptune-sized HAT-P-26b to be ${4.8}_{-4.0}^{+21.5}\times$ solar using the retrieved water abundance from the 1.4 μm feature. This is very different from the retrieved metallicity for GJ 436b from Morley et al. (2016), who found the metallicity to be greater than 106× solar using both the featureless IR transmission spectrum and the day-side spectrum. Large differences in metallicities between these planets may point to differences in their formation, perhaps in formation location or disk properties. In agreement with modeling predictions and population studies, the characterization of individual sub-Jovian exoplanets already indicates that this population is diverse in composition (Fortney et al. 2013; Moses et al. 2013; Venturini et al. 2016; Wolfgang et al. 2016).

  GJ 436 and HAT-P-26b both have STIS spectra that are interpreted as being cloudy and featureless; however, the STIS spectra are strikingly similar, especially regarding an abrupt change in transit depth at about 0.8 μm. This jump is also seen in ground-based data of HAT-P-26b taken at Magellan with LDSS3-C (Stevenson et al. 2016). The transit depth of HAT-P-26b near 0.9 μm is confirmed by WFC/G102 observations from Wakeford et al. (2017). Additionally, ground-based data from Rackham et al. (2017) also found a similar slope at comparable wavelengths, this time in GJ 1214b at Magellan with IMACS. GJ 1214b's optical spectrum is well fit by a model for the effects of unocculted stellar plages on the photosphere of the host star covering 3.2% of the stellar disk with a temperature contrast of ∼350 K. We explore this explanation for the other planets below. Table 5 lists the amplitude, number of transits, host star type, and references for the sub-Jovian planets that show this jump in transit depth. In total, 8 transits of 3 sub-Jovian planets observed from four different instruments from two facilities have shown this increase in transit depth at 0.8 μm.

Table 5.  Comparison of 7000–8000 Å Feature in Warm Neptunes

Planet	Approx. depth of feature (ppm)	Approx. depth of feature ($\delta {R}_{p}/H$)	Number of transits	Host star type	Reference
GJ 436b	250	5.5	2	M2.5V	This work
GJ 1214b	450	2	3	M4.5V	Rackham et al. (2017)
HAT-P-26b	500	5	3	K0V	Stevenson et al. (2016), Wakeford et al. (2017)

Download table as:  ASCIITypeset image

We quantify the significance of this jump in transit depth by comparing ${\chi }^{2}$ fits to a flat line with and without including the data between 0.7 and 0.9 μm. The ${\chi }^{2}$ between a flat line and the GJ 436b and HAT-P-26b STIS data is 55.5 using 17 data points, giving a reduced ${\chi }^{2}$ of 3.47, indicating that the probability that our data could have come from a flat line is less than 0.05%. If we remove all data points between 0.7 and 0.9 μm, we get a ${\chi }^{2}$ of 8.73 using 11 data points, resulting in a reduced ${\chi }^{2}$ of 0.87, indicating that the probability that this data could have come from a flat line is greater than 50%. This simple analysis serves to show that, when taken together, the STIS spectra of these planets are inconsistent with a flat line, primarily due to a jump in transit depth between 0.7 and 0.9 μm. In the event that the error bars are underestimated, we find that they would need to be increased by 35% for the STIS data to be consistent with a flat line. We note, however, that when we include the ground-based and WFC3 data, we find that a flat line is inconsistent with the data, even when excluding the 0.7–0.9 μm region. This makes sense because the GJ 1214b spectrum is far from flat (stellar plages are likely making short-wavelength transit depths lower than long-wavelength transit depths), and HAT-P-26b may have water absorption around 0.9 microns, moving the ground-based and WFC3 data away from a flat line (though this water absorption does not explain the jump in transit depth between 0.7 and 0.9 μm).

The explanation of stellar plages is plausible for GJ 436b. The host star, GJ 436, is an M dwarf like GJ 1214, and the maximum flux difference between the photosphere and a stellar plage with a temperature contrast of 350 K is around 0.75 μm, where we see the jump. In Figure 10, we plot two models of the effect of stellar plages using Equation (11) of Rackham et al. (2017):

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{p}}{{R}_{s}}\right)}_{\lambda ,\mathrm{obs}}=\sqrt{1-\displaystyle \frac{(1-f-{D}_{\lambda }){S}_{o}+{{fS}}_{u}}{(1-f){S}_{o}+{{fS}}_{u}}},\end{eqnarray} \tag{ 3 }$$

where f is the fraction of the stellar disk with stellar plages, ${D}_{\lambda }$ is the transit depth expected if the stellar disk had no stellar plages (i.e., the "true" transit depth), S_u is the flux from a stellar plage, and S_o is the flux typical of the greater photosphere.   We use PHOENIX stellar atmosphere models to model the spectra of the photospheres and the hotter stellar plages.

Zoom In Zoom Out Reset image size

In general, the stellar plage models do not match the jump at 0.8 μm well. Longward of 0.8 μm, transit depth decreases, unlike expectations from the plage models. Similarly, shortward of 0.7 μm, the transit depth does not stay low as in the plage models. The explanation of stellar plages also is not a sufficient explanation for HAT-P-26b, whose host star is of type K0. For temperature contrasts of around 300 K, the effects of stellar plages for a K0 photosphere are flat until about 6500 Å, outside the range of the jump in transit depth. Additionally, the increase in transit depth is rather sudden and does not qualitatively match the trend seen in GJ 1214b that is explained by stellar plages.  

Another explanation could be biases from fixing the LDCs while fitting the transit light curve. However, Stevenson et al. (2016) tested fitting LDCs and found the results to be consistent with LDCs from models. Additionally, Rackham et al. (2017) left one of the quadratic LDCs as a free parameter, further verifying that GJ 1214b's increase in transit depth is not related to biases in limb darkening. This does remain an explanation worth further investigation for the other planets.

A third explanation could be an additional opacity source around 0.8 μm in the atmosphere of these planets. However, in our grid, no model seems to fit this part of the spectrum accurately, and all of the best-fitting models are featureless throughout STIS wavelengths. Any additional opacity source at 0.8 μm would also need to   make sense in the context of the rest of the spectrum, which is either featureless or only shows water absorption. Additionally, the size of such a feature would be unprecedented in sub-Jovian exoplanets. The jump in transit depth at 0.8 μm in HAT-P-26b is about the size of its water absorption feature at 1.4 μm. As Table 5 shows, the jump in transit depth spans up to 5.5 scale heights. Thus we recommend caution before interpreting the 0.8 μm feature as an opacity source. We also note that this effect is not readily apparent in the optical spectrum of hot Jupiters (e.g., Sing et al. 2016).

The wide wavelength coverage and superior light-gathering capability of the James Webb Space Telescope will allow planets like GJ 436b to be studied in more detail. As mentioned above, current JWST GTO plans include GJ 436b eclipse observations in GTO programs of three instruments. These observations will tell us a great deal about the composition and structure of the day-side atmosphere; however, observations of GJ 436b's transmission spectrum are necessary to probe a different region of the planet, namely the terminator. To what degree the day-side atmosphere's composition and structure differ from that of the terminator has important implications for our understanding of the global circulation taking place in these atmospheres (Kataria et al. 2016).

We simulated JWST data using PandExo for three transits observed with the Near Infrared Spectrograph (NIRSpec) G395H grism Bright Object Time Series and Mid-infrared Instrument (MIRI) slitless Low Resolution Spectroscopy with a noise floor of 10 ppm (Batalha et al. 2017). We avoided shorter-wavelength instrument modes that may have trouble with saturation and nonlinearity effects due to the brightness of GJ 436.   The simulated data are shown in Figure 11 with 200 and 10 pixels per bin for NIRSpec and MIRI, respectively. This is comparable to a resolution of R ∼ 20–40 for NIRSpec and R ∼ 10–55 for MIRI. We include four model scenarios from Morley et al. (2017) and described them in Section 4.3 for GJ 436b's transmission spectrum at JWST wavelengths. For the simulated data, we adopted the 100× metallicity scenario with internal heating but no quenching of CH₄, CO, and CO₂ in order to demonstrate JWST's ability to distinguish between model scenarios.

Zoom In Zoom Out Reset image size

All four model scenarios in Figure 11 are consistent with current measurements. However, JWST transit observations can distinguish cloudy and high-metallicity scenarios, quenched and nonquenched scenarios, and different internal heating. The data-point uncertainties approach the imposed noise floor of 10 ppm, but since the spectral features in the transit spectrum have amplitudes much larger than the noise floor, each of the models could be distinguished even if the noise floor were increased.   Morley et al. (2017) found that quenched carbon chemistry with internal heating provided the best fit to the day-side thermal emission spectrum. With transit spectra from JWST, we can test this conclusion at the planet's terminator. Additionally, JWST will be capable of resolving the spectral features of CH₄, CO, and CO₂ in both transmission and thermal emission spectra, a big step forward from the Spitzer photometry used today to infer the presence of these molecules.

It is worth noting that JWST's short-wavelength capabilities end at about 0.7 μm, meaning that space-based optical and UV observations with STIS and ground-based optical and near-UV observations (e.g., the Arizona-CfA-Cátolica Exoplanet Spectroscopy Survey, or ACCESS, and the GTC exoplanet transit spectroscopy survey) will remain a useful and unique way to characterize these exoplanets. However, as we have seen, properly interpreting the optical transmission of an exoplanet requires an accurate understanding of effects from the host star.