ACCESS I. AN OPTICAL TRANSMISSION SPECTRUM OF GJ 1214b REVEALS A HETEROGENEOUS STELLAR PHOTOSPHERE

Following the framework detailed by Jordán et al. (2013), we modeled the observed target light curve $l(t)$ as

$$\begin{eqnarray}&&l(t)={FT}(t,\theta )S(t)\epsilon (t),\end{eqnarray} \tag{ 2 }$$

in which t is time, F is the underlying flux from the target star, $T(t,\theta )$ is the transit signal defined by the vector of transit parameters θ, $S(t)$ is the perturbation signal due to the combination of all systematic variations, and $\epsilon (t)$ is the stochastic noise component. We account for both uncorrelated variations ("white noise") and correlated variations ("red noise") in the light curve with the $\epsilon (t)$ term and modeled them following the wavelet-based method of Carter & Winn (2009)¹⁶ (see Section 4.2.1). We assumed the underlying flux from the target star to be constant, so any actual variations in the star's flux were therefore encapsulated by the noise component. We assumed the perturbation signal to be a linear combination of signals owing to different instrumental and atmospheric effects, which can be represented mathematically as

$$\begin{eqnarray}&&S(t)=\displaystyle \prod _{i=0}^{n}{\alpha }_{i}{s}_{i}(t),\end{eqnarray} \tag{ 3 }$$

in which ${s}_{i}(t)$ represents the different signals and ${\alpha }_{i}$ their respective scaling coefficients. Expressed in logarithmic space, each of these multiplicative signals become additive, so defining $L\equiv \mathrm{log}l(t)$, the base-10 logarithm of the observed target light curve can be written as

$$\begin{eqnarray}&&L=\mathrm{log}F+\mathrm{log}T(t,\theta )+\displaystyle \sum _{i=0}^{n}{\alpha }_{i}{s}_{i}(t)+\mathrm{log}\epsilon (t).\end{eqnarray} \tag{ 4 }$$

The systematic variations represented by the perturbation signal include common variations experienced by the target and comparison spectra due to instrumental and atmospheric effects. Therefore, within this framework, the observed signal from each k comparison star can be modeled as

$$\begin{eqnarray}&&{L}_{k}=\mathrm{log}{F}_{k}+\displaystyle \sum _{i=0}^{n}{\alpha }_{i,k}{s}_{i}(t)+\mathrm{log}{\epsilon }_{k}(t),\end{eqnarray} \tag{ 5 }$$

in which F_k is the baseline flux from the comparison star, ${\alpha }_{k}$ is the set of unique scaling coefficients for the systematic variations as they apply to this comparison star, and ${\epsilon }_{k}(t)$ is the unique noise signal associated with this light curve. Thus the mean-subtracted light curve from each comparison star (${L}_{k}-\mathrm{log}{F}_{k}$) can be used to estimate the perturbation signal $S(t)$, which can then be subtracted from L in Equation (4).

We performed a principal component analysis (PCA) of the mean-subtracted comparison light curves to estimate the independent signals ${s}_{i}(t)$ comprising the perturbation signal $S(t)$. With N comparison stars, we could estimate at most N principal components via PCA. As detailed in Section 4.2.1, we let the scaling coefficients ${\alpha }_{i}$ and the baseline flux from the host star F float as free parameters in the Markov chain Monte Carlo (MCMC) procedure, so we enforced the maximum number of principal components for use in the model to be $M\leqslant N-1$. However, we determined the optimal number of components to be the minimum number that could achieve the same predictive power as the best set of M available components.

We estimated the predictive power of each set of components using a k-fold cross-validation procedure. First, we split the out-of-transit target light curve into 20 segments (or "folds"). For each number of available principal components (from 1 to M), we fit the α scaling coefficients for the available components to the target light curve using 19 "training" folds, and recorded the error (in units of normalized flux) between the target light curve and fit for the single "validation" fold. Repeating this process using each of the 20 folds in turn as a validation fold gave us an estimate and confidence interval of the prediction error for each number of principal components. We identified the number of principal components that produced the smallest median prediction error ("the best set"), and then used the lowest number of principal components with a median prediction error indistinguishable at the $1\sigma$ level from that of the best set.

The final consideration in the PCA-based procedure was to select the comparison stars to ultimately use in generating the principal components. Each comparison light curve provides a noisy estimate of the perturbation signal, so it follows that the combination of comparison stars can be optimized to provide the best estimate with the least noise. Through successive iterations, we found that for each night the brightest comparison stars produced principal components that were most capable of accurately predicting perturbations in the target light curve. Ultimately, we used the brightest 4, 5, and 4 comparison stars to estimate the perturbation signals for Transits 1, 2, and 3, respectively. We found these comparison stars to also provide the best results when utilizing the polynomial-based modeling framework we explored.

We also modeled the observed target light curve following the procedure employed by Bean et al. (2010) in their successful VLT/FORS observations of GJ 1214b. This framework is empirically motivated, as the data show that simply dividing the target light curve by the sum of the comparison light curves removes most of the variations in the out-of-transit flux from the target, leaving only a smoothly varying long-term trend. Like Bean et al. (2010), we find that this trend can be modeled well as a second order polynomial function of time, and attribute it to the color difference between this very red target star and the available comparison stars in the field. In this framework, the target light curve can be expressed as

$$\begin{eqnarray}&&l(t)={FT}(t,\theta )S(t)\displaystyle \sum _{i=0}^{2}{\alpha }_{i}{t}^{i}\epsilon (t),\end{eqnarray} \tag{ 6 }$$

and the comparison light curves as

$$\begin{eqnarray}&&{l}_{k}(t)={F}_{k}S(t){\epsilon }_{k}(t),\end{eqnarray} \tag{ 7 }$$

in which ${\alpha }_{i}$ now describes the polynomial coefficients and the other terms have the same meaning as indicated previously. Dividing by the sum of comparison light curves, the detrended target light curve can be expressed as

$$\begin{eqnarray}&&{l}_{{\det }}(t)=T(t,\theta )\displaystyle \sum _{i=0}^{2}{\alpha }_{i}{t}^{i}{\epsilon }_{c}(t),\end{eqnarray} \tag{ 8 }$$

in which the noise term ${\epsilon }_{c}(t)$ now represents the combined noise from the target and comparison light curves, and the constant F and F_k terms have been subsumed into the ${\alpha }_{0}$ coefficient. We note that this formalism is an approximation, since $S(t)$ will be different for the target and reference stars in real data (as expressed by the ${\alpha }_{i,k}$ coefficients in Equation (5)) and will not divide out exactly.

Both modeling frameworks are tested in the following. We compare their effectiveness in Section 5.1 and report the results of the polynominal-based detrending procedure in this paper.

For both modeling frameworks, the transit parameters were estimated for each night using an MCMC optimization procedure. The transit signal was modeled following the formalism of Mandel & Agol (2002), which accounts for the effect of limb darkening via a quadratic law of the form

$$\begin{eqnarray}&&I(\mu )=I(1)[1-{u}_{1}(1-\mu )-{u}_{2}{(1-\mu )}^{2}],\end{eqnarray} \tag{ 9 }$$

