The Influence of Stellar Contamination on the Interpretation of Near-infrared Transmission Spectra of Sub-Neptune Worlds around M-dwarfs

Transit spectroscopy of planets around small, cool stars offers a promising opportunity to characterize small, cool worlds (e.g., Seager & Sasselov 2000; Swain et al. 2008; Huitson et al. 2012; Gillon et al. 2016, 2017; Luger et al. 2017). An important consideration in characterizing such planets is the nature of the host star (Apai et al. 2018). Transmission spectroscopy of hot Jupiters to terrestrial worlds alike have demonstrated evidence of stellar activity (Pont et al. 2007; Alonso et al. 2008; Zhang et al. 2018; Rackham et al. 2017). M-dwarf hosts, in particular, exhibit higher recorded detections of stellar activity in the form of flares (Hawley & Pettersen 1991; Hawley et al. 2014; Schmidt et al. 2014) and H_α chromospheric emissions (Delfosse et al. 1998) relative to F-, G-, and K-type host stars, indicative of frequent occurrence of photospheric heterogeneities such as stellar spots and faculae. Study of M-dwarf systems are extremely timely, especially with the recent launch of NASA's Transiting Exoplanet Survey Satellite, which will discover thousands of exoplanets around such stars (Muirhead & Veyette 2018; Ballard 2019), providing the need for follow up with future missions like the James Webb Space Telescope (JWST).

Stellar spots that are occulted by a transiting planet contribute to observable changes in the transit light curve. Techniques to model the occulting heterogeneity have been successful to some degree by either excluding the spot-crossing event from the transit fit (Pont et al. 2008; Carter et al. 2011; Narita et al. 2013) or by correcting for epoch-to-epoch variations in the star from relative offsets between observations, therefore accounting for any differences in the planetary spectrum (Zellem et al. 2017). These techniques, however, do not account for homogeneously scattered regions of heterogeneities that are persistent with the entire data in all observations.

Unocculted stellar spots as well as homogeneously covered regions of smaller spots affect the transit baseline in addition to causing discrepancies in the transit depth measurements depending on the fractional coverage in the area of the spots and their temperature contrast relative to the photosphere (Pont et al. 2008; Silva-Valio 2008; Czesla et al. 2009; Wolter et al. 2009; Agol et al. 2010; Berta et al. 2011; Carter et al. 2011; Désert et al. 2011; Sing et al. 2011; Fraine et al. 2014; McCullough et al. 2014; Oshagh et al. 2014; Barstow et al. 2015b; Damasso et al. 2015; Zellem et al. 2015; Rackham et al. 2017). Attempts to understand the nature of unocculted spots have thus far been done with photometric monitoring of the star, i.e., tracking the rotational variability in the modulating brightness on the stellar photosphere translating onto the light curve (Pont et al. 2008, 2013; Berta et al. 2011; Désert et al. 2011; Sing et al. 2011; Knutson et al. 2012; Zellem et al. 2015). Despite reasonable success, efforts with photometric monitoring of unocculted stellar spots have been inadequate in teasing out underlying heterogeneities spread across the entire disk of the star, especially those that are not significant enough to show variability effects on the stellar light curve (Jackson & Jeffries 2012; Rackham et al. 2017, 2018). Studies by Chapman (1987) and Shapiro et al. (2014) have also shown that every differential region of the star's photosphere, whether the transit chord or the entire disk, has its own unique spectrum. Given these challenges, the effect of unocculted stellar spots remains a persistent problem despite corrections from photometric variability monitoring, which at best provide only a lower limit estimate for regions of heterogeneities on the photosphere (Rackham et al. 2018).

Several recent studies (Pont et al. 2013; McCullough et al. 2014; Oshagh et al. 2014; Barstow et al. 2015b, 2015c; Rackham et al. 2017, 2018; Barstow 2018; Bruno et al. 2018; Murgas et al. 2019; Pinhas et al. 2018; Zhang et al. 2018) have demonstrated that unocculted as well as homogeneous spots and faculae on the host star present themselves in a wavelength-dependent manner in the planetary transit data. The wavelength dependence of stellar heterogeneity signatures is particularly influential for transit spectroscopy where these features can often be degenerate with planetary atmosphere features. This necessitates accounting for stellar heterogeneities in order to probe the true nature of a planetary atmosphere. Improper characterization of the effects of stellar surface heterogeneities on the transmission spectrum limits our ability to place precision constraints on basic planetary conditions like temperature structure, composition, and cloud properties. Therefore, it is vital to thoroughly understand the influence of stellar surface heterogeneities on our ability to infer basic planetary properties in the era of high-fidelity observations anticipated from the JWST.

