Phase-curve Pollution of Exoplanet Transmission Spectra

The procedure requires a configuration file to perform the full set of predetermined calculations for a given exoplanet transit or eclipse. The user can choose between several model grids with precalculated stellar spectra, use a blackbody spectrum, or read the stellar spectrum from a file. The choice for exoplanets is between blackbody and user file spectra to model their dayside and nightside emission. In the case of spectra read from user files, it is necessary to specify whether these spectra are to be rescaled by the distance and radius of the source to obtain the observed flux. The input planetary spectra can be given in absolute units or normalized with respect to the stellar spectrum, such as those typically obtained by running forward models of the secondary eclipses with NEMESIS (Irwin et al. 2008), CHIMERA (Line et al. 2013), and Tau-REx II (Waldmann et al. 2015) and III (Al-Refaie et al. 2019).

The user must request one or more instrumental passbands. Some passbands are built into the ExoTETHyS package through text files that report the photon-to-electron conversion efficiencies (PCEs) at given wavelengths. User files with the same format can also be accepted. Optionally, the passbands can be split into multiple wavelength bins, as is usual in exoplanet transit spectroscopy.

Other mandatory input information includes the area of the telescope (${{ \mathcal A }}_{\mathrm{tel}}$), the duration of the simulated observation (Δt_obs), and a set of parameters of the star–planet system. The distance of the system from the Earth (d) and the telescope area can be used to rescale the template spectra in order to match the photon count rates received by the detector. We anticipate here that rescaling by d and ${{ \mathcal A }}_{\mathrm{tel}}$ does not affect the transit depth biases, only the nominal error bars.

The input stellar parameters are the effective temperature (T_*,eff), surface gravity ($\mathrm{log}{g}_{* }$), metallicity ([M/H]_*), and radius (R_*). The stellar models in the grids are identified by T_*,eff, $\mathrm{log}{g}_{* }$, and [M/H]_*. The default values of $\mathrm{log}{g}_{* }=4.5$ and [M/H]_* = 0.0 are implemented if these parameters are missing in the configuration file. The blackbody spectrum is defined by T_*,eff alone. If the stellar spectrum is obtained from a user file, none of these three parameters is necessary, unless T_*,eff is required to calculate the planetary spectra. Here R_* is always mandatory.

The orbital period (P), semimajor axis (a), inclination (i), and planet radius (R_p ) are needed to determine the geometry of the transit. The planet's emission is modeled through two spectra from the dayside and nightside hemispheres, either blackbody or user-supplied spectra. The blackbody spectra can be defined by directly setting the dayside and nightside temperatures (T_day, T_night) or by means of the bond albedo (A_b ) and circulation efficiency (ε).

Each physical parameter, except for temperatures and dimensionless parameters, is defined through two keywords in the configuration file to give its numerical value and the corresponding unit. The astropy.units package is used to interpret the input units and handle operations between physical quantities. Additionally, R_p and a can be expressed in units of R_*, and Δt_obs can be relative to the total transit duration (T₁₋₄).

The output is a text or pickle file containing estimates of transit depth biases and error bars for the requested exoplanet systems and passbands.

The stellar model atmosphere grids consist of one file for each triple of stellar parameters (T_*,eff, $\mathrm{log}{g}_{* }$, [M/H]_*) containing the emergent spectrum at the star surface in units of erg cm⁻² s⁻¹ Å⁻¹. The same files also include spatially resolved intensities, which are used by other subpackages to model limb-darkening. The BOATS algorithm first calculates the spectral photon rates collected by the telescope in units of photon s⁻¹ Å⁻¹. In order to do so, the stellar and planetary spectra are multiplied by ${(R/d)}^{2}{{ \mathcal A }}_{\mathrm{tel}}$, where R is the object's radius. The geometric dilution factor (R/d)² can be set to 1 for a user file spectrum by selecting the no-rescale flux option.

The next step is to compute the electron rates of detectors. This is done by multiplying the photon rates by the PCE and integrating on the passband or wavelength bin. In the following subsections, we will refer to fluxes (symbol F) as equivalent to electron rates in e⁻ s⁻¹, calculated as described here.