in which μ is the cosine of the angle between the stellar surface normal and the line of sight to the observer, and u₁ and u₂ are the limb darkening coefficients. Fixing the limb darkening coefficients has been shown to bias measurements of the transit depth (Espinoza & Jordán 2015). Yet, when left as free parameters, these coefficients are strongly correlated in MCMC retrievals (Pál 2008; Kipping 2013). However, following a rotation onto new principal axes

$$\begin{eqnarray}&&{\omega }_{1}={u}_{1}\cos \phi -{u}_{2}\sin \phi ,\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{\omega }_{2}={u}_{2}\cos \phi +{u}_{1}\sin \phi \end{eqnarray} \tag{ 10b }$$

with $\phi =35\buildrel{\circ}\over{.} 8$ (Kipping 2013), ${\omega }_{1}$ and ${\omega }_{2}$ are essentially uncorrelated, and the first parameter can account for variations induced by the transit geometry, while the second remains constant (Howarth 2011). Therefore we used these rotated coefficients to describe the effect of limb darkening, leaving ${\omega }_{1}$ as a free parameter and fixing ${\omega }_{2}$ to values obtained from a PHOENIX atmospheric model (Husser et al. 2013) for $\mu \geqslant 0.1$, with stellar parameters closest to those identified by Rojas-Ayala et al. (2012): ${T}_{\mathrm{eff}}=3300\,$ K, $\mathrm{log}g=5.0$, and $[M/H]=0.0$. We used a uniform prior on ${\omega }_{1}$, with boundary values set by the triangular sampling method of Kipping (2013).

In each MCMC procedure, we adopted the fixed parameters on the transit model, which included system scale ($a/{R}_{s}$), inclination (i), orbital period (P), eccentricity (e), and argument of periastron (ω), from Kreidberg et al. (2014) to allow for direct comparisons of results. The planet-to-star radius ratio (${R}_{p}/{R}_{s}$), rotated limb darkening coefficient (${\omega }_{1}$), scaling coefficients for principal components or polynomial terms (α), and photometric uncertainty ($\epsilon$) were left as free parameters. The baseline flux (F) was also left as a free parameter in the PCA-based procedure. We placed a Gaussian prior on ${R}_{p}/{R}_{s}$, centered on the median transit depth reported by Kreidberg et al. (2014), with an uncertainty of ${\sigma }_{{R}_{p}/{R}_{s}}=0.01$ to allow the algorithm to thoroughly explore the parameter space, and truncated by the range [0, 1]. We also placed $5\sigma$ Gaussian priors on each α coefficient. In the PCA-based procedure, we determined the mean and standard deviation for the priors using the distributions of the α coefficients obtained for the comparison stars. In the polynomial-based procedure, we performed a bootstrap analysis on the out-of-transit data to determine the mean and standard deviation for the priors. We resampled the out-of-transit data with replacement 1000 times and fit a second order polynomial function of time to each sample via least-squares, recording the fit coefficients. We utilized the mean and standard deviation of the coefficient distributions to determine the $5\sigma$ Gaussian priors for the MCMC. Finally, we placed a $5\sigma$ Gaussian prior on F in the PCA-based procedure using the median and standard error of the out-of-transit flux.

We ran five chains of 130,000 steps and discarded the first 30,000 steps as the burn-in. For each step in the chain, the likelihood of the residuals given the model was calculated via the same likelihood function given by Equation (41) of Carter & Winn (2009), which parameterizes the contributions of uncorrelated and time-correlated noise sources to the observed light curve via the parameters ${\sigma }_{w}$ and ${\sigma }_{r}$, respectively. We placed uniform priors bounded by the interval [0, 1] on both of these parameters.

We combined the results from all chains to determine the posterior values. We evaluated convergence between chains using the Gelman–Rubin statistic $\hat{R}$ (Gelman & Rubin 1992), and considered the chains to be well-mixed if $\hat{R}\leqslant 1.03$ for all parameters. We constructed the posterior distribution by sampling each chain at intervals spaced by $10\times$ the half-life of the autocorrelation in the chain. This wide spacing ensured independent sampling. We adopted the median and 68.27% confidence interval defined by the 15.87th and 84.14th percentiles of the posterior distribution as the final posterior value and uncertainty for each parameter.

We first analyzed the white light curve, using flux from the complete spectra (4500–9260 Å) of GJ 1214 and the comparison stars. We performed the procedure as outlined previously, additionally including the time of mid-transit (t₀) as a free parameter. We placed a uniform prior on t₀ spanning the complete observation.

Following the white light analysis, we repeated the MCMC procedure for each wavelength bin. We allowed the same parameters to float, except t₀, which we fixed to the value obtained from the white light analysis.

In determining the width of wavelength bins for the spectroscopic analysis, we aimed to (1) maximize the number of bins, while (2) collecting enough signal in the bins to discriminate between model transmission spectra with signals on the order of ${\rm{\Delta }}({R}_{p}/{R}_{s})\sim 0.005$. Taking into account our usable wavelength coverage of 4500–9260 Å and the chip gaps for the f/4 and f/2 setups, we determined that dividing the spectra into 14 bins with roughly even signal to be optimal: we use 10 bins with a mean width of $\,200$ Å (174–243 Å) at wavelengths greater than $7250\,\mathring{\rm A}$, where the spectra are brightest, and four bins with a mean width of 616 Å (310–1157 Å) at shorter wavelengths, where the spectra are considerably fainter (Table 2).

Table 2.  Planet-to-Star Radius Ratios

