NEW ANALYSIS INDICATES NO THERMAL INVERSION IN THE ATMOSPHERE OF HD 209458b

As more exoplanets are discovered and studied every year, it becomes necessary to develop a comprehensive understanding of exoplanet characteristics. Characterizing the atmospheres of exoplanets is the first step to probing their interior structures and formation histories. With this in mind, an important question to ask is whether or not an exoplanet possesses a thermal inversion in its atmosphere. If an atmosphere has a large opacity at optical or ultraviolet wavelengths relative to the opacity at thermal infrared wavelengths, then the atmosphere can absorb incident stellar radiation at high altitudes without efficiently radiating this energy back to space. This would warm a region of the upper atmosphere relative to a deeper one, constituting a thermal inversion. A non-inverted atmosphere simply decreases in temperature with increasing altitude. The thermal structure of an atmosphere depends on its opacities, and therefore its composition.

Given their size relative to their host stars, we have so far been most successful at detecting and analyzing the atmospheres of hot Jupiters. There have been notable detections of thermal inversions in several of these exoplanets, including HD 149026b (Harrington et al. 2007), HD 209458b (Burrows et al. 2007; Knutson et al. 2008), and XO-1b (Machalek et al. 2008). The diversity of hot Jupiters in which thermal inversions have been detected has led to concerted efforts to classify hot Jupiters on the basis of their atmospheric profiles, host stars, and chemical abundances (Fortney et al. 2008; Knutson et al. 2010; Madhusudhan 2012).

All thermal inversion detections are based on data from the Spitzer Space Telescope. Until the cryogen was depleted in 2009, Spitzer allowed for the acquisition of photometric data in the 3.6, 4.5, 5.8, and 8.0 μm bandpasses. Data from these wavelength ranges can place constraints on the atmospheric abundances of H₂O, CO, CO₂, and CH₄. By observing in the Spitzer bandpasses as a planet passes behind its host star, in what is termed a secondary eclipse, it is possible to detect light directly from the planet and to determine at which wavelengths there is absorption or emission of a given chemical species. A spectral feature in absorption indicates a monotonically decreasing temperature with altitude, while a spectral feature in emission indicates a thermal inversion.

There has been difficulty in determining the nature of the chemical species that would account for thermal inversions. TiO and VO could exist in the gas phase at high altitudes in the atmospheres of irradiated giant planets (Hubeny et al. 2003), but there are questions as to whether these heavy absorbers could remain at high altitude in the hydrogen-dominated atmospheres of hot Jupiters (Spiegel et al. 2009). Aside from the intense vertical mixing necessary to keep TiO/VO aloft, these absorbers could be depleted by the nightside cold traps of tidally locked hot Jupiters (Showman et al. 2009; Parmentier et al. 2013) or be dissociated by the intense ultraviolet radiation from the nearby host star (Knutson et al. 2010).

One of the first secondary-eclipse observations by Spitzer resulted in the detection of emission features in the spectrum of HD 209458b (Knutson et al. 2008). Previous models invoked a thermal inversion to fit the quoted eclipse depths (Burrows et al. 2007; Madhusudhan & Seager 2009; Line et al. 2014), but all of them rely upon four bandpass-averaged photometric points acquired in 2005, predating many systematic corrections now commonplace for Spitzer observations.

Recently, the thermal inversions of other exoplanets have been called into question. The detection of a thermal inversion in the atmosphere of HD 149026b rested on a single observation at 8.0 μm (Harrington et al. 2007), but was refuted by Stevenson et al. (2012a), who used additional observations at more wavelengths to rule out an inversion. The proposed thermal inversion on the relatively cool hot Jupiter XO-1b may be better explained by a supersolar C/O ratio rather than an optical absorber in the atmosphere (Madhusudhan 2012). Given these findings and advancements in the field of exoplanet atmosphere characterization, we feel it is appropriate to reexamine the original data used to identify the thermal inversion in HD 209458b, as well as to analyze newer secondary-eclipse data that have since become available. We use up-to-date techniques and models to perform a complete, self-consistent analysis of the Spitzer data, with the goal of investigating the thermal structure of HD 209458b's atmosphere.

We investigate HD 209458b secondary-eclipse data acquired with Spitzer Space Telescope's InfraRed Array Camera (IRAC; Fazio et al. 2004). IRAC has two kinds of detectors. The InSb detector arrays observe in the 3.6 and 4.5 μm bandpasses, while the Si:As detector arrays observe in the 5.8 and 8.0 μm bandpasses. Until 2009, when the cryogen supply on board was depleted, this instrument was used to collect data in these four broad photometric bandpasses. The Spitzer Warm Mission continues to observe with the InSb detector arrays in the 3.6 and 4.5 μm bandpasses.