The BOATS calculations rely on a toy model to describe the phase-curve modulations. The observed planet flux is

$$\begin{eqnarray}&&{F}_{p}=\psi {F}_{\mathrm{day}}+(1-\psi ){F}_{\mathrm{night}},\end{eqnarray} \tag{ 1 }$$

where F_day and F_night denote the emergent flux from the two hemispheres of the planet, and ψ is the fraction of the area of the planet's projection showing the dayside. Other assumptions (valid for ExoTETHyS version 2.0.0) are that the planet lies on a circular, edge-on orbit in synchronous rotation with zero obliquity, and the dayside is centered on the substellar point. Under these hypotheses, we have

$$\begin{eqnarray}&&\psi =\displaystyle \frac{1-\cos \varphi }{2},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&\varphi =\displaystyle \frac{2\pi }{P}t\end{eqnarray} \tag{ 3 }$$

is the orbital phase angle, and t is the time from mid-transit. This toy model is identical to that adopted by the ExoNoodle software (Martin-Lagarde et al. 2020a, 2020b) but with one less degree of freedom relative to the centering of the dayside hemisphere with respect to the substellar point. It relies on the simplest hypotheses that enable order-of-magnitude estimates of the possible self- and phase-blend effects, solely based on the easiest-to-measure system parameters.

If Blackbody is the selected grid of planetary models, the two temperatures, T_day and T_night, can also be determined from the atmospheric parameters A_b and ε (Cowan & Agol 2011),

$$\begin{eqnarray}&&\left\{\begin{array}{l}{T}_{\mathrm{day}}={T}_{\mathrm{irr}}{\left(1-{A}_{b}\right)}^{\tfrac{1}{4}}{\left(\displaystyle \frac{2}{3}-\displaystyle \frac{5}{12}\varepsilon \right)}^{\tfrac{1}{4}}\\ {T}_{\mathrm{night}}={T}_{\mathrm{irr}}{\left(1-{A}_{b}\right)}^{\tfrac{1}{4}}{\left(\displaystyle \frac{\varepsilon }{4}\right)}^{\tfrac{1}{4}},\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where the irradiation temperature is

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}={T}_{* ,\mathrm{eff}}\sqrt{\displaystyle \frac{{R}_{* }}{a}}.\end{eqnarray} \tag{ 5 }$$

In general, the dayside flux is the sum of planetary intrinsic emission and reflected starlight,

$$\begin{eqnarray}&&{F}_{\mathrm{day}}={F}_{\mathrm{day}}^{\mathrm{emission}}+{F}_{\mathrm{day}}^{\mathrm{reflection}},\end{eqnarray} \tag{ 6 }$$

where ${F}_{\mathrm{day}}^{\mathrm{emission}}$ is calculated from the corresponding blackbody or user spectrum, and

$$\begin{eqnarray}&&{F}_{\mathrm{day}}^{\mathrm{reflection}}={A}_{b}{\left(\displaystyle \frac{{R}_{p}}{2a}\right)}^{2}{F}_{* },\end{eqnarray} \tag{ 7 }$$

with F_* being the stellar flux. The nightside flux consists of a pure planetary emission component.

The next important step is to calculate the mean flux of the planet during and out of the transit. Given the symmetry of the phase-curve model described in Section 2.3, and assuming that the observation interval is centered on the transit, it is sufficient to consider the half observation after the mid-transit point. The transit event is delimited by the external contact points; the corresponding duration is (Seager & Mallén-Ornelas 2003)

$$\begin{eqnarray}&&{T}_{1-4}=\displaystyle \frac{P}{\pi }\arcsin \left(\displaystyle \frac{\sqrt{{\left(1+\tfrac{{R}_{p}}{{R}_{* }}\right)}^{2}-{\left(\tfrac{a}{{R}_{* }}\cos i\right)}^{2}}}{\tfrac{a}{{R}_{* }}\sin i}\right).\end{eqnarray} \tag{ 8 }$$

