A new statistical method for characterizing the atmospheres of extrasolar planets

TL;DR

This paper introduces a Bayesian statistical method to improve the characterization of exoplanet atmospheres by better accounting for systemic errors and uncertainties in observational data.

Contribution

The paper presents a novel Bayesian approach for exoplanet atmospheric analysis that enhances robustness against data variability and underestimation of errors.

Findings

More accurate probability distributions of atmospheric parameters.

Reduced sensitivity to data variability.

Better reflection of uncertainties in parameter estimates.

Abstract

By detecting light from extrasolar planets,we can measure their compositions and bulk physical properties. The technologies used to make these measurements are still in their infancy, and a lack of self-consistency suggests that previous observations have underestimated their systemic errors.We demonstrate a statistical method, newly applied to exoplanet characterization, which uses a Bayesian formalism to account for underestimated errorbars. We use this method to compare photometry of a substellar companion, GJ 758b, with custom atmospheric models. Our method produces a probability distribution of atmospheric model parameters including temperature, gravity, cloud model (fsed), and chemical abundance for GJ 758b. This distribution is less sensitive to highly variant data, and appropriately reflects a greater uncertainty on parameter fits.

Tables

Table 1. Table 1: References: (1) - Vigan et al. ( 2016 ) (2)- Janson et al. ( 2011 ) , recalibrated for 15.5 pc

Filter	$\lambda (\mu m$)	Absolute Magnitude	Ref.
Y2	1.022	19.19$\pm$0.20	1
Y3	1.076	18.43$\pm$0.10	1
J2	1.190	19.06$\pm$0.25	1
J	1.250	17.60$\pm$0.20	2
J3	1.273	16.83$\pm$0.18	1
CH4s	1.580	17.76$\pm$0.12	2
H2	1.593	16.59$\pm$0.12	1
H	1.650	18.18$\pm$0.20	2
H3	1.667	18.88$\pm$0.42	1
CH4l	1.690	$\ge$18.68	2
K${}_{\text{c}}$	2.098	17.14$\pm$0.20	2
K1	2.110	17.03$\pm$0.21	1
K2	2.251	17.78$\pm$0.35	1
L’	3.776	15.0$\pm$0.1	2
M${}_{\text{s}}$	4.670	$\ge$13.1	2

Equations

$$\textrm{Prob}(\theta|X_{1},X_{2},X_{3}...X_{n})=\int_{p=0.0}^{p=1.0}\prod_{i=1}^{n}\big{[}[p\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{given}})+(1-p)\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{large}})]\cdot\textrm{Prob}(\theta)\big{]}\,dp$$

$$\int\limits_{\textrm{all models}}^{\theta}\int_{p=0.0}^{p=1.0}\prod_{i=1}^{n}\big{[}[p\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{given}})+(1-p)\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{large}})]\cdot\textrm{Prob}(\theta)\big{]}\,dp\,d\theta=1$$

$$\textrm{Prob}(\theta|p,X_{1},X_{2},X_{3}...X_{n})=\prod_{i=1}^{n}\big{[}[p\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{given}})+(1-p)\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{large}})]\cdot\textrm{Prob}(\theta)\big{]}$$

$$\textrm{Prob}(\textrm{ith good})=\frac{\int\limits^{\theta}\int_{p=0.0}^{p=1.0}p\cdot\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{given}})\cdot\prod_{k=1,\ k\neq i}^{n}\big{[}[p\cdot\textrm{Prob}(X_{k}|\theta_{k},\sigma_{\textrm{given}})+(1-p)\cdot\textrm{Prob}(X_{k}|\theta_{k},\sigma_{\textrm{large}})]\cdot\textrm{Prob}(\theta)\big{]}\,dp\,d\theta}{\int\limits^{\theta}\int_{p=0.0}^{p=1.0}\prod_{k=1}^{n}\big{[}[p\cdot\textrm{Prob}(X_{k}|\theta_{k},\sigma_{\textrm{given}})+(1-p)\cdot\textrm{Prob}(X_{k}|\theta_{k},\sigma_{\textrm{large}})]\cdot\textrm{Prob}(\theta)\big{]}\,dp\,d\theta}$$

$$

Full text

††pubyear: 2017††pagerange: A new statistical method for characterizing the atmospheres of extrasolar planets–A.4

1 Introduction

