On the Dayside Atmosphere of WASP-12b

The thousands of exoplanets discovered to date span a wide range of properties and conditions, from small, rocky bodies to hot, Jupiter-like gas giants (Batalha 2014; Winn & Fabrycky 2015). Understanding the compositions of these planets provides real-world tests for atmospheric simulations and formation theories. Characterizing exoplanetary atmospheres requires an observed spectrum. Presently, most exoplanets can only be characterized via transit (when the planet moves in front of its host star, as seen from Earth) measurements. These observations capture starlight that has filtered through the exoplanet's atmosphere at the day–night terminator, imprinting information about its composition. For hot exoplanets, its secondary eclipse (when the planet moves behind the star, as seen from Earth) can be observed, which measures the planet's thermal emission. This provides better data than transits to constrain the atmospheric properties (Deming & Seager 2017).

The inference of atmospheric conditions from observed spectra is known as atmospheric retrieval (Madhusudhan 2018). For exoplanets, atmospheric retrieval involves the proposal of atmospheric models from some prior distribution (e.g., uniform or Gaussian), computation of the theoretical observed spectra, and determination of how well the proposed models explain the observations. Unlike solar system observations which only require least-squares minimization, a Bayesian approach is better suited to estimating the uncertainties of exoplanet retrievals due to the high relative noise levels. The Bayesian sampler explores the parameter space and accepts/rejects new models with some probability based in part on the goodness of fit. The collection of accepted models forms the posterior distribution, which informs the range and relative likelihood of values for the model parameters.

WASP-12b stands out as one of the hottest exoplanets found to date. It orbits an F9V star with a temperature of $6360{\pm }_{140}^{130}$ K at $0.02340{\pm }_{0.00050}^{0.00056}$ au every ∼1.09 days (Hebb et al. 2009; Collins et al. 2017). With a mass of $1.47{\pm }_{0.069}^{0.076}$ M_J and radius of $1.90{\pm }_{0.055}^{0.057}$ R_J , its density of $\mathrm{0.0.266}{\pm }_{0.014}^{0.015}$ g cm⁻³ is less than a quarter of Jupiter's density of 1.326 g cm⁻³. Due to its extreme equilibrium temperature (>2500 K), it is expected to be in thermochemical equilibrium (Moses 2014).

WASP-12b has been the target of numerous observations and analyses since its discovery in 2008. Its secondary eclipse has been observed across the near- and mid-infrared by a variety of instruments, including the Hubble Space Telescope (HST) Wide Field Camera 3 (WFC3), Spitzer Space Telescope Infrared Array Camera (IRAC), Canada–France–Hawaii Telescope (CFHT) Wideband Imaging Camera (WIC), Michigan-Dartmouth-MIT (MDM) Observatory TIFKAM, and Apache Point Observatory (APO) Near-Infrared Camera (López-Morales et al. 2010; Campo et al. 2011; Croll et al. 2011; Zhao et al. 2012; Cowan et al. 2012; Crossfield et al. 2012; Swain et al. 2013; Föhring et al. 2013). Combinations of these data have been used for retrievals of the dayside T(p) pressure–temperature profile and molecular abundances by Madhusudhan et al. (2011), Line et al. (2014), Stevenson et al. (2014), and Oreshenko et al. (2017) to investigate the atmospheric properties of this highly irradiated hot Jupiter. Madhusudhan et al. (2011) uses less data than the others due to the limited data at the time of publication. Further, that retrieval occurred before the discovery of WASP-12's binary M-dwarf companions, which reduces the measured eclipse depths (Bergfors et al. 2011; Bechter et al. 2014). As a result, our paper does not thoroughly compare those results to other investigations.