The phase angles corresponding to the end of the transit and observation are

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Phi }}}_{1-4}}{2}=\displaystyle \frac{\pi }{P}{T}_{1-4}{\rm{a}}{\rm{n}}{\rm{d}}\,\displaystyle \frac{{\rm{\Delta }}{\varphi }_{{\rm{o}}{\rm{b}}{\rm{s}}}}{2}=\displaystyle \frac{\pi }{P}{\rm{\Delta }}{t}_{{\rm{o}}{\rm{b}}{\rm{s}}}.\end{eqnarray} \tag{ 9 }$$

By using Equations (1)–(3), the phase-averaged fluxes are

$$\begin{eqnarray}&&{F}_{p}^{{\rm{i}}{\rm{n}}}=\displaystyle \frac{2}{{{\rm{\Phi }}}_{1-4}}{\int }_{0}^{{\textstyle \tfrac{{{\rm{\Phi }}}_{1-4}}{2}}}{F}_{p}d\phi ={F}_{{\rm{d}}{\rm{a}}{\rm{y}}}{\bar{\psi }}^{{\rm{i}}{\rm{n}}}+{F}_{{\rm{n}}{\rm{i}}{\rm{g}}{\rm{h}}{\rm{t}}}\left(1-{\bar{\psi }}^{{\rm{i}}{\rm{n}}}\right){\rm{a}}{\rm{n}}{\rm{d}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{F}_{p}^{\mathrm{out}} & = & \displaystyle \frac{2}{{\rm{\Delta }}{\varphi }_{\mathrm{obs}}-{{\rm{\Phi }}}_{1-4}}{\int }_{\tfrac{{{\rm{\Phi }}}_{1-4}}{2}}^{\tfrac{{\rm{\Delta }}{\varphi }_{\mathrm{obs}}}{2}}{F}_{p}d\phi \\ & = & \,{F}_{\mathrm{day}}{\overline{\psi }}^{\mathrm{out}}+{F}_{\mathrm{night}}\left(1-{\overline{\psi }}^{\mathrm{out}}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

with

$$\begin{eqnarray}&&{\bar{\psi }}^{{\rm{i}}{\rm{n}}}=\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{\sin \left({\textstyle \tfrac{{{\rm{\Phi }}}_{1-4}}{2}}\right)}{{\textstyle \tfrac{{{\rm{\Phi }}}_{1-4}}{2}}}\right){\rm{a}}{\rm{n}}{\rm{d}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\overline{\psi }}^{\mathrm{out}}=\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{\sin \left(\tfrac{{\rm{\Delta }}{\varphi }_{\mathrm{obs}}}{2}\right)-\sin \left(\tfrac{{{\rm{\Phi }}}_{1-4}}{2}\right)}{\tfrac{{\rm{\Delta }}{\varphi }_{\mathrm{obs}}}{2}-\tfrac{{{\rm{\Phi }}}_{1-4}}{2}}\right).\end{eqnarray} \tag{ 13 }$$

The transit depth is defined as the planet-to-star projected area ratio,

$$\begin{eqnarray}&&{p}^{2}={\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right)}^{2}.\end{eqnarray} \tag{ 14 }$$

Finally, the self- and phase-blend effects are calculated as (Martin-Lagarde et al. 2020b; see their Figure 2)

$$\begin{eqnarray}&&{\left({\rm{\Delta }}{p}^{2}\right)}_{{\rm{s}}{\rm{e}}{\rm{l}}{\rm{f}}{\textstyle \text{-}}{\rm{b}}{\rm{l}}{\rm{e}}{\rm{n}}{\rm{d}}}=-\displaystyle \frac{{F}_{p}^{{\rm{o}}{\rm{u}}{\rm{t}}}}{{F}_{\ast }+{F}_{p}^{{\rm{o}}{\rm{u}}{\rm{t}}}}{p}^{2}{\rm{a}}{\rm{n}}{\rm{d}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\left({\rm{\Delta }}{p}^{2}\right)}_{\mathrm{phase} \mbox{-} \mathrm{blend}}=\displaystyle \frac{{F}_{p}^{\mathrm{out}}-{F}_{p}^{\mathrm{in}}}{{F}_{* }+{F}_{p}^{\mathrm{out}}},\end{eqnarray} \tag{ 16 }$$