The Spitzer data sets that are the focus of this work each include an occultation of HD 209458b as well as out-of-eclipse baselines of the system both before and after the eclipse. Data points outside of the eclipse help to constrain the instrument systematics as well as provide a comparison to the eclipse depth from which we derive constraints on the atmospheric composition.

Spitzer's pipeline processes the raw detector data taken in subarray mode, removing well-understood instrumental signatures, and provides basic calibrated data (BCD) files. In order to observe a sixth K-magnitude star like HD 209458 without saturating the IRAC detectors, the observers worked in subarray mode, capturing sixty-four 32 × 32 pixel frames per BCD file. We then use our Photometry for Orbits, Eclipses, and Transits (POET) pipeline to further reduce the data and remove subtle instrument systematics from the light curve of the secondary eclipse (Stevenson et al. 2012b; Cubillos et al. 2013). At the start of the POET pipeline, we create a mask for bad pixels in the BCD frames and combine it with the bad pixel mask provided by Spitzer. We then determine the two-dimensional (2D) Gaussian center in each frame in order to perform aperture photometry. We test photometric aperture radii in 0.1 pixel increments and determine the best aperture by fitting a common model to the data at each aperture increment. We then compare the resulting standard deviation of the normalized residuals (SDNR). We proceed with the aperture size that produces the lowest SDNR value.

Once we determine the best aperture size, the bulk of the work of producing a clean light curve goes into exploring and modeling the instrument systematics unique to each data set. In all data sets, we use the uniform source equations described by Mandel & Agol (2002) and employ the Levenberg–Marquardt algorithm to find the best-fit results for the free parameters. We fix the parameters of secondary-eclipse duration and ingress/egress times to the calculated values of 3.08 hr and 0.42 hr, respectively. We then implement a differential evolution Markov Chain Monte Carlo (DEMCMC) with 10⁵ steps in each of 10 chains in order to explore correlations in the parameter space and estimate uncertainties (ter Braak 2006).

A key consideration is the intra-pixel effect, by which certain areas of a given pixel are more sensitive to incoming photons, and thus slight variations in stellar position translate into variations in flux that are larger than the secondary-eclipse depth we are looking for (Charbonneau et al. 2005; Knutson et al. 2008). This effect is most notable in the 3.6 and 4.5 μm bandpasses, but it can also arise when working with small aperture sizes in the longer wavelength bandpasses. To model this systematic, we employ a BLISS map, which uses a spline to fit the subpixel sensitivity at high resolution (Stevenson et al. 2012a). We fit this and other systematics simultaneously with the secondary eclipse in order to derive our eclipse depths and uncertainties.

We investigate all publicly available Spitzer data to obtain our results. As outlined in Table 1, these data sets include three eclipses at 3.6 and 4.5 μm (Channels 1 and 2, respectively), one eclipse at 5.8 μm (Channel 3), and two eclipses at 8.0 μm (Channel 4). In Section 2.1, we discuss our treatment of the original 2005 data sets (one eclipse in each channel), which yielded the result of a thermal inversion for Knutson et al. (2008). In Section 2.2, we discuss our treatment of more recently acquired data sets (eclipses in Channels 1, 2, and 4), including our investigation of the presence of time-correlated noise in the data.

2.1. Reanalysis of Previously Published Data

In 2005, Spitzer program 20523 (PI: David Charbonneau) used all four bandpasses for a duration of 8.1 hr to observe a single occultation as well as baseline of HD 209458b (Knutson et al. 2008). To keep from saturating the 3.6 μm detector, the 2005 data sets were taken with 0.1 s exposure times (Table 1). Observing the same eclipse in all four IRAC bandpasses eliminates the concern of variability from one eclipse to another; however, in subarray mode, multiple detector arrays cannot acquire target data simultaneously. All four IRAC detectors collected data throughout the duration of the observation, but the telescope continuously repointed, cycling between the detectors such that only one acquired data of HD 209458b in eclipse at a time. Over the course of one cycle, a single detector collected four BCD files of target data (a "batch") before the telescope repointed. While one detector was acquiring target data, the rest were collecting data from adjacent fields. Taking into account the time needed to repoint, this effectively lowers the duty cycle of each detector to less than 25% compared to typical IRAC data.

