${\mathcal{T}}$-REx. II. RETRIEVAL OF EMISSION SPECTRA

${\mathcal{T}}$-REx is a novel, Bayesian atmospheric spectroscopy retrieval suite designed for extrasolar planets. In Waldmann et al. (2015, hereafter W15), we introduce the overall architecture, data handling, minimization/sampling routines, handling of molecular line lists, Bayesian model selection and the transmission spectroscopy forward model. In this publication, we present the emission spectroscopy forward model as well as the atmospheric TP profile retrieval implemented.

Figure 1 shows a schematic diagram of the ${\mathcal{T}}$-REx architecture for the emission retrieval. The program can be subdivided into five sections.

• 1.  
  Inputs: program inputs such as parameter files, molecular line-lists (ExoMol Tennyson & Yurchenko 2012), (HITEMP Rothman et al. 2010), stellar spectra (PHOENIX, Allard et al. 2012), and spectroscopic observations to be analyzed.
• 2.  
  Model and Data handling: the Central Data Module manages all calls to the TP-Profile and Forward Model Modules and provides a standardized, common interface to different sampling routines implemented.
• 3.  
  Sampling: ${\mathcal{T}}$-REx features three independent sampling routines: Maximum Likelihood estimation (MLE, based on the LM-BFGS minimization in W15), Markov Chain Monte Carlo (MCMC), and Nested Sampling (NS). We refer the reader to W15 for implementation details.
• 4.  
  Post analysis and refinement: This stage differs from W15 by providing two iteration stages for the determination of the TP-Profile, see Section 3.
• 5.  
  Output: the final posterior distributions are analyzed and the final model spectrum, TP-Profile, and mixing ratios are returned.

Zoom In Zoom Out Reset image size

Analytical global-average TP-profiles exist in various flavors ranging from two-stream purely radiative to radiative-convective approximations and global circulation models (GCMs). We refer the reader to the relevant literature for a more in-depth discussion on the various analytic approaches (e.g., Liou 2002; Hubeny et al. 2003; Burrows et al. 2008; Hansen 2008; Madhusudhan & Seager 2009; Showman et al. 2009; Guillot 2010; Pierrehumbert 2010; Heng et al. 2012, 2014; Robinson & Catling 2012). For the Stage 1 approach, ${\mathcal{T}}$-REx features a purely radiative two-stream model as well as a choice of simpler TP-profiles based on purely geometric considerations (see 3.1.1). As described in Section 3.2, the exact form of the TP-profile is less relevant in our case because  the results will be refined in a second stage fitting.

Following Guillot (2010), the mean global temperature profile for a simple radiative downstream–upstream approximation can be expressed as

$$\begin{eqnarray}&&{T}^{4}(\tau )=\displaystyle \frac{3\;{T}_{\mathrm{int}}^{4}}{4}\;\left(\displaystyle \frac{2}{3}+\tau \right)+\displaystyle \frac{3\;{T}_{\mathrm{irr}}^{4}}{4}{\xi }_{{\gamma }_{1}}(\tau ),\end{eqnarray} \tag{ 8 }$$

where T_int is the planet internal heat flux, T_irr is the solar flux at the top of the atmosphere, and

$$\begin{eqnarray}{\xi }_{{\gamma }_{i}} & = & \displaystyle \frac{2}{3}+\displaystyle \frac{2}{3{\gamma }_{i}}\left[1+\left(\displaystyle \frac{{\gamma }_{i}\tau }{2}-1\right){e}^{-{\gamma }_{i}\tau }\right]\\ & & +\displaystyle \frac{2{\gamma }_{i}}{3}\left(1-\displaystyle \frac{{\tau }^{2}}{2}\right){E}_{2}({\gamma }_{i}\tau ),\end{eqnarray} \tag{ 9 }$$

where ${\gamma }_{1}={\kappa }_{\nu }/{\kappa }_{\mathrm{IR}}$ is the ratio of mean opacities in the optical (${\kappa }_{\nu }$) and infra-red (${\kappa }_{\mathrm{IR}}$) and E₂ is the second-order exponential integral given by