The results of the Line et al. (2014), Stevenson et al. (2014), and Oreshenko et al. (2017) retrievals are inconsistent in some respects. When considering CO, CO₂, CH₄, and H₂O, Line et al. (2014) do not find evidence for a high C/O due to high abundances of CO₂ and H₂O, while Stevenson et al. (2014) find a bimodal C/O, both of which have a high abundance of CO₂. However, these analyses find an abundance of CO₂ that is greater than both CO and H₂O, which has been shown to be highly improbable in the atmosphere of a planet like WASP-12b (Madhusudhan 2012; Moses et al. 2013; Heng & Lyons 2016). Line et al. (2014) comment that this is implausible and place an upper limit of 10⁻⁵ on the CO₂ mixing ratio; retrievals under this limit drive the H₂O mixing ratio over 100 parts per million, resulting in a more realistic CO₂ mixing ratio and maintaining a C/O near solar. Stevenson et al. (2014) also mention the implausibility of their retrieved CO₂ abundance and propose the addition of HCN and C₂H₂ into the retrieval model to solve this problem. These species have been shown to exist when C/O > 1 (Madhusudhan 2012; Moses et al. 2013) and have spectral features in Spitzer's IRAC channel 2; this allows the Bayesian sampler to fit the eclipse depths in that channel using the added species. Consequently, they retrieve an abundance of CO₂ that is less than CO, which is a physically plausible result. This C-rich result is more probable than their O-rich result by a factor of 670. They exclude an isothermal model at 7σ significance.

Oreshenko et al. (2017) performed retrievals using the same data as Stevenson et al. (2014) and expands upon that work by including clouds in their model. They find that the cloud compositions are unconstrained, an expected result considering the degeneracy between cloud composition and gas mixing ratios at low spectral resolutions. When considering CO, CO₂, CH₄, and H₂O, they replicate the results of Line et al. (2014) and Stevenson et al. (2014) of an unrealistically high CO₂ mixing ratio. In general, they find that their retrieval results are prior dominated. When assuming Gaussian priors for the C/H and O/H ratios of WASP-12b matching that of its host star (Teske et al. 2014), the resulting C/O is close to solar. However, when increasing the uncertainties on the values, C/O > 1 becomes the favored solution, consistent with Madhusudhan et al. (2011). Including HCN and C₂H₂ in the model results in an HCN mixing ratio of 10⁻²–10⁻¹, which they comment is implausible for reasons discussed in Moses et al. (2013); this mixing ratio is over three orders of magnitude greater than that found in Stevenson et al. (2014). An interesting result of their retrieval models is the wide range of possible T(p) profiles based on the assumed prior, ranging from no to strong inversions. However, despite this wide range of possibilities, the resulting spectra are qualitatively similar.

A curious difference among the retrieved T(p) profiles of these three investigations is that Line et al. (2014) find the temperature of the atmosphere to be almost entirely above 3000 K with the potential for an inversion, while Stevenson et al. (2014) and Oreshenko et al. (2017) find an upper limit of 3000 K with no inversion. Note that when computing opacities using HITRAN databases, the available data to calculate partition functions has an upper limit of 3000 K (Laraia et al. 2011). The CHIMERA code used by Line et al. (2014) probes temperatures above this limit by assuming that cross sections at temperatures >3000 K are equal to those at 3000 K (M. Line, 2018 private communication). In general, cross sections will differ between temperatures >3000 and 3000 K, so this assumption is more likely to give misleading results the greater the deviation from 3000 K. However, extrapolation beyond 3000 K could be even more misleading as there is no guarantee that it will match the true cross sections. This assumption in CHIMERA could lead to a better model fit by allowing the code to explore background emission temperatures greater than 3000 K. Nonetheless, this underscores the need for higher-temperature data to more accurately characterize planets in this temperature regime.

The temperatures found in these retrievals bring attention to another implicit assumption put into these models. Retrieval models consider some set of molecules to fit an observed spectrum. Omitting a molecule that is present in the real object will therefore bias the results: Stevenson et al. (2014) showed that the omission of HCN and C₂H₂ drives up the inferred CO₂ abundance, while including the additional molecules allows for a more realistic fit. At the temperatures retrieved for WASP-12b, both H₂ and H₂O thermally dissociate, forming H (Arcangeli et al. 2018; Kreidberg et al. 2018; Parmentier et al. 2018). Some H gain an electron from ionized metals, forming H⁻. To date, WASP-12b retrieval models have omitted H⁻, which provides an important continuum opacity source from bound–free and free–free transitions (John 1988), as more thoroughly discussed in Parmentier et al. (2018) and Arcangeli et al. (2018).

