A TEMPERATURE AND ABUNDANCE RETRIEVAL METHOD FOR EXOPLANET ATMOSPHERES

The P–T profiles of planetary atmospheres are governed by basic physics. Here, we qualitatively describe the physics and shape of a representative P–T profile starting from the bottom of the atmosphere. We focus on the generality that the temperature structure at a given altitude depends on the opacity at that altitude, along with density and gravity.

In the deepest layers of the planet atmosphere, convection is the dominant energy transport mechanism. The high pressure (equivalently, high density) implies a high opacity, making energy transport by convection a more efficient energy transport mechanism than radiation. For hot Jupiters, the dayside radiative–convective boundaries usually occur at the very bottom layers of the atmospheres, at pressures higher than ∼100 bar. In the layers immediately above this boundary, the optical depth is still high enough that the diffusion approximation of radiative transport holds. The diffusion approximation, and the incident stellar flux predominantly absorbed higher up in the atmosphere, lead to a nearly isothermal temperature structure in the layers immediately above the radiative–convective boundary in hot Jupiters. It is assumed here that the energy sources in the planet interior are weak compared to the stellar irradiation. In the solar system giant planet atmospheres (and cooler exoplanet atmospheres as yet to be observed), the radiative–convective boundary occurs at lower pressures in the planet atmosphere than for hot Jupiters.

Above the isothermal diffusion layer, i.e., at lower pressures, the atmosphere becomes optically thin and the diffusion approximation breaks down. These optically thin layers are at altitudes where thermal inversions can be formed, depending on the level of irradiation from the parent star and the presence of strong absorbing gases or solid particles. Thermal inversions, or "stratospheres," are common to all solar system planets and have recently been determined to exist in several hot Jupiter atmospheres (Burrows et al. 2008). For the solar system giant planets, photochemical haze due to CH₄ is the primary cause of stratospheres, and for Earth it is O₃. For hot Jupiters, TiO and VO may be responsible for the thermal inversion (but see Spiegel et al. 2009), but the identification of the absorbers is still under debate. At still lower pressures, below P ∼10⁻⁵ bar, the optical depths eventually become so low that the layers of the atmosphere are transparent to the incoming and outgoing radiation.

We propose a parametric P–T profile, motivated by physical principles, solar system planet P–T profiles, and 1D "self-consistent" exoplanet P–T profiles generated from model atmosphere calculations reported in the literature. Our "synthetic" P–T profile encapsulates stratospheres, along with the low- and high-pressure regimes of the atmosphere. We will refer to our synthetic P–T profile as the "parametric P–T profile."

The atmospheric altitudes of interest are those where the spectrum is formed. Nominally, we consider this range to be between 10⁻⁵ bar ⩽P ⩽ 100 bar, and we will refer to this pressure range as the "atmosphere" in the rest of the paper. At pressures lower than 10⁻⁵ bar, the atmosphere layers are nearly transparent to the incoming and outgoing radiation at visible and infrared wavelengths. And, above pressures of ∼100 bar (or even above pressures of ∼10 bar), the optical depth is too high for any radiation to escape without being reprocessed. The corresponding layers do not contribute significantly to the emergent spectrum. Figure 1 shows a schematic parametric P–T profile. In our model, we divide the atmosphere into three representative zones or "Layers" as shown in Figure 1. The upper-most layer, Layer 1, in our model profile is a "mesosphere" with no thermal inversions. The middle layer, Layer 2, represents the region where a thermal inversion (a "stratosphere") is possible. And, the bottom-most layer, Layer 3, is the regime where a high optical depth leads to an isothermal temperature structure. Layer 3 is used with hot Jupiters in mind; for cooler atmospheres this layer can be absent, with Layer 2 extending to deeper layers.

Zoom In Zoom Out Reset image size

Our proposed model for the P–T structure in Layers 1 and 2 is a generalized exponential profile of the form

$$\begin{equation} P = P_oe^{\alpha (T - T_o)^\beta }, \end{equation} \tag{ 1 }$$

where, P is the pressure in bars, T is the temperature in K, and P_o, T_o, α, and β are free parameters. For Layer 3, the model profile is given by T = T₃, where T₃ is a free parameter.

Thus, our parametric P–T profile is given by

$$\begin{eqnarray} P_0 < P < P_1 : P = P_{0} e^{\alpha _1(T - T_0)^{\beta _1}}\; &{\rm Layer 1}, \nonumber \\ P_1 < P < P_3 : P = P_2 e^{\alpha _2(T - T_2)^{\beta _2}} \;&{\rm Layer 2}, \\ P > P_3 : T = T_3\; &{\rm Layer 3} \nonumber. \end{eqnarray} \tag{ 2 }$$

In this work, we empirically set β₁ = β₂ = 0.5 (see Section 2.3). The space of P–T profiles is spanned in all generality with the remaining parameters, rendering β redundant as a free parameter. Then, the model profile in Equation (2) has nine parameters, namely, P₀, T₀, α₁, P₁, P₂, T₂, α₂, P₃, and T₃. Two of the parameters can be eliminated based on the two constraints of continuity at the two layer boundaries, i.e., Layers 1–2 and Layers 2–3. And, in the present work, we set P₀ = 10⁻⁵ bar, i.e., at the top of our model atmosphere. Thus, our parametric profile in its complete generality has six free parameters.

Our P–T profile consists of 100 layers in the pressure range of 10⁻⁵–100 bar, encompassing the three zonal "Layers" described above. The 100 layers are uniformly spaced in log(P). For a given pressure, the temperature in that layer is determined from Equation (2), using the form T = T(P). The kinks at the layer boundaries are removed by averaging the profile with a box-car of 10 layers in width.

The parametric P–T profile described in Section 2.2 fits the actual temperature structures of a wide variety of planetary atmospheres. Figure 2 shows the comparison of our model P–T profile with published P–T profiles of several solar system planets and hot Jupiters. For each case, the published profile was fitted with our six-parameter P–T profile using a Levenberg–Marquardt fitting procedure (Levenberg 1944; Marquardt 1963). We were able to find best fits to all the published profiles with the β parameter fixed at 0.5, indicating that β is a redundant parameter. Therefore, we set β₁ = β₂ = 0.5 for all the P–T profiles in this work. In addition, it must be noted that, for the parametric P–T profile to have a thermal inversion, the condition P₂>P₁ must be satisfied, and for one without a thermal inversion, P₁ ⩾ P₂.

Zoom In Zoom Out Reset image size