Bin			${R}_{p}/{R}_{s}$
N	${\lambda }_{\min }$	${\lambda }_{\max }$	Transit 1	Transit 2	Transit 3	Final
1	4500.0	5657.0	${0.1121}_{-0.0020}^{+0.0020}$	${0.1129}_{-0.0021}^{+0.0021}$	${0.1172}_{-0.0027}^{+0.0026}$	${0.1136}_{-0.0013}^{+0.0013}$
2	5657.0	6196.0	${0.1119}_{-0.0018}^{+0.0017}$	${0.1139}_{-0.0018}^{+0.0017}$	${0.1140}_{-0.0021}^{+0.0021}$	${0.1132}_{-0.0011}^{+0.0011}$
3	6196.0	6506.0	${0.1153}_{-0.0015}^{+0.0015}$	${0.1125}_{-0.0025}^{+0.0024}$	${0.1149}_{-0.0013}^{+0.0013}$	${0.1147}_{-0.0009}^{+0.0009}$
4	6596.0	7054.0	${0.1112}_{-0.0015}^{+0.0015}$	${0.1153}_{-0.0025}^{+0.0024}$	${0.1157}_{-0.0019}^{+0.0019}$	${0.1133}_{-0.0010}^{+0.0011}$
5	7254.0	7485.0	${0.1147}_{-0.0012}^{+0.0012}$	${0.1151}_{-0.0016}^{+0.0017}$	${0.1138}_{-0.0008}^{+0.0008}$	${0.1142}_{-0.0006}^{+0.0006}$
6	7485.0	7728.0	${0.1130}_{-0.0013}^{+0.0013}$	${0.1161}_{-0.0013}^{+0.0014}$	${0.1149}_{-0.0012}^{+0.0011}$	${0.1146}_{-0.0007}^{+0.0007}$
7	7728.0	7967.0	${0.1155}_{-0.0013}^{+0.0013}$	${0.1154}_{-0.0011}^{+0.0011}$	${0.1138}_{-0.0011}^{+0.0012}$	${0.1148}_{-0.0007}^{+0.0007}$
8	7967.0	8142.0	${0.1157}_{-0.0016}^{+0.0015}$	${0.1174}_{-0.0012}^{+0.0013}$	${0.1138}_{-0.0010}^{+0.0011}$	${0.1153}_{-0.0007}^{+0.0007}$
9	8142.0	8316.0	${0.1150}_{-0.0019}^{+0.0018}$	${0.1172}_{-0.0019}^{+0.0019}$	${0.1146}_{-0.0010}^{+0.0010}$	${0.1152}_{-0.0008}^{+0.0008}$
10	8316.0	8502.0	${0.1133}_{-0.0013}^{+0.0013}$	${0.1162}_{-0.0016}^{+0.0016}$	${0.1153}_{-0.0011}^{+0.0011}$	${0.1148}_{-0.0008}^{+0.0007}$
11	8502.0	8692.0	${0.1150}_{-0.0014}^{+0.0014}$	${0.1159}_{-0.0016}^{+0.0015}$	${0.1168}_{-0.0014}^{+0.0013}$	${0.1159}_{-0.0008}^{+0.0008}$
12	8692.0	8868.0	${0.1136}_{-0.0014}^{+0.0013}$	${0.1161}_{-0.0018}^{+0.0018}$	${0.1128}_{-0.0013}^{+0.0014}$	${0.1138}_{-0.0008}^{+0.0008}$
13	8868.0	9065.0	${0.1164}_{-0.0017}^{+0.0017}$	${0.1161}_{-0.0017}^{+0.0017}$	${0.1141}_{-0.0017}^{+0.0016}$	${0.1155}_{-0.0010}^{+0.0010}$
14	9065.0	9260.0	${0.1165}_{-0.0023}^{+0.0022}$	${0.1153}_{-0.0016}^{+0.0016}$	${0.1147}_{-0.0014}^{+0.0014}$	${0.1153}_{-0.0010}^{+0.0009}$

Note. For individual transits, we report the medians and 68% confidence intervals on the posterior distributions from the MCMC optimization procedure. The final radius ratios are the weighted means of the three transit measurements (see Section 5.4 for details).

Download table as:  ASCIITypeset image

Along with the limb darkening prescription described previously, we investigated two additional methods for handling the effect of stellar limb darkening on the transit light curve. In the first approach we fixed the limb darkening coefficients ${\omega }_{1}$ and ${\omega }_{2}$ during the MCMC analysis to the values from the PHOENIX atmospheric model (Husser et al. 2013). In the second we utilized the triangular sampling method proposed by Kipping (2013), which involves fitting for new limb darkening coefficients ${q}_{1}={({u}_{1}+{u}_{2})}^{2}$ and ${q}_{2}=0.5{u}_{1}{({u}_{1}+{u}_{2})}^{-1}$ with uniform priors on the interval [0,1] in order to efficiently sample the physically bounded parameter space. The ${R}_{p}/{R}_{s}$ values for the final transmission spectrum did not differ significantly among the three limb darkening prescriptions. With respect to the nominal results, adopting fixed limb darkening coefficients caused the ${R}_{p}/{R}_{s}$ values on average to differ by $5\times {10}^{-4}$ ($0.6\sigma$). Similarly, the triangular sampling method led to ${R}_{p}/{R}_{s}$ values that differed by $3\times {10}^{-4}$ ($0.3\sigma$). Both methods enlarged the $1\sigma$ uncertainties on the transmission spectrum by $0.2 \%$. In this paper we report the ${R}_{p}/{R}_{s}$ values determined via our nominal limb darkening prescription described in Section 4.2.1.

We investigate how heterogeneities across the star's photosphere could affect the observed transmission spectrum of the exoplanet using a basic model to incorporate the effects of a composite photosphere and atmospheric transmission along the planet's limb (hereafter the "CPAT model"). We consider the simplest case in which the emergent spectrum of the star is composed of two distinct components: the spectrum typical of the occulted transit chord S_o (the "occulted" spectrum) and the unocculted spectrum S_u, which is fully outside of the transit chord (Figure 10). It is important to note that the occulted spectrum includes the transit chord but is not limited to it. For example, the planet could transit a region with a spectrum that is typical of $80 \%$ of the photosphere, while another $20 \%$ of the photosphere produces a distinct spectrum that is not probed by the transit chord. Additionally, the unocculted region does not need to be continuous in this model. We note that the planet could transit multiple regions with distinct spectral characteristics (as is the case while crossing starspots; e.g., Sanchis-Ojeda & Winn 2011), or multiple unocculted regions with distinct spectra could exist, though we show here that this simple model is sufficient to reproduce the observations.

Zoom In Zoom Out Reset image size

In cases where the stellar photosphere is composed of two distinct components, the observed transmission spectrum given by the CPAT model is

$$\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{ 11 }$$

in which ${({R}_{p}/{R}_{s})}_{\lambda ,\mathrm{obs}}$ is the observed, wavelength-dependent planet-to-star radius ratio; F is the fraction of the stellar disk covered by the unocculted spectrum S_u; and ${D}_{\lambda }=({R}_{p}/{R}_{s}{)}_{\lambda }^{2}$ is the transit depth expected from the true, wavelength-dependent planet-to-star radius ratio. The numerator within the square root gives the in-transit flux, and the denominator the out-of-transit flux. In the following sections, we apply the CPAT model to three cases: (1) heterogeneous absorbers, (2) generalized temperature heterogeneities, and (3) cool starspots with parameters fixed to values inferred from long-term photometric monitoring of GJ 1214b.

Significant chemical heterogeneities are known to exist for magnetically active hot (B and A-type) main sequence stars that show peculiar chemical abundances (e.g., Pyper 1969; Khokhlova 1985). While no strong correlation was observed between the line profile variations and magnetic field strengths, it has been proposed that the chemical abundance patterns emerge due to anisotropic diffusion of the elements in a strong magnetic field (e.g., Michaud 1970; Urpin 2016). Simultaneous Doppler imaging mapping of chemical heterogeneities and magnetic field geometry argue for highly complex configurations across the stellar disk (e.g., Piskunov & Kochukhov 2002). Similarly, both partially and fully convective mid-M dwarfs have been found to store the bulk of their magnetic flux in small scale components that are non-axisymmetric (Reiners & Basri 2009). Field strengths up to $4\times {10}^{3}$ G have been detected for very active M4.5 dwarfs (Johns-Krull & Valenti 1996) as well as rapidly rotating mid- to late-M dwarfs (Reiners & Basri 2010; Reiners 2012). If indeed anisotropic diffusion is the process that leads to chemically heterogeneous stellar photospheres, then it is expected that this effect will only be important for stars with strong magnetic fields ($\sim {10}^{3}$ G) and without fully convective atmospheres; however, weaker magnetic fields—not able to produce strong enough chemical heterogeneities to lead to varying absorption line profiles—may lead to low-level heterogeneities that could still potentially influence the optical transmission spectrum of transiting exoplanets. The detailed analysis of this effect is beyond the scope of this paper but may be important for future transit spectroscopy of planets orbiting stars with strong magnetic fields, including rapidly rotating mid- to late-M dwarfs.