and the total bias in transit depth is

$$\begin{eqnarray}&&{\rm{\Delta }}{p}^{2}={\left({\rm{\Delta }}{p}^{2}\right)}_{\mathrm{self} \mbox{-} \mathrm{blend}}+{\left({\rm{\Delta }}{p}^{2}\right)}_{\mathrm{phase} \mbox{-} \mathrm{blend}}.\end{eqnarray} \tag{ 17 }$$

In the photon noise limit, the minimum error bar depends on the number of electrons produced in and out of the transit. Ignoring the ingress, egress, and limb-darkening effects, these numbers are

$$\begin{eqnarray}&&{N}^{{\rm{i}}{\rm{n}}}=\left({F}_{\ast }(1-{p}^{2})+{F}_{p}^{{\rm{i}}{\rm{n}}}\right){T}_{1-4}{\rm{a}}{\rm{n}}{\rm{d}}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{N}^{\mathrm{out}}=\left({F}_{* }+{F}_{p}^{\mathrm{out}}\right)\left({\rm{\Delta }}{t}_{\mathrm{obs}}-{T}_{1-4}\right).\end{eqnarray} \tag{ 19 }$$

The nominal error bar is

$$\begin{eqnarray}&&{\sigma }_{{p}^{2}}=\left(1-{p}^{2}\right)\sqrt{\displaystyle \frac{1}{{N}^{\mathrm{in}}}+\displaystyle \frac{1}{{N}^{\mathrm{out}}}}.\end{eqnarray} \tag{ 20 }$$

We note that the actual error bar obtained by fitting the light-curve data can be larger because of the parameter correlations and additional noise components.

We took the transmission spectra observed with HST/WFC3 as starting points. In particular, we reanalyzed the low-resolution spectra at 1.1–1.7 μm reported by Tsiaras et al. (2018) using the novel Tau-REx III retrieval algorithm (Al-Refaie et al. 2019). As the narrow wavelength range limits the ability to detect molecules other than H₂O, we adopted atmospheric chemical equilibrium (ACE) models (Agúndez et al. 2012; Venot et al. 2012) instead of the most popular free chemistry to include more absorbing species. The molecular absorption cross sections were computed by the ExoMol group (Chubb et al. 2021). We also included collisionally induced absorption (CIA) by H₂–H₂ and H₂–He (Gordon et al. 2017) and Mie scattering (Lee et al. 2013). We assumed an isothermal temperature profile.

The free parameters in our retrievals were the planet radius at 10 bars, atmospheric temperature, Mie scattering cloud particle size, mixing ratio, and composition. The carbon-to-oxygen (C/O) ratio and metallicity were fixed to solar values. The best-fit solutions obtained with these settings have the same Bayesian evidence as those reported by Tsiaras et al. (2018) to within ${\rm{\Delta }}\mathrm{log}E\leqslant 2$, despite using only five free parameters instead of 12 and 10 for WASP 12 b and WASP 43 b, respectively.

We calculated forward model spectra that extend the best-fit solutions over a broader wavelength range, then bin-averaged. In this way, we obtained the transmission spectra associated with three JWST instruments using the same observing modes planned for the Transiting Exoplanet Community ERS program (Bean et al. 2018):

• 1.  
  the Near-InfraRed Imager and Slitless Spectrograph (NIRISS) with the Single-Object Slitless Spectroscopy (SOSS) mode,
• 2.  
  the Near-InfraRed Spectrograph (NIRSpec) with the G235H and G395H gratings, and
• 3.  
  the MIRI with the Low-Resolution Spectroscopy (LRS) slitless mode.

Table 1 reports the nominal wavelength ranges for the three instruments and the size of the spectroscopic bins. Figure 3 shows the corresponding PCEs. The PCE for each observing setup is the combined throughput from the telescope, instrument optics, detector efficiency, quantum efficiency, and filter throughput (if applicable). We adopted the reference data files from the Pandeia Exposure Time Calculator (Pontoppidan et al. 2016) and precalculated response functions for the MIRI LRS slitless (Kendrew et al. 2015) and NIRISS SOSS (Maszkiewicz 2017). The resulting PCEs are among the passbands available within the ExoTETHyS package.