Due to the constant repointing of the telescope and the time lapses between batches of target data, systematic errors grow and information is lost, and the observational strategy employed in the 2005 HD 209458b campaign is no longer used. Since our goal is to provide a complete picture of the atmosphere of HD 209458b, we perform a thorough analysis of this data set and compare it to subsequent data sets taken in the 3.6, 4.5, and 8.0 μm bandpasses. The 2005 data set provides our only constraints on HD 209458b's atmosphere at 5.8 μm.

In analyzing this original data set, Knutson et al. (2008) found a systematic ramp that occurred in the first BCD file of every batch of target data taken in a given bandpass. They therefore decided to remove this first BCD file from each batch, depleting the remaining usable data in each bandpass by an additional 25%. We apply a different approach in order to retain a maximal amount of data. By separately analyzing the first, second, third, and fourth BCD files of every batch of data in a given bandpass, we are able to investigate the systematic noise that may occur over the course of a single pointing of Spitzer. In the 3.6, 4.5, and 8.0 μm bandpasses, the first BCD files differ significantly from the others with regard to the system flux, or baseline; there is little deviation in the first BCD file of the 5.8 μm bandpass. In each bandpass, we perform a joint fit across the 4 BCD files of every batch, allowing only the system flux to vary from one BCD file to the next, while the rest of the model and ramp parameters are shared. With this method we are able to keep all frames available in each bandpass, while still accounting for changes in baseline flux that arise due to the constant repointing of Spitzer during observation.

The noise in the 2005 data set at 5.8 μm is comparable in size to the secondary-eclipse signal, and our ramp models have difficulty finding the ingress and egress points. We take advantage of the fact that the four data sets obtained in this observational campaign target the same eclipse. We perform a joint fit between all four bandpasses to determine a weighted average eclipse time. We then fix this eclipse time for the 5.8 μm bandpass data in order to produce a light-curve model.

As a test, we perform a separate analysis of the 2005 data set with the parameters outlined by Knutson et al. (2008). We employ similar systematic model components and aperture sizes, and we eliminate the first BCD file from every batch in a given bandpass. Even with this similar setup of parameters, our new treatment of the data, especially with regard to our handling of systematic errors, recovers discrepant values from those reported by Knutson et al. (2008).

The difference between our best results for the 2005 data set and those reported by Knutson et al. (2008) is over 1σ at 3.6 μm, over 2σ at 4.5 μm, over 3σ at 5.8 μm, and almost 2σ at 8.0 μm. Figure 1 and Table 2 illustrate these discrepancies in eclipse depth. These differences may be due in part to the fact that we received our data from Spitzer pipeline version S18.18.0, while Knutson et al. (2008) received it from version S13.0. Moreover, our data reduction techniques employ current methodology, especially with regard to mapping intra-pixel sensitivities in the 3.6 and 4.5 μm bandpasses.

Zoom In Zoom Out Reset image size

Table 2. Comparison of Derived Eclipse Depths from Spitzer Program 20523 (PI: Charbonneau)

Label	Wavelength	Aperture	Ramp	Intra-Pixel	Eclipse Depth
	(μm)	(pixels)	Modelc	Sensitivity Map	(%)
Channel 1 a	3.6	2.3	Quadratic (1)	BLISS	0.113 ± 0.010
Channel 1 (K08)b	3.6	5.0	Quadratic	quadratic	0.094 ± 0.009
Channel 2 a	4.5	2.1	Linear (2)	BLISS	0.167 ± 0.014
Channel 2 (K08)b	4.5	5.0	Quadratic	quadratic	0.213 ± 0.015
Channel 3 a	5.8	2.2	Quadratic (1)	None	0.134 ± 0.035
Channel 3 (K08)b	5.8	3.5	Quadratic of ln	None	0.301 ± 0.040
Channel 4 a	8.0	2.8	Quadratic (1)	None	0.303 ± 0.023
Channel 4 (K08)b	8.0	3.5	Quadratic of ln	None	0.240 ± 0.026

Notes. ^aValues and functions derived from this work, using the same data set, obtained in 2005, as Knutson et al. (2008). ^bValues and functions quoted from Knutson et al. (2008). ^cRamp equations used in this work are designated by number and can be referenced in Section 2.1.1.

Download table as:  ASCIITypeset image

2.1.1. The 3.6 and 4.5 μm Bandpasses