We use Equation (11) to model the effects of a heterogeneous distribution of absorbers in the stellar photosphere on the observed transmission spectrum (the "CPAT-absorber model"). We utilize PHOENIX model stellar spectra (Husser et al. 2013) with ${T}_{\mathrm{eff}}=3300$ K and $\mathrm{log}g=5.0$ to generate the S_o and S_u spectra for our model. As a proxy for the strength of absorbers, we employ models with a range of metallicities as defined by [Fe/H]. Higher metallicity models demonstrate deeper spectral features due to the larger absorber abundances, while lower metallicity models possess relatively muted absorption features. We do not consider alpha-enhanced or depleted models, so [Fe/H] is synonymous with the overall metallicity Z. The PHOENIX model grid includes spectra for metallicities of $-4.0\leqslant [\mathrm{Fe}/{\rm{H}}]\leqslant -2.0$ in steps of 1.0 and $-2.0\leqslant [\mathrm{Fe}/{\rm{H}}]\leqslant +1.0$ in steps of 0.5. We include each of these models in our analysis.

We employ an MCMC approach to find the best-fitting CPAT-absorber model. We conduct an MCMC optimization for each of the exoplanet atmosphere models in the grids of cloudy and hazy models described in Section 6.2. At each step in the Markov chain, the measured ${R}_{p}/{R}_{s}$ values from the joint data set are compared with the CPAT-absorber model produced by the combination of the heterogeneous stellar photosphere and the input exoplanet atmosphere. The spectra S_o and S_u are generated using the values for ${[\mathrm{Fe}/{\rm{H}}]}_{o}$ and ${[\mathrm{Fe}/{\rm{H}}]}_{u}$ (the metallicities of the occulted and unocculted spectra, respectively) and linearly interpolating between the closest spectra in the PHOENIX model grid. Following the results of Doppler imaging studies of M dwarfs, we adopt $F=0.032$, which is the mean spot filling factor Barnes et al. (2015) found for the M4.5 dwarf GJ 791.2A over its rotation period. We fix the metallicity of the occulted spectrum to the best-fit value for the metallicity of GJ 1214, ${[\mathrm{Fe}/{\rm{H}}]}_{o}=0.20$ (Rojas-Ayala et al. 2012). The free parameters in the model are the metallicity contrast ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$, which determines the metallicity of the unocculted spectrum relative to that of the occulted spectrum, and the uniform offset applied to the exoplanet's model transmission spectrum, ${({R}_{p}/{R}_{s})}_{o}$. We place a uniform prior on ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$ to allow the algorithm to fully explore the parameter space, with interval boundaries $[-4.2,+0.8]$ to keep ${[\mathrm{Fe}/{\rm{H}}]}_{u}$ within the PHOENIX model range $[-4.0,+1.0]$. While allowing the metallicity of 3.2% of the stellar disk to vary will change the mean metallicity of the photosphere, we note that the mean will never vary more than $1\sigma$ from the measured metallicity, $[\mathrm{Fe}/{\rm{H}}]=0.20\pm 0.17$ (Rojas-Ayala et al. 2012). We place a Gaussian prior on ${({R}_{p}/{R}_{s})}_{o}$ using the mean and standard deviation of the residuals found by fitting the joint data set to the exoplanet's model transmission spectrum (Section 6.2) and bounded on the interval $[-1,+1]$. We run five chains of 10⁵ steps with an additional 10⁴ steps discarded as the burn-in. We consider the chains to be well-mixed if the Gelman–Rubin statistic $\hat{R}\leqslant 1.03$ for all parameters.

