Understanding the Effects of Systematics in Exoplanetary Atmospheric Retrievals

To study how non-Gaussian errors can bias the retrieval results, we perform retrievals on simulated data generated with and without correlation in the noise. Here, by using the same forward model to create the synthetic spectrum and to retrieve it, we remove model dependencies as much as possible and are able to examine the bias due to noise in isolation. To obtain the statistics of retrievals, we use the following procedure (also shown schematically in Figure 2).

• 1.  
  Choose the input parameters of a planet, such as radius, temperature, and metallicity.
• 2.  
  Run the forward platon (Zhang et al. 2019) model to produce the unpolluted spectrum of the atmosphere.
• 3.  
  Bin the full-resolution unpolluted spectrum to the chosen wavelength bins.
• 4.  
  Choose a noise model (independent or correlated) and noise parameters that simulate the noise of a real observation.
• 5.  
  Sample a noise realization of the binned spectrum using the noise model.
• 6.  
  Perform retrieval analysis on the simulated data.
• 7.  
  Repeat steps 5 and 6 a sufficient number of times such that the retrieved parameters can be combined to generate reliable statistics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (5 panels)

Tarball of figure set images

We use platon (version 3.0; Zhang et al. 2019), an open-source retrieval tool, as the forward and retrieval model for the transmission spectra. platon has the advantage of being extremely fast for a retrieval code ( ≲ 30 minutes per run), which is suitable for our application as we perform hundreds of retrievals with randomly sampled noise realizations. We perform this process only on transmission spectroscopy, as the geometry allows for the assumption of an isothermal atmosphere, greatly reducing the number of free parameters in our model as well as the computation time per run.

We repeat the above procedure for five cases of observation: a clear hot Jupiter, a clear hot Jupiter with offsets between data sets, a cloudy hot Jupiter, a clear hot Jupiter at higher (James Webb Telescope (JWST)-like) precision, and a warm Neptune. In what follows, we describe the forward model parameters, the noise model, and the retrieval setup. A summary of the input parameters and the assumed priors for each set of retrievals is provided in Table 1.

Table 1. Summary of Equilibrium Chemistry Retrievals and the Input Parameters and Priors Used

		Clear HJ			Clear HJ with Offset
Parameter	Name	Truth Value	Prior		Truth Value	Prior
Mp	Planetary mass (MJ)	0.73	15%		0.73	15%
Rp	Planetary radius (RJ)	1.42	15%		1.42	15%
T	Temperature (K)	1130	[300, 1500]		1130	[300, 1500]
C/O	Carbon-to-oxygen ratio	0.53	[0.2, 2.0]		0.53	[0.2, 2.0]
$\mathrm{log}a$	Log-scattering factor	0	[−0.3, 0.3]		0	[−0.3, 0.3]
γ	Scattering slope	4.3	[2.0, 5.5]		4.3	[2.0, 5.5]
$\mathrm{log}Z$	Log-metallicity	0.3	[−1, 3]		0.3	[−1, 3]
$\mathrm{log}{P}_{\mathrm{cloud}}$	Cloud-top pressure (log(Pa))	4.7	[2.5, 6.5]		4.7	[2.5, 6.5]
Offset 1	STIS offset (ppm)	⋯	⋯		−100	[−300, 300]
Offset 2	Spitzer offset (ppm)	⋯	⋯		150	[−300, 300]
Error multiple	⋯	Unity	[0.5, 2.0]		Unity	[0.5, 2.0]
Other parameters
Rs	Stellar radius (R⊙)	1.19			1.19
Ts	Stellar temperature (K)	6090			6090
Data error	(ppm)	75			75
# of runs		440			660

		Cloudy HJ			Clear HJ, High-precision			Warm Neptune