$$\begin{eqnarray}&&{E}_{n+1}(z)=\displaystyle \frac{1}{n}[{e}^{-z}-{{zE}}_{n}(z)].\end{eqnarray} \tag{ 10 }$$

We note that similar parameterizations exist in the literature (e.g., Equation (18), Robinson & Catling 2012; Pierrehumbert 2010). We also include the variation by Line et al. (2012) and Parmentier & Guillot (2014) including two optical opacity sources ${\kappa }_{{\nu }_{1}}$ and ${\kappa }_{{\nu }_{2}}$ and a weighting factor between optical opacities (left as a free parameter) α,

$$\begin{eqnarray}{T}^{4}(\tau ) & = & \displaystyle \frac{3\;{T}_{\mathrm{int}}^{4}}{4}\left(\displaystyle \frac{2}{3}+\tau \right)+\displaystyle \frac{3\;{T}_{\mathrm{irr}}^{4}}{4}(1-\alpha ){\xi }_{{\gamma }_{1}}(\tau )\\ & & +\displaystyle \frac{3\;{T}_{\mathrm{irr}}^{4}}{4}\alpha {\xi }_{{\gamma }_{2}}(\tau ).\end{eqnarray} \tag{ 11 }$$

The temperature as a function of opacity τ can be mapped to a pressure grid by assuming the following relation

$$\begin{eqnarray}&&\tau =\displaystyle \frac{{\kappa }_{\mathrm{IR}}P}{g}.\end{eqnarray} \tag{ 12 }$$

3.1.1. Other TP-profiles

In ${\mathcal{T}}$-REx, we have worked toward a hybrid solution of the above mentioned approaches in which the final solution is not bound by a parametric model, but does not inherit the potentially poor convergence properties (and susceptibility to noise) of a layer-by-layer approach. Here we proceed in two stages: (1) we perform a parametric model retrieval and (2) we take the retrieval solution to "guide" a second layer-by-layer retrieval, relaxing the parametric model constraint and thereby fine-tuning the TP-profile.

3.2.1. Stage 1

In Stage 1, we adopt a classical parametric model retrieval using one of the TP-profile parameterizations implemented in ${\mathcal{T}}$-REx. The idea is to constrain an initial best-fit of the model and while a realistic model (i.e., one well suited to describe the underlying TP-profile) is preferable, it is not a hard requirement. We now fit the solution using either of ${\mathcal{T}}$-REx's sampling routines (MLE, MCMC, NS). The error on the sampled parametric model parameters is converted to one σ lower and upper temperature bounds for a layer-by-layer model. This is done using a numerical model permutation analysis of the Stage 1 parameters.

We then calculate the following distance matrix

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{ij}}^{2}=| {\hat{T}}_{i}-{\hat{T}}_{j}{| }^{2}+{({\sigma }_{i}+{\sigma }_{j})}^{2},\end{eqnarray} \tag{ 13 }$$

where $\hat{T}$ is the maximum likelihood temperature estimator of the parametric model fit for the ith and jth atmospheric layer and σ is the respective error. Equation (13) captures the layer to layer temperature variations in the TP-profile and is hence conceptually similar to the Jacobian of the profile. We normalize Δ_ij in terms of the minimal and maximal temperature variations found in the TP-profile,

$$\begin{eqnarray}&&{{\mathcal{C}}}_{{ij}}=1.0-\displaystyle \frac{{{\rm{\Delta }}}_{{ij}}-\mathrm{argmin}\;{\rm{\Delta }}}{\mathrm{argmax}\;{\rm{\Delta }}}\end{eqnarray} \tag{ 14 }$$

which can be understood as a positively defined temperature correlation matrix with layers most similar in temperature featuring the highest correlation and layers most dissimilar resulting in a very low correlation. An example of ${\mathcal{C}}$ for a given TP-profile can be found in Figure 2.

Zoom In Zoom Out Reset image size

3.2.2. Stage 2

In the second stage of our TP-profile retrieval, we relax the original solution found in Stage 1 and, by that, "fine-tune" the TP-profile. We construct a second correlation matrix imposing an exponential correlation length across pressure levels (Rodgers 1976, 2000)