The results of the CPAT-absorber model fitting for the full grids of cloudy and hazy transmission spectra are illustrated in Figure 11. As can be seen by comparing Figures 8 and 11, the best-fitting parameters for the exoplanet's atmosphere are not substantially changed by incorporating the effect of a heterogeneous stellar photosphere, though the CPAT-absorber models can provide better fits to the data. The complete results from the fitting procedure are provided in Tables 6 and 7, for the equilibrium cloud and photochemical haze models, respectively. The values for the free parameters, ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$ and ${({R}_{p}/{R}_{s})}_{o}$, that we report there and quote as follows are the median and 68% confidence intervals from the MCMC optimization procedure. In the case of equilibrium clouds, we find the best-fitting model to have a very high metallicity and thick clouds in the exoplanet's atmosphere and a relatively low metallicity in the unocculted region of the star's photosphere ($1000\times$ solar, ${f}_{\mathrm{sed}}=0.1$, ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]=-{1.58}_{-0.37}^{+0.28}$, ${({R}_{p}/{R}_{s})}_{o}={8594}_{-26}^{+28}$ ppm, ${\chi }^{2}=54.7$, $\mathrm{dof}=36$, ${\chi }_{\mathrm{red}}^{2}=1.52$). Considering photochemical haze models, we find the best-fitting model to include the same mode particle size and haze-forming efficiency as when considering only the exoplanetary contribution to the transmission spectrum (Section 6.2) and a low metallicity for the unocculted region similar to that found for the CPAT-absorber cloud model (r = 0.3 μm, ${f}_{\mathrm{haze}}=10 \%$, ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]=-{1.31}_{-0.39}^{+0.35}$, ${({R}_{p}/{R}_{s})}_{o}=-{5758}_{-26}^{+27}$ ppm, ${\chi }^{2}=53.3$, $\mathrm{dof}=36$, ${\chi }_{\mathrm{red}}^{2}\ =1.48$). Figure 12 shows the posterior distributions for the free parameters in each of the best-fitting cloud and haze models. A model using a constant, achromatic value of ${R}_{p}/{R}_{s}$ for the exoplanet's transmission spectrum together with a CPAT-absorber model for the star provides a fit to the data comparable to that of the best-fitting CPAT-absorber cloud and haze models (${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]=-{1.02}_{-0.24}^{+0.33}$, ${R}_{p}/{R}_{s}=0.11604\pm 3\times {10}^{-5}$, χ² = 55.8, dof = 36, ${\chi }_{\mathrm{red}}^{2}=1.55$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6.  Results from MCMC Fits to CPAT-absorber Models for Cloudy Atmospheres

Metallicity	fseda	${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$	${({R}_{p}/{R}_{s})}_{o}$	${\chi }^{2}$	dof	${\chi }_{\mathrm{red}}^{2}$
(× Solar)			(ppm)
100	0.01	$-{3.38}_{-0.29}^{+0.55}$	$-{1364}_{-96}^{+155}$	119.7	36	3.32
100	0.1	$-{2.23}_{-0.15}^{+0.14}$	$-{32}_{-35}^{+42}$	354.6	36	9.85
100	1	$-{1.72}_{-0.14}^{+0.12}$	${1956}_{-25}^{+27}$	3914	36	109
100	cloud-free	$-{1.53}_{-0.18}^{+0.11}$	${2392}_{-25}^{+26}$	6863	36	191
150	0.01	$-{2.80}_{-0.15}^{+0.18}$	${601}_{-55}^{+64}$	87.68	36	2.44
150	0.1	$-{2.19}_{-0.17}^{+0.16}$	${1638}_{-33}^{+42}$	227.9	36	6.33
150	1	$-{1.68}_{-0.14}^{+0.13}$	${3288}_{-26}^{+25}$	2580	36	71.7
150	cloud-free	$-{1.46}_{-0.23}^{+0.13}$	${3646}_{-25}^{+25}$	4590	36	127
200	0.01	$-{2.71}_{-0.14}^{+0.15}$	${1990}_{-54}^{+51}$	73.65	36	2.05
200	0.1	$-{2.14}_{-0.19}^{+0.18}$	${2900}_{-31}^{+40}$	163.5	36	4.54
200	1	$-{1.65}_{-0.16}^{+0.14}$	${4340}_{-25}^{+26}$	1834	36	50.9
200	cloud-free	$-{1.42}_{-0.16}^{+0.21}$	${4645}_{-25}^{+25}$	3306	36	91.8
250	0.01	$-{2.65}_{-0.15}^{+0.15}$	${3075}_{-52}^{+52}$	66.64	36	1.85
250	0.1	$-{2.09}_{-0.21}^{+0.19}$	${3895}_{-31}^{+37}$	127.1	36	3.53
250	1	$-{1.63}_{-0.16}^{+0.15}$	${5180}_{-24}^{+27}$	1378	36	38.3
250	cloud-free	$-{1.41}_{-0.16}^{+0.21}$	${5450}_{-25}^{+26}$	2508	36	69.7
300	0.01	$-{2.60}_{-0.15}^{+0.15}$	${3940}_{-52}^{+51}$	62.5	36	1.74
300	0.1	$-{2.03}_{-0.23}^{+0.20}$	${4688}_{-30}^{+36}$	105.6	36	2.93
300	1	$-{1.62}_{-0.18}^{+0.15}$	${5854}_{-26}^{+25}$	1074	36	29.8
300	cloud-free	$-{1.41}_{-0.17}^{+0.21}$	${6098}_{-24}^{+26}$	1981	36	55
1000	0.01	$-{1.98}_{-0.54}^{+0.23}$	${8231}_{-37}^{+54}$	61.43	36	1.71
1000	0.1	$-{1.58}_{-0.37}^{+0.28}$	${8594}_{-26}^{+28}$	54.67	36	1.52
1000	1	$-{1.54}_{-0.24}^{+0.20}$	${9205}_{-24}^{+26}$	212.9	36	5.91
1000	cloud-free	$-{1.39}_{-0.22}^{+0.21}$	${9339}_{-26}^{+25}$	404.9	36	11.2

Note. The free parameters in the fitting procedure, ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$ and ${({R}_{p}/{R}_{s})}_{o}$, are described in Section 6.3.2. We report the median and 68% confidence intervals from the MCMC optimization procedure for each free parameter. With 38 data points and 2 fitted parameters, we assume 36 DOF when calculating the ${\chi }_{\mathrm{red}}^{2}$ statistic.

^aCloud-free models do not have an f_sed value.

Download table as:  ASCIITypeset image

Table 7.  Results from MCMC Fits to CPAT-absorber Models for Hazy Atmospheres

r	fhaze	${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$	${({R}_{p}/{R}_{s})}_{o}$	${\chi }^{2}$	dof	${\chi }_{\mathrm{red}}^{2}$
(μm)	(%)		(ppm)
0.01	1	$-{4.16}_{-0.04}^{+0.02}$	$-{4007}_{-28}^{+28}$	512.4	36	14.2
0.01	3	$-{3.93}_{-0.17}^{+0.15}$	$-{4626}_{-54}^{+49}$	175.6	36	4.88
0.01	10	$-{2.57}_{-0.15}^{+0.15}$	$-{4064}_{-54}^{+52}$	76.34	36	2.12
0.01	30	$-{1.65}_{-0.43}^{+0.41}$	$-{3201}_{-28}^{+31}$	61.98	36	1.72
0.03	1	$-{4.16}_{-0.04}^{+0.03}$	$-{4777}_{-31}^{+27}$	333.8	36	9.27
0.03	3	$-{3.67}_{-0.17}^{+0.18}$	$-{4743}_{-56}^{+53}$	215	36	5.97
0.03	10	$-{2.47}_{-0.18}^{+0.15}$	$-{3879}_{-60}^{+54}$	74.38	36	2.07
0.03	30	$-{1.61}_{-0.44}^{+0.40}$	$-{2928}_{-27}^{+31}$	61.51	36	1.71
0.1	1	$-{3.55}_{-0.18}^{+0.18}$	$-{5320}_{-57}^{+55}$	184.2	36	5.12
0.1	3	$-{2.77}_{-0.15}^{+0.15}$	$-{5538}_{-52}^{+53}$	115.4	36	3.21
0.1	10	$-{1.77}_{-0.35}^{+0.45}$	$-{3871}_{-28}^{+35}$	61.18	36	1.7
0.1	30	$-{1.35}_{-0.44}^{+0.35}$	$-{2965}_{-26}^{+28}$	56.68	36	1.57
0.3	1	$-{1.69}_{-0.17}^{+0.15}$	$-{3184}_{-26}^{+26}$	1079	36	30
0.3	3	$-{1.76}_{-0.29}^{+0.26}$	$-{6049}_{-27}^{+29}$	65.6	36	1.82
0.3	10	$-{1.31}_{-0.39}^{+0.35}$	$-{5758}_{-26}^{+27}$	53.34	36	1.48
0.3	30	$-{1.24}_{-0.33}^{+0.36}$	$-{4403}_{-26}^{+26}$	54.21	36	1.51
1	1	$-{1.61}_{-0.12}^{+0.10}$	$-{635}_{-25}^{+25}$	5518	36	153
1	3	$-{1.54}_{-0.19}^{+0.15}$	$-{2976}_{-25}^{+26}$	1691	36	47
1	10	$-{1.29}_{-0.28}^{+0.26}$	$-{6518}_{-25}^{+26}$	90.08	36	2.5
1	30	$-{1.21}_{-0.31}^{+0.32}$	$-{7258}_{-26}^{+26}$	53.73	36	1.49

Note. The free parameters in the fitting procedure, ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]$ and ${({R}_{p}/{R}_{s})}_{o}$, are described in Section 6.3.2. We report the median and 68% confidence intervals from the MCMC optimization procedure for each free parameter. With 38 data points and 2 fitted parameters, we assume 36 DOF when calculating the ${\chi }_{\mathrm{red}}^{2}$ statistic.

Download table as:  ASCIITypeset image

Figure 13 shows the best-fitting equilibrium cloud and photochemical haze CPAT-absorber models. Additionally shown is the constant-${R}_{p}/{R}_{s}$ CPAT-absorber model, which illustrates the effect of a heterogenous distribution of absorbers in the stellar photosphere on an otherwise flat transmission spectrum. It shows that variations in the strength of absorbers across the stellar disk can produce large deviations from a flat spectrum in the optical, while simultaneously preserving a flat transmission spectrum in the near-infrared.

Zoom In Zoom Out Reset image size

For both the cloud and haze cases, the best fits to the data using the CPAT-absorber model require 3.2% of the unocculted stellar disk to possess weaker absorption features than the transit chord. In effect, the unocculted stellar disk is brighter in optical absorption bands than the region of the photosphere typified by the transit chord, which decreases the observed ${R}_{p}/{R}_{s}$ within those absorption bands. This allows the CPAT-absorber model to provide better fits to the data than the flat line model or models for atmospheric transmission alone considered in Section 6.2, but still does not fully explain the observed spectrum.

We also investigate the contribution of stellar temperature heterogeneities to the observed transmission spectrum using the CPAT model (the "CPAT-temperature model"). We employ PHOENIX models with a surface gravity ($\mathrm{log}g=5.0$) and a metallicity ($[\mathrm{Fe}/{\rm{H}}]=0.0$), the closest PHOENIX grid values to those identified for GJ 1214 by Rojas-Ayala et al. (2012), and temperatures in the range $2700\,{\rm{K}}\leqslant {T}_{\mathrm{eff}}\leqslant 4700\,{\rm{K}}$. We conduct the MCMC optimization procedure as described in Section 6.3.2, linearly interpolating in temperature space within the PHOENIX model grid to generate the S_o and S_u spectra based on the effective temperatures of the occulted and unocculted regions (T_o and T_u, respectively) at each step. We fix ${T}_{o}=3252$ K using the effective temperature of GJ 1214 (Anglada-Escudé et al. 2013) and $F=0.032$ following the results of Barnes et al. (2015). The free parameters in the model are the temperature contrast ${\rm{\Delta }}T$, which determines T_u relative to T_o, and ${({R}_{p}/{R}_{s})}_{o}$. We place a Gaussian prior on ${\rm{\Delta }}T$ centered on 0 K, with an uncertainty equal to 10% of the T_eff of GJ 1214 to thoroughly explore the parameter space and truncate it on the range $[-552\,{\rm{K}},+1448\,{\rm{K}}]$, to keep T_u within the PHOENIX model range we consider $[2700\,{\rm{K}},4700\,{\rm{K}}]$. As with the CPAT-absorber model, we conduct an MCMC optimization for each of the exoplanet atmosphere models in the grids of cloudy and hazy models described in Section 6.2.

Figure 14 illustrates the goodness-of-fit for the full grids of cloudy and hazy transmission spectra modulated by the CPAT-temperature model. As with the CPAT-absorber model, the best-fitting parameters for the exoplanet's atmosphere are similar to those found when considering the contribution of the exoplanet atmosphere alone to the transmission spectrum (Figure 8). However, allowing for temperature heterogeneities provides better fits to the data than either the exoplanet-alone or CPAT-absorber models. The complete results for the CPAT-temperature model fitting for the equilibrium cloud and photochemical haze model grids are provided in Tables 8 and 9, respectively. The best-fitting cloud model includes an exoplanet atmosphere with the highest metallicity and thickest clouds of the models we considered and a temperature contrast of ${\rm{\Delta }}T\simeq +0.1{T}_{\mathrm{eff}}$ for the unocculted region of the star's photosphere ($1000\times$ solar, ${f}_{\mathrm{sed}}=0.01$, ${\rm{\Delta }}T\ ={354}_{-47}^{+46}$ K, ${({R}_{p}/{R}_{s})}_{o}={9054}_{-94}^{+110}$ ppm, ${\chi }^{2}=45.2$, $\mathrm{dof}=36$, ${\chi }_{\mathrm{red}}^{2}=1.25$). The best-fitting haze model includes a smaller mode particle size and the same haze-forming efficiency as the best-fitting exoplanet-only model and ${\rm{\Delta }}T\simeq +0.1{T}_{\mathrm{eff}}$ as well (r = 0.1 μm, ${f}_{\mathrm{haze}}=10 \%$, ${\rm{\Delta }}T={354}_{-46}^{+46}$ K, ${({R}_{p}/{R}_{s})}_{o}\ =-{3064}_{-102}^{+101}$ ppm, ${\chi }^{2}=40.5$, $\mathrm{dof}=36$, ${\chi }_{\mathrm{red}}^{2}=1.13$). This model is the best-fitting of all those we considered in this analysis. The posterior distributions for the free parameters in each of the best-fitting cloud and haze CPAT-temperature models are shown in Figure 15. The ${\rm{\Delta }}T$ and ${({R}_{p}/{R}_{s})}_{o}$ parameters are positively correlated due to the fact that larger temperature contrasts depress model ${R}_{p}/{R}_{s}$ values at all wavelengths and thus require larger offsets to bring models in line with observations. A model using a constant ${R}_{p}/{R}_{s}$ value for the exoplanet's transmission spectrum together with a CPAT-temperature model for the star provides a fit to the data comparable to that of the best-fitting CPAT-temperature cloud model but not as good as that of the best-fitting CPAT-temperature haze model (${\rm{\Delta }}T={336}_{-45}^{+54}$ K, ${R}_{p}/{R}_{s}=0.11680\pm 1\times {10}^{-4}$, χ² = 45.6, dof = 36, ${\chi }_{\mathrm{red}}^{2}\ =1.27$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 8.  Results from MCMC Fits to CPAT-temperature Models for Cloudy Atmospheres

Metallicity	fseda	${\rm{\Delta }}T$	${({R}_{p}/{R}_{s})}_{o}$	${\chi }^{2}$	dof	${\chi }_{\mathrm{red}}^{2}$
(×Solar)		(K)	(ppm)
100	0.01	${417}_{-41}^{+39}$	$-{45}_{-89}^{+91}$	142.9	36	3.97
100	0.1	${328}_{-43}^{+43}$	${759}_{-90}^{+99}$	368.9	36	10.2
100	1	${344}_{-32}^{+30}$	${2735}_{-68}^{+73}$	3887	36	108
100	cloud-free	${389}_{-40}^{+29}$	${3266}_{-83}^{+75}$	6785	36	188
150	0.01	${408}_{-41}^{+43}$	${1734}_{-90}^{+98}$	102	36	2.83
150	0.1	${329}_{-45}^{+42}$	${2424}_{-96}^{+95}$	235.7	36	6.55
150	1	${331}_{-34}^{+34}$	${4039}_{-71}^{+82}$	2563	36	71.2
150	cloud-free	${359}_{-28}^{+30}$	${4457}_{-66}^{+69}$	4555	36	127
200	0.01	${399}_{-41}^{+43}$	${3073}_{-92}^{+97}$	80.93	36	2.25
200	0.1	${333}_{-45}^{+43}$	${3689}_{-97}^{+95}$	166.3	36	4.62
200	1	${324}_{-35}^{+39}$	${5076}_{-81}^{+84}$	1827	36	50.8
200	cloud-free	${347}_{-30}^{+29}$	${5426}_{-66}^{+69}$	3274	36	90.9
250	0.01	${392}_{-42}^{+44}$	${4126}_{-98}^{+95}$	69.26	36	1.92
250	0.1	${332}_{-44}^{+44}$	${4676}_{-98}^{+95}$	128.5	36	3.57
250	1	${317}_{-38}^{+40}$	${5903}_{-84}^{+87}$	1368	36	38
250	cloud-free	${338}_{-32}^{+32}$	${6213}_{-72}^{+72}$	2489	36	69.1
300	0.01	${384}_{-42}^{+45}$	${4957}_{-93}^{+102}$	61.85	36	1.72
300	0.1	${332}_{-42}^{+46}$	${5465}_{-99}^{+94}$	104.1	36	2.89
300	1	${315}_{-39}^{+41}$	${6572}_{-87}^{+91}$	1070	36	29.7
300	cloud-free	${329}_{-33}^{+35}$	${6842}_{-76}^{+78}$	1967	36	54.6
1000	0.01	${354}_{-47}^{+46}$	${9054}_{-94}^{+110}$	45.15	36	1.25
1000	0.1	${324}_{-47}^{+47}$	${9331}_{-103}^{+103}$	46.13	36	1.28
1000	1	${309}_{-42}^{+47}$	${9910}_{-96}^{+99}$	210.9	36	5.86
1000	cloud-free	${298}_{-42}^{+47}$	${10019}_{-92}^{+100}$	404.8	36	11.2

Note. The free parameters in the fitting procedure, ${\rm{\Delta }}T$ and ${({R}_{p}/{R}_{s})}_{o}$, are described in Section 6.3.3. We report the median and 68% confidence intervals from the MCMC optimization procedure for each free parameter. With 38 data points and 2 fitted parameters, we assume 36 DOF when calculating the ${\chi }_{\mathrm{red}}^{2}$ statistic.

^aCloud-free models do not have an f_sed value.

Download table as:  ASCIITypeset image

Table 9.  Results from MCMC Fits to CPAT-temperature Models for Hazy Atmospheres

r	fhaze	${\rm{\Delta }}T$	${({R}_{p}/{R}_{s})}_{o}$	${\chi }^{2}$	dof	${\chi }_{\mathrm{red}}^{2}$
(μm)	(%)	(K)	(ppm)
0.01	1	${630}_{-29}^{+31}$	$-{1974}_{-73}^{+82}$	673	36	18.7
0.01	3	${546}_{-36}^{+36}$	$-{2860}_{-89}^{+88}$	255	36	7.08
0.01	10	${404}_{-38}^{+46}$	$-{3014}_{-91}^{+98}$	55.96	36	1.55
0.01	30	${345}_{-45}^{+48}$	$-{2416}_{-98}^{+106}$	42.93	36	1.19
0.03	1	${616}_{-29}^{+34}$	$-{2780}_{-79}^{+83}$	542.8	36	15.1
0.03	3	${527}_{-37}^{+39}$	$-{3089}_{-93}^{+87}$	252.2	36	7.01
0.03	10	${392}_{-41}^{+46}$	$-{2883}_{-89}^{+108}$	51.3	36	1.43
0.03	30	${344}_{-47}^{+46}$	$-{2145}_{-103}^{+103}$	42.46	36	1.18
0.1	1	${487}_{-38}^{+38}$	$-{3791}_{-87}^{+91}$	227.4	36	6.32
0.1	3	${445}_{-39}^{+40}$	$-{4334}_{-91}^{+90}$	96.09	36	2.67
0.1	10	${354}_{-46}^{+46}$	$-{3064}_{-102}^{+101}$	40.52	36	1.13
0.1	30	${329}_{-45}^{+51}$	$-{2222}_{-101}^{+107}$	42.54	36	1.18
0.3	1	${318}_{-37}^{+41}$	$-{2459}_{-82}^{+89}$	1071	36	29.7
0.3	3	${325}_{-44}^{+47}$	$-{5304}_{-102}^{+97}$	53.34	36	1.48
0.3	10	${322}_{-48}^{+48}$	$-{5029}_{-102}^{+107}$	41.24	36	1.15
0.3	30	${314}_{-45}^{+52}$	$-{3691}_{-104}^{+105}$	43.98	36	1.22
1	1	${358}_{-28}^{+27}$	${174}_{-68}^{+62}$	5470	36	152
1	3	${297}_{-39}^{+44}$	$-{2298}_{-87}^{+93}$	1689	36	46.9
1	10	${272}_{-51}^{+51}$	$-{5893}_{-109}^{+107}$	91.88	36	2.55
1	30	${302}_{-47}^{+51}$	$-{6572}_{-102}^{+110}$	47.19	36	1.31

Note. The free parameters in the fitting procedure, ${\rm{\Delta }}T$ and ${({R}_{p}/{R}_{s})}_{o}$, are described in Section 6.3.3. We report the median and 68% confidence intervals from the MCMC optimization procedure for each free parameter. With 38 data points and 2 fitted parameters, we assume 36 DOF when calculating the ${\chi }_{\mathrm{red}}^{2}$ statistic.

Download table as:  ASCIITypeset image

The transmission spectra for the best-fitting equilibrium cloud and photochemical haze CPAT-temperature models are shown in Figure 16. The constant ${R}_{p}/{R}_{s}$ model, shown as a dashed black line, illustrates the effect of a temperature heterogeneity in the stellar photosphere on an otherwise flat transmission spectrum. Its similarity in the optical to the cloud and haze models owes to the fact that the best-fitting exoplanet transmission spectra are essentially flat in the optical, and the observed variations are due to features imprinted by the heterogeneous stellar photosphere. In each of these cases, a region of the unocculted stellar disk is brighter than that occulted by the transiting planet, which effectively decreases the observed ${R}_{p}/{R}_{s}$ during the transit. The effect is chromatic, producing the largest change in ${R}_{p}/{R}_{s}$ in the optical, where the difference in emergent flux between the stellar spectral models is most pronounced.

Zoom In Zoom Out Reset image size

Cool, unocculted spots, on the other hand, would increase the observed ${R}_{p}/{R}_{s}$ values. We use the CPAT-temperature model to investigate how cool, unocculted starspots reported in the literature could affect the transmission spectrum. Long-term monitoring shows that GJ 1214 demonstrates a $1 \%$ peak-to-peak variability in the MEarth (Nutzman & Charbonneau 2008) bandpass ($715\,{\rm{nm}}\lt \lambda \lt 1000\,{\rm{nm}}$) on a timescale that is an integer multiple of 53 days (Berta et al. 2011), which has been attributed to rotational modulation of cool starspots (Charbonneau et al. 2009; Berta et al. 2011; Fraine et al. 2013). We calculate the spot-covering fraction implied by the variability using PHOENIX model spectra with $\mathrm{log}g=5.0$ and $[\mathrm{Fe}/{\rm{H}}]=0.0$. Assuming the spots are 10% cooler than the unspotted surface, following results from Doppler imaging of M dwarfs (Barnes et al. 2015), we utilize models with $T=3252$ K and $T=2927$ K for the unspotted and spotted surfaces, respectively, by linearly interpolating in temperature space between models in the PHOENIX grid. Integrating over the MEarth bandpass,¹⁷ we find that a spot-covering fraction $F=0.03$ of the stellar disk can reproduce the reported variability in the MEarth bandpass. The effect of such a spot configuration on the observed transmission spectrum is shown in gray on Figure 16. In the near-infrared, the scale of the effect is smaller than the reported measurement uncertainties. The most pronounced change is present in the optical, where the observed transmission spectrum is increased as much as ${\rm{\Delta }}({R}_{p}/{R}_{s})=8\times {10}^{-4}$ above the near-infrared spectrum. The optical transmission spectrum would decrease by a similar amount if these starspots were instead occulted by the transiting exoplanet (Pont et al. 2013). However, the observed decrease in the optical transmission spectrum is $\sim 3\times$ larger than what would be caused by these spots. Additionally, no brightening events due to occulted starspots are evident in our data (Figure 3).

In summary, we find that cool starspots cannot account for the decreased optical ${R}_{p}/{R}_{s}$ values, though unocculted bright regions of the photosphere with ${\rm{\Delta }}T\simeq +0.1{T}_{\mathrm{eff}}$ can decrease a flat optical transmission spectrum to the observed values.

The results of our modeling efforts indicate that the deviations from a flat optical transmission spectrum for GJ 1214b could be introduced by heterogeneities in the stellar photosphere, either in terms of temperature or the distribution of absorbers (Figures 13 and 16). However, the near-infrared transmission spectrum is largely unaffected by the stellar photosphere. The combination of the optical and near-infrared data, then, allows us to probe the stellar and planetary contributions to the transmission spectrum simultaneously.

The optical spectrum of a M4.5 dwarf like GJ 1214 is largely driven by opacity from TiO molecular bands (Morgan et al. 1943). Using the CPAT-absorber model framework, we find that an unocculted region of the stellar photosphere with relatively weak absorption in these bands could imprint spectral features on an otherwise flat exoplanetary transmission spectrum. However, the best-fitting CPAT-absorber model requires a metallicity contrast of ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]=-1.31$, corresponding to a depletion in the abundance of absorbers by a factor of 20 in the unocculted region. Efforts to measure a latitudinal dependence of the solar spectrum have found elemental abundances to be within 0.005 dex across latitudes (Kiselman et al. 2011). Moreover, as stars with masses less than 0.35 M_⊙ (corresponding to spectral types M3.5 and later) are fully convective (Chabrier & Baraffe 1997), a viable mechanism for maintaining significant sustained abundance differences in the stellar photosphere of a M4.5 dwarf like GJ 1214 is not immediately apparent.

Temperature heterogeneities, however, are known to exist in stellar photospheres. We find a temperature heterogeneity could produce the observed discrepancy in ${R}_{p}/{R}_{s}$ between the optical and near-infrared if $\sim 3.2 \%$ of the unocculted stellar disk is ∼350 K hotter than the remaining photosphere. GJ 1214 is known to host cool starspots from starspot-crossing events detected in transit light curves (Carter et al. 2011; Kreidberg et al. 2014). Mid- to late-M dwarfs have also been found to have abundant polar spots from Doppler imaging. Barnes et al. (2015) found the M4.5 dwarf GJ 791.2A to have a mean spot coverage of 3.2% and a maximum spot coverage of 82.3% during its rotation, assuming spots are 300 K cooler than the photosphere. On the Sun, spots are accompanied by bright plage, which include faculae with temperature contrasts of 300–500 K with their surroundings (Topka et al. 1997). Brightening from these faculae overpowers spot darkening, causing total solar irradiance to increase during solar maxima (Fröhlich & Lean 1998; Meunier et al. 2010). Faculae cover roughly 0.36% of the sky-projected solar disk during periods of low solar activity and 3% during high activity (Shapiro et al. 2014). Additionally, faculae are much more common than spots on the Sun, covering $100\times$ more disk-area during periods of low activity, and $10\times$ more area during periods of high activity (Shapiro et al. 2014).

Like the Sun, old, slowly rotating FGK stars are known to have activity cycles dominated by faculae, unlike their younger counterparts, which are spot-dominated (Radick et al. 1983; Lockwood et al. 2007). With an age of 3–10 Gyr (Charbonneau et al. 2009), the photometric variability of GJ 1214 may be faculae-dominated as well. Our data are most consistent with unocculted faculae in the photosphere of GJ 1214 producing the observed offset in the optical transmission spectrum of GJ 1214b.

The results of our CPAT-absorber and CPAT-temperature modeling efforts demonstrate that, in principle, a heterogeneous stellar photosphere can provide a transmission spectroscopy signal larger than that introduced by the exoplanetary atmosphere in the optical, adding another degeneracy to the modeling and interpretation of exoplanet spectra. At longer wavelengths, however, the effect is less pronounced for a star with the T_eff of GJ 1214, and differences between the atmospheric models are more apparent. We model the effect out to the maximum wavelength of the PHOENIX stellar models (5.5 μm), and find the largest differences between transmission spectra for cloudy and hazy atmospheres at 2–5 μm for both CPAT-absorber and CPAT-temperature cases. This would be a promising region to target in transmission spectroscopy studies with the James Webb Space Telescope (JWST), in order to distinguish between cloudy and hazy atmospheres for this planet. Optical spectra like those presented here can complement JWST investigations of GJ 1214 and other targets by constraining the contribution of heterogeneous stellar photospheres to observed transmission spectra.

This analysis suggests that facular brightening may contribute more to transmission spectra of exoplanets than has been previously recognized. If common for M dwarfs, this photospheric heterogeneity could complicate optical transmission spectroscopy studies of exoplanets around small stars, such as GJ 1132b (Berta-Thompson et al. 2015) and TRAPPIST-1b, c, and d (Gillon et al. 2016), including searches for Rayleigh scattering. While unocculted starspots can mimic the transmission signature of scattering in an exoplanetary atmosphere (e.g., McCullough et al. 2014), facular brightening has the potential to mask a scattering signature by reducing the Rayleigh scattering slope in a transmission spectrum. However, the method we provide here can be used to take unocculted faculae into account if high SNR transmission spectra are obtained. Additionally, the same exoplanets would provide a unique opportunity to gain spatial information about M dwarf surfaces if the exoplanetary and photospheric contributions can be uniquely identified. The Transiting Exoplanet Survey Satellite (TESS), which will utilize a long-pass filter from 0.6 to 1.0 μm, will measure transit depths for potentially hundreds of planets around M dwarfs; comparing the TESS transit depths with follow-up measurements in the near-infrared could be one way of probing for spectral features introduced by a heterogeneous stellar photosphere. Along a different avenue of research, future work can investigate the range of stellar temperatures and metallicities for which this effect will be important.