Parameter	Truth Value	Prior		Truth Value	Prior		Truth Value	Prior
Mp	0.73	15%		0.73	15%		0.0736	15%
Rp	1.42	15%		1.42	15%		0.395	15%
T	1130	[300, 1500]		1130	[300, 1500]		700	[300, 1500]
C/O	0.53	[0.2, 2.0]		0.53	[0.2, 2.0]		0.53	[0.2, 2.0]
$\mathrm{log}a$	0	[−3.0, 3.0]		0	[−0.3, 0.3]		0	[−2, 2]
γ	4.3	[2.0, 5.5]		4.3	[2.0, 5.5]		4.3	[, ]
$\mathrm{log}Z$	0.3	[−1, 3]		0.3	[−1, 3]		0.3	[−1, 3]
$\mathrm{log}{P}_{\mathrm{cloud}}$	3.0	[0.0, 6.5]		4.7	[2.5, 6.5]		5.0	[2.5, 6.5]
Offset1	⋯	⋯		⋯	⋯		⋯	⋯
Offset2	⋯	⋯		⋯	⋯		⋯	⋯
Error multiple	Unity	[0.5, 2.0]		Unity	[0.5, 2.0]		Unity	[0.5, 2.0]
Other parameters
Rs	1.19			1.19			0.683
Ts	6090			6090			4780
Data error	75			10			75
# of runs	660			400			400

Note. The stellar parameters shown are not retrieved for. For the mass and radius, the priors show the width of the Gaussian prior, relative to the input value. For other parameters, the priors are uniform in the interval. For log values, the priors are log-uniform in the interval.

Download table as:  ASCIITypeset image

To best imitate retrievals on actual observations, we choose input parameters and spectrum bins similar to HST and Spitzer observations of the canonical hot Jupiter HD 209458b and exo-Neptune GJ436b which have reliable data and have been studied in the context of retrievals (Evans et al. 2015; Deming et al. 2013; Knutson et al. 2007; Deming et al. 2011; Dittmann et al. 2009; Lothringer et al. 2018). We adopt the measured mass, radius, and temperature of each planet, and set the log-metallicity to 0.3 and carbon-to-oxygen ratio to the solar value of 0.53 (Asplund et al. 2009).

We also choose to include cloud and hazes in our model. Such aerosols have been found to be ubiquitous in exoplanetary atmospheres (e.g., Kreidberg et al. 2014; Sing et al. 2016) and produce a spectral signature that can be degenerate with other parameters (Deming & Seager 2017; Marley et al. 2014). platon accounts for clouds and hazes using a parametric model. The user can specify a gray cloud-top pressure, the atmosphere absorption below which is truncated, and an amplitude and slope in the optical end of the spectrum to account for a non-Rayleigh slope caused by Mie scattering. In summary, the aerosol opacity κ_aer is given as

$$\begin{eqnarray*}&&\begin{array}{l}P\gt {P}_{\mathrm{cloud}}:{\kappa }_{\mathrm{aer}}=\infty \\ P\lt {P}_{\mathrm{cloud}}:{\kappa }_{\mathrm{aer}}\propto a{\lambda }^{-\gamma },\end{array}\end{eqnarray*}$$

where a = 1 and γ = 4 corresponds to Rayleigh scattering from the gaseous atmosphere.

For our cloud-free simulations, we choose a low-altitude cloud-top pressure of 0.5 bar. For our cloudy simulations, we choose 0.1 mbar such that clouds obscure roughly half of the spectrum while preserving some molecular features. Similarly, we adopt a nearly Rayleigh slope of 4.3 and an amplitude of 1, indicating no excess Rayleigh scattering from the aerosol particles. We stress that since we use the same forward model in the retrieval to isolate the effects of noise from model systematics, the specific choices in parameters are not of great importance as long as they can be correctly retrieved and as long as we select a set of forward models that span a representative set of exoplanet atmospheres.

We choose a wavelength range that spans observations from the HST spectrographs most commonly used for exoplanet atmospheric investigations—specifically the Space Telescope Imaging Spectrograph (STIS) and WFC3—as well as photometric observations from the Spitzer Space Telescope. In the STIS wavelength range we follow the bins of Knutson et al. (2007), in the WFC3 wavelength we choose 33 equal sized bins between 1.01 μm and 1.64 μm (Gennaro 2018), and for Spitzer we include observations in the photometric bands of the IRAC instrument at 3.6 μm and 4.5 μm. The resulting wavelength bins are comparable to a complete panchromatic data set from current space-based observations. The resulting simulated spectra are shown in Figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To simulate the observed data, we sample multiple noise instances centered around the unpolluted spectrum, treating the simulated unpolluted spectra as a multivariate Gaussian distribution with the unpolluted depths as the mean and the reported error at each bin as the width. In addition, we adopt a nondiagonal covariance matrix with an exponential kernel to simulate correlated noise, such that the matrix element that correlates the wavelength bin at λ_i and λ_j is given by