For the solar system planets (top left panel of Figure 2), published profiles show detailed temperature structures obtained via several direct measurements and high-resolution temperature retrieval methods (Lindal et al. 1981, 1985, 1987, 1990; COESA 1976; http://atmos.nmsu.edu/planetary_datasets/). For hot Jupiters, on the other hand, observations are limited. Consequently, the published profiles for hot Jupiters (Figure 2) are obtained from self-consistent 1D models reported in the literature attempting to explain the available data on each system. As is evident from Figure 2, the parametric P–T profile provides good fits to all the known profiles. We emphasize that, while the fits alone are not a test of our retrieval method, they show that our parametric P–T profile is generally enough to fit the wide range of planetary atmospheres represented in Figure 2.

We also tested the emergent spectrum obtained with our parametric P–T profile by comparing it with a self-consistent model. In this regard, we used a P–T profile from a self-consistent atmosphere model (Seager et al. 2005) to generate an emergent spectrum for HD 209458b. The emergent spectrum calculated by our model (for the same molecular composition) agrees with that of the self-consistent model.

The fundamental significance of our parametric P–T profile is the ability to fit a disparate set of planetary atmosphere structures, with very different atmospheric conditions (Figure 2). In particular, the published hot Jupiter P–T profiles were obtained from several different modeling schemes. The different planetary atmospheres all conform to a basic mathematical model because of the common physics underlying the P–T profiles. We anticipate that our parametric P–T profile will open up a new avenue in retrieving P–T structure of exoplanetary atmospheres by enabling an efficient exploration of the temperature structures and compositions allowed by the data.

An important point concerns the adaptability of the model proposed in this work. The introduction of the isothermal Layer 3 in Equation (2) is motivated by the fact that hot Jupiters, which are the prime focus of this work, have convective–radiative boundaries that are deep in the atmosphere, and the presence of an isothermal Layer 3 is physically plausible (as is evident from self-consistent calculations in literature). However, the model can be adapted to cases of cooler atmospheres where the convective–radiative boundary occurs at lower pressures than for hot Jupiters, and where there may not be an isothermal layer. The appropriate model P–T profile in such cases would be one with only the two upper layers (Layers 1 and 2). As can be seen from Figure 2, model fits to the solar system planets belong to the category of models with only Layers 1 and 2.

We calculate emergent spectra using a 1D line-by-line radiative transfer code, assuming LTE. Our model atmosphere consists of 100 layers in the pressure range between 10⁻⁵ bar and 100 bar, and uniformly spaced in log(P). The key aspect of our model is the parameterization of the P–T profile and the chemical composition. In each layer, the parametric P–T profile (described in Section 2.2) determines the temperature and pressure. And, the density is given by the ideal gas equation of state. The radial distance is determined by the requirement of hydrostatic equilibrium. And, the equation of radiative transfer, assuming LTE and no scattering, is then solved along with the adopted molecular abundances (see Section 3.2 below) to determine the emergent flux. The geometry and the radiative transfer formulation are modified appropriately for calculating transmission spectra.

Our model directly computes the thermal emission spectrum and indirectly accounts for stellar heating and scattered radiation. In our model, we assume that the effects of heating by stellar radiation are included in the allowed "shapes" of the parametric P–T profile. As described in Section 2.3, our parametric P–T profile fits the P–T profiles obtained from self-consistent models that directly include absorbed stellar radiation, including profiles with thermal inversions.

One reason the approach of indirectly including stellar irradiation is valid is that the visible and infrared opacities are largely decoupled (Marley et al. 2002). For example, Na, K, TiO, VO can have strong absorption at visible wavelengths in brown dwarfs and exoplanets, yet are mostly negligible absorbers in the infrared. H₂O dominates infrared absorption, and CH₄ and CO are also spectroscopically active in the IR. Another reason the indirect account for stellar irradiation via a parametric P–T profile is justified is the uncertain chemistry (via unknown optical absorbers that in some cases cause thermal inversions) and physics that is hard to account for directly in a self-consistent 1D model. We note that it is not known whether the common assumption of radiative equilibrium in self-consistent models, while one possible way of including stellar irradiation, holds for hot Jupiter atmospheres.

We include molecules of hydrogen (H₂), water vapor (H₂O), carbon monoxide (CO), carbon dioxide (CO₂), methane (CH₄), and ammonia (NH₃). Our H₂O, CH₄, CO, and NH₃ molecular line data are from Freedman et al. (2008), and references therein. Our CO₂ data are from R. S. Freedman (2009, private communication) and Rothman et al. (2005). And, we obtain the H₂–H₂ collision-induced opacities from Borysow et al. (1997) and Borysow (2002).

Molecular abundances in exoplanet atmospheres may depart from chemical equilibrium (Liang et al. 2003; Cooper & Showman 2006). Our code, therefore, has parametric prescriptions to allow for non-equilibrium quantities of each molecule considered. We begin by calculating a fiducial concentration of each molecule in a given layer, based on the assumption of chemical equilibrium. For chemical equilibrium we use the analytic expressions in Burrows & Sharp (1999), and assume solar abundances for the elements. The molecules we include are H₂O, CO, CH₄, CO₂, and NH₃. The molecule CO₂ is thought to originate from photochemistry and no simple expression for equilibrium composition is available; we therefore choose an arbitrary fiducial concentration. We denote the concentration of each molecule by its mixing ratio; i.e., number fraction with respect to molecular hydrogen.

The concentration of each molecule in a given layer is calculated by multiplying the corresponding fiducial concentration by a constant factor specific to each molecule. The constant factor corresponding to each molecule is a parameter of the model. Therefore, corresponding to the four prominent molecules, H₂O, CO, CH₄, and CO₂, in our model, we have four parameters: $f_{\rm {H_2O}}$, f_CO, $f_{\rm {CH_4}}$, and $f_{\rm {CO_2}}$. For instance, $f_{\rm {H_2O}}$ is the ratio of the concentration of H₂O to the fiducial concentration of H₂O in each layer. Thus, although the concentration of H₂O is different in each layer, it maintains a constant ratio ($f_{\rm {H_2O}}$) with the fiducial concentration in that layer. By varying over a range in these four parameters, we are essentially varying over a wide range of molecular mixing ratios and elemental abundances.

For NH₃, we restrict the concentration to the fiducial equilibrium concentration, i.e., $f_{\rm {NH_3}} = 1$. This is because, NH₃ has limited and weak spectral features in the Spitzer and HST bands under consideration in this work. Additionally, NH₃ is not expected to be abundant in enhanced quantities at such high temperatures as are seen in hot Jupiters.

Clouds are not included in our model because HD 209458b and HD 189733b likely do not have clouds that scatter or emit at thermal infrared wavelengths. Several reasons have been proposed in literature justifying the use of cloud-free models for these planets: finding weak effects of clouds on the P–T profile (Fortney et al. 2006, 2008); fast sedimentation rates in radiative atmospheres (Barman et al. 2005); clouds forming too deep in the planetary atmosphere to affect the emergent spectrum (Fortney et al. 2006); spectral features in transmission dominated by molecular absorption over cloud particles (Showman et al. 2009, and references therein); haze absorption at blue wavelengths on HD 189733b drops off rapidly with increasing wavelength and should be negligible at IR wavelengths (Pont et al. 2008). That said, clouds can be included in our model as another free parameter, via a wavelength-dependent opacity, and we plan to make this extension in the future.

A planetary atmosphere model must satisfy the fundamental constraint of energy conservation (here called energy balance). We consider a 1D plane-parallel atmosphere with normal incidence. While using a plane-parallel model, it is imperative to weight the incident stellar flux appropriately in order to represent an average dayside atmosphere.

The average stellar flux incident on the planetary dayside is given by F_s = fF_⋆ where, F_⋆ is the wavelength integrated stellar flux at the substellar point of the planet. f is a geometric factor by which the stellar flux at the substellar point must be weighted so as to represent an average 1D plane-parallel incidence. A value of ${\it f} = \frac{1}{2}$ represents uniform distribution of stellar flux on the planet dayside; it comes from F_⋆ × (πR²_p)/(2πR²_p), i.e., the ratio of the planetary surface projected perpendicular to the incident stellar flux to the actual area of the planetary surface receiving the flux. However, it has been shown that just before secondary eclipse the dayside average flux should be biased toward a higher value of ${\it f} = \frac{2}{3}$ (Burrows et al. 2008). Therefore, in the present work, we adopt a value of ${\it f} = \frac{2}{3}$.

In the context of our 1D atmosphere, the constraint of energy balance means that the wavelength integrated emergent flux at the top of the planetary atmosphere matches the wavelength integrated incident stellar flux, after accounting for the Bond albedo (A_B) and possible redistribution of energy onto the night side (f_r). Our prescription for the day–night redistribution is similar to that of the P_n prescription of Burrows et al. (2006). We define f_r as the fraction of input stellar flux that is redistributed to the night side. The input stellar flux is given by F_in = (1 − A_B)F_s, i.e., the incident stellar flux (F_s), less the reflected flux (A_BF_s). Therefore, the energy advected to the night side is given by f_rF_in. And hence, the amount of flux absorbed on the dayside is given by F_in,day = (1 − f_r)F_in, which equals (1 − A_B)(1 − f_r)F_s. It is assumed that the intrinsic flux from the planet interior is negligible compared to the incident stellar irradiation.

Thus, the energy balance requirement on the dayside spectrum is given as

$$\begin{equation} F_p = (1-A_B)(1-f_r)F_s, \end{equation} \tag{ 3 }$$

where, F_p is the wavelength integrated emergent flux calculated from the model spectrum. We use a Kurucz model for the stellar spectrum (Castelli & Kurucz 2004, ³). We adopt the stellar and planetary parameters from Torres et al. (2008).

We use energy balance, Equation (3), in two ways. The first, as described above, is to ensure that our model atmosphere satisfies energy conservation. The constraint on the model comes from the fact that A_B and f_r are bounded in the range [0, 1]. Thus, models for which η = (1 − A_B)(1 − f_r) is greater than unity are discarded on grounds of energy balance. Second, for models which satisfy energy balance, Equation (3) gives us an estimate of η for a given model. If an A_B is assumed, this gives a constraint on f_r, and vice versa.

The parameter space of our model consists of N = N_profile + N_molec free parameters, where N_profile is the number of parameters in the P–T profile, and N_molec is the number of parameters pertaining to the chemical composition. As described in Section 2, our parametric P–T profile has N_profile = 6, the six parameters being T₀, α₁, P₁, P₂, α₂, and P₃. And, corresponding to the four prominent molecules, H₂O, CO, CH₄, and CO₂, in our model, we have N_molec = 4 parameters: $f_{\rm {H_2O}}$, f_CO, $f_{\rm {CH_4}}$, and $f_{\rm {CO_2}}$ (see Section 3.2). Thus, our model has 10 free parameters.

Given a set of observations of a planet, for example, a dayside spectrum, our goal is to determine the regions of parameter space constrained by the data, and to identify the degeneracies therein. We investigate models with and without thermal inversions. To this effect, we run a large number of models on a grid in the parameter space of each scenario. Each parameter in either scenario has a finite range of values to explore, and we have empirically determined the appropriate step size in each parameter. A grid of reasonable resolution in this 10-parameter space can approach a grid-size that is computationally prohibitive.

Even before running the grid of models, however, the parameter space of allowed models can be constrained to some extent based on the data and energy balance. First, the model planet spectrum is bounded by two blackbody spectra, corresponding to the lowest and highest temperatures in the atmosphere. Correspondingly, the maximum and minimum temperatures, T_max and T_min, of the P–T profile are partly constrained by the maximum and minimum blackbody temperatures admissible by the data, and allowing a margin accounting for the fact that observations are available only at limited wavelengths. Second, an upper limit can be placed on the temperature at the base of our model atmosphere (T₃ in the P–T profile) based on energy balance. The maximum of T₃ is given by the effective temperature of the planet dayside assuming zero albedo and zero redistribution of the energy to the night side, and allowing for a factor close to unity (∼1–1.5) to take into account the changes in the blackbody flux due to line features. Therefore, given the data, the space of P–T profiles can be constrained to some extent even before calculating the spectra, using the above conditions.

The final grid is decided by running several coarse grids, and obtaining an empirical understanding of the parameter space. At the end, a nominal grid in 10 parameters, with or without inversions in the P–T profiles, consists of ∼10⁷ models. We thus computed a total of ∼10⁷ models for each planet under consideration. We ran the models on ∼100 2 GHz processors on a Beowulf cluster, at the MIT Kavli Institute for Astrophysics and Space Research, Cambridge, MA.

Broadband observations of exoplanet atmospheres are currently limited to the six channels of broadband infrared photometry using Spitzer. Additionally, a few narrowband spectrophotometry have also been reported using the HST NICMOS G206 grism (1.4–2.6 μm), and the Spitzer IRS spectrograph (5–14 μm). The number of broadband observations available in a single data set are typically smaller than the number of parameters in our model. Also, given that the parameters in our model are highly correlated and the degeneracies not well understood, the nature of existing data does not allow spectral fitting in the formal sense.

Consequently, our approach in this work is to compute a large number of models on a pre-defined grid in the parameter space. For each model, we calculate a "goodness-of-fit" with the data using a statistic defined by a weighted mean squared error, given by

$$\begin{equation} \xi ^2 = \frac{1}{N_{\rm obs}} \sum _{i = 1}^{N_{\rm obs}} \bigg (\frac{f_{i,\rm model} - f_{i,\rm obs}}{\sigma _{i,\rm obs}} \bigg)^2. \end{equation} \tag{ 4 }$$

Here, f_i is the flux observable and N_obs is the number of observed data points. f_i,obs and f_i,model are the observed f_i and the corresponding model f_i, respectively. For secondary eclipse spectra, f_i is the planet-star flux contrast, whereas, for transmission spectra, f_i is the transmission. σ_i,obs is the 1σ measurement uncertainty in the observation. In order to obtain the model fluxes in the same wavelength bins as the observation, we integrate the model spectra with the transmission functions of the instruments. Where spectrophotometry is available, we additionally convolve our model spectra with the instrument point-spread function.

The goodness-of-fit (ξ²), as defined in Equation (4), is similar in formulation to the conventional definition of the reduced χ² statistic, if N_obs is replaced by the number of degrees of freedom. In the current context, we refrain from using the χ² statistic to asses our model fits because our number of model parameters (N) are typically more than the number of available broadband data points, leading to negative degrees of freedom. A value of χ² with negative degrees of freedom lacks meaning in that we cannot relate it to a confidence level. Therefore, in using ξ² as our statistic of choice, we are essentially calculating χ² per data point as opposed to χ² per degree of freedom. We discuss this further in Section 6.2.

HD 189733b has the best data set for any hot Jupiter known to date, including both signal-to-noise and wavelength coverage. Several remarkable observations of the planet dayside have been made by determining the planet-star flux ratios just before secondary eclipse. Grillmair et al. (2008) reported an infrared spectrum of the dayside in the 5–14 μm range, obtained with the Spitzer IRS instrument. Charbonneau et al. (2008) presented broadband photometry in the six Spitzer channels (3.6, 4.5, 5.8, 8, 16, and 24 μm).⁴ Separate photometric observations in the Spitzer 8 μm and 24 μm channels were also reported by Knutson et al. (2007, 2009a), respectively. Additionally, an upper limit was placed on the 2.2 μm flux contrast by Barnes et al. (2007). More recently, Swain et al. (2009a) observed the planet dayside with HST NICMOS spectrophotometry in the 1.5–2.5 μm range.

Together, these data sets place important constraints on the dayside atmosphere of the planet. The high signal-to-noise ratio (S/N) IRS spectrum reported by Grillmair et al. (2008), obtained with 120 hr of integration time, is the best exoplanet spectrum known. This spectrum forms the primary data set for our present study of HD 189733b. In addition, we also use the photometric observations of Charbonneau et al. (2008), and the spectrophotometric observations of Swain et al. (2009a) to place complementary constraints.⁵ We use the two independent observations of Knutson et al. (2007, 2009a), at 8 μm and 24 μm, only as a guide for the models, and exclude them from our fits. This is because we place constraints from each individual data set separately, for reasons explained below, and we cannot place any constraints from a single data point from each of Knutson et al. (2007, 2009a). However, it must be noted that the 8 μm and 24 μm values of Knutson et al. (2007, 2009a) agree with the IRS spectrum and the Charbonneau et al. (2008) observations, respectively, which are being used in the current work. We present the data in Figure 3, along with model fits.

Zoom In Zoom Out Reset image size

We analyze each data set separately. The reason we use a separate treatment is that some individual observations at similar wavelengths differ from each other. Specifically, the 8 μm Spitzer Infrared Array Camera (IRAC) data point from Charbonneau et al. (2008) is ∼2σ above the IRS spectrum at the same wavelength (Grillmair et al. 2008). Also, 8 μm IRAC observations taken at separate times differ by over 2σ (Charbonneau et al. 2008). From a model interpretation viewpoint, a 2σ spread means that almost no constraints on molecular abundances can be drawn (based on the reported uncertainties in the data). In addition, the very strong evidence of CO₂ that is visibly obvious around 2 μm in the HST/NICMOS data set is not so strong in the 4.5 μm and 16 μm Spitzer channels. In this work, we are unable to find models that fit all the three data sets of HD 189733b simultaneously within the ξ² = 2 level. The data and model results below may well indicate atmospheric variability in this planet. See Section 6.3 for further discussion about possible atmospheric variability.

In what follows, we describe the constraints on the atmospheric properties of HD 189733b placed by each data set. Ideally, we would like to report constraints at the ξ² = 1 surface, i.e., only those models which, on an average, fit the observations to within the 1σ error bars. However, we adopt a more conservative approach and lay emphasis on the results at the ξ² = 2 surface as well, allowing for models fitting the observations to within ∼1.4σ on an average. The latter choice is somewhat arbitrary, and the informed reader could choose to consider other ξ² surfaces; we show contours of ξ² in the different atmospheric properties.

4.1.1. Molecular Abundances

The constraints on the abundances due to the different data sets are shown in Figure 4. The figure shows contours of ξ² in the parameter space of our model atmosphere. The contours show minimum ξ² surfaces. At each point in a given two-parameter space, several values of ξ² are possible because of degeneracies with other parameters in the parameter space. In other words, the same point in a pair of parameters can have multiple values of ξ² corresponding to evaluations at various values in the remaining parameters. For the contours in Figure 4, we consider the minimum possible value of ξ² at each point in a two-parameter space. Consequently, each colored surface shows the region of parameter space which allows a minimum possible ξ² corresponding to that color. The mixing ratios of species are shown as ratios by number.

Zoom In Zoom Out Reset image size

Spitzer IRS spectrum. We begin with the IRS spectrum (Grillmair et al. 2008) because this data set has the broadest wavelength coverage and the best S/N. Also, the IRS data set has more data points than the number of parameters in our model. Our most significant new finding is the existence of models fitting the IRS data set that are totally consistent with efficient redistribution of energy from the planet dayside to the night side (see Figure 4 for the range of allowed values). This addresses a significant problem in previously published models for dayside spectra (Grillmair et al. 2008) not being able to explain the efficient redistribution found in the phase curves of Knutson et al. (2007, 2009a; see Charbonneau et al. 2008 and Grillmair et al. 2008; but also see Barman 2008). With the IRS data, we also confirm the absence of a thermal inversion in HD 189733b.

Despite the high S/N, the IRS spectrum allows only weak constraints on the molecular abundances. The reason being that the IRS spectral range has molecular absorption features of only CH₄ and H₂O; even for these molecules, their features at other wavelengths remain completely unconstrained. The significant features are of CH₄ (7.6 μm), and H₂O (6.3 μm). For example, the ξ² = 2 surface constrains the mixing ratio of water vapor to 10⁻⁶–0.1, and that of CH₄ to be less than 10⁻² (but cf. Section 6.4).

Spitzer broadband photometry. Our model fits to the Spitzer broadband data (Charbonneau et al. 2008) require less efficient circulation of incident stellar energy compared to the best model fits to the Spitzer IRS spectrum. Specifically, the value for redistribution in our framework is f_r ≲ 0.25, and this agrees with previous results (Barman 2008; Charbonneau et al. 2008) which find that model fits to the broadband data require inefficient day–night redistribution of energy. The discrepancy in redistribution arises from different planet-star contrast ratios. The 8 μm IRAC point of Charbonneau et al. (2008) is ≳2σ different from the IRS spectrum, and from the 8 μm IRAC measurement of Knutson et al. (2007). Similar to the IRS spectrum, the Spitzer broadband photometry shows no sign of a thermal inversion.

The Spitzer broadband photometry enables significant constraints on molecular abundances. Again, we take the ξ² = 2 surface, corresponding to an average difference of 1.4σ between the model and data. We find the lower limit of the H₂O mixing ratio to be 10⁻⁶, and the upper limit of the CH₄ mixing ratio to be 10⁻⁵ (Figure 4). At the ξ² = 2 surface there are no constraints on the presence or abundances of CO and CO₂. At the ξ² = 1 surface, the best-fit models constrain the molecular mixing ratios as follows: 10⁻⁵⩽ H₂O ⩽10⁻³; CH₄ ⩽ 2 × 10⁻⁶; 7 × 10⁻⁷⩽ CO₂ ⩽ 7 × 10⁻⁵; CO is unconstrained because it is degenerate with CO₂. The identification of CO₂, at the ξ² = 1 surface, is intriguing, and this is the first time CO₂ is identified with Spitzer photometry. Our model also finds the C/O ratio to be 7 × 10⁻³ ⩽ C/O ⩽ 1 for both the ξ² = 1 and ξ² = 2 surfaces.

HST/NICMOS spectrophotometry. The HST/NICMOS data (Swain et al. 2009a) are useful because the wavelength range includes significant features of H₂O, CH₄, CO, and CO₂ (Figure 3). We can therefore constrain the abundances even at the ξ² = 2 surface (Figure 4). We find mixing ratios of: H₂O ∼10⁻⁴; CH₄ ⩽ 6 × 10⁻⁶; 2 × 10⁻⁴⩽ CO ⩽2 × 10⁻²; and CO₂ ∼ 7 × 10⁻⁴. It is remarkable that at the ξ² = 2 surface we can constrain the abundances so stringently. The exception is CH₄, where the upper limit indicates that methane is not present in large quantities on the dayside of HD 189733b. The C/O ratio is between 0.5 and 1. Our inferred abundances are consistent with those reported in Swain et al. (2009a), with the major exception of CO₂. While we find CO₂ mixing ratios of 7 × 10⁻⁴, Swain et al. (2009a) find 10⁻⁷⩽ CO₂ ⩽ 10⁻⁶. At present we have no explanation for this discrepancy.

A mixing ratio of 7 × 10⁻⁴ is seemingly high from a theoretical standpoint; photochemistry which is thought to be the primary source of CO₂ in hot Jupiters allows mixing ratios of CO₂ up to ∼10⁻⁶, assuming solar metallicities (Liang et al. 2003; Zahnle et al. 2009). One possible explanation of a large inferred CO₂ mixing might be a high mixing ratio concentrated in only a few layers of the atmosphere. The data can be fit equally well with the same mixing ratio (7 × 10⁻⁴) of CO₂ present only between 10⁻³ and 0.3 bar causing an average mixing ratio of 10⁻⁶ over the entire atmosphere.

4.1.2. Pressure–Temperature Profiles

Our results include constraints on the 1D averaged P–T profile. We explore the parameter space with a grid of ∼10⁴P–T profiles, shown in Figure 5. For P–T profiles with no thermal inversions, the parameters that influence the spectrum most are the two pressures (P₁ and P₃), which contribute the temperature differential in Layer 2 of the atmosphere, and the corresponding slope (α₂). Additionally, P₃ governs the base flux, corresponding to the blackbody flux of temperature T₃, on which the absorption features are imprinted.

Zoom In Zoom Out Reset image size

The constraints on the parameters of the P–T profiles (Section 2.2) from the three data sets are shown in Figure 4. The IRS spectrum (Grillmair et al. 2008) constrains the P–T profiles toward larger values of α₂, leading to less steep temperature gradients. For the same base temperature (T₃) of the planet dayside, the maximum temperature gradients allowed by the broadband photometry and the HST spectrum are always higher than that allowed by the IRS data set. The ξ² = 2 surface for the IRS data set constrains α₂ to be greater than 0.2, while the constraints due to the broadband photometry and the HST data are 0.18–0.35 and 0.15–0.27, respectively. The same surface constrains the mean temperature gradient (<dT/dlog P>) in Layer 2 to 100–400 K, 180–440 K, and 290–470 K for the IRS, the broadband photometry, and the HST data sets, respectively. The ξ² = 2 surface also shows that the IRS data allow for lower base temperatures as compared to the other data sets. The skin temperature (T₀) at the top of the atmosphere (at P ∼ 10⁻⁵ bar) is relatively unconstrained because of its degeneracy with the temperature gradient, and because of the low contribution of spectral features at that pressure.

4.1.3. Albedo and Energy Redistribution

The transmission spectrum of HD189733b was presented by Swain et al. (2008b), who reported the detection of CH₄ and H₂O in the 1.4–2.5 μm range covered by the HST/NICMOS spectrophotometry. Transmission data are complementary to secondary eclipse data, because a transmission spectrum probes the part of the planet atmosphere in the vicinity of the limb, whereas secondary eclipse data probe the planetary dayside. In addition, a transmission spectrum reveals line features of the species in absorption, imprinted on the stellar flux, with the self-emission of the planet contributing negligibly to the observed flux.

We computed ∼7 × 10⁵ models to place constraints on allowed limb P–T profiles and quantify the allowed abundances by the ξ² surfaces. In order to compute the ξ², our models were convolved with the instrument point-spread function, integrated over the grism transmission function, and finally binned to the same bins as the data (as described in Swain et al. 2008b). We find that 2 of 18 binned data points (at 1.888 μm and 2.175 μm) as reported by Swain et al. (2008b) have unusually high values, and amount to an unusually large contribution to the net ξ². In the framework of our models, we do not find any spectral features that might be able to explain those values. Consequently, we leave these two points out of our present analysis. Figure 6 shows the model spectra and data.

Zoom In Zoom Out Reset image size

4.2.1. Molecular Abundances

Our results confirm the presence of H₂O and CH₄ in the limb, as also reported by Swain et al. (2008b). At the ξ² = 2 surface, as shown in Figure 7, the H₂O mixing ratio is constrained to 5 × 10⁻⁴–0.1, and the CH₄ mixing ratio lies between 10⁻⁵ and 0.3. These ranges of mixing ratios for H₂O and CH₄ are consistent with the mixing ratios found for the best-fit model reported by Swain et al. (2008b). Although, these constraints are roughly consistent with the constraints on the same molecules on the dayside inferred in Section 4.1, the high mixing ratios of CH₄ allowed is appalling at first glance. This hints at a region of high CH₄ concentration in the limb, possibly on the night side. The effect of other possible molecules, CO, NH₃, and CO₂ on the spectrum is minimal. Although each of these molecules have some spectral features in the observed wavelength range, the upper limits on their compositions placed by the data are not statistically significant.

Zoom In Zoom Out Reset image size

4.2.2. Pressure–Temperature Profiles

Figure 8 shows the P–T profiles explored for the HD 189733b transmission data set. The best-fit models with ξ² ⩽ 2 are shown in red. Figure 7 shows the ξ² surface in the parameters space of the P–T profiles. The best-fit models are consistent with an atmosphere with no thermal inversions. The α₂ parameter is constrained by the ξ² = 2 surface to be greater than ∼0.15, and the temperature T₃ at the limb is constrained to lie between 1500 K and 2250 K. The mean temperature gradient (<dT/dlog P>) in Layer 2 constrained by the ξ² = 2 surface lies between 130 K and 550 K.

Zoom In Zoom Out Reset image size

Our results indicate the presence of H₂O, CO, CH₄, and CO₂ in the atmosphere of HD 209458b. As will be discussed in Section 5.2, the allowed P–T profiles show the existence of a deep thermal inversion in the atmosphere of HD 209458b. Consequently, all the line features of the molecules are seen as emission features, rather than absorption features as in the case of HD 189733b. The constraints on the concentrations of all the molecules are shown in Figure 10.

Zoom In Zoom Out Reset image size

The ξ² = 2 surface places an upper limit on the mixing ratio of H₂O at 10⁻⁴, with the ξ² = 1 surface allowing mixing ratios up to 10⁻⁵. The limits on H₂O are governed by the band features of the molecule in the 5.8 μm and 24 μm channels. While the upper limit on the mixing ratio of H₂O is well constrained, the constraint on the lower limit is rather weak, allowing values as low as 10⁻⁸ at ξ² = 1. The stringent upper limit on H₂O comes primarily from the low planet-star flux contrast observed in the 24 μm channel. Similarly, the high flux contrast observed in the 5.8 μm channel should, in principle, yield a stringent lower limit. However, the weak lower limit is because of the large observational uncertainty on the flux contrast in the 5.8 μm channel. The 1σ error bar on the observed 5.8 μm flux contrast is the largest of uncertainties in all the channels, and 1.5 times larger than the error bar on the 24 μm flux contrast.

The ξ² = 2 constraint on the mixing ratio of CH₄ is 10⁻⁸–0.04, based on band features in the 3.6 μm and 8 μm channels. Although, the features in the 3.6 μm channel are degenerate with some features of H₂O in the same channel, the exclusive CH₄ features in the 8 μm channel break the degeneracy. The ξ² = 2 surface places an absolute lower limit on CO at 4 × 10⁻⁵, and requires a mixing ratio of CO₂ between 2 × 10⁻⁹ and 7 × 10⁻⁶. The need for CO₂ arises from the high contrast observed between the 4.5 μm and the 3.6 μm channels. The 4.5 μm channel could in principle be explained by CO alone. However, such a proposition would require large amounts of CO, close to mixing ratios of 1. CO₂ provides additional absorption in the same 4.5 μm channel, with a relatively reasonable abundance. In addition, the amount of CO₂ is also constrained by the 16 μm channel—too high of a CO₂ abundance would not agree with the data because CO₂ has a strong absorption feature in the 16 μm channel. Thus, the 4.5 μm and 16 μm channels together constrain the amounts of CO₂, thereby also constraining the CO feature in the 4.5 μm channel. If the 16 μm data did not exist, there would have been a large degeneracy in the abundances of CO and CO₂.

If we consider the ξ² = 1 surface, instead of the ξ² = 2 surface, our results place tighter constraints on all the molecules. At this surface, the constraints on the mixing ratios are: 10⁻⁸⩽ H₂O ⩽10⁻⁵; 4 × 10⁻⁸⩽ CH₄ ⩽ 0.03; CO ⩾4 × 10⁻⁴; and 4 × 10⁻⁹⩽ CO₂ ⩽ 7 × 10⁻⁸. However, in considering the ξ² = 1 surface, we might be overinterpreting the data.

We find definitive detections of CH₄, CO, and CO₂. Even at the ξ² = 4 surface, the data require non-zero mixing ratios for the carbon-bearing molecules. Furthermore, the results suggest a high C/O ratio (see Seager et al. 2005). Although the model fits give C/O ratio close to 1 or higher, our molecular abundance grid resolution is not high enough to put a precise lower limit to the C/O ratio. As seen from Figure 10, the allowed models reach values as high as ∼60 for ξ² = 1, and ∼600 for ξ² = 2. As discussed previously, and in Section 6.1, high molecular concentrations and high C/O ratios could be localized effects overrepresented by the 1D averaged retrieval.

The most interesting aspect about the atmosphere of HD 209458b is the presence of a thermal inversion on the planet dayside, as previous studies have pointed out (e.g., Burrows et al. 2008). While modeling atmospheric spectra, one often encounters a degeneracy between temperature structure and molecular compositions. However, by exploring a large number of P–T profiles we find that the dayside observations for HD 209458b cannot be explained without a thermal inversion, for any chemical composition, at the ξ² < 2 level. On the other hand, our results place stringent constraints on the 1D averaged structure of the required thermal inversion.

Our results presented here explore ∼7000 P–T profiles with thermal inversions. Figure 11 shows all the P–T profiles with ξ² ⩽ 5. The ranges in the P–T parameters explored are shown in Figure 10. In the presence of a thermal inversion, the spectra are dominated by the location and extent of the inversion layer, quantified by the parameters P₂, T₂, and α₂. We find the models with ξ² ⩽ 2 to have an inversion layer in the lower atmosphere, at pressures (P₂) above ∼0.03 bar, with the best-fit models preferring the deeper layers, i.e., higher P₂. We find that the ξ² ⩽ 1 region extends to inversion layers as deep as P₂ ∼ 1 bar. However, we limit P₂ to a value of 0.5 based on the physical consideration that, at pressure ≳ 1 bar collisional opacities are likely to dominate over the molecular opacities in the optical which are expected to be the cause of thermal inversions. And, therefore, we do not expect inversions at such high pressures.

Zoom In Zoom Out Reset image size

The allowed models (Figure 10) cover a wide range in the base temperature (T₃) owing to the fact that the spectrum is dominated by the maximum temperature differential (T₁ − T₂) attained in the stratosphere, and is less sensitive to T₃ by itself. The stratospheric temperatures of all the models with ξ² ⩽ 2 lie between 800 K and 2200 K. The temperature gradient in the inversion layer is governed by the α₂ parameter. As can be seen from Figure 10, the data restrict the ξ² ⩽ 2 models to have α₂ between 0.05 and 0.25, with the ξ² ⩽ 1 models preferring values between 0.05 and 0.15. These values of α₂ translate into mean temperature gradients (<dT/dlog P>) in the inversion layer to be 230–1480 K for ξ² ⩽ 2, and 400–1480 K for ξ² ⩽ 1. Consequently, large thermal inversions spanning more than a few dex in pressure, and hence with less steeper thermal gradients, are ruled out by the data. This is also evident from Figure 11. Finally, ξ² is not very sensitive to the parameters in Layer 1 of the P–T profile, since the dominant spectral features in this planet arise from the inversion layer.

Our results also show the constraints on the albedo and day–night redistribution. Figure 10 shows the ξ² surface in the space of T_eff versus η, where, η = (1 − A_B)(1 − f_r); A_B is the Bond albedo, and f_r is the fraction of input stellar flux redistributed to the night side. As described in Section 4.1.3, using η one can estimate the constraint on A_B given f_r, or vice versa. For HD 209458b, the ξ² = 2 surface in Figure 10 constrains η in the range 0.24–0.81, with the ξ² = 1 surface in the range 0.26–0.72. If we assume an A_B of 0.24 (2σ upper limit, scaled from the 1σ estimate of Rowe et al. 2008), the constraints on f_r for the ξ² = 2 and ξ² = 1 surfaces are 0.0–0.68 and 0.05–0.66, respectively. This range of values allows for very efficient redistribution in the atmosphere of HD 209458b, although values of f_r>0.5 seem physically unrealistic.

We have constrained the P–T profile and abundances on the planetary dayside. But which part of the atmosphere is really being sampled by the observations? Here, we will discuss the contribution to the net emergent flux from the surface elements at the top of the atmosphere (see Knutson et al. 2009a for a discussion of radial contribution functions). We begin with the simple picture of radiative equilibrium, with no atmospheric circulation, in a homogeneous and isotropically radiating atmosphere. In this case, the amount of flux from a given surface element is proportional to the amount of incident flux from the star, where the incident stellar flux drops off as cos θ; θ is the angle away from the substellar point. In this picture, the planetary infrared flux is greatest at the substellar point, where the planet is being heated the most. The area contribution to the net planet flux, however, increases as sin θ. The bigger annuli correspond to larger areas on the projected planetary disk as viewed by the observer. The lower flux emitting regions (the outer regions of the disk, with lower incident stellar intensity), therefore, have more area than the higher flux regions. There is a sweet spot at 45 deg where the contribution to the total observed flux peaks, with a distribution of sin(2θ). In other words, the total flux is an average over a wide region on the projected planetary disk.

A more realistic scenario is one not limited to radiative equilibrium; indeed atmospheric circulation models show dramatic differences of emergent flux across the planetary disk caused by strong hydrodynamic flows (see, e.g., Showman et al. 2009; Langton & Laughlin 2008; Cho et al. 2003; Dobbs-Dixon & Lin 2008). Without knowing the atmospheric circulation and the resulting surface pattern, can we say anything about which part of the dayside planetary disk our P–T profiles and abundances preferentially sample? Given a retrieved abundance for a molecule, we do not know if this is from a localized concentration or from a global abundance, or something in between. For instance, we have found evidence for a high CO₂ abundance in HD 189733 (see also Swain et al. 2009a). But is this just a local effect of transported photochemical products? Or is the whole dayside evenly enriched in CO₂? If we are not measuring a true average, we have to be careful about how we interpret the retrieved CO₂ abundance.

A thorough examination of the contribution of different parts of the planetary disk has to take into account the emergent surface flux profiles that are more realistic than the assumption of homogeneity and radiative equilibrium. For example, 3D models of hot Jupiter atmospheres by Showman et al. (2009) show a much slower than cos θ drop off from the substellar point caused by hydrodynamic flows. And at the planetary photosphere the "substellar hot spot" is blown downwind, to 30 deg; a detail that would not be apparent in the assumption of an isotropic and homogenous atmosphere. Moreover, the Spitzer 8 micron measurements by Knutson et al. (2007) show a hot spot at lower separation from the substellar point (and also show a cold spot on the same hemisphere). For instance, if CO₂ were predominantly formed in localized regions and blown downwind, it could land at geometrically favored regions in terms of the 1D averaged spectrum. Alternatively, vertically localized concentrations of CO₂ at favorable pressure levels could give rise to strong spectral signatures, similar to those formed by having the same concentration of CO₂ over the entire atmosphere. In these cases, a high value of CO₂ may not be representative of the entire dayside.

In the current work, we have described most of our results at the ξ² = 2 level. Our formulation for ξ², as described in Section 3.5, would be equivalent to the conventional definition of reduced χ², if N_obs were to be replaced by the number of degrees of freedom. Normally, one aims for fits at the reduced χ² = 1 level (we have also provided results for fits at the level of ξ² = 1). A χ², however, is devoid of meaning when the number of degrees of freedom is less than zero—in that we cannot relate χ² to a confidence level. In the absence of confidence levels, we are considering ξ² as a weighted mean squared error. In this framework, a ξ² = 1 means that, on average, the model deviates from the mean values of the data by 1σ, and ξ² = 2 means that the deviation is $\sqrt{2}$ = 1.4σ. In this work, the number of data points in the Spitzer broadband photometry is less than the number of free model parameters, meaning less than zero degrees of freedom. In both the Spitzer IRS and the HST/NICMOS data sets, there are, in principle, more data points than the number of free model parameters. But both these data sets sample a narrow wavelength range, with a limited number of absorption features. In addition, the degeneracies in our model parameter space are poorly understood, if at all. Without a robust parameter-fitting algorithm, we cannot assign confidence levels even with the HST and IRS data sets. Which error surface should be used to encompass the best-fit models, given that we have no direct correspondence to confidence levels? There is no conclusive answer and, ideally, one might like to consider only those models which fit the data to within ξ² = 1, or even ξ² = 0 indicating perfect fits. However, we have used a conservative choice of the ξ² = 2 surface for interpretation of our results, allowing for models fitting the observations to within ∼1.4σ of each data point, on an average.

One approach to accurately consider the uncertainty on the data is as follows. One could take tens of thousands of realizations of each data set, uniformly sampled over the error distribution. Fitting each realized data set would mean running millions of models for each case, leading up to billions of models per data set per planet. Such an approach is computationally prohibitive at the present time.

There is evidence that the dayside atmosphere of HD 189733b is variable. One highly relevant example to our interpretation is the 8 μm Spitzer IRAC measurement of HD 189733b reported by Knutson et al. (2007) versus that reported by Charbonneau et al. (2008) in the same channel at a different time. The flux ratios differ by ∼2.3σ. Also significant is the disagreement in the planet-star contrast levels between the HD 189733b IRS spectrum and the IRAC broadband photometry at 8 μm, but not at 5.8 μm (Grillmair et al. 2008).

Atmospheric variability is also suggested by the very strong CO₂ feature in the HD 189733 HST/NICMOS data set and the lack of a noticeably strong absorption feature in the 4.5 μm and 16 μm Spitzer channels. It is well known that the CO₂ molecule has an extremely strong absorption feature at 15 μm; yet there is no indication of strong CO₂ absorption observed in the 16 μm channel. Whether this data variability is due to intrinsic variability of the planet atmosphere or instead is due to the different data reduction techniques is under study (Desert et al. 2009; Beaulieu et al. 2009).

Turning to our model results, we have been unable to find fits to all three data sets of HD 189733b at the ξ² ⩽ 2 level, in the large range of parameter space we have explored. There may be befitting P–T profiles and compositions that our grid did not cover, but we think this is unlikely, mainly because the data sets show noticeable differences even without running the models, e.g., the observed features of CO₂. We plan to implement a more exhaustive exploration of the P–T parameter space by extending our work with a new optimization technique in the future.

Our approach allowed us to run a large ensemble of models in the parameter space. Nevertheless, future calculations are needed to estimate the covariances between the model parameters. Several model parameters are degenerate; for example, molecular abundances of species that share absorption features in the same photometric channel (e.g., CO and CO₂). Also well-known is the degeneracy between the P–T profile and the molecular abundances. Understanding the covariances between the parameters might help break the degeneracies in the future.

We chose the ∼10⁴P–T profiles based on physical arguments and preliminary fits to the data. Even though having a grid of ∼10⁴P–T profiles in the six-parameter P–T space, and ∼10⁷ models overall, means that our grid is coarse, we believe it covers most physically plausible profiles. Nonetheless, a more complete exploration requires a more sophisticated approach, and may find P–T profiles and the corresponding abundance constraints that may be slightly different from those reported here.

Our model does not include clouds or scattering. Although cloudless models are justified for hot Jupiters (see Section 3.2), clouds are needed for extending our model to cooler planetary atmospheres, such as those of habitable zone super Earths and cooler giant planets. Our model also does not include photochemistry, but instead varies over different molecular abundances. Once a set of best-fit models are found by our approach, the molecular abundances can be used to guide photochemical models; conversely, photochemical models can be used to check the plausibility of the best-fit molecular abundances. And, although we have used the latest available high-temperature opacities, further improvements are sought, particularly for CH₄ and CO₂.

Our model presented in this work is suitable only for interpreting the infrared spectra of planetary atmospheres; it was motivated by existing observations. We do not include the sources of opacity and scattering that are required to model the visible part of the spectra. However, our model can be extended to visible wavelengths, and we aim to discuss it in a future work. Constraints placed by visible observations will help to further narrow down the model parameter space already constrained by IR observations.

We have not included the intrinsic flux from the planet interior, because, for hot Jupiters stellar irradiation is the dominant energy source. However, extensions of our model to brown dwarf atmospheres or young, or massive, Jupiters will require inclusion of the interior energy. And, although irradiation is naturally included in the form of the parametric P–T profile, our model gives no information about the nature of the absorber. Nevertheless, given the best-fit P–T profile from our approach, one could quantify the range of parameters of a potential absorber.

We have presented a radical departure from the standard exoplanet atmosphere models. Our temperature and abundance retrieval method does not use the assumption of radiative + convective equilibrium in each layer to determine the P–T profile. By searching over millions of models, we are extracting the P–T profile from the data.

We have not necessarily suggested to replace the so-called self-consistent forward models, those models that include radiative transfer, hydrostatic equilibrium, radiative + convective equilibrium, day–night energy redistribution, non-equilibrium chemistry, and cloud formation. Instead of a stand-alone model, our method can be used as a starting point for narrowing the potential parameter space to run forward models in. We anticipate that this temperature retrieval method will enable modelers to run enough forward models to find quantitative constraints from the data, instead of running only a few representative models.

Atmospheric temperature and abundance retrieval is not new. Solar system planet atmosphere studies use temperature retrieval methods, but in the present context there is one major difference. Exquisite data for solar system planet atmospheres means that a fiducial P–T profile can be derived. The temperature retrieval process for a solar system planet atmosphere therefore involves perturbing the fiducial profile. For exoplanets, there is no starting point to derive a fiducial model because of the limited data. Building on methods used for solar system planets, Swain et al. (2008b, 2009b) have used a temperature retrieval method using published P–T profiles from other models, and perturbing the profiles to find fits to data (also see Sing et al. 2008). In this work, we report a new approach for atmospheric retrieval of exoplanets which allows us to explore a wide region of the parameter space, and to possibly motivate a new generation of retrieval methods tailored to exoplanet atmospheres.