By comparing the measured spectral energy distributions of extrasolar planets to atmospheric models, it is possible to constrain exoplanet bulk physical properties such as temperature, cloudiness, and metallicity. Exoplanet atmospheric characterization necessarily pushes instruments beyond their design capabilities, so systematic effects often dominate over atmospheric signatures. Underestimated errors and erroneous measurements have lead to incorrect inferences of exoplanet properties (Burrows, 2014). Evidence of underestimated errors exist. For example, Hansen et al. (2014) note discrepancies in multi-epoch Spitzer photometry of transiting exoplanets. Yet more measurements at greater resolution are needed for accurate characterization (Croll et al., 2010; Lupu et al., 2016; Line et al., 2016). A rigorous method to evaluate and recognize erroneous measurements is necessary to provide accurate interpretations of exoplanet properties.

In this paper, we study the properties of a direcly-imaged companion, GJ 758b, using a Bayesian technique that simultaneously evaluates the correctness of individual photometry points while fitting the data to a suite of atmospheric models. GJ 758b is a T dwarf companion separated by 30 AU from a Solar type star with an age of 5-9 Gyr (Thalmann et al., 2009; Janson et al., 2011). This object has been characterized before, however, previous work has demonstrated that its photometry is difficult to fit with atmospheric models or spectra of analog objects (Vigan et al., 2016; Nilsson et al., 2017).

In Section 2 we provide literature photometry data, information on the models we use, and a basic model fit. In Section 3 we explain why the photometry on GJ 758b gives us reason to believe that some of the data is erroneous. In Section 4 we describe our Bayesian approach to evaluate the quality of each data point, which is based on similar work to estimate the Hubble constant in the context of discrepant data (Press, 1997). We compare GJ 758b’s derived atmospheric properties with and without the Bayesian approach.

2 Data and Model Fitting

2.1 Photometry

We compiled literature data on GJ 758b from three sources: the VLT planet finder, SPHERE (Vigan et al., 2016), Gemini North using NIRI, and Keck observatories using NIRC2 (Janson et al., 2011). Vigan et al. (2016) provide absolute magnitudes calibrated for a distance of 15.76 pc, while Janson et al. (2011) provide absolute magnitudes calibrated for 15.5 pc, which we recalibrated for 15.76 pc to maintain consistency with Vigan et al. (2016). All absolute magnitudes are tabulated in Table 1. The distance change accounts for a difference of 0.02 magnitudes, which in the cases of L’, and Ms, does not change the magnitude within the known significant figures. We used filter curves from the VLT, Gemini North, and Keck websites - and used cryogenically measured filter transmission curves when available to calculate effective wavelengths and model photometry for section 2.3.

2.2 Models

We calculate model spectra assuming radiative–convective and chemical equilibrium, and including sulfide and salt condensates as described in Morley et al. (2012, 2014a). Various updates to the chemistry and opacities will be described in detail in an upcoming paper (Marley et al. in prep.) Briefly, the methane line lists have been updated using Yurchenko & Tennyson (2014) and the alkali line lists have been updated to use the results from Allard et al. (2005). Chemical equilibrium calculations are based upon previous thermochemical models (Lodders & Fegley, 2002; Visscher et al., 2006; Visscher, 2012), and have been revised and extended to include higher metallicities.

The model grid parameters we vary include temperature, gravity, metallicity, and cloud sedimentation parameter fsed (see Ackerman & Marley (2001)). Our model grid contains 660 models, comprising temperatures of 450, 475, 500, 525, 550, 575, 600, 625, 650, 675, 700 K, surface gravities, of 30, 100, 300 and 1000 m/s2, metallicities of [M/H]=0, 0.5, and 1.0, and cloud sedimentation efficiency of fsed=1, 2, 3, 5 and cloud-free.

We linearly interpolate within this grid in order to better sample the error distribution of our parameters, increasing our effective grid size from 660 models to 17000 interpolated models. We have shown in Skemer et al. (2016) that this interpolation minimally effects the results for similar ranges of parameters.

2.3 Model Fitting

To determine how well a model fits, we estimate model absolute magnitudes using the model SEDs and the filter transmission curves described in Section 2.1. We compare the model magnitudes to the measured magnitudes in Table 1 to analyze the goodness of fit per model.

We used a Vega model from Rieke et al. (2008) which we smoothed by 0.3%, as recommended, before using it to generate zero point fluxes.