$$\begin{eqnarray}&&{K}_{{ij}}={\epsilon }_{{ij}}{\sigma }_{i}{\sigma }_{j}\exp \left(-\displaystyle \frac{| {\lambda }_{i}-{\lambda }_{j}| }{l}\right),\end{eqnarray} \tag{ 1 }$$

where σ_i is the reported error at wavelength λ_i , and _ij is 1 for a pair of points from the same data set and 0 otherwise. We choose the scaling factor l to be the distance to the neighboring bin. We select this covariance matrix in particular because it allows for the best-fit spectrum to the WFC3 observations of HD 97658b in Guo et al. (2020) to yield a reduced chi-squared of ∼1, as opposed to 2.5 when the errors are construed to be independent and Gaussian. We choose a random noise error of 75 ppm for all instruments, which represents a moderate quality data for STIS and WFC3 but is better than average for typical Spitzer observations (Sing et al. 2016; Garhart et al. 2020). We also assume a uniform transmission across the wavelength range of each Spitzer filter in both the forward and retrieval models. In practice, we find that the two broadband Spitzer points provide little constraint on the retrieved parameters, and using the same error for all instruments does not give undue importance to the Spitzer points. For the high-precision hot Jupiter case, we use 10 ppm errors to match the best current data quality (Colón et al. 2020).

We show a few randomly selected noise realizations in Figure 5 for the Gaussian and correlated noise models. It is discernible from the bottom row of panels that the correlated noise has slightly redder residuals compared to the Gaussian noise; that is, there are less zero-crossings as the neighboring points are correlated. We also stress that overall this is a rather subtle effect; without knowing the unpolluted depths a priori to produce the residuals, from the spectrum alone it is hardly obvious that there is a distinction between the two.

Zoom In Zoom Out Reset image size

Additionally, to examine the effects of including the offset between data sets as a retrieved parameter, we create spectra with and without a fixed offset between the data sets. Namely, we add a vertical offset of −100 ppm to the STIS points and +150 ppm to Spitzer points, holding the WFC3 points constant. The specific amount of offset is a somewhat arbitrary choice. The relevant heuristic is that the offset should be exactly retrieved in the absence of degeneracy.

Using platon, setting up the retrieval involves choosing the priors and use of the statistical sampling method. We choose the priors to be comparable to a real retrieval analysis. platon accepts either a Gaussian prior or a uniform prior in a user-specified interval. We set Gaussian priors for the planet radius and mass, as these are often constrained via other methods such as radial velocity and transit measurements prior to observing the transmission spectrum. The Gaussian prior is centered around the input value and with a standard error of 15%. This precision is comparable to or slightly overestimates the typical uncertainty in the measurement of mass via the radial velocity method and provides sufficient tolerance for the retrieved value to deviate from the input value, if necessary. For the instrumental offsets, to ensure that the prior is broad enough, we choose a uniform prior offset to be two- and three-fold of the offsets. For all other parameters we opt for uniform priors that are as wide and uninformative as possible and adopt the full parameter range supported.

For now, we only choose to do retrievals with equilibrium chemistry, where the composition of the atmosphere at a given temperature and pressure is dictated by the the global elemental abundances set by metallicity and the carbon-to-oxygen ratio. While disequilibrium chemistry is indeed expected for planets below T_eq ≤ 1200K, most of its effects take place below the altitude that is typically probed by transmission spectroscopy and have no easily discernible effect on the spectrum at the data precision of current instruments. The actual discrepancy in the relevant pressure range (∼1 mbar) is smaller than the uncertainties we can currently obtain (e.g., Line et al. 2011; Miller-Ricci Kempton et al. 2012; Fortney et al. 2020).

