Retrieving Exoplanet Atmospheres Using Planetary Infrared Excess: Prospects for the Night Side of WASP-43 b and Other Hot Jupiters

We developed a modified version of the CaltecH Inverse ModEling and Retrieval Algorithms (CHIMERA) code (Line et al. 2013b, 2014b) to operate using the PIE technique. CHIMERA is a flexible, open-source, and well-established retrieval framework used to infer the atmospheric composition of exoplanets from emission and transmission spectroscopy measurements in the optical to mid-infrared for a wide variety of planet types (terrestrial to ultra-hot Jupiter) given arbitrary molecular abundances (constant with altitude or variable with altitude), thermal structures, and cloud/aerosol properties (droplet sizes, single scatter albedo/asymmetry parameter, vertical distribution; Mai & Line 2019). The core radiative transfer scheme accounts for full multiple scattering of outgoing thermal radiation utilizing the methods outlined by Toon et al. (1989). We incorporate correlated-K ("resort-rebin"; Amundsen et al. 2017) opacities by generating pre-computed, high-resolution cross section grids sourced from a variety of databases including ExoMol (Tennyson et al. 2016; Tennyson & Yurchenko 2018), utilizing the EXOCROSS routine (Yurchenko et al. 2018), and HITRAN/HITEMP (Rothman et al. 2010), via the HAPI routine (Kochanov et al. 2016). The transmission and emission routines have been validated on numerous occasions against other codes (Line et al. 2013a, 2013b; Morley et al. 2015) and used to analyze brown-dwarf observations (Line et al. 2014a, 2015, 2017).

The most significant change required to implement PIE into the standard CHIMERA emission spectrum forward model is to use the absolute flux from the system, F_sys. This differs from the standard emission spectroscopy approach, which models the wavelength-dependent secondary eclipse depth as the planet-to-star flux ratio, ΔF = F_p /F_s . In this case, we define the system flux as the sum of all major flux sources in the exoplanetary system:

$$\begin{eqnarray}&&{F}_{\mathrm{sys}}={F}_{s}+{F}_{p}+{F}_{\epsilon },\end{eqnarray} \tag{ 1 }$$

where F_s is the stellar flux, F_p is the planet flux, and F is a catch-all term to account for any residual astrophysical flux source in the field (such as an additional, cooler planet or exozodiacal dust). For simplicity, we set F = 0 and focus on joint modeling of the stellar and planetary flux. We explore the ramifications of F = F_ez due to exozodiacal light in Section 3.4.

2.1.1. Planetary-flux Model

We use the standard CHIMERA emission-spectrum forward model presented by Line et al. (2013b, 2014b) and Batalha et al. (2018) with no modifications to the physics and chemistry. We employ the free chemistry (henceforth "free") forward model in our exploration of the PIE technique, which allows for the abundances of the gases in the state vector to take on any (physical, even if chemically implausible) values and assumes evenly mixed volume-mixing ratio vertical profiles. As such, we directly fit for the (log) gas abundances. Although we performed tests using the "chemically consistent" CHIMERA forward model (e.g., Kreidberg et al. 2015), which assumes thermochemical-equilibrium gas-vertical profiles, our results were generally consistent with those acquired with the "free" model and have been omitted for clarity and brevity.

For the temperature-pressure (TP) profile, we use the three-parameter, analytic radiative equilibrium model from Guillot (2010; see also Parmentier & Guillot 2014; Line et al. 2013b; Mai & Line 2019). The irradiation temperature T_irr, the infrared opacity $\mathrm{log}({\kappa }_{\mathrm{IR}})$, and the single channel visible to IR opacity ratio $\mathrm{log}({g}_{1})$ are treated as free parameters, and the internal temperature T_int is fixed at 200 K.

Bringing together all of the physical, atmospheric, and thermal parameters that define our planet model, the planet flux is given as,