We assume a radius for GJ 758b of 1.05 $R_{J}$ based on gas-giant planet evolution models which are well constrained to within 5% by the age, mass, and separation of the system (Fortney et al., 2007). When not using a prior on the radius, it tended towards unphysically small radii without significantly improving the fit.

We discarded models that were inconsistent with 3$\sigma$ upperlimits in the CH4l and Ms filters (and otherwise we did not include these upperlimits in our new statistical method in Section 4).

Using a standard Bayesian fit for Gaussian errors ($P\propto\frac{1}{\sigma}e^{-\frac{{\chi}^{2}}{2}}$) normalized over all models, we found best fit parameters which we detail in Figure 2, and a best fit model shown in Figure 1.

Notably, no model was able to simultaneously fit the CH4s and H2 photometry, which overlap in wavelength as shown in Figure 1. The best fit model reached $2.6\sigma$ for CH4s and $2.5\sigma$ for H2. Models which fit CH4s better fit H2 worse, and vice versa. Also, no model was able to fit K2: the best models had K2 magnitudes $>8\sigma$ off. The best fit from the interpolated grid of models, had a reduced chi-square of 7.01.

3 Bad models or bad data?

3.0.1 H2 and CH4s

The CH4s and H2 filters overlap in wavelength, but report different flux measurements. We were not able to simultaneously fit the CH4s and H2 photometry within $2.5\sigma$ with self consistent atmospheric models, and a model that fits would require a very strange shape.

To demonstrate the problem, we constructed an arbitrary box model (free parameters: cut-on wavelength, cut-off wavelength and height) that best fits the CH4s and H2 photometry (see Figure 3). This box could fit both points within 1 $\sigma$, due to its tuned shape, which sits on a low part of the H2 filter curve, and a high part of the CH4s filter curve. However, the box model is far too discrepant from the self-consistent atmosphere models to be plausible. The box model function has a reduced chi square of 0.89 and the atmosphere model, using only H2 and CH4s, has a reduced chi square of 6.88. The atmospheric model was chosen for the best simultaneous CH4s and H2 fit, with parameters $T=625K,g=300m/s^{2},f_{\rm sed}=5.0,[M/H]=0.0$ and a radius R=$0.85R_{J}$ chosen for best fit. Despite allowing the radius to vary to unphysically low amounts, we were still unable to find a model to fit both points simultaneously.

Since no realistic model could fit both H2 and CH4s, we conclude that these data points are inconsistent with themselves.

3.0.2 Flux in K2 band

We were unable to find any physical models that could fit K2 within 3$\sigma$. Flux is not expected at this wavelength (Morley et al., 2014b). The spectrum in the model spectra at the wavelengths of the K2 filter ( 2.2-2.3 micron) is dominated by methane opacity. For the model spectra to be brighter in K2 to match the observed data, methane would need to be present at lower abundances than predicted. Disequilibrium carbon chemistry, which can decrease the methane abundance, is observed in brown dwarfs and exoplanets of similar temperatures (Hubeny & Burrows (2007); Saumon et al. (2006); Konopacky et al. (2013)). However, for GJ 758b, decreasing the methane abundance to match the K2 photometry would also change other regions of the model spectrum. For example, the observed flux within the H3 filter ( 1.6 micron), which is also dominated by methane opacity, would significantly increase. By this argument, we conclude that the the observed K2 photometry is not consistent with the rest of the spectrum, unless there is an unexpectedly substantial problem with the methane line list at near-infrared wavelengths.

Assuming the K2 photometry is incorrect, including K2 in a statistical analysis will artificially constrain which models are preferred. Removing this photometry point, or increasing its formal error, would increase the quality of our atmospheric fitting analysis.

K2 photometry could be correct and the models wrong, but the result of directly using highly variant photometry prevents productive model fitting of the rest of the data.

4 Press and Kochanek’s Method

To address the problem of data that is impossible to fit, we implemented Press and Kochanek’s method (Press, 1997), which considers the possibility of underestimated errorbars.

4.1 Background

Press and Kochanek’s method was originally used to determine the Hubble constant, using varied and contradictary data, where some select data points with small errors weighted the averaged value for the constant disproportionately compared to other measurements, despite being rather far from most other measurements. The method was introduced in order to address this problem rigorously, without discarding data entirely or arbitrarily.