platon supports either Markov chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) or nested sampling methods for the posterior estimation (Feroz & Hobson 2008). We note that while both are statistically robust and widely used, we observe a minor discrepancy in the resulting posteriors between the two methods, in which the posteriors estimated by MCMC tend to be slightly broader than those estimated by nested sampling. This does not pose a major issue for this work inasmuch as we are concerned with biases due to data idiosyncrasies, and we consistently use one algorithm across our analyses. Nevertheless, we draw attention to this point as it requires extra scrutiny when comparing retrieval results. We choose nested sampling as it is known to perform better in estimating multimodal or oddly shaped distributions.

Here we present the effects of correlation in data on the retrieval overall. To do this, we must reduce a retrieval result into a simpler metric. In retrievals on actual data this is commonly done by using the reduced chi-squared between the observed data and the best-fit or median spectrum. Here, as we begin from a known simulated ground truth, we can also compare the retrieved posterior directly to the input values to measure the accuracy of the retrieval by using a probability integral transform (PIT), which is the cumulative distribution function evaluated at the input value. To do this, in each retrieved posterior distribution, we draw an isolikelihood contour of the input parameters and sum the relative weights contained within, producing a confidence interval between 0 (the input was the most likely sample) and 1 (the least likely). The distribution of the PIT values should follow a uniform ${ \mathcal U }(0,1)$ were the retrievals accurate.

Our main finding is that, on average, correlation in the data allows for overfitting the spectrum, thereby weakening the overall accuracy of the retrieval. In other words, the best-fit spectrum is more likely to achieve a reduced chi-squared lower than unity in the presence of correlated noise, while simultaneously the retrieved posterior distribution rules out the input with higher significance.

In Figure 6, we present the comparison between Gaussian and correlated noise for the clear hot Jupiter case. For both retrievals with Gaussian (black) and correlated noise (blue), we show the distributions of: χ² between the unpolluted spectrum and the data instances, the reduced chi-squared between the best-fit spectrum and the data instances, and the PIT values showing the accuracy of retrieval, as described above. A few observations can be made:

• 1.  
  The fit ${\chi }_{{\rm{r}}}^{2}$, or the goodness of fit of the retrieval, on synthetic data is on average skewed to better than unity in the presence of correlated noise. In other words, it is more likely that the retrieval will overfit the data with forward models.
• 2.  
  The accuracy of the retrieval, shown as the PIT value, on the other hand, is worse in the presence of correlated noise. We also show the cumulative distribution in Figure 7 to demonstrate the worsening of the accuracy. The Kolmogorov–Smirnov (K-S) statistic between the two cases is 0.19.
• 3.  
  For the normal noise case, even in the absence of the correlated noise, the retrieval accuracy is close to but slightly worse than the expected uniform distribution. This minor discrepancy is likely due to the fact that the retrieved posterior distribution is already non-Gaussian for a few parameters (discussed in Section 3.2), as well as due to degeneracy between certain parameters.
• 4.  
  While the goodness of fit is on average better in the presence of correlated noise, it is not the case that the distributions of fit ${\chi }_{{\rm{r}}}^{2}$ are so discrepant that one can deduce the presence of correlated noise from the goodness of fit alone. That is, a given overfit spectrum can plausibly be construed as either a consequence of correlated noise or as an unlucky instance of Gaussian noise that happens to lie at the tail of the ${\chi }_{{\rm{r}}}^{2}$ distribution. As such, we stress that the effect of correlated noise is manifested statistically, and no individual value of ${\chi }_{{\rm{r}}}^{2}$, good or bad, is uniquely a diagnostic of correlated noise in a single retrieval instance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given that the correlated noise degrades the overall accuracy of retrievals, it is necessary to then look at which parameters are affected. We do this by marginalizing the posterior distribution over each parameter, as is typically done in retrieval analyses. The effect on each marginalized distribution can be twofold—the mean can be biased, the estimated error can be affected, or both. Either a shift in the mean away from the input parameter or an underestimation of the error can worsen the accuracy of the retrieval. As such, we examine the retrieved mean and the retrieved error separately for each parameter.