To understand these degeneracies and the impact of not accounting for unocculted star spots on atmospheric inference, we leverage the utility of atmospheric retrieval. In Section 2, we briefly explore the effect of stellar heterogeneities on the shape of a planetary transmission spectrum analytically, providing intuition for degeneracies that might occur in interpretations of transit spectra, followed by a description of our atmospheric retrieval model. In Section 3 we investigate the level of bias that arises while inferring planetary properties from JWST transmission spectra when ignoring corrections for stellar heterogeneities. Additionally, we address the question of whether a priori knowledge of stellar heterogeneity, given the current state of M-dwarf stellar models, could improve precision on the inferred planetary characteristics from JWST transmission spectra. In this section, we also investigate the potential biases arising directly from the stellar models when inferring planetary properties. Section 4 presents a discussion on various sources of bias that can affect appropriate characterization of the transmission spectrum in addition to the stellar contamination correction, followed by a brief summary of our findings. In sum, our work demonstrates that, regardless of the coverage area or temperature contrast of a given stellar heterogeneity feature, one can infer unbiased planetary properties retrieved from contaminated transmission spectra provided an appropriate correction for stellar heterogeneity has been incorporated; however, these are conditional on our confidence of stellar atmospheric models being accurate representations of the true photosphere.

Here we briefly build our intuition for how the effects of stellar contamination can alter the shape of the planetary transmission spectrum. The transmission spectrum corrected for stellar surface heterogeneities is given by

$$\begin{eqnarray}&&{\alpha }_{\lambda ,c}={\epsilon }_{\lambda }{\alpha }_{\lambda }\end{eqnarray} \tag{ 1 }$$

where α_λ is the native planetary transit spectrum and _λ is the contamination factor (McCullough et al. 2014; Zellem et al. 2017; Rackham et al. 2018, see their Equations (1) and (2)) given by

$$\begin{eqnarray}&&{\epsilon }_{\lambda }=\displaystyle \frac{1}{\left(1-{f}_{s}\left(1-{Q}_{\lambda ,s}\right)\right)}\end{eqnarray} \tag{ 2 }$$

where Q_λ,s is the spot-to-photospheric stellar flux ratio (F_λ,s/F_λ,p) and f_s is the fractional coverage of spots. The shape of the transmission spectrum (when corrected for stellar heterogeneities) is given by the derivative of Equation (1) combined with the expressions given in Line & Parmentier (2016; see their Equations (6)–(9) for a gas+gray cloud system):

$$\begin{eqnarray}&&\displaystyle \frac{d{\alpha }_{\lambda ,c}}{d\lambda }={\epsilon }_{\lambda }\displaystyle \frac{2{R}_{p}}{{R}_{* }^{2}}H\left(\displaystyle \frac{1}{1+{\gamma }_{\lambda }}\right)\displaystyle \frac{d\mathrm{ln}({\sigma }_{\lambda })}{d\lambda }+\displaystyle \frac{d{\epsilon }_{\lambda }}{d\lambda }\left(\displaystyle \frac{{{Rp}}^{2}}{{R}_{* }^{2}}+\displaystyle \frac{2{R}_{p}{z}_{\lambda }}{{R}_{* }^{2}}\right)\end{eqnarray} \tag{ 3 }$$

where R_p and R_* are the radius of the planet and star respectively, z_λ is the wavelength-dependent sharp occulting disk radius, H is the planetary atmospheric scale height, σ_λ is the absorption cross-section for a given gas, and γ_λ is the gas-to-gray cloud opacity ratio.

The shape of the stellar contamination spectrum, which is simply the wavelength-dependent derivative of the contamination factor (Equation (2)), is given by

$$\begin{eqnarray}&&\displaystyle \frac{d{\epsilon }_{\lambda }}{d\lambda }=\displaystyle \frac{\left({f}_{s}\tfrac{{{dQ}}_{\lambda ,s}}{d\lambda }\right)}{{\left(1-{f}_{s}+{f}_{s}{Q}_{\lambda ,s}\right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

Simplistically assuming that the stellar fluxes behave like blackbodies and also that the temperature contrast between the photosphere and spot is small, consistent with Doppler imaging observations for M-dwarfs (Strassmeier 2010; Andersen & Korhonen 2015), i.e., T_p ≈ T_s, therefore we can write

$$\begin{eqnarray}&&{Q}_{\lambda ,s}\approx 1+\displaystyle \frac{{hc}}{\lambda {k}_{B}}\left(\displaystyle \frac{1}{{T}_{p}}-\displaystyle \frac{1}{{T}_{s}}\right)\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dQ}}_{\lambda ,s}}{d\lambda }\approx -\displaystyle \frac{{hc}}{{\lambda }^{2}{k}_{B}}\left(\displaystyle \frac{1}{{T}_{p}}-\displaystyle \frac{1}{{T}_{s}}\right).\end{eqnarray} \tag{ 6 }$$

At longer wavelengths ${{dQ}}_{\lambda ,s}/d\lambda$ (as well as Equation (4)) tends to zero with _λ approaching a constant "stretching" term, thus nulling the second term in Equation (3). From the surviving first term, we can see clear degeneracies between the stellar contamination _λ and the scale height H which both act like "stretching" terms on the opacity structure; e.g., an unaccounted stretch due to the spot could be compensated with scale height, leading to a bias in the unknown terminator temperature and/or molecular weight.

At shorter wavelengths, on the other hand, both terms in Equation (3) survive, resulting in a complex shape interaction. The scaling behavior from the first term persists with a power-law-like behavior from the second term, yielding additional degeneracies with z_λ—the wavelength-dependent radius alongside the factors that contribute to the shape of the stellar contamination spectrum, $\tfrac{{{dQ}}_{\lambda }}{d\lambda }$ i.e., the temperature contrast between T_p and T_s and fractional coverage in area of spots, f_s.