To do this, Press (1997) defined a probability function where (1) individual data points have an unspecified probability of being “correct” or “incorrect”, (2) data that are “correct” are assumed to have Gaussian-distributed errors as defined by the data’s error bars, and (3) data that are “incorrect” are assumed to have Gaussian-distributed errors that are increased by some factor, which is defined at the outset of the analysis. By marginalizing over the probability of each data point being “correct” or “incorrect”, the technique is able to rigorously reject erroneous data. A mathematical representation of this method is described in the appendix, Appendix A.

In (Press, 1997)’s original implementation, the implicit assumption was that there is a single value of the Hubble constant that does not vary between measurements. Discrepancies between measurements are then internally inconsistent. When applying the method to SED fitting, the analogous assumption is that one of our model SEDs is correct and that the astrophysical source does not vary between measurements. The assumption that one’s model grid contains a “correct” model is implicit in all Bayesian model fitting, however we acknowledge that it is an additional assumption not considered by Press 1997. In the case of GJ 758, the additional assumption is justified by multiple measurements at different wavelengths, which together, conflict with models so badly that they cannot all be correct. Although our use of the method is different, the mathematical formalism is similar, as described in Appendix A.

Overall, the method serves to desensitize the analysis to data that is hard to fit, which would ordinarily skew the result.

4.2 Application to GJ 758b

Given the difficulty of fitting CH4s and H2 simultaneously, as well as fitting K2, we consider the possibility that the errorbars for our data set from Table 1 are underestimated by a factor of 2 (note that increasing this fixed parameter has a negligable effect on the final results as shown below). All photometry points are treated the same—for each photometry point, we assume, with a probability that is a free parameter, that the error bar might be underestimated by a factor of 2. We evaluate the multi-dimensional posterior probability on a grid that includes the atmosphere model parameters, and the probability that each photometry point is correct or incorrect.

We plot the marginalized posterior probabilities of this method in Figure 4. We note that the variance in the parameter PDF is larger than without the method, which suggests that we have removed the constraints of incorrect photometry. Wall & Jenkins (2012) provides a probability of correctness per data point, which we plot in Figure 5. This formula is detailed in Appendix A.4, but in short it measures the consistency of each data point compared to the rest of the set. The best fit model using this technique ended up being the same as the best fit model from Section 2.3 (shown in Figure 1).

4.3 Method with prior

As shown in Figure 5 the method’s estimate that each individual photometry point is correct is rather low: most photometry points have a $\sim$30% chance of being correct, while the outlier points are assigned lower probabilities. If we want to take the optimistic viewpoint that most observations result in correct error estimates, we can apply a prior that gives the observer the benefit of the doubt. The prior fixes the probability that an individual photometry point is correct, which requires more substantial improvements to fits for the increased errorbar to become relevant. The mathematical formalism is presented in Appendix A.3.

4.4 Method with varied parameters

We try using a variety of priors, and changing the factor which “incorrect” errors are increased by, to see how important method parameters are in changing the atmospheric parameter fit distributions. We plot our temperature fits and uncertainties in Figure 6 for different paremeters used in the method. Note that the least aggressive use of the method (assume that 90% of the data are correct) still has a considerable effect on the temperature probability distribution.

4.5 Results for GJ 758b

The results of our analysis suggest that GJ 758 b’s derived parameter errors increase when some photometry errors are assumed to be underestimated. With a prior of 90% correctness on the data (i.e. 1 or 2 incorrect data points) - the method decisively concludes that K2 has an underestimated errorbar, and that between CH4s and H2, CH4s is more likely to have an underestimated errorbar.

The use of a prior on the correctness of each data point affects our results as shown in Figure 6. In short the variance increases and the best fit changes marginally with any use of the method, but a more aggressive use (low prior or high measurement variance), within reason, does not have significantly different results from a less aggressive (high prior or low measurement variance) use of the method.

The difference in parameter fit demonstrated in Figure 6 indicates the change in temperature fit as a result of the method. The fit shifts up by 10K on the temperature when a larger (reasonable) variance factor is used. This varies from 1 to 2 $\sigma$ of the temperature fit, which may not be important for temperature but could dramatically effect the fit for other parameters, especially when an atmospheric retrieval method is used for generating models. Poor constraints from highly variant or inconsistent data ought to lead to poor fits on parameters, and we find that this statistical method can accomplish this for photometry with underestimated uncertainties.