The distribution of the retrieved means is shown in Figure 8. For the case in which noise is independent (black), the retrieved means form clean normal distributions around the input values (red) for most parameters, as expected from the central limit theorem. The two exceptions are the C/O ratio and the cloud-top pressure. This is most likely due to the fact that the retrieved distributions for these parameters are not Gaussian in the first place. For cloud-top pressure, the retrieved distribution is at best a flat distribution with a lower bound, ruling out a cloudy atmosphere as per the clear atmosphere in the input used. For the C/O ratio, we suspect that the distribution is skewed due to the increasing influence the parameter has over its range. That is, as one sweeps through the C/O ratio, the spectrum changes more rapidly over the range above the solar value of 0.53, and thus the means are naturally skewed to values lower than the solar value where there is a greater density of near-consistent solutions.

Zoom In Zoom Out Reset image size

The presence of correlated noise has a few interesting effects on the retrieved means. First off, the error multiple parameter is biased to less than unity. This intuitively follows from the global result that correlated noise allows for overfitting of the spectrum, tricking the error multiple parameter to believe that the error bars are overestimated. This means that the error multiple parameter is more likely to behave pathologically in a situation where one may expect it to be useful, such as if the reported error bars truly were underestimated due to unknown and unaccounted systematics. The retrieval instead selects a less-than-unity value of the error multiple, incorrectly implying that the data precision is better than initially reported. This is possible if the domain of input parameters and the forward model can still reproduce the spectrum polluted with systematics.

In the correlated noise case, the retrieved means generally show a wider distribution to varying degrees for each parameter. Specifically, the radius, mass, and temperature are the most affected, while the effect is the least pronounced for metallicity and the cloud-top pressure. This result may be explained by considering the wavelength scale the former three parameters have on the spectrum. Mass and temperature affect the scale height of the atmosphere, which affects the overall vertical extent of the transmission spectrum. The radius affects the baseline transit depth as well as the scale height. These are global parameters in the sense that the transit depths in all of the bins are affected together. As such, a wavelength-dependent correlation can bias these parameters. On the other hand, metallicity, while it also affects the scale height (via the mean molecular weight), directly controls the individual transit depths. This has a more local effect in that it changes the actual shape of the spectrum.

The distribution of the retrieved errors is shown in Figure 9. The effect of the correlated noise is clearly visible for all of the parameters in that the retrieved error bars show a tendency to be underestimated. For instance, the retrieved error on log-metallicity is on average underestimated by ∼0.2 dex. While this disparity is smaller still than typical constraints, it is worth bearing in mind as this is a statistical effect; the spread over the retrieved error is by itself broad enough that the actual effect of a given instance can be much larger than this value. Additionally, when JWST allows for the precise measurement of metallicity, this level of uncertainty may not be negligible when one considers analyzing archival data simultaneously. The same consideration applies to other parameters. As such, in this context we suggest that the retrieved constraints for parameters, in the face of the potential for correlated noise, are best understood as lower limits.

Zoom In Zoom Out Reset image size

In this section, we present our results for planet scenarios other than the baseline clear hot Jupiter case to understand the sensitivities of our results to various system and data set parameters. We show histograms of the retrieved mean and retrieved standard error for each parameter in figure sets for the remaining planetary scenarios (the hot Jupiter with offsets, cloudy hot Jupiter, high-precision hot Jupiter, and warm Neptune). We generally find that the main results stated so far hold true for all cases: correlated noise causes both overfitting in ${\chi }_{{\rm{r}}}^{2}$ and worsening of the accuracy of retrieval (i.e., larger PIT values). This point is summarized in Figure 10, in which we show the medians of the ${\chi }_{{\rm{r}}}^{2}$ and PIT value distributions for each planet realization, i.e., distilling down the results of Figure 6 and the like to values quantifying the peak and the spread.

Zoom In Zoom Out Reset image size