$$\begin{eqnarray}&&{S}_{{ij}}={({S}_{{ii}}{S}_{{jj}})}^{1/2}\mathrm{exp}\left(-\displaystyle \frac{| \mathrm{ln}({P}_{i}/{P}_{j})| }{c}\right),\end{eqnarray} \tag{ 15 }$$

where c is the correlation length in terms of atmospheric scale heights. We set c = 3.0 by default, but it can otherwise be user defined. Rodgers (2000) gives an advisable range between 1.0 and 8.0, with Irwin et al. (2008) using a default of c = 1.5 and Line et al. (2013) c = 7.0. Larger values of c correspond to a stronger regularization of the TP-profile (i.e., smoothing). A stronger regularization can be useful in poorly constrained data sets (either low S/N, sparsely sampled, or both). Equation (15) is identical to that used in the layer-by-layer approach by Irwin et al. (2008). We now construct a hybrid correlation matrix by combining Equations (14) and (15)

$$\begin{eqnarray}&&{{\mathcal{D}}}_{{ij}}(\alpha )=\alpha {{\mathcal{C}}}_{{ij}}+(1-\alpha ){S}_{{ij}},\end{eqnarray} \tag{ 16 }$$

where α is a scaling factor ranging from zero to one. Figure 3 shows instances of ${\mathcal{D}}(\alpha )$ at varying values of α. We now correlate our layer-by-layer TP-profile with ${\mathcal{D}}(\alpha )$ whereby leaving α as free parameter to be fitted. By allowing α to vary, we can dynamically relax the parametric model solution from Stage 1, $\alpha =1.0$, to an unconstrained solution, $\alpha =0.0$. The advantage of such an approach is twofold: (1) since the Stage 1 fitting should have already achieved a solution close to the real maximum likelihood, convergence in the second stage toward the volume of lowest χ² is significantly faster and (2) ${\mathcal{T}}$-REx can freely decide to relax the Stage 1 solution should this be favored by the data. Practically, this happens frequently at later stages of the fitting. Whereas ${\mathcal{T}}$-REx initially converges toward the Stage 1 solution (i.e., $\alpha \to 1$), at later stages of the fitting the code begins to reject the parametric model (i.e., $\alpha \to 0$) as it "fine-tunes" the original solution. This "tuning" can achieve significantly higher maximum likelihoods than the Stage 1 fitting alone.

Zoom In Zoom Out Reset image size

In the approach explained in Section 3.2, we keep an even sampling of atmospheric layers in log(P) for Stage 2. For well sampled and high S/N data, this approach is adequate. However, for coarsely sampled and/or poor S/N data, it is often advisable to reduce the number of free parameters to a minimum to aid convergence. In these cases, we can utilize the sparsity of the Stage 1 TP-profile solution to devise a nonlinear sampling of the exoplanetary atmosphere. We base our compression algorithm on the fact that only changes in the TP gradient need to be modeled, i.e., an isothermal TP-profile can be perfectly described by a single temperature parameter whereas areas of non-isothermal structure need more parameters to capture variation. The compression algorithm uses the Stage 1 correlation matrix ${\mathcal{C}}$ to only retain TP-profile layers corresponding to a $\gt 2\%$ change in the TP gradient with respect to the previously retained layer. We graphically depict this in Figure 4. In addition to TP-profile layers sampled this way, we also include an extra sampling layer whenever no change in thermal gradient was detected for $\gt 10$ layers. The majority of the TP-profile variations should be captured by the Stage 1 retrieval and Stage 2 is "fine-tuning" this solution. The inclusion of a coarse sampling in isothermal regimes does allow the Stage 2 retrieval to deviate from any Stage 1 solution, should this be supported by a higher global likelihood. Using the nonlinear sampling, we can reduce a 100 layer atmospheric model to typically 15–25 free parameters.

Zoom In Zoom Out Reset image size