Zoom In Zoom Out Reset image size

Additionally, we considered another set of reference transmission spectra with identical parameters, except for C/O = 1.5.

We used ExoTETHyS.BOATS to calculate the bias in transit depth on the same wavelength bins as the transmission spectra. Table 2 reports the system parameters. Table 3 contains the dayside and nightside temperatures based on current estimates from observations. We modeled the emission from the two planetary hemispheres in the following ways:

• 1.  
  assuming blackbody spectra with T_day and T_night from Table 3 (most likely configurations) and
• 2.  
  calculating the emission spectra with Tau-REx III, assuming ACE with solar C/O ratio (and C/O = 1.5) and metallicity, CIA, Rayleigh scattering, and temperature profiles as in Guillot (2010) with the same T_day and T_night.

We then added the spectroscopic bias in three forms to the pure transmission spectra to obtain the biased spectra:

• 1.  
  total bias resulting from the self- and phase-blend effects;
• 2.  
  total bias minus average offset, such that the sum of the biases on the wavelength bins is zero; and
• 3.  
  self-blend bias only.

The case of total bias minus average offset may occur if a constant correction is applied at all wavelengths. For example, when the spectral light curves are divided by the white one (e.g., Tsiaras et al. 2016), if the bias in transit depth is removed for the white light curve, the spectral values will also be shifted, but they remain biased. The self-blend bias is what remains if the phase-blend effect is removed from the data for each wavelength (Martin-Lagarde et al. 2020b). Figure 4 shows the pure transmission and biased spectra considered for WASP 12 b and WASP 43 b. The biased spectra appear to have an offset and a different slope than the pure transmission ones in the NIRISS passband. The effect on the slope is smaller at the longer wavelengths, and it is negligible in the MIRI passband. The bias for WASP 12 b is much larger than the 1σ error bars in most configurations, while for WASP 43 b, it is within about 1σ.

Zoom In Zoom Out Reset image size

We performed atmospheric retrievals on the pure and biased transmission spectra using Tau-REx III. The free fitting parameters were the planet radius, the isothermal temperature, three Mie scattering parameters, and the abundances of all of the molecules used by the ACE module to compute the transmission spectra. We take the approximation of uniform chemical distribution with height for the retrievals, neglecting the pressure dependence of the input models.

In all cases, the retrieval algorithm found a solution that provides an excellent match to the simulated spectrum, albeit with different than input parameters. Figure 5 compares the atmospheric parameters retrieved from pure and biased transmission spectra with their underlying values and/or intervals from the input models. The spectrum bias can affect the H₂O abundance measured on WASP 12 b by orders of magnitude, the temperature by ∼500 K, and the planet radius by 2%–3%. In some cases, the measured discrepancies may include the effect of parameter degeneracies, which are discussed below. The largest effects occur with JWST/NIRISS, for which the biased spectra have a significantly different slope than the pure transmission ones. We obtained similar trends using blackbody or ACE spectra for planetary emission, but it is possible to have larger differences in the presence of more pronounced molecular features. The constant offsets change only the apparent radius of the planet, not the other atmospheric parameters.

Zoom In Zoom Out Reset image size

Unsurprisingly, the parameters obtained for WASP 43 b are not significantly affected by the modest spectroscopic bias. If the phase-blend effect is corrected, the self-blend bias is negligible for both WASP 12 b and WASP 43 b.

We note that, in some cases, the retrieved parameters do not match the underlying values, even in the absence of spectrum bias. These discrepancies can only be partially explained by the approximation of constant gas profiles adopted in the retrievals (Changeat et al. 2019). The error bars obtained from atmospheric retrievals can sometimes be underestimated, not accounting for the full parameter degeneracies (Rocchetto et al. 2016). Performing retrievals over the broader wavelength range covered by the three instruments should help to significantly reduce the degeneracies between the planet's radius and atmospheric parameters (Chapman et al. 2017; Welbanks & Madhusudhan 2019). A discussion of the robustness of the retrievals themselves is beyond the scope of this paper.