Figure 1 shows example contamination spectra for both spots and faculae with various temperature contrasts and covering fractions (e.g., similar to Rackham et al. 2017 and Zhang et al. 2018) for realistic stellar spectra (using PHOENIX models from Allard et al. 2003, 2007 and Allard & Freytag 2009). Consistent with the preceding analytical result, with realistic spectra we also see that increasing the spot-covering fraction f_s from a small 0.9% to a large 12% value produces a scaling effect at longer wavelengths. Moreover, at shorter wavelengths, we see a power-law (or slope) effect, similar to what may be expected from scattering in the sub-Neptune atmosphere, as also demonstrated for the case of hot Jupiter HD 189733b (McCullough et al. 2014).

Zoom In Zoom Out Reset image size

Now, armed with a first-order intuition for the effect of stellar contamination on the transmission spectrum, we can explore numerically, through retrievals, the influence of spot contamination.

In this section we aim to explore the role of spot contamination by addressing the following questions.

• 1.  
  Will there be a bias in JWST results if we do not account for stellar contamination in atmospheric retrieval models? How can we rigorously explore this bias to decode degeneracies that arise in transmission spectra?
• 2.  
  Does a priori knowledge of stellar contamination parameters reduce bias in our interpretations?

We first describe our retrieval setup, followed by a numerical exploration addressing each of the aforementioned questions.

3.1. Model Description

We use the planetary transmission/retrieval tool CHIMERA¹ (Line et al. 2013; Swain et al. 2014; Kreidberg et al. 2014, 2015; Mai & Line 2019) to simulate and retrieve upon JWST-like spectra from 0.6 to 5 μm (e.g., NIRISS, NIRCam, or NIRSpec). The specific model used here is from Mai & Line (2019) but with modifications to accommodate for the transit light source effect (Apai et al. 2018; Rackham et al. 2018) arising from unocculted starspots (e.g., Equation (2)). The relevant model variables are given in Table 1. Briefly, the model² computes a limb transmission spectrum via the relations in Brown et al. (2001) and Tinetti et al. (2012). We include correlated-K opacities (derived from the line-by-line cross-sections of Freedman et al. 2008, 2014) for H₂–H₂/He collision-induced absorption, H₂O, CH₄, CO, CO₂, NH₃, Na, K, TiO, VO, C₂H₂, HCN, H₂S, and FeH. The retrieval uses the "chemically consistent" approach (Kreidberg et al. 2015, 2018) given a metallicity and C/O assuming equilibrium gas+condensate chemistry (computed using the NASA Chemical Equilibrium with Applications; Gordon & McBride 1996), computed along the three-parameter analytic T–P profile (Guillot 2010; Line et al. 2013). Hazes and clouds are parameterized with a simple power law (Etangs et al. 2008) and gray uniform opacity, respectively.

Table 1.  Model Parameters, Nominal "Truth" Values, and Uniform Prior Ranges

Parameter	Description	Initial Value	Status	Prior Range
log(κir)a	Thermal profile gray IR opacity (cm2 g−1)	0.306	Fixed	N/A
log(γv)a	Thermal profile Vis/IR opacity	−2.02	Fixed	N/A
Rpb	Adopted 10 bar planetary radius (Rjup)	0.239	Fixed	N/A
Mpb	Planetary mass (Mjup)	0.0204	Fixed	N/A
R${}_{* }$ b	Stellar radius (R⊙)	0.211	Fixed	N/A
Tirra	Irradiation temperature (K)	794.0	Retrieved	[300, 3000]
[M/H]	Planetary metallicity	1.61 (40×) or 2.47 (300×)	Retrieved	[−2, 3]
log(C/O)	C-to-O ratio	−0.72	Retrieved	[−2, 2]
×Rp	10 bar radius scaling	1.0	Retrieved	[0.5, 1.5]
log(${\sigma }_{0}$)c	Haze cross-section amplitude	−10.0	Retrieved	[−15, 2]
βc	Haze scattering slope	4.	Retrieved	[0, 6]
log(κc)	Well-mixed gray cloud opacity	−35.5 (clear) or −29 (cloudy)	Retrieved	[−40, −20]
Tp	Stellar photosphere temperature (K)	3300	Retrieved	[2000, 4000]
Ts	Stellar spot temperature (K)	2838	Retrieved	[2000, 4000]
fs	Fractional coverage of spots	Case specificd	Retrieved	[0, 1]
Tf	Stellar faculae temperature (K)	3400	Retrieved	[2000, 4000]
ff	Fractional coverage of faculae	Case specificd	Retrieved	[0, 1]

Notes.

^aThree-parameter model for thermal profile (Guillot 2010). ^bGJ1214/b, from exoplanets.org (Han et al. 2014). ^cFrom Etangs et al. (2008). Haze cross-section given by ${\sigma }_{\lambda }={\sigma }_{0}{({\lambda }_{0}/\lambda )}^{\beta }$ w ith λ₀ = 0.43 μm and σ₀ given in units relative to H₂ Rayleigh scattering (2 × 10⁻²⁷ cm²). ^dInitial values are specific to cases of spot- and faculae-covering fractions, listed in Table 2, also given in Rackham et al. (2018) and Zhang et al. (2018).