$$\begin{eqnarray}&&{F}_{p}={F}_{p}({\theta }_{p})\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\theta }_{p}=\{{R}_{p},{T}_{\mathrm{irr}},\mathrm{log}({\kappa }_{\mathrm{IR}}),\mathrm{log}({g}_{1}),\mathrm{log}({{\rm{H}}}_{2}{\rm{O}})\\ \times \mathrm{log}({\mathrm{CH}}_{4}),\mathrm{log}(\mathrm{CO}),\mathrm{log}({\mathrm{CO}}_{2}),\mathrm{log}({\mathrm{NH}}_{3})\}.\end{array}\end{eqnarray} \tag{ 3 }$$

2.1.2. Stellar-flux Model

We use Phoenix stellar models as the basis for our stellar-flux model (Allard et al. 2003, 2007, 2012). We access the stellar grids using the pysynphot code (STScI Development Team 2013) and linearly interpolate in effective temperature T_eff, surface gravity $\mathrm{log}(g)$, and metallicity $\mathrm{log}([M/H])$ between computed stellar models. We scale the stellar flux using the stellar distance, d, and stellar radius, R_s , so as to reproduce stellar-flux observations at Earth. This results in a five-parameter stellar-flux model that is parameterized via the following state vector,

$$\begin{eqnarray}&&{F}_{s}={F}_{s}({\theta }_{s})\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\theta }_{s}=\{{T}_{\mathrm{eff}},\mathrm{log}(g),\mathrm{log}([M/H]),d,{R}_{s}\}.\end{eqnarray} \tag{ 5 }$$

2.1.3. Inverse Model

We use PandExo to simulate observations of WASP-43b with JWST (Batalha et al. 2017, 2018). PandExo uses the official JWST exposure-time calculator (ETC), Pandeia (Pontoppidan et al. 2016), to perform noise modeling for exoplanet observations. In this work, we focus exclusively on observations in the 0.6–12 μm wavelength range using NIRISS SOSS (Doyon et al. 2012), NIRSpec G395 (Ferruit et al. 2012), and MIRI LRS (Kendrew et al. 2015) because these instruments provide complete wavelength coverage over the stated range without saturating. We bin our simulated JWST data to a constant resolving power of ${ \mathcal R }=100$ for our retrieval analyses.

To study the efficacy of the PIE technique, we simulated high-resolution model spectra using the forward model described in Section 2.1 before implementing noise from PandExo. Furthermore, we do not add random Gaussian jitter to the data. Although these assumptions yield idealized observations that, by design, will be capable of being well fit by our nominal model, this exercise is common practice in the exoplanet retrieval literature (see Feng et al. 2018, for more details) and allows us to characterize bias and degeneracies in the retrieval procedure, which are of keen interest in this study. Although we calculate the precision of our synthetic spectra using only PandExo noise calculations, we considered the effect of absolute flux calibration accuracy in the Appendix as it may complicate the study of broad wavelength spectra constructed using multiple instruments with offsets in absolute flux between them.

The first two columns of Table 1 provide the default choice of model parameters that were used to simulate the JWST observations. We used Gaia DR2 data for the stellar parameters (Gaia Collaboration et al. 2018; Bailer-Jones et al. 2018). For the planetary parameters, we used gas abundances and a TP profile that are consistent with the current estimates from retrievals of WASP-43b's transmission spectrum (Kreidberg et al. 2014). However, one notable departure from this convention is that we set the irradiation temperature (T_irr) to 1000 K to be representative of the night side of this hot Jupiter (Beatty et al. 2019). The third column of Table 1 shows the Bayesian prior probability distributions used throughout this investigation, unless otherwise stated.

Table 1. Retrieved Constraints on WASP-43b using Multiple JWST Instruments to Simulate the PIE Technique

			NIRISS SOSS	NIRSpec G395	MIRI LRS	All	All