The 3.6 and 4.5 μm bandpasses exhibit significant intra-pixel variability. In both bandpasses, we have few good frames relative to subsequent data sets (Table 1), meaning that our measured eclipse depths for these data sets have relatively high uncertainties. The BLISS map employed in this analysis relies on having enough photons received by a given subpixel region to map its sensitivity. Looking at a 2D histogram of the data in Figure 2, it is clear that there are not enough photons clustered on the same pixel to do this effectively, so the map becomes highly flexible with large uncertainty. We attempted to use a quadratic function in both the x and y positions to model intra-pixel sensitivity, but while this method can fit the overall ramp, in many cases it was unable to detect the occultation.

Zoom In Zoom Out Reset image size

At 3.6 μm, we find that the best aperture has a radius of 2.3 pixels with a time-dependent quadratic function of

$$\begin{equation} R(t) = 1 + r_{2} (t - 0.5) + r_{3} (t - 0.5)^2 \end{equation} \tag{ 1 }$$

as a best fit for the ramp in the data, where t is time in units of phase, and r₂ and r₃ are free parameters. With these parameters, we determine an eclipse depth of 0.113% ± 0.010% (Figure 1; Table 2). At 4.5 μm, we use an aperture radius of 2.1 pixels and employ a linear function

$$\begin{equation} R(t) = 1 + r_{2} (t - 0.5) \end{equation} \tag{ 2 }$$

as a best fit for the ramp. Here, we determine an eclipse depth of 0.167% ± 0.014% (Figure 1; Table 2).

Given the difficulties in developing an effective pixel map for these bandpasses, we turn to other data sets taken in 2010 and 2011 to further explore the parameter space and to obtain more accurate eclipse depths and uncertainties (Section 2.2).

2.1.2. The 5.8 and 8.0 μm Bandpasses

The 5.8 and 8.0 μm bandpasses do not exhibit intra-pixel sensitivity. We attempted to use the BLISS map in both bandpasses since the small aperture sizes we use can lead to greater pixelation effects. In neither bandpass were we able to find a subpixel bin size at which the BLISS map outperformed a nearest-neighbor interpolation (Stevenson et al. 2012a). We ultimately achieved consistent results without using the BLISS map, and we do not include it in our final analysis of these bandpasses.

At 5.8 and 8.0 μm, we use aperture sizes of 2.2 and 2.8 pixels, respectively, and a quadratic ramp to model the data. At 5.8 μm, we achieve an eclipse depth of 0.134% ± 0.035%, and at 8.0 μm, we achieve an eclipse depth of 0.303% ± 0.023% (Table 2; Figure 1). While we have no further data at 5.8 μm against which to compare our results, we do look at a subsequent data set in the 8.0 μm bandpass, taken in 2007 (Section 2.2).

2.2. Analysis of Newer Data Sets

In 2007, Spitzer program 40280 (PI: Heather Knutson) took a half-orbit phase curve of HD 209458b in the 8.0 μm bandpass (Table 1). In 2010 and 2011, Spitzer acquired data of HD 209458b as part of a campaign to capture a complete phase curve, this time in the 4.5 and 3.6 μm bandpasses, respectively (Spitzer program 60021; PI: Heather Knutson). From these observations, we extract one occultation of HD 209458b in the 8.0 μm bandpass and two consecutive occultations of HD 209458b in each of the 3.6 and 4.5 μm bandpasses; however, only one eclipse in each of the lower two bandpasses is "clean," while the other exhibits abundant noise (Figure 3). It is unclear whether the abundant noise in these eclipses is a result of short-timescale stellar activity or unknown instrument systematics. Given that HD 209458 is a relatively inactive star (Knutson et al. 2010), we determine that these particular data sets are plagued by time-correlated noise that does not manifest predictably.

Zoom In Zoom Out Reset image size

We clip the phase-curve data around the secondary eclipses to prevent the periodic flux variations from affecting the measured secondary-eclipse depths. Clipping the phase curve to just the secondary eclipses plus a few hours of baseline means that we can use simple linear or quadratic equations to fit the temporal systematics in the data.

In performing aperture photometry, we again choose aperture sizes that produce the lowest SDNR values. The best apertures tend to be small (⩽3.6 pixel radii), but smaller apertures are more susceptible to flux variations from imprecise stellar centering. We therefore use Time-series Image Denoising (TIDe) to remove high-frequency jitter in the stellar centering while performing photometry on the unfiltered images (Stevenson et al. 2012b). The temporal continuity of the data and short exposure times for each frame allow us to employ this method. The 3.6 μm bandpass data were captured with 0.1 s exposures, and we see a decrease in SDNR when TIDe is used. The 4.5 and 8.0 μm bandpass data were captured with 0.4 s exposures, and thus we see only a slight decrease in SDNR and almost no change in eclipse depth when we apply TIDe.