WASP-76b (West et al. 2013) is a very hot-Jupiter orbiting a late F7 star. It is highly inflated at ${1.83}_{-0.04}^{+0.06}$ R_Jup, $0.92\pm$0.03 M_Jup, and ${T}_{\mathrm{equ}}\sim 2200$ K. We take its bulk and orbital properties and generate a simulated observation at a resolution of 300, a spectral range of 0.5–20 μm, and 100 ppm errors. We set the atmospheric composition to $1\times {10}^{-4}$ H₂O, $1\times {10}^{-5}$ CO, and $1\times {10}^{-5}$ CO₂. The input TP-profile is shown in Figure 7 (black line). It is important to note that here (as well as in Section 4.2) the input TP-profile was generated using a script external to ${\mathcal{T}}$-REx with no relation to either Stage 1 parametrization. This provides an adequate test-bed for the Stage 1 fitting to accurately retrieve an arbitrary atmospheric profile.

Zoom In Zoom Out Reset image size

As described in Section 3.2, we retrieve the TP-profile and abundances in two stages. For computational efficiency reasons, here we only compute the NS solutions (which are also the most accurate, see W15). Tests were run with both maximum-likelihood retrievals and MCMC retrievals and solutions between all sets of solutions are in good agreement.

The NS Stage 1 solution returns four local likelihood maxima of which the global maximum was selected. Figures 6 (top), 7 (left), and Table 1 show the retrieved spectrum, TP-profile, and retrieved abundances respectively. The Stage 1 TP-profile mostly captures the input TP-profile but shows unrealistic bumps and wiggles as well as unrealistic distributions of the one sigma error bar. These are artefacts of the parameterization in Section 3.1. The three local maxima TP-profiles shown in Figure 7 underestimate the bulk temperature of the planet, driving the retrieval to assume unrealistically low abundances of molecular trace gases. In this example, this is found to be a numerical effect that disappears by increasing the sampling grid density of the NS. Nonetheless, the potential degeneracy between TP-profile and molecular abundances in atmospheric retrievals is worth noting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As described earlier, we computed the TP-profile covariance, ${\mathcal{C}}$ (Figure 5), and tuned the Stage 1 retrieval by relaxing the parametric solution. Figures 6 (bottom), 7 (right), and Table 1 show the Stage 2 retrieval results. Inspecting Figure 6, both Stage 1 and Stage 2 retrievals fit the input spectrum sufficiently well but Stage 2 comes significantly closer to the "true" TP-profile and trace gas abundances. Figure 8 shows the normalized emission contribution functions of both retrievals and their difference. This shows the planetary emission mainly emanating in the non-isothermal regions (up to the tropopause) of the TP-profile, as expected. However, Stage 2 emissions originate significantly lower (by nearly an order of magnitude) in the atmosphere (blue band in the bottom panel).

Zoom In Zoom Out Reset image size

We simulated an emission spectrum of the hot-super-Earth 55 Cnc e (Fischer et al. 2008). We use trace gas compositions of $1\times {10}^{-4}$ H₂O, $1\times {10}^{-5}$ CO, and $1\times {10}^{-5}$ CO₂ and a sharp TP-profile shown in Figure 10. Because the previous WASP-76 example is currently unrealistically optimistic, given the combination of high S/N, moderate resolution (R ∼ 300) and a very broad wavelength coverage, we opt for a more realistic example here. We calculated the expected resolution and S/N for 100 co-added eclipses obtained by a one-meter-class transiting spectroscopy space-mission (e.g., the ESA M4 proposal ARIEL) over a wavelength range spanning 2–8 μm. Using EChOSim, an end-to-end simulator for transit spectroscopy space-missions (Pascale et al. 2014; Waldmann & Pascale 2014) developed for the ESA M3 EChO proposal (Tinetti et al. 2012), we calculated realistic error bars for this hot super-Earth, shown in Figure 9.

Zoom In Zoom Out Reset image size