Download table as:  ASCIITypeset image

We use a grid of PHOENIX (Allard et al. 2003, 2007; Allard & Freytag 2009) stellar models (from the STScI pysynphot routine; STScI Development Team 2013) in effective temperature to construct the photospheric as well as spot/facular region spectrum (fixing the stellar metallicity to solar and log(g) to 5.0). Within the forward model/retrieval, the photosphere/spot spectra are interpolated to the appropriate effective temperature using the Scipy package griddata (Virtanen et al. 2019).

With this modeling framework, we investigate the influence of stellar contamination on a generic M-dwarf-orbiting sub-Neptune (a clear 40× solar, and a cloudy 300× solar atmosphere; Table 1) observed with a JWST-like platform from 0.6 to 5 μm (broken up between NIRISS-like: 0.6–2.5 μm and NIRCam/NIRSpec G395-like: 2.4–5.0 μm.). In all cases the resolving power (R) is assumed at 100 and, for simplicity, we assume a wavelength-independent uncertainty of 30 ppm per bin, loosely motivated by reasonable JWST performance expectations (e.g., Greene et al. 2016). We do not randomize the data-points so as to remove any bias due to arbitrary noise instantiations (e.g., Feng et al. 2018 and Mai & Line 2019). Parameter estimation and model selection (used to assess the need for spot-correction) are performed with the PyMultiNest (Feroz et al. 2009; Buchner et al. 2014) nested sampling routine (Skilling et al. 2006) with uniform priors, given in Table 1.

Figure 2 illustrates the influence of unocculted spot-contamination on sub-Neptune transmission spectra. Consistent with recent studies (Pinhas et al. 2018; Rackham et al. 2018), there can be a significant shape change in the transmission spectrum, where the wavelength-dependent stellar contamination parameters could mimic atmospheric effects on the planet's spectrum (e.g., molecular features and haze-like power-law slopes).

Zoom In Zoom Out Reset image size

3.2. To What Degree Is There Bias if Not Accounting for Unocculted Star Spots?

3.2.1. JWST NIRISS-like Bandpass

We choose spot/faculae-covering fractions and contrasts to be consistent with 1% I-band rotational variability for mid-to-late M-dwarfs (Newton et al. 2016; Rackham et al. 2018). Specifically, We have assumed T_s = 0.86 × T_p K (or C_s = 0.86) and T_f = T_p + 100 K (or C_f = (T_p+100)/T_p) (Gondoin 2008; Afram & Berdyugina 2015; Rackham et al. 2018). We also explore a range of planetary atmospheric conditions (e.g., clear, cloudy, high metallicity (300× solar)/molecular weight) given in Table 2 and shown in Figure 2.

Table 2.  Bayes Factor between the Contamination-included Model and the Atmosphere-only Model

Stellar Contamination	Atmosphere Type	Bayes Factor (ln B)
NIRISS-like
fs = 0.9%	Clear	10
	Cloudy	2.5
	High metallicity + cloudy	0.23
fs = 12%	Clear	3191
	Cloudy	2653
	High metallicity + cloudy	1036
fs = 20%	Clear	10,634
	Cloudy	6085
	High metallicity + cloudy	2479
fs = 0.5%, ff = 4.6%	Clear	14,634
	Cloudy	14,628
	High metallicity + cloudy	6389
fs = 14%, ff = 63%	Clear	1,575,792
	Cloudy	1,543,800
	High metallicity + cloudy	66,478
NIRCam-like
fs = 14%, ff = 63%	Cloudy	63,295

Notes. Bayes Factor between an atmosphere-only model and a model that accounts for the spot contamination under the several atmosphere scenarios (Table 1) and a NIRISS-like (30 ppm, R = 100, 0.6–2.5 μm) and NIRCam-like (2.4–5 μm) observational setup.

Download table as:  ASCIITypeset image

For each scenario in Table 2 (as observed from 0.6 to 2.5 μm, R = 100, with 30 ppm precision) we perform a series of retrievals: first, a control case that assumes a homogeneous stellar photosphere (no spot contamination correction, retrieves only the atmospheric properties T_irr, [M/H], log(C/O), log(σ₀), β, and log(κ_c)) and, second, accounting for the spot contamination (e.g., additionally retrieves for T_p, T_s, f_s, T_f, and f_f ). Two figures of merit are used when assessing the level of "bias" in the inferred planetary properties. The first (Figures 3 and 4) is a qualitative metric describing the bias by simply measuring the actual deviation of the posterior probability distributions (medians of the distributions) from the truth values. We only report the bias in terms of the width of the posterior probability distributions (i.e., shift of the histogram in units of σ width relative to the retrieved precision) for a select few cases with distributions lying fully within the chosen prior range. Instances where the posterior histograms converge partly on the edge of the priors skew the true level of bias by falsely portraying additional width to the histograms; therefore, this is purely a "qualitative metric" to measure the bias in the inferences (e.g., see Figure 3, right column for high metallicity and cloudy planetary atmosphere with f_s = 12% for the carbon-to-oxygen ratio: log(C/O)). The second metric (discussed below) is a quantitative measure of the level of bias in the inferred planetary parameters, accounted for by computing the Bayes factor (Table 2, last column) between the atmosphere-only and the atmosphere+spot model retrieval results for all cases with varying levels of stellar heterogeneities.