Due to the continuity and equal spacing of the data, we are able to better estimate our uncertainties by taking into account the contribution of unknown time correlations in the data, or red noise (Carter & Winn 2009). Investigating the red noise does not change our eclipse depths, but it does ensure that our uncertainties are large enough to account for these time correlations.

The red noise follows a power spectral density expressed by the equation 1/f^γ (Carter & Winn 2009). We explore all available wavelets in the Python PyWavelets package and choose the best one, based on SDNR and Bayesian information criterion values, to model the time-correlated noise in the data. We determine the coefficients for white noise, red noise, and γ by incorporating these parameters into the DEMCMC used to fit the eclipse parameters, temporal systematics, and BLISS map. A result of γ = 1 implies that the noise in the data is made up of equal parts uncorrelated white noise (γ = 0) and time-correlated red noise (γ = 2). Assuming a γ of one, as is done in the example provided by Carter & Winn (2009), can result in an under- or overestimation of the uncertainties, depending on the amount of correlated noise. Freeing γ constrains the relative amounts of white and red noise present in the data and allows us to more accurately predict uncertainties for our eclipse depths.

As an experiment, we fix γ to several values between 0.0 and 2.0 for the Channel 1 (2011a) data set. We find an increase in uncertainty with increasing γ. Had we fixed γ to one, we would have overestimated our eclipse depth uncertainty for this data set by 30%. From looking at Table 3, it is clear that the "noisy" eclipses have γ values greater than one, implying that there is more red noise than white noise in these data sets. In the case of the Channel 4 (2007) data set, the γ parameter is driven so low that we conclude that there is no discernible time-correlated noise in the data. We plot the normalized rms residuals versus bin size and verify that they follow the predicted standard error for Gaussian noise (Pont et al. 2006).

Table 3. Derived Eclipse Depths from Spitzer Programs 40280, 60021 (PI: Knutson)

Label	Wavelength	Aperture Size	Ramp	Intra-Pixel	γ	Eclipse Depth
	(μm)	(pixels)	Model	Sensitivity Map		(%)
Channel 1 (2011a)	3.6	2.3	Quadratic	BLISS	0.69	0.119 ± 0.007
Channel 1 (2011b)	3.6	2.6	Quadratic	BLISS	1.70	0.105 ± 0.011
Channel 2 (2010a)	4.5	2.8	Linear	BLISS	1.83	0.132 ± 0.005
Channel 2 (2010b)	4.5	3.0	Quadratic	BLISS	0.95	0.123 ± 0.006
Channel 4 (2007)	8.0	3.6	Quadratic	BLISS	⋅⋅⋅	0.215 ± 0.008

Download table as:  ASCIITypeset image

Once we perform photometry, we fit the eclipse parameters, temporal systematics, BLISS map, white noise, red noise, and gamma value simultaneously to determine our final results and uncertainties. At 3.6 μm, we use an aperture of 2.3 pixels and a quadratic ramp to model the data, and we achieve an eclipse depth of 0.119% ± 0.007%. At 4.5 μm, we use an aperture of 3.0 pixels and a linear ramp to achieve an eclipse depth of 0.123% ± 0.006%. We are able to achieve a good fit and comparable eclipse depths to within 1σ in the 8.0 μm channel, with or without a BLISS map (Figure 4). We choose to include the BLISS map, and for our final analysis, we use an aperture of 3.6 pixels and a quadratic ramp to achieve an eclipse depth of 0.215% ± 0.008% (Table 3).

Zoom In Zoom Out Reset image size

As a final check for our measured eclipse-depth uncertainties, we inject 30 fake transit signals into out-of-eclipse baseline regions of the Channel 1 (2011) phase curve, prior to applying any of our models or the BLISS map. Upon successfully retrieving all of the light curves, we find the mean and distribution of eclipse depths to be consistent with our reported best-fit depth and uncertainty.

In order to characterize the atmosphere of HD 209458b, we perform a Bayesian retrieval analysis of photometric data resulting from our best light curves (Table 4). We choose to perform our retrieval on only these results because they are least affected by correlated noise or, in the case of the 5.8 μm channel, because it is the only available light curve. We do not include our other light curves from the 2005 data set because they were captured in a suboptimal observation mode and have large uncertainties.