Parameters	Defaults	Priors	(1 hr)	(1 hr)	(1 hr)	(1 hr)	(4 hr)
Teff	4305.0	${ \mathcal N }(4305,174.75)$	±0.020	±0.065	±0.16	±0.018	±0.0086
$\mathrm{log}([M/H])$	−0.05	${ \mathcal N }(-0.05,0.17)$	±0.000090	±0.00080	±0.0018	±0.000069	±0.000034
$\mathrm{log}(g)$	4.646	${ \mathcal N }(4.646,0.052)$	±0.000054	±0.00031	±0.00069	±0.000050	±0.000025
Rs	0.6629	${ \mathcal N }(0.6629,0.0554)$	±0.0025	±0.0024	±0.0025	±0.0025	±0.0025
d	86.7467	${ \mathcal N }(86.75,0.33)$	±0.33	±0.32	±0.32	±0.32	±0.33
Rp	0.93	${ \mathcal N }(0.93,0.08)$	±0.080	±0.079	±0.068	±0.026	±0.014
Tirr	1000	${ \mathcal U }(300,3000)$	±100	±105	±303	±41	±20
$\mathrm{log}({\kappa }_{\mathrm{IR}})$	−1.5	${ \mathcal U }(-3,0)$	±0.89	±0.44	±0.78	±0.18	±0.090
$\mathrm{log}({g}_{1})$	−0.7	${ \mathcal U }(-3,-1)$	±1.0	±0.83	±1.1	±0.18	±0.044
H2O	−3.0	${ \mathcal U }(-15,0)$	±5.1	±1.7	±2.3	±0.39	±0.16
CH4	−8.0	${ \mathcal U }(-15,0)$	±4.9	±3.6	±5.1	±3.1	±3.0
CO	−7.0	${ \mathcal U }(-15,0)$	±5.0	±4.7	±5.1	±3.7	±3.3
CO2	−8.0	${ \mathcal U }(-15,0)$	±4.6	±3.5	±5.0	±2.9	±2.7
NH3	−10.0	${ \mathcal U }(-15,0)$	±5.1	±3.8	±5.3	±3.5	±3.2

Note. ${ \mathcal U }$ denotes a uniform prior probability distribution described in terms of the lower and upper bounds, and ${ \mathcal N }$ denotes a normal prior described by the mean and standard deviation.

Download table as:  ASCIITypeset image

Our initial retrieval used simulated JWST data for 4 hr out-of-transit observations with NIRISS SOSS, NIRSpec G395, and MIRI LRS to achieve continuous wavelength coverage between 0.6 and 12 μm. Now we examine the result of using individual instruments. Table 1 shows our retrieval results where each column lists 1σ constraints on each retrieved model parameter (rows) for multiple experiments using different JWST instrument modes (columns). Note that the final column on the right shows results for 4 hr with all three instruments and corresponds to the results shown in Figure 1 and Figure 2, while the other columns correspond to 1 hr exposures.

As expected, retrievals on the full wavelength range provide the best constraints on the stellar and planetary parameters, but the performance of individual instruments varies from parameter to parameter. For example, we find that shorter wavelength observations that resolve the peak of the stellar SED offer the best constraints on the stellar parameters, with the tightest constraints on the effective stellar temperature (T_eff) provided by NIRISS, then NIRSpec, then MIRI. Similarly, the planet's 1000 K night-side emission peaks at 2.9 μm; therefore, its irradiation temperature (T_irr) is best constrained by either NIRISS or NIRSpec, while MIRI offers a 3× worse constraint. However, the planet radius (R_p ) is essentially unconstrained by the individual use of NIRISS or NIRSpec (where the posteriors are consistent with the priors), but shows an improved posterior constraint using only MIRI data. Despite the lack of radius constraint from using NIRISS or NIRSpec on their own, the planet radius is significantly better constrained using data from all three instruments. This result is consistent with the presence of a correlation between planet radius and irradiation temperature shown in Figure 2, which we will revisit in Section 3.2 for non-transiting planets. Regarding constraints on H₂O, NIRSpec G395 provides the best individual constraint of ±1.7 dex, MIRI LRS provides the next best constraint of ±2.3 dex, and NIRISS SOSS is much worse at ±5.1 dex. These abundance results that favor the use of NIRSpec reflect a combination of the rising sensitivity to the planet thermal emission with wavelength and the specific wavelengths that exhibit H₂O absorption features.