To further quantify this point, we perform a K-S test for the goodness of fit and the accuracy-of-retrieval metrics to measure the discrepancy between the results of Gaussian and correlated noise. In Table 2, we show the K-S statistics, D, for the fit ${\chi }_{{\rm{r}}}^{2}$ and the PIT value representing the accuracy of retrieval. We find that the clear hot Jupiter happens to be the best case for the smallest discrepancy of accuracy of retrieval between Gaussian and correlated noise, and that other cases generally result in further discrepancy between results with Gaussian and correlated noise.

Table 2. K-S Statistic of ${\chi }_{r}^{2}$ and PIT Values

	${\chi }_{{\rm{r}}}^{2}$	PIT value
Case	G-C	${ \mathcal U }$-G	${ \mathcal U }$-C	G-C
Clear hot Jupiter	0.28	0.19	0.60	0.47
Clear hot Jupiter w/offset	0.36	0.19	0.67	0.58
Cloudy hot Jupiter	0.27	0.12	0.54	0.48
Clear hot Jupiter, high precision	0.28	0.04	0.48	0.47
Warm Neptune	0.35	0.16	0.58	0.52

Note. The K-S statistic, D, measures the maximum vertical discrepancy between the cumulative distributions (see Figure 7), of the goodness of fit, and the retrieval accuracy for each planet scenario. The first column shows the two-sample D between the distributions of fit ${\chi }_{{\rm{r}}}^{2}$ for the Gaussian and correlated noise. The next two columns contain the D between the expected uniform distribution and the PIT value distributions from retrievals with Gaussian and correlated noise, respectively. The final column shows the D between the two distributions. In all cases, the discrepancy is due to overfitting and worsening of retrievals.

Download table as:  ASCIITypeset image

3.3.1. Bias in the Non-Rayleigh Scattering Slope

3.3.2. Retrieving Offsets

3.3.3. Effects of Clouds

3.3.4. Effects of Higher Precision

We find that in the hot Jupiter retrieval with high-precision (10 ppm) data, the broad conclusions again still hold. Correlated noise leads to an underestimation of retrieved uncertainty for all parameters. Compared to other cases however, correlated noise does not shift the retrieved means as much, which is to be expected since every noise instance only has a minor deviation from the unpolluted spectrum (specifically a factor of 7.5 times smaller than in our baseline case), even with correlation.

Comparing the high-precision case to the baseline case with 75 ppm errors, we find that, naturally, both the estimated means are retrieved closer to the input values and the retrieved parameter uncertainties are concurrently smaller. Interestingly, the uncertainties shrink more than the means approach the input values; consequently, in the high-precision case, the retrieval more readily rules out the input. This is shown in Figure 11, in which the retrieved means are normalized by their retrieved uncertainty to show the number of standard errors within which the input value for each parameter is retrieved. The high-precision case (dashed line) actually has more retrievals further from the input value when normalized. We find an identical trend for the retrievals with correlated noise.

Zoom In Zoom Out Reset image size

Additionally, we find the PIT value distribution for the retrievals with Gaussian noise much more successfully follows the ideally expected uniform distribution, with a K-S statistic of D = 0.04 (see Table 2 and Figure 12). Given that the accuracy of the retrievals for all other cases (with 75 ppm error bars instead) are at least somewhat discrepant from the expected distribution even for the Gaussian noise retrievals, this result suggests that, even for hot Jupiters, 75 ppm error bars are too large to assume a priori that the retrieved posterior will follow a multivariate Gaussian. This has implications for parameter estimation methods that need this assumption of a Gaussian posterior, such as optimal estimation or some of recent machine learning-based retrievals (Line et al. 2013; Cobb et al. 2019). These methods require that the data uncertainty is small enough such that the forward model behaves linearly over the parameter uncertainties. Retrievals using these methods must be trusted only when the data has an exceptional signal-to-noise ratio.

Zoom In Zoom Out Reset image size

A major compounding issue is that, when retrieving on real data, model assumptions and unknowns contribute to and act as systematic errors in addition to the data systematics themselves. In short, bad data are degenerate with bad models. In our study we generated the synthetic observations using the identical forward model as that used in the retrieval in order to minimize any model-dependent effects and to isolate the effects of data systematics. In interpreting real data, the fact that our forward models are a simplified incomplete representation of complex atmospheric phenomena will act as a source of systematic error that will remain pervasive, even if the observed data were perfect and free from their own systematics. We therefore remain open to the possibility that the observed examples of potential systematic noise in the data are in fact due to unaccounted for obscure physics.