Similarly to Section 4.1, Figures 9–11 show the Stage 1 and 2 retrieved spectra, TP-profiles, and Stage 1 TP-profile covariance respectively. Table 2 shows that Stage 2 retrieval converges at a significantly higher global Evidence and presents an improvement in the accuracy of abundances retrieved as well as the TP-profile retrieved. Figure 12 shows the absolute difference between the Stage 1 and Stage 2 model fits at a spectral resolution of 1000. Here the discrepancies between both spectral fits are of the order of $2\times {10}^{-5}$ or less. This is not very significant in terms of a naive χ² fit to relatively low S/N data, but does significantly drive the retrieval of the TP-profile and trace gas abundances. This is reflected in the Bayesian Evidence measured for Stage 1 and Stage 2 models. Figure 13 shows the contribution functions for both retrievals. Here we have positive and negative temperature differences (see bottom left plot), resulting in a more complex contribution differential (bottom left plot) than for WASP-76b. It should be noted that this contribution differential is highly wavelength dependent and illustrates the sensitivity of varying wavelength ranges on the TP-profile retrievability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we test ${\mathcal{T}}$-REx on the emission spectrum of the hot-Jupiter HD 189733b. Its emission spectrum is among the most complete and best studied (e.g., Madhusudhan & Seager 2009; Swain et al. 2009; Lee et al. 2011; Line et al. 2012, 2014). We have compiled data sets ranging from the NIR (1.4 μm) to the MIR (24μm), namely: Hubble Space Telescope/NICMOS (Swain et al. 2009), Spitzer/IRAC 3.6, 4.5 μm (Knutson et al. 2012), 8.0 μm (Agol et al. 2010), Spitzer/IRS spectroscopy (Grillmair et al. 2008), Spitzer/IRS 16 and 24 μm photometry (Charbonneau et al. 2008). As in previous studies, we consider the active trace gasses H₂O, CH₄, CO, and CO₂ and otherwise assume a hydrogen dominated atmosphere.

Figure 14 shows the retrieved emission models with retrieved abundances and log(Evidence) reported in Table 3. Figures 15 and 16 report the retrieved Stage 1 and Stage 2 TP-profiles and associated Stage 1 TP-covariance, respectively. Figures 17 and 18 are the marginalized and conditional posteriors for the active trace gasses considered. Figure 19 illustrates the difference between Stage 1 and Stage 2 emission spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As described above, the Stage 1 retrieval was first performed using the parmetric TP-profile by Guillot (2010) with the additional optical opacity terms proposed by Line et al. (2012). Here, all parameters of the TP model were allowed to vary. We modeled the atmosphere over 25 scale heights sampled onto a 150 pressure-layer grid. Following this initial retrieval, the TP-covariance was computed resulting in a predominantly isothermal atmosphere down to the lowest ∼20 layers (∼0.1 bar pressure). The obtained TP-profile is well constrained but shows systematics (e.g., at 5×10⁻⁴ bar) typical of this type of parametrization. The posterior distributions of the active trace gases, Figure 17 shows good constrains on abundances. From Figure 17, we can see the retrieved values for methane to only constitute an upper limit.

The Stage 2 hybrid model contained 22 free parameters on its nonlinearly sampled TP-profile grid. All parameters were allowed to vary freely between fully constrained ($\alpha =1.0$ in Equation (16)) and unconstrained scenarios. Figure 15 (right) shows the Stage 2 TP-profile with $\alpha =0.54$, indicating a significant shift away from the parametric solution of ${\text{}}{Stage}\;1$. The Stage 2 model features a lower tropopause pressure along with reduced water and carbon-dioxide abundances to achieve a very significant increase in Bayesian Evidence, see Table 3 and the discussion in Section 5. The Stage 2 TP-profile error bounds are slightly larger. Given a significantly higher Evidence for Stage 2, it is indicative of Stage 1 TP-profile errors to be underestimated by the parameterization used. The posteriors distributions of Stage 2 now show a convergence of methane beyond the upper-bound constrained of Stage 1.

Compared to Table 3 in Line et al. (2014), we obtain lower abundances of CH₄ and CO_2, but higher values for CO. Whereas these values differ from previous analyses, we find them significantly closer to the chemical model predictions of HD 189733b (e.g., Line et al. 2010; Venot et al. 2012; Moses et al. 2013). Note also that Stage 1 results exhibit the same trend of lower CH₄ and CO₂ abundances. We find two aspects by which our analysis differs from others in the literature: (1) a finer and complete sampling of correlated likelihoods, in particular, compared to maximum likelihood methods that are often less efficient with sparsely sampled data; (2) the completeness of molecular line-lists used. This is particularly true for the C/O ratio determination using CH₄ and CO/CO₂ as tracers. In this work, we use the new YT10to10 CH₄ line-list (Yurchenko & Tennyson 2014).