To examine how marginalizing over uncertainties in the stellar spectrum decreases the precision on retrieved planetary parameters when using the PIE technique, we ran a series of retrievals with the stellar flux fixed at the true input value. This simulates the traditional secondary eclipse retrieval approach. We kept the exposure time per instrument fixed at 4 hr and simulated individual cases for NIRISS, NIRSpec, and MIRI. Using a fixed stellar flux resulted in tighter constraints on, in particular, the planet's irradiation temperature and radius, with the largest effect seen toward longer wavelengths. When marginalizing over the stellar fit using the PIE technique we observed a 1.0×, 2.4×, and 3.5× decrease in the precision on the retrieved irradiation temperature using NIRISS, NIRSpec, and MIRI, respectively, and a 1.1×, 2.3×, and 2.5× decrease in the precision on the retrieved planet radius. We observed no decrease in the precision of molecular abundances. Thus, long wavelength PIE observations in the MIR that lack the NIR reference wavelength range to constrain the stellar spectrum suffer the greatest loss in precision of retrieved planetary information relative to traditional secondary eclipse measurements. One potential caveat that must be considered when applying the PIE technique to absolute fluxes is that JWST's absolute flux calibration accuracy may limit the ability to use the PIE technique on panchromatic spectra stitched together from observations using different JWST instruments with non-negligible offsets caused by their absolute flux calibration. We investigated this in detail in the Appendix and found that if JWST achieves its nominal requirement of 2% absolute flux calibration accuracy, offset-correction parameters can be used within the PIE retrieval to enable the analysis of spectra obtained with different instruments. Using 2% offsets between instruments, we find that the retrieved molecular abundances in the planetary atmosphere are completely unaffected and a relatively small loss in precision (the 1σ uncertainties grow by <50%) is incurred for the stellar radius, planetary radius, and irradiation temperature relative to the baseline results using all three instruments without calibration offsets (shown in Table 1). Even with 2% instrument offsets to simulate absolute flux calibration, our results indicate that using all three JWST instruments—NIRISS SOSS, NIRSpec G395, and MIRI LRS to observe near-continuous wavelength coverage between 0.6 and 12 μm still outperforms any single instrument for night-side characterization of WASP-43b with the PIE technique.

By using spectral information to separate planets from stars, the PIE technique offers an opportunity to study the thermal emission from non-transiting planets. For the purposes of this work, we assume that non-transiting planets have known orbits but an unknown radius, like those detected from radial-velocity (RV) surveys. However, we remain focused on a WASP-43b-like exoplanet. Instead of using the Gaussian prior on planet radius shown in Table 1 from transit observations, we use an uninformative uniform prior given by ${ \mathcal U }(0.5,1.5)$.

Figure 3 shows 1D and 2D marginalized posterior distributions for WASP-43b night-side PIE retrievals assuming no informative prior information about the planet radius. The 1D marginal posteriors show that, if only one instrument is used, planet temperature can be better constrained with shorter wavelength observations and planet radius can be better constrained at longer wavelengths. Both the planet radius and temperature can be much better constrained by combining data from two or three instruments. Using only two instruments, the irradiation temperature is best constrained with combined data from NIRISS and NIRSpec, while the planet radius is best constrained with NIRISS and MIRI. Intriguingly, the 2D posterior projection between irradiation temperature and planet radius shows how their covariance changes with each instrument.

Zoom In Zoom Out Reset image size