For the same reasons discussed above, this will adversely affect high S/N observations in particular, in which the fine (and the not-so-fine) details of the model become discernible. There has recently been a growing body of work that studies the biases incurred by model assumptions and parameterization. To list a few examples for demonstration, MacDonald et al. (2020) performed one-dimensional retrievals on three-dimensional synthetic spectra to show that the retrieval biases the limb temperature to few hundred Kelvins cooler than the actual day–night mean temperature. Lacy & Burrows (2020) extended this study to cloudy atmospheres, finding that the presence of aerosols exacerbate the biases induced by three-dimensional effects when not accounted for. Changeat et al. (2019) found that using a vertically constant chemical abundance profile may no longer be sufficient to fully capture signatures of disequilibrium processes in the spectrum. Perhaps most inconspicuously, Barstow et al. (2020) compared and performed cross-retrieval between retrieval codes developed by three groups and found that JWST-quality data is now sensitive to rather rudimentary model unknowns such as the line lists used to generate the opacities and the precision of fundamental constants used.

In our framing of describing biases, these results generally effect shifts in the estimated means of the retrieved parameters. Incorporating our conclusion that correlated noise generally leads to an underestimated error of the retrieved parameters means that biases due to model limitation now strike with a stronger statistical significance. Furthermore, the wavelength-dependent effect of both missing physics and systematics now leave the possibility for degenerate interpretations.

To better understand the significance of our results, it would be useful to consider the different sources of how systematics can arise and evaluate their prevalence especially in the context of future space missions such as JWST. While other sources of systematics are possible, such as starspots (Rackham et al. 2018; Bruno et al. 2020), inaccurate orbital parameters, or time-dependent telluric contamination in ground-based observations, the key source of systematics that we will discuss here is instrumental. However, we remind the reader that our formulation of correlated noise can be generalized to any effect that results in wavelength-dependent correlation.

While instrumental systematics are expected to be ubiquitous to some extent, the exact magnitude of their effect in generating wavelength-correlated noise has not been fully understood. Yet a handful of observations exist that hint at the existence of such systematics. Colón et al. (2020) argued that, with current facilities, these systematics are visible at the highest level of precision (∼15 ppm), inferred from an unusual behavior of the residuals in the H₂O band. In the spectrum of HD 97658b in Guo et al. (2020), our motivating example in Section 1, we inferred from the inability to fit the data as well as the fact that no obvious physics were missing that there must be some systematics present, even at a lower precision (∼25 ppm). These examples indicate that some wavelength-correlated noise must be present, unless there is unaccounted for physics in the retrieval model.

There is some reason to surmise that these systematics are more prevalent (or, at least, more noticeable) in the case of bright host stars. A brighter host star allows for a better S/N and higher precision and thereby naturally makes the presence of these systematics more conspicuous compared to a lower S/N data. Additionally, even at the same data quality, a brighter host star requires fewer stackings of observations; for a dimmer host star, by contrast, the number of stackings required to achieve the same S/N naturally averages out any nonrepeatable correlated noise. Finally, as alluded to in Section 1, given that instrumental systematics can also behave differently with bright sources, it is not out of the question that there is a separate effect at play here beyond the S/N which may persist through multiple observations in a repeatable fashion.

Comparing the WFC3 spectra of GJ 1214b (Kreidberg et al. 2014) and of HD 97658b (Guo et al. 2020) illustrates this point. While both planets are comparable sub-Neptunes with featureless spectra and have a similar level of precision, the spectrum of the latter displays an unusual upward trend in transit depth in the redder end and other wavelength-correlated residuals throughout the WFC3 bandpass. The relevant difference here may be the host star brightness (9.8 versus 6.2 in the J-band magnitude, respectively). The spectrum of GJ 1214b is the combination of stacking 15 transits, whereas that of HD 97658b has 4. Further, the spectrum of HD 97685b presented in Knutson et al. (2014), which only had the first two visits, shows the most obvious possible example of correlated noise due to systematics.