In this paper, we perform retrievals using the Bayesian Atmospheric Radiative Transfer (BART) code (Blecic et al. 2022; Cubillos et al. 2022; Harrington et al. 2022) matching the setups of Line et al. (2014) and Stevenson et al. (2014) to investigate the discrepancies in their results, and we present a new, physically motivated model that includes additional species not considered in previous investigations. Section 2 describes the BART code, and Section 3 discusses the setup for each of the nine models. In Section 4, we discuss our results in the context of previous analyses of WASP-12b's dayside atmosphere. Finally, we draw conclusions from our findings in Section 5.

Our retrieval code, BART (Blecic et al. 2022; Cubillos et al. 2022; Harrington et al. 2022), pairs the Transit radiative-transfer code (Rojo 2006) with Multi-Core Markov chain Monte Carlo (MCcubed, Cubillos et al. 2017), a Bayesian framework. The user specifies a parameter space to be explored for some model parameters (e.g., the T(p) profile and molecular abundances). Other inputs include the observational data and its type (e.g., transit or eclipse depths) as well as the instrument filters associated with each data point. For each proposed atmospheric model, the theoretical spectrum is calculated at a high resolution, binned according to the filters, and compared to the observational data. The abundance profiles begin from a user-specified atmosphere (e.g., uniform profiles, or thermochemical equilibrium for a certain T(p) profile), and MCcubed scales the abundances of the molecules being fit. For credible regions estimated from the posterior, BART computes the steps per effective independent sample (SPEIS) and effective sample size (ESS) to estimate the uncertainty in a given credible region, as detailed in Harrington et al. (2022).

To investigate the previously published retrievals of WASP-12b, we replicate their setups to the best of BART's ability, and we expand upon those setups to delve deeper into the nature of the discrepancies in results. Our nine retrieval models are that of