Here we briefly consider the effect of varying the exposure time. We repeated our previous retrieval that used data from 1 hr exposures with all three JWST instrument modes, now using 0.25 hr and 4 hr exposures to halve and double the signal-to-noise ratio (SNR) of the data, respectively, compared to the fiducial case. Table 1 provides the constraints for the 1 hr and 4 hr cases. With the exception of prior-dominated parameters (e.g., the stellar radius and distance) and unconstrained parameters with non-Gaussian posteriors (e.g., CH₄, CO, etc.), our exposure-time experiment yielded results that scaled with SNR exactly as expected. For example, the 1σ constraint on the planet radius dropped (rose) from ±0.026 R_J to ±0.014 R_J (±0.051 R_J ), approximately a factor of 2, for a 2× increase (decrease) in data SNR. We note that these tests do not include the effect of a noise floor in the observations, which may cause wavelength-dependent deviations in the scaling of SNR with the square root of the exposure time that could propagate through the parameter inference. The effects of various noise floors on PIE retrievals will be explored in an upcoming paper. Our 4 hr NIRISS SOSS spectrum reaches a maximum precision of 17 ppm at 1.4 μm with ${ \mathcal R }=100$. While the exact wavelength coverage needed to precisely constrain the planet parameters will change with planet temperature, for the non-transiting exoplanet case we recommend maximizing the wavelength coverage instead of repeated observations with a specific instrument mode.

A key concern for the PIE technique is the presence of exozodiacal dust in exoplanetary systems and how the thermal emission from such dust may be incorrectly attributed to planetary thermal emission and could propagate biases into the retrieval of atmospheric parameters. To estimate the magnitude of this effect within the context of our previous findings, we conducted a series of retrieval simulations with exozodiacal emission injected into the observed spectrum. We varied the intensity of the exozodiacal signal injected in the simulated data, but initially we did not attempt to retrieve the exozodiacal signal; instead we seek to characterize the bias that manifests in the inference if exozodiacal dust is not accounted for by the retrieval model.

We modeled exozodiacal flux (emitted and reflected) in the WASP-43 system using Zodipic (Kuchner 2012). We assumed the following properties for WASP-43, R_* = 0.663 R_⊙, L_* = 0.136 L_⊙, T_* = 4305 K, d=86.7467 pc, $\mathrm{log}g$ = 4.492, and solar metallicity. By default, the dust in the disk is considered to scatter isotropically and scattered light is only considered when λ < 4.2 μm. The inner spherical dust radius corresponds to a dust temperature of 1500 K, which is appropriate for removal from the system via sublimation in the absence of a planet. In reality, a planet could change the morphology of the disk and cause dust to collect at various mean-motion resonances and clear out annuli (Stark et al. 2013), but we assume the simple case of an isolated and smooth disk out to 3.28 au. We modeled the spectral contribution from 0.3 μm to 30 μm assuming a pixel size of 1 mas.

For this exozodiacal sensitivity study, we considered only the case with the highest quality data, which used all three JWST instruments to span 0.6–12 μm and assumed a 4 hr exposure time for each instrument mode. We tested adding 1 zodi, consistent with a solar system level of zodiacal dust, 3 zodis, consistent with the median zodi level for Sun-like stars (Ertel et al. 2020 found the median, $m={3}_{-3}^{+6}$), and 10 zodis.

We find that neglecting exozodiacal emission in PIE retrievals can bias the retrieved planet radius and, to a lesser extent, the planet temperature for systems with high exozodiacal levels. The left panel of Figure 4 shows the retrieved 1D and 2D marginalized posteriors for the planet radius and irradiation temperature for a set of simulations with varying exozodiacal contamination in the observed spectrum. For the case without exozodiacal contamination (N_ez = 0), we retrieve posteriors that are centered at the input value used to generate the JWST spectra. Increasing the exozodiacal contamination in the simulated system biases the retrieved planet radius to higher values in order to account for the additional system flux. This radius bias corresponds to 0.00558 R_J (∼ 400 km) per zodi in the WASP-43 system. Similarly, the retrieved irradiation temperature shifts to lower values with increasing zodi, consistent with the anticorrelation seen between radius and temperature. We find that the bias in radius can be statistically significant with 1.1σ and 4.2σ radius increases for N_ez = 3 and N_ez = 10, respectively, but the corresponding bias in irradiation temperature is negligible (<1σ for N_ez = 10) for this JWST test case. We note that for our high exozodiacal case with N_ez = 10, the stellar gravity $\mathrm{log}(g)$, TP profile parameter $\mathrm{log}({g}_{1})$, and H₂O abundance are all marginally biased by about 1σ.