These types of systematics may be even more pernicious for future high S/N observations from JWST for a few reasons. First, as the noise floor is lower, correlated noise will be relatively more prominent even if it actually manifests at weaker levels. Then one can no longer reliably assume that the observed noise is strictly photon-dominated. This requires an additional step of modeling out now wavelength-dependent systematics during the data reduction, which is necessarily (although perhaps not practically) incomplete. Given that, secondly, the high S/N per transit means that stacking will be unnecessary for most targets. As such, nonrepeating systematics do not get averaged out. Third, we have demonstrated that higher precision leads to biases of stronger significance. This is true even in the absence of systematics in the sense that a retrieval will be more sensitive to the observational instance. Fourthly, we can predict that our understanding of the characteristics of JWST instruments and their appropriate data reduction tools will be only partially correct, at least during the initial few cycles of JWST before practical experience accumulates. Finally, as planets around bright host stars allow for achieving high S/N, they will make attractive targets for JWST observation. However, if the above intuition that bright host stars can exacerbate instrumental systematics is true, it adds another dimension to consider when selecting targets for observation, in addition to the S/N.

Given that this is the case, it would be worthwhile to put the above heuristic that bright stars bring about correlated noise to a more formal test. This can be accomplished with JWST if the correlation strength can actually be measured, and by marginalizing over the magnitude of the host star to obtain a trend. While this would not comprise a main scientific objective of any program, correlation strength is a parameter we can try to measure for all observations, so this is a test we can perform at no extra cost in observation time.

platon originally supports equilibrium chemistry retrievals only. Using this method the molecular abundance of each species at a given temperature and pressure is set by the metallicity and carbon-to-oxygen ratio. A popular alternative method to constrain the chemistry is to use free retrievals, in which the abundances of each species are allowed to vary independently. This accounts for any nonequilibrium chemistry effects in the atmosphere, brought on by vertical mixing or photochemical interactions. To establish how wavelength-dependent systematics or outliers can bias certain species, we implement free retrievals in platon and perform some basic tests in addition to the suite of retrievals already presented that used equilibrium chemistry models.

platon already naively supports inputting custom chemical profiles to its forward models, but only accepts equilibrium chemistry parameters during retrievals. We extend its capability by allowing it to accept custom chemical abundances during retrievals as well. As such, all other details regarding how the spectrum is calculated during the retrieval remain exactly the same as the original implementation in platon. We assume that each species has a vertically fixed mixing ratio. This is a reasonable approximation over the pressure range probed by transit spectroscopy at the present data quality, and most current one-dimensional free retrieval codes parameterize the composition using this assumption (Caldas et al. 2019). JWST-quality data may merit a more complex prescription, such as a two-part vertical abundance profile (Changeat et al. 2019), but for now we do not consider these possibilities. Additionally, for these tests we remove the simulated Spitzer data points from our synthetic spectra to limit the wavelength range, thereby reducing the choice of necessary species to be included. We assume an H/He-dominated background atmosphere. We include H₂O, Na, K, TiO, and VO as they are the primary detectable species in the remaining wavelength and temperature range; we are mostly interested in how differently species with broad absorption features (e.g., H₂O) versus species with narrow ones (e.g., Na) can be biased.

We show a free retrieval on two observational instances of the same baseline hot Jupiter as described in Section 2, where we did not add correlated noise or any instrumental offsets. Figure 17 shows the two random input data realizations and the best-fit spectra, and Figure 18 shows the retrieved posteriors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The point of interest here is how the retrieved abundances compare for species with broad spectral features (H₂O, 1–2 μm) and with a narrow feature (Na, 0.58 μm). As one may expect, for both of the two random spectra, the retrieved posterior distribution for water shows a tighter constraint than that for Na, which is dictated almost entirely by one data point. Consequently, between the two posteriors as well, the retrieved distributions for Na show little overlap, ruling out each other by ≥ 2σ. These are simply two randomly drawn samples, but it demonstrates the point that measurements of water abundance or metallicity are more robust compared to that of say Na or K abundances, because of the impact the former parameters have across a broader wavelength range.