For our retrieval analysis, we use the CHIMERA suite (Line et al. 2013; Line & Yung 2013; Line et al. 2014). Briefly, CHIMERA uses three retrieval approaches to determine the allowed range of temperature profiles and abundances that are consistent with the data. We summarize the results from the DEMCMC approach, which is the most comprehensive of the three CHIMERA algorithms. We use a parameterized temperature profile based off of an analytic gray radiative equilibrium solution (e.g., Guillot 2010; Robinson & Catling 2012) and four molecular absorbers, H₂O, CH₄, CO, and CO₂. Further details on the opacity databases and atmospheric parameterization can be found in Line et al. (2013).

We compare our results to those derived by Line et al. (2014), who performed a retrieval analysis of the eclipse depths reported by Knutson et al. (2008). Figure 5 displays the median and the 1σ and 2σ spreads in the spectral fits to our best secondary-eclipse depth measurements (red) as compared to those from Line et al. (2014, blue). Using the Knutson et al. (2008) eclipse depths, Madhusudhan & Seager (2009) and Line et al. (2014) were able to confirm the presence of a thermal inversion; however, when applying the retrieval analysis to our eclipse depths, we find no evidence for a thermal inversion at the pressure regions probed by our observations. We stress, though, that we have little sensitivity at altitudes above the ∼10 m bar level, so it is possible that a weak inversion may persist at these high altitudes with minimal impact on our data.

Zoom In Zoom Out Reset image size

We find that volumetric mixing ratio for water is well bounded from ∼7 × 10⁻⁶ to 5 × 10⁻⁴. This is in stark contrast to the water abundance found with the inversion by Line et al. (2014), where they were only able to achieve an upper limit of 1 × 10⁻⁶ for water. The inversion solution derived from the Knutson et al. (2008) eclipse depths primarily results from the high degree of flux in the 4.5 and 5.8 μm bandpasses. These two photometric bandpasses overlap with both CO and CO₂ absorption. The abundance of CO, which absorbs more strongly at 5.8 μm than does CO₂, must be driven high in order to match these points. The discrepancy between the two solutions is due to the need to suppress the water abundance in the inversion model in order to prevent strong emission in the longer wavelength bandpasses.

We also find an upper limit on the volumetric mixing ratio of methane of ∼1 × 10⁻⁷ from this analysis, consistent with Line et al. (2014). Our results suggest a much lower abundance of CO than reported by Line et al. (2014), again because in the inversion scenario the CO abundance must be driven to an unphysically high abundance in order to produce strong emission at 5.8 μm. All of our volumetric mixing ratios are consistent, to within 1σ, with a solar-composition thermochemical equilibrium atmosphere.

Zellem et al. (2014) analyzed the full phase curve of HD 209458b taken in the 4.5 μm bandpass in 2010 (Spitzer program 60021; PI: Heather Knutson). Although discrepant to our 4.5 μm final result by 1.7σ, their secondary-eclipse depth measurement, reported from a full phase-curve analysis, is consistent with our interpretation of a non-inverted atmosphere.

This analysis of historical Spitzer data serves the dual purpose of exploring the systematics in the data using up-to-date techniques, as well as providing an interpretation of HD 209458b's atmosphere that is consistent with current methods. HD 209458b is known to be the prototypical exoplanet atmospheric thermal inversions; however, the results of our analysis do not corroborate this claim. Our best-fit atmospheric models do not require any species in emission at the pressures probed by the observations in order to explain our measured eclipse depths. This suggests a temperature profile that decreases with increasing altitude.

At the time of writing, there are several new observational campaigns of HD 209458b with warm Spitzer (program 90186, PI: Kamen Todorov; program 10103, PI: Nikole Lewis). We look forward to comparing the results of these new observations with the ones we have derived from the historical data. Ultimately, we look to future spectroscopic observations, from Wide Field Camera 3 on the Hubble Space Telescope or from the heavily anticipated James Webb Space Telescope, to provide data that can better constrain our models, and thereby definitively characterize the thermal profile and atmospheric composition of HD 209458b.

We thank Dorian Abbot for his thoughtful comments. We also thank contributors to SciPy, Matplotlib, the Python Programming Language, and the free and open-source community. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration (NASA). Funding for this work has been provided by NASA grant NNX13AJ16G. J.L.B. acknowledges support from the Alfred P. Sloan Foundation. K.B.S. performed work in part during a Sagan Fellowship supported by NASA and administered by the NASA Exoplanet Science Institute (NExScI).