4.6 Implications of Exoplanet Variability

We considered the possibility that GJ 758b exhibits variability that effectively causes an additional, unaccounted error term. Radigan et al. (2014) measure T-dwarf variabilities up to 1.6% in the near-infrared due to rotation. However, this is far below the $\sim$10-20% errors on photometric points used in our analysis. For isolated objects, variability may contribute to “under-estimated” error bars. For exoplanet imaging in the high-contrast regime, derived errors tend to be larger than the underlying variability (Apai et al., 2016).

5 Conclusion

To characterize exoplanet atmospheres, we must be able to fit them with atmospheric models. When we have reason to consider some of our errorbars to be underestimated, we can use Bayesian statistics to handle erroneous data quantitatively.

We implemented a statistical method to consider some of our errorbars as being larger than initially estimated on photometry for substellar companion GJ 758b, in order to address inconsistency in the data. Using our method consistently increased the standard deviation on the parameter fit of our atmospheric models, making it less sensitive to a small number of highly variant data points. Hansen et al. (2014) has suggested that a key to characterizing exoplanets is translating underestimated uncertainties to uncertainties on model parameter fits, and we find our result to be consistent with this goal. In the case of the photometry used in this analysis, our method also rigorously put into question our data points in the CH4s and K2 bands.

We conclude that the problem of characterizing exoplanet atmospheres with photometry that has underestimated errorbars can be addressed by increasing the errorbars rigorously using this method. We emphasize its use for low signal-to-noise photometry-based characterizations, including transit photometry, especially for new instruments such as James Webb Space Telescope, where systemic error may not be immediately understood. We suggest that this method can also be used to make claims about the consistency of previous characterizations, as well as general trends of underestimated errorbars, as is necessitated by Burrows (2014) and Hansen et al. (2014).

6 Acknowledgments

We thank the referee for providing their very helpful report. We also thank the UC Santa Cruz LAMAT program as well as the UC Santa Cruz Other Worlds Lab for their generous funding for this project.

Appendix A Appendix: Press and Kochanek statistics

A.1 Composite probability formula

We’ve adapted Press and Kochanek’s formula from Wall & Jenkins (2012) for use in characterizing exoplanet atmospheres. As for regular model fitting, each measurement $X_{i}$ is compared to a model value $\theta_{i}$ which is computed using the filter curve used for $X_{i}$ and the model $\theta$ (as well as a Vega model). These values are compared using a Gaussian distribution and their given uncertainty, and the probability of the measurement being good is $\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{given}})$.

The product of all $\textrm{Prob}(X_{i}|\theta_{i},\sigma_{\textrm{given}})$ gives a likelihood function for a $\theta$’s goodness, and this allows us to find the best fitting model. The probability of each model being correct given some measurements and errorbars follows in Equation 1.

$\textrm{Prob}(\theta)$ represents the Bayesian normalization condition that one model must be correct. In this formula, $p$ iterates over the probability that the data are correct. This formula assumes no knowledge about the correctness of the observer.

Because the probability functions are added, the larger errorbar term will dominate only if it significantly improves the fit.

A.2 Normalizing the probability function

For calculating the probability of each data point being correct, we must normalize each data point to each model. The resultant normalization condition is Equation 2.

Therefore each term in the product must be multiplied by some constant contained within $\textrm{Prob}(\theta)$ which is the $\frac{1}{n}th$ root of the unnormalized summed product, before integration. In this method, the only free parameter is the relation between $\sigma_{\textrm{given}}$ and $\sigma_{\textrm{large}}$.

A.3 Adding a prior on p

We can add a prior on $p$ by instead evaluating the composite probability formula for p = prior. This leads to Equation 3 for the probability of a model being correct.

This results in the same normalization condition as before, except without the probability integration.

A.4 Point correctness formula

Press and Kochanek’s formula for the correctness of a given data point is adapted to Equation 4.

Simply put, this is the ratio of goodness of fit when the original probability is used as opposed to when the composite probability is used for the $ith$ data point in particular. It will be higher when using a larger errorbar does not result in a better fit, and lower when it does. Therefore, the method weighs how well a data point fits, for each model, against the fit of every other data point in the set, to compose the probability that it is correct. This measures the consistency of a data point with all other data points, using the models as a metric.