• 1.  
  Line et al. (2014 null case,
• 2.  
  Line et al. (2014) ellipsoidal case,
• 3.  
  Line et al. (2014) null case plus the data from López-Morales et al. (2010) corrected by Crossfield et al. (2012),
• 4.  
  Stevenson et al. (2014) case without HCN and C₂H₂,
• 5.  
  Stevenson et al. (2014) case with HCN and C₂H₂,
• 6.  
  Stevenson et al. (2014) without HCN and C₂H₂ with the assumptions of CHIMERA about cross sections,
• 7.  
  Stevenson et al. (2014) with HCN and C₂H₂ mixing ratios fixed to their reported C-rich best-fit values,
• 8.  
  Model 3, with H⁻,
• 9.  
  Model 8, with HCN and C₂H₂,
• 10.  
  Model 6, with H⁻,
• 11.  
  Model 10, with HCN and C₂H₂,
• 12.  
  Model 3, with NH₃, HCN, C₂H₂, TiO, and H⁻, and
• 13.  
  Model 6, with the molecules of Model 12.

At the time of writing, BART does not have a realistic cloud model, so we do not try to replicate the Oreshenko et al. (2017) result that finds cloud composition to be unconstrained. Table 1 summarizes the setup of each model regarding data sources and molecules besides H, H₂, and He.

For this investigation, BART's only restrictions on the atmospheric models are (1) the sum of molecular abundances must equal 1 for each layer, and (2) the ratio of H₂ to He is held constant by adjusting their abundances to satisfy condition 1. For models that do not use CHIMERA's assumption about cross sections at T > 3000 K, BART also enforces that the temperature of the atmosphere must remain within the line-list limits.

The atmospheric models consist of 100 log-spaced layers spanning 10⁻⁸–100 bar. For radiative-transfer calculations, only layers above where the optical depth reaches ≥10 are considered. We assume uniform abundances for the species present, consistent with previous publications. Each model has a free parameter for the abundance of each opacity-contributing species, as well as five free parameters for the T(p) profile (the Planck mean infrared opacity, the ratios of the Planck mean visible and infrared opacities for two streams, the partition between the two streams, and a general parameter for albedo/emissivity/energy recirculation; Line et al. 2013). The free parameters for the partition between streams and albedo/emissivity/recirculation have uniform priors; all other free parameters have log-uniform priors. We do not vary the H⁻ or e⁻ abundances because previous publications indicate the atmosphere is nearly isothermal at ∼3000 K in the regions with sensitivity, and the abundances change by only a factor of ∼2 for a change as large as 200 K. H⁻ and e⁻ are fixed to an abundance of 10⁻⁹ and 10⁻⁶, respectively, which is roughly consistent with thermochemical equilibrium at 3000 K (NASA Chemical Equilibrium with Applications code; Gordon & McBride 1994) at pressures probed by the observations. A wide range of values are allowed for each free parameter without consideration of physical plausibility. For this investigation, we use the DEMCzs sampling algorithm of ter Braak & Vrugt (2008) because we found that the DEMC algorithm of ter Braak (2006) occasionally leads to rogue chains that do not converge. Since DEMCzs only considers the goodness of fit of each proposed model, it can explore both realistic and unrealistic solutions. The initial samples of parameters for the DEMCzs algorithm are drawn randomly from a uniform distribution; most of the parameters are the logarithm of the true parameter, so those true parameters are randomly sampled from a log-uniform space.

We include HITEMP opacities for CO, CO₂, and H₂O (Rothman et al. 2010), and HITRAN opacities for NH₃, C₂H₂, and HCN (Rothman et al. 2013). Models 1–11 use the Rothman et al. (2013) CH₄ line list for consistency with Line et al. (2014). while Models 12 and 13 use the new CH₄ HITEMP line list (Hargreaves et al. 2020). TiO opacities are sourced from Schwenke (1998). We include H₂-H₂ and H₂-He collision-induced absorptions (Richard et al. 2012) as well as H⁻ bound–free and free–free absorption (John 1988), where appropriate.

Note that Equation (3) of John (1988) does not lead to the correct bound–free opacity values necessary to reproduce Table 1 of John (1988); a factor of 10 for the bound–free cross sections is required to obtain agreement. We refer the reader to our compendium for detailed proof. It is unclear whether the given equations or table provide the correct opacities. Considering the numerous fitted constants used in the equations, we have chosen to assume that the table is correct, as it appears to lead to agreement with the H⁻ opacities plotted in Figure 1 of Arcangeli et al. (2018) and Figure 4 of Parmentier et al. (2018). Additionally, by assuming the greater of the two possibilities, we can assess an upper limit on the impact of H⁻ when retrieving WASP-12b's atmospheric properties.

The accepted T(p) profiles with 1σ and 2σ regions, normalized contribution functions, best-fit spectrum, and 1D marginalized posteriors are shown for Models 12 and 13 in Figures 1 and 2, respectively. The electronic supplement provides additional figures: 2D pairwise posteriors and trace plots for Models 12 and 13, as well as the corresponding set of six plots for Models 1–11. Table 2 contains the best-fit values and 68.27% interval for the retrieved log abundances of each molecule for all 13 models, the best-fit values and 68.27% interval reported by Line et al. (2014) for the null and ellipsoidal cases, the best-fit values reported by Stevenson et al. (2014) for the C-rich and O-rich cases, and the lower/upper limits on abundances shown in the top of Figure 3 of Oreshenko et al. (2017). We also report the best-fit and 68.27% credible region for C/O for each case. We have excluded extreme outliers (C/O > 1000) from the density estimation, as they represent a small percentage of the total models and cause problems for the density estimation algorithm, and we do not consider HCN or C₂H₂ when calculating C/O due to the lack of evidence to support their inclusion in the model. Values of "..." indicates that the model does not contain that molecule. Table 3 lists the SPEIS, ESS, and associated uncertainty in the 68.27% credible region for each model considered. We also provide a compendium with the data and commands necessary to reproduce this work; the link is at the end of the text.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (13 panels)

Tarball of figure set images

For Model 1 (Figure 1, figure set panel 1.1), BART's results generally agree with that of the null case of Line et al. (2014). Our best-fit abundances have CO₂ > CO, as Line et al. (2014) find, which, as mentioned previously here and in their paper, is implausible. However, there is large uncertainty in this result, as indicated by the nearly flat posterior of both CO and CO₂. Like Line et al. (2014), we find that CO and CH₄ are unconstrained, according to the flat posteriors, but we also find that CO₂ is unconstrained. Except for CH₄, the best-fit values fall within the 68.27% region reported by Line et al. (2014). The T(p) profile 1σ regions overlap, with both best-fit profiles favoring an inversion.

For Model 2 (Figure 1, figure set panel 1.2), BART's results differ in some respects from the ellipsoidal case of Line et al. (2014). Most notably, the T(p) profile 1σ region shows a non-inverted atmosphere where the upper atmosphere is <3000 K, whereas Line et al. (2014) find an inverted atmosphere with the upper atmosphere >3000 K. However, the normalized contribution functions indicate minimal sensitivity below a pressure of 10⁻⁶ bars; in the best-constrained region (0.1–10 bars), the retrieved T(p) profiles agree. For abundances, the best-fit values disagree, but, except for CO₂, the 68.27% intervals overlap. However, the case of Line et al. (2014) where CO₂ has an upper limit agrees more closely with our retrieved interval. It is unclear why BART does not find the high CO₂ mode found by Line et al. (2014).

The results for Model 3 (Figure 1, panel 1.3) in general agree with Model 1. CO, CO₂, and CH₄ are similarly unconstrained, while the 68.27% region for H₂O overlaps (Table 2). A notable difference is that Model 3 finds a lower minimum for that region. This is likely from the additional data point of López-Morales et al. (2010) providing an additional constraint on the background emission.

BART's results for Model 4 (Figure 1, panel 1.4) in many ways match that of Model 2. The T(p) profile 1σ regions closely match, and the marginalized posteriors exhibit many similarities. Similar to Model 2, BART does not find the unrealistically high CO₂ abundance reported by Stevenson et al. (2014) when retrieving with this setup (Table 2). Since the parameter space for CO₂ extended to a log mixing ratio of −1, these models must have been discarded by BART's sampler. It is uncertain whether this difference can be attributed to the sampling algorithm, the opacity sources, or some other difference in the retrieval algorithm. Nevertheless, the best-fit T(p) profiles of Stevenson et al. (2014) fall within BART's reported 1σ region where there is sensitivity. Further, their C-rich best-fit values are within BART's 68.27% regions for all but CO. Our 68.27% interval for H₂O rules out their O-rich result, which features an H₂O abundance of 10^−3.3.

The results for Model 5 (Figure 1, panel 1.5) mostly agree with the findings of Stevenson et al. (2014). The C-rich best-fit values for CO₂, CH₄, H₂O, and C₂H₂ reported by them fall within BART's 68.27% regions (Table 2). Our 68.27% interval for H₂O similarly rules out their O-rich result. BART finds HCN to be unconstrained. The retrieved T(p) profiles closely match the results of Model 4 and are therefore consistent with the retrieved T(p) profiles of Stevenson et al. (2014). The best-fit values and 68.27% region found for Model 5 closely matches the results of Model 4, which does not include HCN or C₂H₂. Thus, the inclusion of these additional molecules only minorly affects the retrieved result.

Model 6 (Figure 1, panel. 1.6) examines the effect of allowing T(p) profiles with T > 3000 K on Model 4. BART's results closely match those of Models 4 and 5. The T(p) profile 1σ regions generally agree, with a general upper limit of ∼3000 K. The marginalized posteriors are similar, and the 68.27% regions for the molecular abundances closely overlap. This importantly demonstrates that the eclipse data, not the model assumptions, are driving the result.

For Model 7 (Figure 1, panel 1.7), BART's results generally agree with those of Model 5. Even with the HCN and C₂H₂ abundances fixed to the C-rich best-fit values reported by Stevenson et al. (2014), BART does not favor their reported CO abundance. We suspect that this difference can be attributed to differences in line lists, though we were unable to ascertain the line lists used by Stevenson et al. (2014) via private communication to test this hypothesis.

Models 8, 9, and 12 (Figure 1, panels 1.8, 1.9, and 1.12), which are equivalent to Model 3 but with additional opacity sources, bear qualitatively similar results to Model 3. While the 68% regions may differ slightly, these are largely due to numerical differences in the density estimation of the marginalized posteriors. Visually, the posteriors are quite similar; CO, CO₂, and CH₄ have flat marginalized posteriors, while H₂O has a preference for log abundances greater than −6. In Model 12, the flat posteriors of HCN, C₂H₂, and NH₃ demonstrate that the data are insensitive to their inclusion. While the log abundance TiO is loosely constrained to be −7.5 ± 2.5, consistent with thermochemical equilibrium for the retrieved T(p) profile in the regions probed by the observations, the posterior allows for TiO to be absent from the atmosphere within 2σ.

Similarly, Models 10, 11, and 13 (Figure 2; see also figure set panels 1.10, 1.11, and 1.13), which are equivalent to Model 6 except with additional opacity sources, yield results similar to Model 6. The resulting thermal profiles agree closely. As with Models 8, 9, and 12, the 68% regions differ slightly due to the density estimation algorithm, though the marginalized posteriors are qualitatively similar to Model 6. That is, CO tends to favor smaller values, CO₂ is tightly constrained within an order of magnitude from 10⁻⁵, CH₄ features a constraint of about 10⁻³ but could also be absent from the atmosphere (the posterior has a non-negligible tail), and H₂O has an upper limit of around 10⁻⁴. For Model 13, BART finds upper limits of 10⁻² for NH₃ and 10⁻⁶ for TiO.

Zoom In Zoom Out Reset image size

As Models 8 and 10, 9 and 11, and 12 and 13 bear identical setups/assumptions aside from the eclipse depths, the differences between them are thus solely attributable to the data. Both data sets favor a ∼3000 K atmosphere in the regions probed by the observations, as evidenced by the normalized contribution functions. However, their retrieved abundances, and by extension the inferred C/O, are incompatible. The data set of Models 8, 9, and 12 yields evidence of a high water abundance and possibly TiO in thermochemical equilibrium, with no evidence of other molecules, while the data set of Models 10, 11, and 13 constrain CO₂ and CH₄, with upper limits for CO and H₂O. None of the data sets considered offer meaningful constraints on HCN, C₂H₂, or NH₃. Despite the theoretical expectation that H⁻ plays an important role in the atmosphere of hot Jupiters like WASP-12b (Parmentier et al. 2018), its inclusion does not make a significant difference in the retrieval results (Figure 3).

Zoom In Zoom Out Reset image size

Our results are broadly consistent with those of Oreshenko et al. (2017), except for HCN (Table 2). In general, the models using the Line et al. (2014) data favor an atmosphere that is isothermal or has an inversion, whereas the models using the Stevenson et al. (2014) data favor an atmosphere that has no inversion. Lothringer et al. (2018) show that thermal inversions are likely for this planet class even without VO or TiO. While the results using the Line et al. (2014) null data set are consistent with this finding, it is important to consider that these results also favor a high H₂O abundance, which is inconsistent with the expected thermal dissociation of H₂O. While the Line et al. (2014) ellipsoidal case favors a low H₂O abundance, it does not show strong evidence of a thermal inversion and more closely matches the results of Stevenson et al. (2014)

Over the reported 68.27% credible regions, the C/O for these models can take a wide range of values. However, they do not indicate that C/O ≫ 1. Rather, it is due to a combination of reasons. For one, the DEMCzs sampler is free to explore the parameter space without regard for C/O. For example, in the case of Model 4, the best-fit C/O value is 52.8, which is notably different than the 68.27% region; this is a product of a high best-fit abundance of CH₄ and low best-fit abundance of H₂O, which provides a better statistical fit than models with more reasonable C/O values. Additionally, thermal dissociation of H₂O would lead to the formation of O (Parmentier et al. 2018), causing an apparent increase in the measurable C/O but not the true C/O. There are also other oxygen-bearing species not considered here, particularly condensates at the limb (Wakeford et al. 2017), that would contribute to C/O, if included in the model. However, Oreshenko et al. (2017) demonstrated that the cloud composition is degenerate with gas mixing ratios; higher-resolution data is necessary to explore a model with various species of condensates.

Typically, when considering multiple models, the Bayes factor of each model is compared to choose the best model. In this investigation, however, this would be erroneous: the 13 models presented do not use the same data sets, and they are not competing for the best model. Rather, the models demonstrate that the retrieval results are data driven and independent of the model selected. Consequently, we do not compute the Bayes factor as it would be a misleading metric. We emphasize that the results show that the previous retrieval analyses of Line et al. (2014) and Stevenson et al. (2014) are consistent, when considering the data set used in each investigation. Follow-up observations are required to determine which data set, if either, represents the true nature of WASP-12b.

In general, we are able to reproduce the published results of Line et al. (2014), Stevenson et al. (2014), and Oreshenko et al. (2017) using BART. We confirm the finding of Stevenson et al. (2014) that excludes an isothermal profile when mimicking their setup. By following the model assumption of Line et al. (2014) allowing temperatures above 3000 K with the Stevenson et al. (2014) data, the range of possible T(p) profiles expands to include inverted profiles but still favors a non-inverted atmosphere. Note that an inverted profile is expected for ultra-hot Jupiters like WASP-12b (Lothringer et al. 2018).

We find that current data does not support the inclusion of HCN, C₂H₂, H⁻, e⁻, NH₃, and TiO, as they do not significantly affect the posterior. As new telescopes in the near future provide higher quality data, these molecules should be reconsidered, as in this investigation, to determine if they inform the results.

Some aspects were unable to be reproduced, namely, the unrealistically high CO₂ abundances reported in the Line et al. (2014) ellipsoidal case and the Stevenson et al. (2014) case without HCN or C₂H₂, and the unrealistically high HCN abundance reported by Oreshenko et al. (2017). While we suspect this discrepancy is due to differences in the retrieval model, further investigation is necessary to definitively determine the origin.

We have demonstrated that differences in eclipse depth data sets primarily drive the differences between the results of Line et al. (2014) and Stevenson et al. (2014), with more subtle differences likely attributable to the retrieval model and input data sources (e.g., line lists). Many of the eclipse depths come from the same set of observations but are analyzed using different reduction pipelines. Our study shows that subtle differences in the reduction pipelines used (e.g., the binning of WFC3 spectra into discrete channels) can drive radically different results. This emphasizes the need for standard data sets to be used for benchmarking photometry and spectroscopy reduction pipelines.

The conflicting results of previous publications highlight the complexities of retrieval modeling and the importance of clearly communicating model assumptions. This will be especially important as retrieval models become more sophisticated with the introduction of more complicated techniques such as 3D modeling and machine learning (Márquez-Neila et al. 2018; Zingales & Waldmann 2018; Waldmann & Griffith 2019; Cobb et al. 2019). Transparency allows for published analyses to be easily reproduced by others and encourages quicker resolution of conflicting results, which can lead to better work in the field.

Future retrieval studies should use multiple binnings of the same data to explore whether the retrievals are consistent across different binnings of the data, including the unbinned data. This requires that data analyses publish unbinned spectra to support future retrieval studies. Retrieval results dependent on the binning, such as those shown here, indicate the need for higher quality spectroscopic data. A comprehensive retrieval analysis of WASP-12b with additional data from future flagship observatories, such as the James Webb Space Telescope, will provide deeper insight into the nature of this extreme exoplanet.

The Reproducible Research Compendium for this work is available for download. ¹

We thank Michael Line and Julianne Moses for helpful discussions during the preparation of this manuscript. We also thank the anonymous referee for valuable comments that improved the quality of this manuscript. We thank contributors to SciPy, Matplotlib, and the Python programming language, the free and open-source community, and the NASA Astrophysics Data System for software and services. Part of this work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This work was supported by NASA Planetary Atmospheres grant NNX12AI69G, NASA Astrophysics Data Analysis Program grant NNX13AF38G, and NASA Fellowship Activity under NASA Grant 80NSSC20K0682.