Zoom In Zoom Out Reset image size

In all cases (see Table 2)³ with a JWST NIRISS-like transmission spectrum (covering 0.6–2.5 μm, with 30 ppm precision) of a generic sub-Neptune orbiting an M-dwarf, we demonstrate that there is significant "qualitative bias" in the retrieved atmospheric properties when not correcting for stellar contamination. For instance, in the case with a spot-covering fraction f_s = 0.9%, an uncorrected transmission spectrum for a clear atmosphere sub-Neptune yields a bias in the retrieved planetary irradiation temperature (T_irr), bulk metallicity ([M/H]), and the bulk carbon-to-oxygen ratio (log(C/O)) of 1.6σ, 13.7σ, and 95.5σ respectively (illustrated in Figures 3 and 4). Moreover, in Figure 4, a correlation between the planetary radius (fiducial radius at 10 bar × R_p) and the stellar spot-covering fraction (f_s) emerges from the retrieval 2D posterior histograms. Consistent with the previously illustrated analytical result (Equation (3)), the wavelength-dependent stellar contamination term is degenerate with both the planetary atmospheric scale height H as well as the wavelength-dependent sharp occulting disk radius z_λ, altogether affecting the shape of the transmission spectrum. We also see this effect in Figure 2, where increasing the level of stellar photospheric heterogeneity (a linear combination of the spot-covering fraction and the temperature contrast between the spot and the immaculate photosphere) modifies the shape of the transmission spectrum and hence the derived planetary radius significantly, beyond anticipated spectral precisions (15–120 ppm). Focusing, however, only on 30 ppm precision as an example case, the trend of increasing bias with the first metric (deviance from truth) persists for a combination of stellar heterogeneities, ranging from 0.5% to 20% in spot-covering fraction and 4.6%–63% for faculae-covering fraction, as illustrated in Figure 3.

Zoom In Zoom Out Reset image size

We also find that, for the smallest heterogeneity fractions, muting the transit spectral features with clouds and high metallicity reduces the effective "qualitative bias" as seen in the retrieved planetary parameters (see Figure 3, top right panel). For f_s = 0.9% in particular, we find that there is a bias of 0.48, 0.55, and 0.01σ for the inferred irradiation temperature, planetary metallicity, and carbon-to-oxygen ratio respectively, for an uncorrected transmission spectrum. With increasing stellar heterogeneity covering fractions, however (ranging from f_s = 12% to 20% alone and for f_s = 0.5% along with f_f = 4.6% up to 63%), the qualitative increase in the bias (deviation from the truth values) of the posterior probability distributions for cloudy and high-metallicity atmospheres (although slightly reduced) are consistent with those of a clear sub-Neptune atmosphere. Only for small spot-covering fractions do we find that the level of bias falls well within 1σ of the true value, and this effect is persistent through a range of spectral precisions (15–120 ppm; see Figure 5, left column for f_s = 0.6%). This is primarily because, as the transit feature signal-to-noise decreases (whether due to clouds and/or high mean molecular weight), thereby increasing degeneracy in the transmission spectrum, this results in larger parameter uncertainties. These uncertainties therefore seemingly reduce the effective bias under our first figure of merit, in turn motivating the question: at what degree of stellar heterogeneity in fractional coverage does the bias in retrievals matter?

Zoom In Zoom Out Reset image size

To address this question, we start by qualifying our second figure of merit, the Bayes factor, which provides a quantitative measure of the need to include as "nuisance" parameters a stellar correction scheme (last column in Table 2, for 30 ppm precision). Specifically, we compute the log of the Bayes factor (ln B; Cornish & Littenberg 2007) which is a metric permitting quantitative comparison of two models with different complexity (e.g., with and without a spot-correction parameterization). For all cases in Table 2, we highlight two points illustrating the behavior of the Bayes factor: (1) ln B increases with an increase in contamination covering fraction and (2) it decreases with increasing cloudiness/molecular weight in the planetary transmission spectrum. The first point is consistent for all planet atmosphere scenarios, as seen previously with our first figure of merit (Figure 3). For instance, a clear atmosphere scenario with increasing spot-covering fraction f_s from 0.9 to 20% increases ln B from 10 to 10,634, with the entire range demonstrating strong evidence for preferring the model including stellar contamination as per the Jeffrey's scale (Trotta 2008). The increasing trend in ln B is similar through all atmospheric scenarios including for cases with cloudy and high-metallicity atmospheres, where the ln B ranges from 0.23 to 2479 for f_s increasing from 0.9% to 20% (see Table 2). This pattern is consistent for combinations of spot- and faculae-covering fractions both for clear as well as cloudy and high-metallicity atmospheres, thereby quantitatively indicating the degree to which the model including stellar contamination is preferred. The second point shows a decrease with increasing feature muting (due to clouds/metallicity). For the entire range of stellar contamination cases in Table 2 we see this behavior, especially for the case of a small spot-covering fraction of 0.9% (under the cloudy and high-metallicity atmosphere), which yields a ln B = 0.23, demonstrating inconclusive evidence favoring either model (true for all ln B values on the Jeffrey's scale <1, Trotta 2008). Despite the decreasing trend of ln B with increasing atmospheric degeneracies, for other spot- and spot+faculae-covering fraction combinations, the Bayes factors strongly prefer the model including stellar corrections to prevent bias in the inferences, i.e., ln B > 5. From this trend, we can say that, regardless of the planetary atmosphere properties explored here, a system with a host star of spot-covering fraction >1% must include a parameterization for stellar activity to mitigate the bias atmospheric inferences.