Zoom In Zoom Out Reset image size

As a final test of exozodiacal contamination, we incorporated a simple exozodiacal emission model into our PIE retrieval in order to explore the ability to mitigate the aforementioned biases. In this case, Equation (1) becomes

$$\begin{eqnarray}&&{F}_{\mathrm{sys}}={F}_{s}+{F}_{p}+{F}_{\mathrm{ez}},\end{eqnarray} \tag{ 6 }$$

where F_ez—the exozodiacal spectrum—replaces F . The exozodiacal light is assumed for simplicity to be dependent only on the level of exozodiacal dust in the system, N_ez, which linearly scales the fiducial spectrum. We then performed the same retrieval experiment as above with N_ez = 10 contamination, except with N_ez as an additional free parameter.

The right panel of Figure 4 shows the results of our exozodiacal retrieval for a subset of the free parameters that are most tied to the effects of exozodiacal emission. The zodi level in the system is constrained to N_ez = 10.26 ± 2.84, which is consistent with our injected value. The covariance between planet radius and exozodiacal level shows a clear anticorrelation, but no significant correlation is seen between zodi and irradiation temperature. Including N_ez as a free parameter successfully mitigates the radius bias; however, this comes at the cost of additional radius uncertainty due to the marginalization over the zodi uncertainty propagated through the radius–zodi covariance. Compared to the case with no exozodiacal light (N_ez = 0), the planet radius uncertainty in this scenario grows by a factor of ∼1.5. Although we achieve a slightly better fit to the spectrum when we include N_ez as a free parameter, using the BIC for model selection marginally favors the model without N_ez (ΔBIC = 1.7) due to the penalty of adding another dimension to the inference. This suggests that it may be difficult, in practice, to confidently diagnose the presence of exozodiacal contamination and mitigate potential biases, but systems with even greater zodi or a different system configuration may not have this issue.

Although we have demonstrated that the PIE technique may be a viable strategy for thermal studies of transiting and non-transiting hot exoplanet atmospheres, more work will be needed to expand on our findings and to relax assumptions that we have made (e.g., static stellar and planetary-flux sources, single-planet chemical composition, no systematic noise, etc.). Most critically, we have assumed a simple stellar spectral model that does not account for stellar granulation, spots, or faculae, which will all add complexity to the stellar spectrum inference problem that must be solved for any planet information to be inferred. A more sophisticated stellar model with multiple spectral components for cool spots and hot faculae may cause degeneracies in the stellar fit that could decrease the precision on inferred stellar parameters from those found in this work. It remains unclear if these stellar complexities will be degenerate with the planetary spectrum and parameters since their temperatures differ sufficiently well to be separated. Additionally, we have considered only a WASP-43b-like planet with a ∼ 1000 K night side. The PIE technique will only increase in difficulty toward smaller and cooler planets that emit less flux. More work is needed across a broad range in planet and star temperatures to identify the optimal and plausible bounds for the PIE technique with JWST and other future mission concepts.

Finally, our exploration into exozodiacal contamination highlights how additional unresolved flux sources may confuse and bias the PIE technique. We showed that increasing the complexity of our PIE retrieval model to account for exozodiacal light did mitigate the bias, but in the same stroke decreased the precision of our inferred planet radius. This serves as a cautionary warning if too many unresolved sources and/or yet unidentified complicating effects enter the picture: the fidelity of planetary information could be lost to these degeneracies. These questions and others require additional study in order to advance the maturity of the PIE technique. Ultimately, our work has demonstrated that the out-of-transit baseline observations that are already planned with JWST for many hot Jupiters will provide an opportune test bed for night-side atmospheric characterization using the PIE technique.