We now focus on quantifying the bias as a function of spectral precision, under the more likely scenario of a cloudy, high-metallicity atmosphere (Fortney et al. 2013). We do this by generating a sparse grid of retrievals—varying spot-covering fractions on the first dimension, ranging from 0.6% to 12% and spectral precision on the second dimension ranging from 15 to 120 ppm, consistent with expectations for the JWST based on achieved Hubble Space Telescope/Wide Field Camera 3 precisions (e.g., down to 15 ppm; Line et al. 2016). For each point on the grid, we perform two sets of retrievals: the first with the atmosphere-only model and the second with stellar corrections included. Our main findings from this exercise are: for the entire range of precisions, an increase in fractional coverage in area of spots (f_s) on the stellar photosphere increases the significance with which the model with stellar corrections is favored (ln B). Particularly for spot-covering fractions f_s > 1%, we find that the Bayes factor ranges from moderately favoring the model with stellar contamination included (e.g., case of 120 ppm precision at f_s = 3% in Figure 6 is in the moderate evidence regime as per the Jeffrey's scale) to strongly favoring the model for all other sensitivities. This trend provides strong evidence for the spectral presence of stellar heterogeneities in the transmission spectrum throughout these levels of precision.

Zoom In Zoom Out Reset image size

For cases of f_s < 1% and for all sensitivities >30 ppm on our grid (see Figure 6), we only see weak to moderate evidence for stellar contamination. This is consistent with the reduction in our "qualitative bias" metric (see Figure 3, top panel, right column and Figure 5, left column) as demonstrated previously, due to increasing the level of degeneracies in the spectrum for cloudy/high-metallicity scenarios, especially at small spot-covering fractions. For the case of 15 ppm precision, however, the uncertainties on the transmission spectrum are significantly small enough that even a region of heterogeneity as small as f_s = 0.6% shows a very high ln B, providing strong spectral evidence for the contamination, although this is not apparent from our "qualitative" bias metric as seen in the top left panel of Figure 5 for f_s = 0.6%. Overall, we see that for spot-covering fractions >1%, for our defined grid of precision levels, increasing f_s leads to an increase in the "quantitative bias" (ln B), providing strong evidence for stellar contamination; however, this increase is penalized by the level of uncertainties on the JWST transmission spectrum, thereby reducing the apparent "qualitative bias" (deviance from the truth) as we see in the right column of Figure 5.

3.2.2. JWST NIRCam/NIRSPEC G395-like Bandpass

In Section 2 we showed that stellar contamination has a stretching effect on the shape of the transmission spectrum at longer wavelengths, and was therefore degenerate with the scale height. In this section we explore the role of stellar contamination over the JWST NIRCam/NIRSPEC G395-like bandpass (2.4–5.0 μm, and again assuming a constant 30 ppm precision) for a cloudy sub-Neptune atmosphere. The 3–5 μm region is critical for determining the carbon inventory in planetary atmospheres since large CO, CO₂, and CH₄ features are present throughout. We explore a scenario with a 14% spot and 63% faculae fraction, under the same temperature contrast as above.

We find that when not including the contamination model parameters there is a deviation of 21.7, 197, and 16.9σ respectively from the truth values, for the metallicity, C/O, and irradiation temperature along with a log Bayes factor value of 63,295, which falls under the regime of strong evidence (Trotta 2008) for stellar contamination (Table 2, Figures 7 and 8). From the 2D posterior probability distributions (Figure 8), we also find noteable correlations that emerge between the spot-covering fraction (f_s) and the gray cloud opacity (log(κ_c)) as well as a correlation between log(C/O) and the spot-covering fraction. These degeneracies are consistent with our qualitative findings in Section 2 (Equation (3)), between cloud opacity and metallicity. We also note in Figure 7 that the stellar contamination spectrum is a major contributor to the variations in the shape of the transmission spectrum within this bandpass. These results suggest that stellar contamination can have a similarly large influence on longer wavelengths due to the presence of stellar molecular spectral features (our Figure 8; see also Rackham et al. 2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. How Precisely Do We Need to Know the Stellar/Spot Properties to Remove Any Bias?

The natural next issue to arise in light of the influence of spot contamination on retrieved atmospheric properties is determining the level of required a priori spot contamination knowledge to mitigate the bias. In this section, we investigate whether having any prior knowledge of the level of stellar contamination on the host star's photosphere would mitigate atmospheric retrieval bias.

First, we again simulate a cloudy, high-metallicity sub-Neptune atmosphere with a 12% region of unocculted spot as observed under the JWST–NIRISS-like bandpass with 30 ppm precision. We choose this spot-covering fraction as it falls under the regime that is sensitive to detect stellar contamination bias (f_s > 1%) as demonstrated in Section 3.2 and Figure 6. For this exercise, we perform retrievals using two modeling approaches: (1) joint fit model, which simultaneously fits the stellar spectrum (assuming a linear combination of the the photosphere and spot spectra) and transmission spectrum for the contamination parameters under the assumption that the out-of-transit stellar photosphere has similar properties as the in-transit photosphere (e.g., Zhang et al. 2018), and (2) corrected transmission spectrum model, which only fits for the contamination parameters from the transmission spectrum itself (as in the above sections). We stress that we implicitly operate under the assumption that the stellar model is an accurate representation of an actual stellar spectrum, e.g., we assume no systematic "bias" in stellar model fits (however, we explore this assumption in Section 3.4 and Figure 10). We also assume that the absolute flux calibration uncertainty of a stellar spectrum is incorporated into our assumed photometric precision, which may not necessarily be true.

Zoom In Zoom Out Reset image size

The posterior probability distributions of f_s (fractional coverage of the spot region on the photosphere), T_s (temperature of the spot region), and T_p (temperature of the photosphere) retrieved from the disk-integrated spectrum under both approaches are presented in Figure 9, top panel. We see here that the precision in the inferred stellar contamination parameters improves by several orders of magnitude for f_s, T_s, and T_p.

However, utilizing the stellar spectrum itself (the joint fit model) shows no statistically significant improvement in the atmospheric parameter precisions (i.e., we cannot reject the null hypothesis that both samples are drawn from the same population by performing a two-sample Kolmogorov–Smirnov test (Smirnov 1939); see Figure 9, bottom panel). This important result suggests, at least in this example, that simply including a parameterization for the contamination spectrum, and marginalizing over said parameters within a retrieval, is enough to "remedy" the stellar inhomogeneity problem; no need for a priori stellar contamination knowledge (e.g., as in Pinhas et al. 2018). This result, however, may only be true if the stellar models are accurate representations of the true photosphere, which we explore in the next section.

3.4. What if the Stellar Models are Inaccurate?

We have shown that stellar spot/faculae contamination will strongly influence our ability to retrieve atmospheric properties from JWST transmission spectra of sub-Neptune-like planets orbiting M-dwarfs. We find that, for a range of plausible near-infrared spectrophotometric precisions, there is significant bias when marginalizing over planetary parameters from the transmission spectrum if not correcting for these contamination sources appropriately. The bias appears to be less significant with decreased transmission feature signal-to-noise (e.g., muting due to clouds/metallicity). We also demonstrate that our ability to adequately marginalize over the stellar contamination parameters is only as good as the degree to which stellar photosphere models match reality.

Within the context of our work, there is very little knowledge of physically reasonable coverage fractions of stellar heterogeneities as well as their distributions and occurrence frequencies throughout the photosphere of an M-dwarf in the current literature. Priors for spot-covering fractions we used were derived from I-band variability of 1% (Newton et al. 2016; Rackham et al. 2018) or an order of 0.05 or less from variability in G and R photometric bands (Rockenfeller et al. 2006a, 2006b). Spectroscopic methods to measure TiO lines emerging from spot regions of M-dwarfs have produced dramatically different estimates for spot-covering fractions when compared to methods such as Doppler imaging, varying from ∼50% to ∼10% (O'Neal et al. 1998; Barnes et al. 2011). A significant source of uncertainty with photometric variability measurements is also the assumption that starspot patterns vary as a function of rotation (Browning et al. 2010; Barnes et al. 2011); however, there is only a weak correlation between rotation and activity for M-dwarfs according to Reiners & Basri (2010). Moreover, uncertainties associated with the frequency of the occurrence of homogeneously scattered stellar spots provide difficulties in performing multi-wavelength stitching of non-simultaneous observations (e.g., Barstow et al. 2015a and Bruno et al. 2018) as we could potentially expect changes in the stellar contamination spectra in a time-dependent manner.

Several uncertainties arise in determining accurate temperature contrast ratios between stellar photospheric and heterogeneity regions, as these have previously only been modeled assuming single spots and not for uniformly covered spot regions (Barnes et al. 2011; Afram & Berdyugina 2015). However, Barnes et al. (2011) show a clear relation between stellar activity induced jitter in radial velocity observations as a function of wavelength, especially at low-temperature contrasts between the stellar photosphere and spot regions. Leveraging this relation could potentially provide one avenue for acquiring observational constraints on temperature contrasts between the photosphere and heterogeneity regions of M-dwarfs.

Numerous methods exist to estimate M-dwarf photospheric inhomogeneaities; however, it remains a challenging problem. The aim of this work is not to debate the presence of stellar contamination but to approach the problem from a "marginalization-in-retrievals" perspective. By doing so we found the following.

• 1.  
  There are degeneracies between photospheric contamination properties and transiting planet atmosphere properties within the transmission spectra.In Section 2 we gained a qualitative intuition for the expected degeneracies. From the relations in Section 2 it is apparent that stellar contamination _λ and planetary scale height H are degenerate in that both result in a "stretching" effect in the spectrum at longer wavelengths. Ignoring stellar contamination therefore will lead to incorrect estimates of the inferred terminator temperature or molecular weight of the planetary atmosphere. At shorter wavelengths, in particular, additional degeneracies emerge between the stellar contamination temperature contrast C_s, covering fraction in area f_s, and the wavelength-dependent radius z_λ. Therefore, stellar contamination parameters can alter the apparent opacity structure of the planetary spectrum in a wavelength-dependent manner.
• 2.  
  There can be significant biases in the near-infrared transmission spectra for sub-Neptune planets orbiting M-dwarfs for a range of stellar contamination scenarios.We showed that biases can exist in the atmospheric properties of interest, here the metallicity, C/O, and scale height temperature (see Figures 3 and 4). Failure to account for contamination within a retrieval would result in incorrect atmospheric parameter estimates, by several sigma, depending on the degree of the heterogeneities and the planetary feature signal-to-noise (e.g., muted versus non-muted atmospheres). We also showed through model comparison (the Bayes factor) that the evidence favoring the inclusion of a contamination parameterization increased with increasing spectrophotometric precision and contamination fractions larger than 1%. We also identified over the 2.5–5 μm spectral window notable degeneracies between the spot-covering fraction and the atmospheric C/O as well as gray cloud opacity, and also a degeneracy between the stellar spot temperature and the atmospheric metallicity.
• 3.  
  There is no substantial improvement in the atmospheric parameter constraints if we have a priori knowledge of the spot contamination characteristics.Assuming the spot parameterization and stellar models are accurate representations of the true stellar photosphere, simultaneously retrieving the contamination parameters from the stellar spectrum directly (e.g., the joint fit model) does not provide a statistically significant improvement in the precision on atmospheric parameters. Simply including stellar contamination parameterization as a set of nuisance parameters within the transmission spectrum retrieval alone is sufficient.
• 4.  
  An incorrect parameterization/representation of the contaminated stellar photosphere can result in a substantial bias on the atmospheric properties.Mimicking the differences between model and measured stellar spectra, we find that, despite including corrections for stellar contamination (see Figure 10), retrievals would result in significantly biased atmospheric properties. This is evident for the regime that is sensitive to detecting bias within the retrieval framework for planetary parameters, i.e., f_s > 1%. Therefore, our ability to correct for stellar contamination on the transmission spectrum is only as good as the accuracy of the stellar models. Understanding gaps within the stellar atmosphere models would be an important step toward alleviating this issue. Characterization of exoplanet atmospheres is a rewarding but challenging venture. Upcoming precision facilities/instrumentation will greatly improve our understanding of exoplanet atmospheres, but will also bring numerous new challenges. This work has only focused on the impact of a single "degeneracy" on our ability to infer atmospheric properties (namely the transit light source effect, e.g., Apai et al. 2018; Rackham et al. 2018). This is but one in a large realm of possible retrieval degeneracies/complications of varying degrees that include, but are not limited to, the spectral resolution linked bias (Deming & Sheppard 2017), treatment of opacities (Baudino et al. 2017), cloud model assumptions (e.g., Helling et al. 2008a, 2008b, 2016; Mai & Line 2019), choice of abundance parameterizations (Kreidberg et al. 2015; Kirkpatrick et al. 2016), and three-dimensional effects (Feng et al. 2016; Line & Parmentier 2016; Blecic et al. 2017). Diagnosing the impact and improving our understanding of these biases/degeneracies is critical to our determination of the true nature of extra-solar planetary atmospheres.

We thank our referee for their extremely thorough and timely feedback on our manuscript. We thank Jennifer Patience, Patrick Young, Evgenya Shkolnik, Mark Swain, Jeff Valenti, and Nestor Espinoza for meaningful discussions about this project. M.R.L. & A.R.I. acknowledge support from NASA Exoplanet Research Program award NNX17AB56G and HST-GO-14793 and 15255. A.R.I. also acknowledges support from the ASU SESE Summer Exploration Graduate Research Fellowship. The authors acknowledge Research Computing at Arizona State University for providing HPC and storage resources that have contributed to the research results reported within this paper. URL: https://cores.research.asu.edu/research-computing/about-rc. This research has made use of the Exoplanet Orbit Database, the Exoplanet Data Explorer at exoplanets.org (Han et al. 2014), and the pysynphot database, URL https://pysynphot.readthedocs.io/en/latest/. A.R.I. would also like to thank Edward Buie II and Sean Czarnecki for proofreading the manuscript. This paper has also made use of the pygtc plotting routine for displaying posterior probability distribution corner plots (Bocquet & Carter 2016).

Supplementary information: Figures showing posterior probability distributions of all simulation scenarios discussed in this paper are available to the reader here: doi:10.5281/zenodo.3605726.

Facility: JWST. -

Software: PyMultiNest (Buchner et al. 2014), pygtc (Bocquet & Carter 2016), pysynphot (STScI Development Team 2013).