Global Chemistry and Thermal Structure Models for the Hot Jupiter WASP-43b and Predictions for JWST

ATMO is a 1D/2D atmospheric model that solves the radiative/convective equilibrium with and without irradiation from an host star. It has been used for the study of brown dwarfs and directly imaged exoplanets (Tremblin et al. 2015, 2016, 2017a; Leggett et al. 2016, 2017), and also for the study of irradiated exoplanets (Drummond et al. 2016; Evans et al. 2017; Wakeford et al. 2017). The gas opacity is computed by using the correlated-k method (Lacis & Oinas 1991; Amundsen et al. 2014, 2017), including the following species in this study: H₂–H₂, H₂–He, H₂O, CO, CO₂, CH₄, NH₃, K, Na, Li, Rb, Cs, from the high-temperature ExoMol (Tennyson & Yurchenko 2012) and HITEMP (Rothman et al. 2010) line list databases. We use 32 frequency bins between 0.2 and 320 μm with 15 k-coefficients per bin. The chemistry is solved either at equilibrium or out of equilibrium by a consistent coupling with the chemical kinetic network of Venot et al. (2012). Recently, 1D-ATMO has been benchmarked against Exo-REM and petitCODE (Baudino et al. 2017).

In this study, we have used 2D-ATMO, an extension of 1D-ATMO (Tremblin et al. 2017b) that takes into account the circulation induced by the irradiation from the host star at the equator of the planet. We have taken a Kurucz spectrum (Castelli & Kurucz 2003) for WASP-43 with a radius of 0.667 R_⊙, an effective temperature of 4500 K, and a gravity of log(g) = 4.5. The magnitude of the zonal wind is imposed at the substellar point at 4 km s⁻¹ and is computed accordingly to the momentum conservation law in the rest of the equatorial plane. The vertical mass flux is assumed to be proportional to the meridional mass flux, with a proportionality constant 1/α; the wind is therefore purely longitudinal and meridional if $\alpha \to \infty$, or purely longitudinal and vertical for $\alpha \to 0$. As in Tremblin et al. (2017b), a relatively low value of α drives the vertical advection of entropy/potential temperature in the deep atmosphere that can produce a hot interior, which can explain the inflated radii of hot Jupiters. A high value of α will produce a "cold" deep interior, as in the standard 1D models. In this study, we have used two values of α, 10 and 10⁴ to explore these two limits. The simulation with α = 10⁴ should be more representative of WASP-43b, since the planet is not highly inflated.

SPARC/MITgcm couples a state-of-the-art, nongray, radiative transfer code with the MITgcm (Showman et al. 2009). The MITgcm solves the primitive equations on a cubed-sphere grid (Adcroft et al. 2004). It is coupled to the nongray radiative transfer scheme based on the plane-parallel radiative transfer code of Marley & McKay (1999). The stellar irradiation incident on WASP-43b is computed with a Phoenix model (Hauschildt et al. 1999). The opacities we use are described in Freedman et al. (2008), including more recent updates (Freedman et al. 2014), and the molecular abundances are calculated assuming local chemical equilibrium (Visscher et al. 2010). In the 3D simulation, the radiative transfer calculations are performed on 11 frequency bins ranging from 0.26 to 300 μm, with eight k-coefficients per bin statistically representing the complex line-by-line opacities. For the purpose of calculating the spectra, the final SPARC/MITgcm thermal structure is post-processed with the same radiative transfer code, but using a higher spectral resolution of 196 spectral bins (Fortney et al. 2006). We initialize the code with the analytical planet-averaged pressure-temperature profile of Parmentier et al. (2015a), run the simulation for 300 days, and average all physical quantities over the last 100 days of simulation.

Our baseline model is a solar-composition, cloudless model. We also performed simulations including the presence of radiatively active clouds and radiatively passive clouds following the method outlined in Parmentier et al. (2016). Finally, models assuming a constant [CO]/[CH₄] ratio were also performed, following the method outlined in Steinrueck et al. (2019). These latter models simulate out-of-equilibrium, transport-induced quenching of CO and CH₄.

To address the variations of atmospheric chemical composition with altitude and longitude, we used both 1D and 2D chemical kinetics models. We describe these two codes here.

3.3.1. One-dimensional Chemical Kinetics Model

The thermo/photochemical model developed by Venot et al. (2012) is a full 1D time-dependent model. This model takes into account a detailed chemical kinetics and the out-of-equilibrium processes of photodissociation and vertical mixing (eddy and molecular diffusion). The atmospheric composition is computed for a fixed thermal profile divided in discrete layers, solving the continuity equations for each species until steady state is reached. No flux of species is imposed at the boundaries of the atmosphere. In this study, we used the C₀–C₂ chemical kinetic network, which contains ∼2000 reactions describing the kinetics of 105 species made of H, C, O, and N, with up to two carbon atoms. This chemical scheme has been developed in close collaboration with specialists in combustion and validated experimentally on wide ranges of pressure and temperature as a whole (i.e., not only each reaction individually) leading to a high reliability. Since we consider both direction (forward and reverse) for each reaction, in the absence of out-of-equilibrium processes, thermochemical equilibrium is achieved kinetically. To these ∼2000 reactions, 55 photodissociations have been added to model the interaction of incoming UV flux with molecules. These reactions are not reversed, of course, due to the disequilibrium irreversible nature of this process. A complete description of the model and the chemical scheme can be found in Venot et al. (2012). The kinetic model has been applied to several exoplanetary atmospheres (Venot et al. 2012, 2013, 2014, 2016; Agúndez et al. 2014b; Rocchetto et al. 2016; Tsiaras et al. 2016), as well as the deep atmospheres of Saturn (Mousis et al. 2014), Uranus (Cavalié et al. 2014, 2017), and Neptune (Cavalié et al. 2017).

Using this 1D chemical kinetics model, we determined the chemical composition of the atmosphere of WASP-43b at different longitudes around the equator. The vertical temperature structure as a function of longitude is taken from the 2D radiative-transfer models (Section 3.1) utilizing α = 10⁴, since the planet is not highly inflated. We extended these profiles to higher altitudes using extrapolation, and the resulting profiles are shown in Figure 3. Each longitude has been computed until steady state is reached, independently from each other, starting from the thermochemical equilibrium corresponding to the thermal structure. The stellar zenith angle varied with longitude. Given the characteristics of the host star, we estimated the incident stellar spectrum using the same Kurucz (Castelli & Kurucz 2003) stellar model as in 2D-ATMO for wavelengths ≥2200 Å, an IUE spectrum of GL 15B scaled by a factor of 10 for the 1250–2200 Å region, and the solar maximum spectrum of Woods & Rottman (2002) scaled by a factor of 1.36 (see Czesla et al. 2013) for wavelengths less than 1250 Å. Vertical transport operates through eddy and molecular diffusion, with an assumed eddy diffusion coefficient K_zz profile that varies as K_zz (cm² s⁻¹) = 10⁷[P(bar)]^−0.65 throughout the bulk of the planet's infrared photosphere, independent of longitude, except that K_zz approaches a constant-with-altitude value of 10¹⁰ cm² s⁻¹ at high altitudes and at pressures greater than 300 bar (Figure 2). As a large uncertainty resides for this parameter, we chose to use an expression similar to the one adopted for HD 189733b based on GCM passive tracer transport (Agúndez et al. 2014a). This method gives values for the vertical mixing much lower (up to 1000 times) than the classical root-mean-square method (Agúndez et al. 2014a; Charnay et al. 2015). A correct estimation of this parameter is crucial, as it determines the quenching level—and thus the molecular abundances of some fundamental species (e.g., Miguel & Kaltenegger 2014; Venot et al. 2014; Tsai et al. 2017). Solar proportions of elemental abundances are assumed (Lodders 2010), with a depletion of 20% in oxygen due to the sequestration in silicates and metals (Lodders 2004). No TiO/VO has been included.

Zoom In Zoom Out Reset image size

3.3.2. Two-dimensional Chemical Kinetics Model

Based on the procedure of Agúndez et al. (2014a), for this study we developed a "pseudo-2D" chemical model to track how the atmospheric composition of WASP-43b would vary as a function of altitude and longitude—and hence orbital phase. We assume that longitudes are not isolated from each other (as in the 1D chemical model), but rather are connected through a strong zonal jet. Tidal interactions between host stars and gas giant planets are expected to circularize the planet orbit and synchronize the rotation period of the planet (Lubow et al. 1997; Guillot & Showman 2002). The timescale for this to happen is much shorter than the stellar lifetime when the planet orbit is shorter than 10 days, such as for WASP-43b (see Figure 1 of Parmentier et al. 2015b). With one hemisphere constantly facing the star, the unequal stellar forcing produces strong zonal winds that transport heat and chemical constituents across the dayside and into the nightside of the planet, and back again (e.g., Showman et al. 2010; see also the GCM results in Section 4.4). This zonal transport provides "diurnal" variation in temperatures and stellar irradiation from the point of view of a parcel of gas being transported by the winds. If chemical equilibrium were to prevail in WASP-43b's atmosphere, the composition would be very different on the colder nightside in comparison to the hotter dayside. However, if the horizontal transport by the zonal winds is faster than the rate at which chemical constituents can be chemically converted to other constituents, then the composition can be "quenched" and remain more uniform with longitude (e.g., Cooper & Showman 2006; Agúndez et al. 2014a). Our pseudo-2D model tracks the time variation in atmospheric composition as a function of altitude and longitude, as an atmospheric column at low latitudes experiences this pseudo-rotation.

Our pseudo-2D thermo/photochemical kinetics model uses the Caltech/JPL KINETICS code (Allen et al. 1981) as its base, modified for exoplanets as described in Moses et al. (2011, 2013b, 2013a, 2016), along with the pseudo-2D procedure described in Agúndez et al. (2014a). The model contains 1870 reactions (i.e., 900 forward and reverse pairs) involving 130 different carbon-, oxygen-, nitrogen-, and hydrogen-bearing species whose rate coefficients have been reversed using the thermodynamic principle of reversibility (see Visscher & Moses 2011). Photolysis reactions are included and have not been reversed. The reaction list is taken from the GJ 436b model of Moses et al. (2013a), and includes H, C, O, N species. Molecules with up to six carbon atoms are included, but the possible chemical production and loss pathways in the model become less complete as the molecule becomes heavier. Note that the chemical reactions list used in the pseudo-2D model slightly differs from that of the 1D chemical model. The departures between these two chemical schemes and their implications on the calculated abundances have already been addressed in several studies (Venot et al. 2012; Moses 2014; Wang et al. 2016): depending on the scheme used, the quench level and resulting quenched abundances of some species are somewhat different. For instance, for HD 209458b-like planets, the quenched mixing ratio of CH₄ will be twice as large with Moses et al.'s (2013a) scheme as with Venot et al.'s (2012). However, the goal of our study here is to qualitatively compare the expected longitudinal variation obtained with the 1D and pseudo-2D chemical models in order to evaluate the effect of horizontal circulation on the global composition and transport/eclipse observations, and to determine a rough CO/CH₄ ratio to be used in the 3D GCM. Comparing results between the 1D and pseudo-2D chemical models described here adequately addresses these goals. Other model inputs, including the boundary conditions and vertical diffusion coefficients adopted in the model, are identical to the 1D chemical kinetics model described above.

Following Agúndez et al. (2014a), our pseudo-2D approach is to solve the 1D continuity equations for a vertical column of gas at the equator as it rotates around the planet with a constant average low-latitude jet speed of 4.6 km s⁻¹ (based on the GCM results described in Section 4.4). As it has been discussed in Agúndez et al. (2012, 2014a), assuming a uniform zonal wind is a good approximation for an equatorial region dominated by a super-rotating jet (Kataria et al. 2015). This approximation might be inadequate for latitudes toward the poles, when the circulation regime is more complex. Both the stellar zenith angle and atmospheric temperatures vary with time as the gas column rotates through different longitudes, both affecting the atmospheric chemistry. The vertical temperature structure as a function of longitude is fixed, and is the same as that used in the 1D chemical kinetics model (Figure 3). The eddy diffusion profile (K_zz coefficient profile) is also similar to that of the 1D chemical kinetics model (Figure 2). The planet is divided into different longitude regions, and the system of differential equations making up the continuity equations is integrated over the amount of time it would take a parcel of gas at the equator to be transported from one discretized longitude to the next discretized longitude. At that point, the mixing ratios as a function of pressure at the end of the first longitude calculation are fed in as initial conditions to the next longitude calculation, with its new thermal structure and an incident UV flux that depends on the new zenith angle. The temporal evolution of this equatorial column of gas is followed for 20 full planetary "rotations" to provide sufficient time for the species produced photochemically at high altitudes to be transported down through the atmosphere to deeper regions where thermochemical equilibrium dominates. At that point, the "daily" longitude variations were consistent from one pseudo-rotation to the next. Agúndez et al. (2012) discuss the various advantages of beginning the pseudo-2D calculations at the hottest dayside conditions. Based on their discussion, we use the results of a 1D thermo/photochemical kinetics model for conditions at the substellar point (longitude 0°) that has been run long enough to reach steady state as our initial conditions for the pseudo-2D model.

Zoom In Zoom Out Reset image size

We use the Community Aerosol and Radiation Model for Atmospheres (CARMA; Turco et al. 1979; Toon et al. 1988; Jacobson & Turco 1994; Ackerman et al. 1995) to investigate the vertical and longitudinal distribution of clouds in the atmosphere of WASP-43b. While we do not include the resulting distributions in our GCM simulations, we will be able to extract insights into the effect of longitudinal temperature variations on cloud distributions. CARMA is a 1D cloud microphysics model that generates binned size distributions of aerosol particles as a function of altitude (pressure) in an atmospheric column by explicitly computing and balancing the rates of cloud particle nucleation, growth by condensation and coagulation, loss by evaporation, and transport by sedimentation, advection, and diffusion. This sets CARMA apart from simpler cloud condensation models (e.g., Fegley & Lodders 1994; Ackerman & Marley 2001), which assume cloud formation as soon as the condensate vapor saturates. The equations CARMA solves to evaluate the rates of these processes are presented in the Appendix of Gao et al. (2018). CARMA has been applied to aerosol processes across the solar system (Colaprete et al. 1999; Barth & Toon 2003; Bardeen et al. 2008; Gao et al. 2014, 2017), and has recently been used to simulate Al₂O₃, TiO₂, MgSiO₃ (enstatite), KCl, and ZnS clouds on exoplanets and brown dwarfs (Gao et al. 2018; Gao & Benneke 2018; Powell et al. 2018).

For this work, we include additional condensates that have been hypothesized to dominate the condensate mass in exoplanet atmospheres, including Mg₂SiO₄ (forsterite), Fe, Cr, MnS, and Na₂S (Lodders 1999; Visscher et al. 2006; Helling et al. 2008; Visscher et al. 2010). Additional condensates are possible, as shown in grain chemistry models (e.g., Helling & Woitke 2006), but we do not treat them here, as it would be computationally prohibitive. In addition, to reduce the number of different condensates, we assume that forsterite is the primary silicate condensate rather than modeling both forsterite and enstatite clouds. This is based on the argument that a rising parcel of vapor would see forsterite condense first, due to it having higher condensation temperatures than enstatite; this depletes Mg and SiO, such that the enstatite cloud that forms above the forsterite cloud should have significantly lower mass. We discuss the implications of this assumption in the Section 4.3. As with the treatment of enstatite in Powell et al. (2018), we assume forsterite and Fe clouds form by heterogeneously nucleating on TiO₂ seeds, as direct nucleation of these two species from vapor is slow (Helling & Woitke 2006; Gao et al. 2018). All other condensates are assumed to nucleate homogeneously. The saturation vapor pressures of Cr, MnS, and Na₂S are taken from Morley et al. (2012). The surface energy of Cr is calculated from the Eötvös rule, while for MnS and Na₂S, we assume the same surface energy as that of KCl. The size distribution for each condensate species is calculated separately, and so a distinct size distribution exists for each species.

As in the 1D and pseudo-2D chemical kinetics models, we use fixed pressure–temperature profiles described in Section 3.1 for our background atmosphere. All planetary parameters used are the same as those of the other models presented here, to ensure consistency. In the cloud microphysics model, we use a K_zz profile very similar to the one used in the chemical kinetics models (Figure 2), except we set a minimum K_zz of 10⁷ cm² s⁻¹. This change only affects pressures >1 bar, where the chemical kinetics model K_zz is <10⁷ cm² s⁻¹. We also reduce the high K_zz at pressures >300 bar to our minimum value. This was necessary to reduce model run time and numerical instabilities. An atmospheric column is simulated at each longitude, independently of the others, under the assumption that microphysical timescales are short compared to horizontal transport timescales, though this may not be the case for all pressure levels and particle sizes (Powell et al. 2018). For each column, we investigate distinct clouds composed of Al₂O₃, TiO₂, Mg₂SiO₄, Fe, Cr, MnS, Na₂S, ZnS, and KCl, though which clouds actually form depends on which species is supersaturated and their nucleation rates. We assume solar abundances for the limiting elements of these clouds, which are, in the same order, Al, TiO₂, Mg, Fe, Cr, Mn, Na, Zn, and KCl.

Importantly, the vertical, longitudinal, and particle size distributions computed by this model are not used to generate synthetic observations, as will be presented later in this work. This is due to the uncertainties in the material properties of some of the condensates (e.g., surface energies of MnS and Na₂S) and the way exoplanet clouds form, whether through homogeneous nucleation, heterogeneous nucleation on some foreign condensation nuclei, or grain chemistry (Helling & Woitke 2006). Instead, results from CARMA will be helpful for informing general GCM and retrieval studies, due to its ability to compute the relative abundances of different cloud species in the atmosphere of WASP-43b, thus indicating the species that affect the observations the most. Simplified cloud models can then be used to explore the parameter space around these results.

WASP-43b is the primary target for the "MIRI Phase Curve" observation that will be carried out as part of the Transiting Exoplanet JWST ERS Program. The goal is to observe a full orbit of WASP-43b, including two eclipses and one transit in the wavelength range 5–12 μm at a resolution R ∼ 100 with MIRI LRS (Low Resolution Spectroscopy) in slitless mode (Kendrew et al. 2015). In that program, the planetary emission spectra as a function of longitude will be measured and relevant atmospheric properties retrieved. We simulate the expected outcomes of this observation using the PandExo²¹ software program (Batalha et al. 2017). PandExo is a noise simulator specifically designed for transiting exoplanet observations with JWST and HST, and includes all observatory-supported time-series spectroscopy modes.

The input parameters for the star and planet are those indicated in Table 1. The stellar spectrum is obtained from the NextGen (Hauschildt et al. 1999) grid interpolated at the T_eff and log(g) of WASP-43, and is the same as used in the 3D SPARC/MITgcm. We consider a range of planetary emission spectra derived from the 3D SPARC/MITgcm model described in Section 3.2, with or without clouds, assuming thermochemical equilibrium or a quenched [CH₄]/[CO] ratio. These simulations are performed with inputs similar to those used for the JWST ERS Program proposal (PIs: N. Batalha, J. Bean, K. Stevenson; Bean et al. 2018): the radiative transfer models, the star and planet parameters and input spectra, and the observation parameters are the same (we simulate a broader range of planetary spectra here).

The planetary spectra are calculated from the emission integrated over the visible hemisphere. For this work, we simulate spectra with a spacing of 20° in longitude and use them as inputs for PandExo. We consider that we observe each longitude during one eighteenth of the orbital period (1.08 hr), and we use a baseline of twice the eclipse duration because we will observe two eclipses (2.32 hr). In practice, the longitude 0° will be in-eclipse, so we may have to split the orbit slightly differently, but for these simulations we treat this longitude as the other ones. The resolution and instrumental parameters are those of MIRI LRS. The wavelength range goes up to ∼14 μm, but we consider only the 5–12 μm range because the efficiency of LRS decreases significantly beyond 12 μm. We use a saturation level of 80% of the full well. The details of the noise modeling can be found in Batalha et al. (2017).

To retrieve the atmospheric properties of WASP-43b, we use two models: TauREx²² (Waldmann et al. 2015a, 2015b; Rocchetto et al. 2016) and the Python Radiative Transfer in a Bayesian framework (Pyrat Bay,²³ P. E. Cubillos et al. 2020, in preparation, J. Blecic et al. 2020a, 2020b, in preparation). Both, TauREx and Pyrat Bay are open-source retrieval frameworks that compute radiative-transfer spectra and fit planetary atmospheric models to a given set of observations. The atmospheric models consist of parameterized 1D profiles of the temperature and species abundances as a function of pressure, with atomic, molecular, collision-induced, Rayleigh, and cloud opacities. We decided to use two codes that do not use the same retrieval methods, in order to compare the results obtained and highlight the eventual biases that could emerge. We present hereafter the two codes.

3.6.1. TauREx

3.6.2. Pyrat Bay

The pressure/temperature structure of the atmosphere can be constrained with 2D and 3D models. A comparison is given in Figure 4 between the 3D SPARC/MITgcm model and the 2D-ATMO model with two different interior cases: the hottest one (α = 10) and the coldest one (α = 10⁴). For the two α values, the upper atmosphere is similar. In this part (pressures less than 1 bar), the agreement between the 2D and the 3D models is quite remarkable since there is no tuning of the 2D model on the 3D model, apart from the choice of horizontal wind speed. Since the chemistry and radiative transfer models are independent between the two codes, this agreement is also a sign that convergence between different GCMs and 2D steady-state circulation models can be reached for the pressure/temperature structure at and above the photosphere. As explained in Section 3.1, for pressures greater than 1 bar, the different α values lead to different temperatures in the deep atmosphere. Because the GCM is not fully converged at pressures larger than 10 bars, it likely produces spurious variations of the deep pressure–temperature profile. The shape of the deep flow structure and its influence in the upper atmospheric dynamics is still the subject of active research (Mayne et al. 2014, 2019; Carone et al. 2019) and out of the scope of this paper. Given that both the GCM and the 2D model produce very similar thermal structure in the observable atmosphere, we decided to use the outputs of the 2D model as inputs for the chemical and cloud formation models. We opted for the cold interior model (α = 10⁴), based on the premise that the planet does not appear to be highly inflated. However, a more thorough investigation of the deep thermal structure on the observable cloud properties (e.g., Powell et al. 2018) will be needed in the future to interpret the observations.

Zoom In Zoom Out Reset image size

We study the chemical composition of WASP-43b at different longitudes with our 1D and pseudo-2D models. The results are presented in Figure 5, together with the abundances at thermochemical equilibrium, corresponding to the same longitude-variable thermal structure. In the 1D model, all the longitudes have been computed independently, assuming thermochemical composition as initial composition at each longitude. The different vertical columns do not interact with each other. In the pseudo-2D model, the longitudes are not independent, instead interacting through horizontal circulation. As we explained in Section 3.3, the steady-state composition of the substellar point is given as an initial condition to the adjacent longitude. There, the evolution of chemical composition is calculated over the amount of time necessary for a parcel of gas to reach the next longitude, and so on.

Zoom In Zoom Out Reset image size

In all models, as the temperature is identical at each longitude for pressures greater than ∼10³ mbar, we find that species have also the same abundances, corresponding to the thermochemical equilibrium values. The composition varies with longitude above this region, more or less depending on species. For pressures lower than ∼10³ mbar, many of the atmospheric constituents would vary significantly with longitude if the atmosphere remained in thermochemical equilibrium throughout. Particularly noteworthy is the fact that CH₄ would be the dominant carbon-bearing constituent at high altitudes on the colder nightside in thermochemical equilibrium, while CH₄ would virtually disappear from the dayside and CO would become the dominant carbon-bearing constituent at all altitudes. The 1D kinetic model predicts that vertical quenching will reduce this variation, but there are still several orders of magnitude differences between the abundances of the dayside and those of the nightside. In contrast, the pseudo-2D model predicts much less variation with longitude, particularly in the 0.1–1000 mbar region that is probed at infrared wavelengths. The CO that forms on the hot dayside cannot be chemically converted to CH₄ quickly enough on the nightside, before the atmospheric parcels are carried by the zonal winds back to the dayside. These results confirm the findings by Cooper & Showman (2006), Agúndez et al. (2014a), Mendonça et al. (2018b), Drummond et al. (2018a), and Drummond et al. (2018b) that both vertical and horizontal chemical quenching are important in hot Jupiter atmospheres.

Another species whose abundance predicted by our kinetic model is very different from what is expected by thermochemical equilibrium is HCN. At thermochemical equilibrium, this species has the same abundance profile at each longitude, and its abundance decreases with increasing altitude. The 1D kinetic model predicts that this species will be quenched at around 10² mbar, leading to a higher abundance than what is predicted by thermochemical equilibrium on the nightside, and even a higher abundance on the dayside thanks to a photochemical production. Note that the quenching pressure we determine with our model is, of course, highly dependent on the K_zz profile we assume. At 10² mbar, K_ZZ is about 4.5 × 10⁷ cm² s⁻¹. As we said in Section 3.3, this parameter is rather uncertain and could vary by several orders of magnitude, e.g., typically 10⁶–10¹² cm² s⁻¹ between Parmentier et al. (2013) and Agúndez et al. (2014a). Consequently, with these extreme values, the pressure-level quenching of HCN could vary between 10 and 10⁵ mbar. In contrast to what has been found with the 1D kinetic model, our pseudo-2D model indicates that the abundance of HCN on the nightside will remain very high and close to that of the dayside, thanks to the horizontal circulation, in agreement with Agúndez et al. (2014a). Such a high abundance might be detectable via high-resolution spectroscopic observations in the near-infrared coupled to a robust detrending method (Hawker et al. 2018; Cabot et al. 2019). On JWST/MIRI observations, HCN could eventually appear in the 7–8 μm band, albeit spectra will probably be dominated by water absorption in this region, given the important abundance of H₂O in the atmosphere of WASP-43b (Rocchetto et al. 2016).

Similarly to Agúndez et al. (2014a), we find that, in addition to vertical quenching due to eddy diffusion, the horizontal circulation leads to horizontal quenching of chemical species. Globally, the atmosphere of WASP-43b has a chemical composition homogenized with longitude to that of the dayside. This is particularly true for pressures larger than 1 mbar, while variations of abundances between the day and nightside still remain at lower pressures.

In summary, the pseudo-2D model suggests that CH₄ would be a relatively minor constituent on WASP-43b at all longitudes, that photochemically produced HCN will be more abundant than CH₄ in the infrared photosphere of WASP-43b at all longitudes, and that the key spectrally active species H₂O and CO will not vary much with longitude on WASP-43b. Benzene (C₆H₆) is a proxy for photochemical hazes in the pseudo-2D model, and the strong increase in the benzene abundance at nighttime longitudes suggests that refractory hydrocarbon hazes could potentially be produced at night from radicals produced during the daylight hours (e.g., Miller-Ricci Kempton et al. 2012; Morley et al. 2013, 2015). Note that a recent experimental study demonstrates that refractory organic aerosols can be formed in hot exoplanet atmospheres with a C/O ratio higher than solar (Fleury et al. 2019).

Based on these chemical models, and because the variation of CH₄ with longitude could be observed with MIRI, we ran GCMs assuming chemical equilibrium and a fixed [CH₄]/[CO] ratio of 0.001, which is representative of the 2D chemical model in the 0.1–1000 mbar region.

We use the CARMA model to determine the physical and chemical properties of the clouds along different vertical columns at the equator of the planet. Assuming Mie-scattering particles, we calculate the cloud optical depth profile for each longitude from 0.35 to 30 μm. Figure 6 (top) shows the pressure levels at which the total cloud column optical depth (taking into account all cloud species) equals unity, and reveals large differences between the day- and nightsides. Specifically, the dayside temperature profile is such that most of the forsterite (Mg₂SiO₄) is cold-trapped below 100 bars, and although the forsterite condensation curve crosses the temperature profile again at lower pressures, the abundance of Mg there is sufficiently low so as to prevent optically thick clouds from forming (Figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the nightside, sufficiently low temperatures allow for the condensation of optically thick MnS and Na₂S clouds, such that the optical depth = 1 pressure level is above 0.1 bar blueward of 7 μm. As the typical particle sizes of these clouds are 1 to a few μm (Figure 6, bottom), they become optically thin at longer wavelengths, allowing for forsterite clouds to become visible, as shown by the 10 μm silicate feature. Note that this forsterite cloud is not the cold-trapped cloud at 100 bars. Instead, because the forsterite condensation curve crosses the temperature profile a second time, at higher pressures here than on the dayside, there is sufficient Mg to produce an optically thick "upper" cloud even after accounting for cold trapping. The shape of the effective particle radius profiles shown in Figure 6 (bottom) also reveals this transition in cloud composition with longitude, as particle size tends to increase toward the cloud base due to available condensate vapor supply and size-sorting by lofting and sedimentation. For example, while dayside profiles are largely smooth, corresponding to the dominance of the forsterite cloud, the nightside profiles feature an MnS cloud deck above 1 bar sitting atop the forsterite cloud below (Figure 7).

Our results suggest that whether forsterite or enstatite is considered the primary silicate condensate could strongly impact the dayside cloud opacity. As enstatite condenses at lower temperatures (Figure 3), it would not form a deep cloud at pressures >100 bars on the dayside, like forsterite. This lack of cold trapping may result in an optically thick cloud at lower pressures. This is in contrast to the nightside, where forsterite only dominates the cloud opacity at long wavelengths. We therefore expect that, since the cloud base of forsterite is only ∼50% higher in pressure than enstatite, forsterite and enstatite will have similar effects on the nightside spectra. Whether the forsterite or enstatite clouds are cold-trapped in the deep atmospheric layers depends on the microphysical behavior of the cloud (studied here), as well as the strength of the vertical mixing and the temperature in the deep atmosphere (Powell et al. 2018). Thorngren et al. (2019) recently predicted a connection between planet equilibrium temperature and their intrinsic flux, suggesting that cold traps on certain hot Jupiters may not exist due to high temperatures in the deep atmosphere. Determining the cloud chemical composition in the nightside of WASP-43b through our JWST/MIRI phase curve observation will provide insights into the presence of a deep cold trap, and thus test the predictions from Thorngren et al. (2019).

Our work decouples cloud microphysics from the radiation field and dynamics of the rest of the atmosphere, and thus we cannot treat cloud radiative feedback or cloud advection. Fully coupled 3D models that include cloud microphysics in the form of grain chemistry have been applied to other individual exoplanets in the past, including HD 189733b and HD 209458b, which have temperatures similar to that of WASP-43b (Helling et al. 2016; Lee et al. 2016; Lines et al. 2018a, 2018b). These works show that the mean particle radii vary between 1 and 100 μm between 0.1 and 100 bars, and that the composition of mixed cloud particles is dominated by enstatite, forsterite, iron, SiO, and SiO₂, with forsterite being more abundant than enstatite at most longitudes. This is similar to our results, though we do not consider SiO and SiO₂ in our model, while they do not consider sulfide clouds in theirs. Advection tends to smooth out cloud composition differences, which we do not capture in our work. One other major difference between the grain chemistry models and our model is the high abundance of small particles at low pressures in grain chemistry models, stemming from high nucleation rates at low pressures. In contrast, nucleation rates are the highest at the cloud base in our model (Gao et al. 2018), owing to the high atmospheric density there, and so we lack a low-pressure small particle population.

To summarize, our results show that, if silicates primarily form forsterite clouds, then the dayside of WASP-43b should be cloudless down to 100 bars, while the nightside cloud opacity should be dominated by MnS and Na₂S clouds shortward of 7 μm, and forsterite clouds at longer wavelengths. Cloud particle sizes on the nightside at the pressure levels where clouds become opaque are on the order of 1 to a few μm. On the other hand, if silicates primarily form enstatite clouds, then the dayside should be cloudier at pressures <100 bars, while the nightside cloud opacity would remain dominated by the sulfide clouds.

We use our 3D model to calculate the thermal structure of WASP-43b, assuming different chemical composition (thermochemical equilibrium and disequilibrium) and cloudy conditions (clear, MnS, and MgSiO₃). The temperature structure and CH₄ abundances for the cloudless chemical equilibrium and disequilibrium simulations are shown in Figure 8. From these models, we calculate the corresponding emission spectra for the dayside and nightside.

Zoom In Zoom Out Reset image size

The thermal structure of our cloudless, chemical equilibrium SPARC/MITgcm simulations are very similar to the one presented in Kataria et al. (2015), where the reader can find a thorough description of the atmospheric flows. While Kataria et al. (2015) focused on the effect of TiO, metallicity, and drag, we hereafter discuss the role of disequilibrium chemistry and clouds in shaping the nightside spectrum of the planet.

In the case of quenched carbon chemistry ([CH₄]/[CO] = 0.001), the dayside is slightly cooler and the nightside is slightly warmer at a given pressure level than for our chemical equilibrium case. However, the differences in the spectra seen in Figure 9 are mainly due to change in the opacities rather than changes in the thermal structure. On the dayside, where the [CH₄]/[CO] ratio is small at chemical equilibrium, our quenched and chemical equilibrium simulations are indistinguishable. On the nightside, the quenching removes the CH₄ absorption bands between 3 and 4 μm, and the ones between 7 and 9 μm are weakened in the quenched scenario, leading to a signature detectable by JWST/MIRI. Note that our GCM simulations approximate the [CH₄]/[CO] ratio to be constant throughout the atmosphere (horizontally and vertically) for computational reasons, while the pseudo-2D simulations in Section 4.2, as well as 3D simulations of WASP-43b with a simplified chemistry scheme (Mendonça et al. 2018b) find that the methane abundance is homogenized horizontally but decreases with increasing altitude. However, Steinrueck et al. (2019) found that the effect of disequilibrium chemistry on the thermal structure and phase curve is qualitatively similar for different constant [CH₄]/[CO] ratios as long as CO is the dominant carbon-bearing species. Therefore, it is likely that the effect of a horizontally homogenized [CH₄]/[CO] ratio that decreases with altitude is also qualitatively similar. Thus, our quenched simulation still provides a valuable estimate of the effects of disequilibrium chemistry.

Zoom In Zoom Out Reset image size

The cloudless simulations were also post-processed with cloud opacities. The post-processing allows for a quick estimate of the strength of potential signature of cloud properties in the emission spectrum, without the need to run additional, time-consuming, global circulation models. In Section 4.3, we found that the nightside of WASP-43b could be dominated by MnS, Na₂S, MgSiO₃, and/or Mg₂SiO₄. Following Parmentier et al. (2016), we explore two possible cloud compositions: MnS and MgSiO₃. Forsterite and enstatite have very similar opacities and condensation curves, ergo we chose to include just one of these silicate species. As we found with our microphysical cloud model (Section 4.3), the atmosphere of WASP-43b is cool enough for MgSiO₃ clouds to cover the whole planet, affecting both the dayside and the nightside of the planet. Conversely, MnS clouds can only form on the cooler nightside, and thus only affect the nightside's spectrum. Both MnS and MgSiO₃ clouds are able to sufficiently dim the nightside emission spectrum blueward of 5 μm in order to match the HST and the Spitzer observations. As shown in Figure 9, in all our models, the nightside flux remains observable with JWST/MIRI, even when the thermal emission is extremely small shortward of 5 μm. MnS and MgSiO₃ cloud composition could be distinguished spectrally by our JWST/MIRI phase curve observation through the observation of the 10 μm absorption band seen in the red models of Figure 9 (see also Wakeford & Sing 2015).

The effect of the cloud particle size is explored in Figure 10. Assuming that the formation of MnS clouds is limited by the available amount of manganese in a solar-composition atmosphere, the MnS clouds could be either transparent or optically thick in the JWST/MIRI bandpass depending on the size of their particles. Conversely, MgSiO₃, if present, should always be optically thick in the MIRI bandpass.

Zoom In Zoom Out Reset image size

The radiative feedback effect of the clouds in a hot Jupiter is a subject of intense research. The amplitude and spatial distribution of the cloud heating is extremely dependent on the cloud model used (e.g., Lee et al. 2016; Roman & Rauscher 2017, 2019; Lines et al. 2018b, 2019). In Figure 11, we show the spectrum resulting from a global circulation model incorporating the radiative feedback of MnS clouds. The clouds are opaque up to mbar pressures on the planet nightside and produce a strong greenhouse effect, leading to a warmer nightside—and thus a higher nightside flux. Horizontal heat transport from nightside to dayside changes the dayside thermal structure via thermal inversion, leading to a dayside spectrum dominated by emission features. Qualitatively, any nightside clouds should increase the nightside opacity and warm the atmosphere. The lower the pressure of the cloud photosphere, the higher the greenhouse effect of the clouds should be. As a consequence, the dayside photosphere should warm up through heat transport from the nightside to the dayside. The lower the photospheric pressure on the nightside, the larger the warming effect of the clouds on the dayside photosphere. The example shown here assumes the highest possible cloud, and is therefore likely to overestimate the effect of the nightside cloud on the dayside thermal structure. A deeper cloud will likely have a smaller impact on the dayside spectrum. Overall, because the amplitude and spatial distribution of the cloud heating is extremely dependent on the cloud model used (e.g., Lee et al. 2016; Roman & Rauscher 2017; Lines et al. 2018b, 2019), we decided to focus the reminder of the paper on the post-processed case by comparing the spectral effect of the clouds for a given thermal structure.

Zoom In Zoom Out Reset image size

The (cloud-free) quenched simulations are not a good match to existing HST and Spitzer nightside observations. Most likely, this is because the effect of nightside clouds dominates over the effect of disequilibrium chemistry at the wavelengths of existing observations. This situation is qualitatively similar to what Steinrueck et al. (2019) find for HD 189733b. However, quenched chemistry might still be important on the nightside of WASP-43b. We note that, during the referee process of this paper, Morello et al. (2019) published a new reduction of the Spitzer observations of WASP-43b. The resulting nightside fluxes lay between those of Stevenson et al. (2017) and Mendonça et al. (2018a).

We also note that a new 3D circulation model elicits a dynamical regime for WASP-43b that is different from that of SPARC/MITgcm. Carone et al. (2019) propose that the presence of very deep wind jets down to 700 bar leads to an interruption of super-rotation, and hence also an interruption of dayside-to-nightside heat transfer. In contrast to that, our SPARC/MITgcm model does not display very deep wind jets; it has uninterrupted super-rotation, and hence an efficient dayside-to-nightside heat transport. Thus, the Carone et al. (2019) model yields nightsides that are colder by several 100 K, compared to our model; see also Figure 7 in Carone et al. (2019), for a direct comparison.

For the nightside cloud-free retrievals, both Pyrat Bay and TauREx reproduce the simulated data well (Figure 13, top panel), and they obtain similar results for the retrieved temperature profile and the molecular abundances (Figure 13, bottom panels). Note, however, that the temperatures and abundances of the individual cells in the input 3D model span wide ranges, which vary with latitude, longitude, and pressure (see e.g., Figure 8). Since the output emission spectrum is not a linear transformation of the input temperature, one must consider the averaged values of the input model as guidelines rather than a strict measure of accuracy.

Zoom In Zoom Out Reset image size

Given the properties of the system, the bulk of the emission from WASP-43b comes from the 0.1–1 bar pressure range. At these altitudes, both codes closely follow the hemisphere-averaged profile of the underlying model, and well-fit the noninverted slope of the temperature profile.

Water is the dominant absorber across the MIRI waveband. Its ubiquitous absorption at all wavelengths shapes the emission spectrum. Both retrievals constrain well the water abundance, aided by the NIRSpec prior (Gaussian for Pyrat Bay; uniform for TauREx).

Methane has its strongest absorption band between 7 and 9 μm. Consequently, the larger methane abundance for the equilibrium case over the quenched case produces a markedly lower emission at these wavelengths (more methane concentration leads to stronger absorption, which leads to higher photosphere at lower temperatures, which leads to lower emission). Both retrievals are able to distinguish between these two cases, producing a precise methane constraint for the equilibrium case: Pyrat Bay obtains a median with 68% HPD (highest posterior density) of $\mathrm{log}({\mathrm{CH}}_{4})=-5.06\pm 0.25$, whereas TauREx obtains $\mathrm{log}({\mathrm{CH}}_{4})=-4.82\pm 0.89$. The lower concentration of methane in the quenched case leads to wider posteriors, and the retrieval is only able to provide an upper limit on the methane concentration.

Carbon monoxide only has a strong band at the shorter edge of the observed spectra (5 μm), and thus its abundance is harder to constrain. In the equilibrium case, both retrievals set an upper limit on the CO abundance. In Figure 13, the upper limits of the credible intervals are nearly two orders of magnitude below the averaged CO abundance at 0.3 bar. However, in the input model, CO decreases rapidly with altitude over the probed pressures—from ∼4 × 10⁻⁴ to ∼4 × 10⁻⁵ between 1 and 0.1 bar. Since the retrievals assume a constant-with-altitude profile, the retrieval models require a lower CO abundance to produce the same signal of the input model, which might explain the underestimated retrieved CO values. For the quenched case, both CO posterior distributions peak slightly below the averaged input.

Table 2 and Figure 13 compare the retrieved abundances and the associated uncertainties for the two cloud-free cases (equilibrium and quenched chemistry). Both codes (based on different methods) produce qualitatively similar results for each molecule. The posterior credible intervals are consistent when the molecule is well-constrained, whereas they differ by up to two dex when finding upper limits. These results show that both methods are robust. In light of our results, we expect that in a cloud-free case, we will be able to distinguish between the equilibrium and quenched scenarios using MIRI phase curve observations of WASP-43b, thanks to the methane absorption band seen between 7 and 9 μm.

We ran Pyrat Bay and TauREx nightside cloudy retrievals on the 1 μm MnS and MgSiO₃ JWST/MIRI simulated data sets (Figure 12, middle-right and bottom panels). Although both codes retrieved similar best-fit spectra, temperature profiles, and contribution functions, Pyrat Bay's TSC cloud model was more successful in retrieving the input condensate particle size, cloud number density, and the location of the cloud deck for the MgSiO₃ clouds.

To investigate MnS and MgSiO₃ clouds, TauREx ran "free" retrieval scenarios, while Pyrat Bay ran both "free" and "self-consistent" equilibrium-chemistry retrieval. We chose to include the "self-consistent" scenario, as the clouds in the input synthetic models were post-processed with cloud opacities using the cloudless GCM simulations, in the same way as we add clouds in retrieval (see Section 4.4). For the MgSiO₃ clouds, both Pyrat Bay and TauREx also ran nonconstrained and constrained temperature-profile cases within "free" retrieval, as the data led the parameter exploration to nonphysical solutions.

For Pyrat Bay "free" retrieval, we used the same setup as in the cloud-free retrieval case (see Sections 3.6 and 6.1) and retrieved the same temperature and pressure parameters (κ, γ₁, γ₂, α, and β; see Line et al. 2013), with the addition of five new free parameters describing the cloud characteristics (J. Blecic et al. 2020a, in preparation): the cloud extent, Δi, the cloud profile shape, H_c, the condensate mole fraction, ${\mathrm{log}}_{10}(q* )$, the particle-size distribution, ${\mathrm{log}}_{10}({r}_{\mathrm{eff}})$, and the gas number fraction just below the cloud deck, ${\mathrm{log}}_{10}({X}_{c})$. For the "self-consistent" scenario, we retrieved the temperature–pressure parameters together with the aforementioned cloudy parameters, excluding the chemical species parameters. The abundances of the chemical species were produced using the initial implementation of the RATE code (Cubillos et al. 2019), based on Heng & Tsai five-species solution. We calculated the volume mixing ratio of H₂O, CH₄, CO, CO₂, and C₂H₂ species given a T − P profile. TauREx uses the same setup as in Sections 3.6 and 6.1, in addition to the cloud prescription and free parameters described in Bohren & Huffman (1983), Appendix A. The free parameters of the model are the condensate mole fraction, log₁₀(q*), and the peak of the particle size distribution ${\mathrm{log}}_{10}({r}_{\mathrm{eff}})$. For this paper, we use the particle cloud distribution described in Sharp & Burrows (2007).

Prior to the Pyrat Bay analysis, we ran several different MCMC settings to fully explore the phase space and get the best constraints on the cloud parameters. For all of the tested cases, the posterior histograms of some of the cloud parameters were flat and had no correlations with other parameters, implying that the data are uninformative about them. Thus, we fixed and excluded them from the exploration (Δi was fixed to 80, H_c to 0.75, and ${\mathrm{log}}_{10}({X}_{c})$ to 0.0), allowing only the condensate effective particle size, ${\mathrm{log}}_{10}({r}_{\mathrm{eff}})$, and the condensate mole fraction, ${\mathrm{log}}_{10}(q* )$, to be free. This leaves us with the same parameters as TauREx. For the temperature profile in the MnS case, the solution found with TauREx is degenerated between γ₁ and γ₂. We fix this by disabling the contribution of the second visible opacity (α = 0).

6.2.1. MnS Clouds

In the MnS case (Figure 14), both Pyrat Bay and TauREx well-fit the input spectrum model. While the retrieved temperature profiles are hotter than the hemispheric averaged profile, both profiles converge to similar values at the top of the atmosphere where the contribution functions are located. H₂O abundance and cloud mixing ratio are well-constrained and closely match the input value for both codes. As in the cloud-free retrieval, Pyrat Bay, used a Gaussian prior on the H₂O abundances (log₁₀(H₂O) = −3.52 ± 0.3), while TauREx used a uniform prior (1 × 10⁻¹–1 × 10⁻⁵). Pyrat Bay retrieved a H₂O abundance to log₁₀(H₂O) =−3.48 ± 0.3, i.e., basically the prior value, and TauREx retrieved log₁₀(H₂O) = −4.15 ± 0.71. The retrieved cloud particle sizes have a slightly lower value than the input: Pyrat Bay retrieves ${\mathrm{log}}_{10}({r}_{\mathrm{eff}})=-0.58\pm 0.14\,\mu {\rm{m}}$, while TauREx retrieves ${\mathrm{log}}_{10}({r}_{\mathrm{eff}})=-1.7\pm 0.33\,\mu {\rm{m}}$. Both codes produced no constraints on the CO abundance, while CH₄ abundances could not be constrained in the TauREx run and produced a lower limit (log₁₀(CH₄) > −2.4) in the Pyrat Bay run. The retrieved cloud parameters and calculated chemical abundances of the Pyrat Bay "self-consistent" retrieval nicely match those of the "free" retrieval. A summary of the retrieved abundances is given in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Comparison of the Retrieved Abundances Obtained by TauREx and Pyrat Bay for the MnS Clouds Cases

Case	Molecule	TauREx	Pyrat Bay
	log(H2O)	−4.15 ± 0.71	−3.48 ± 0.30
Free	log(CO)	<−3	unconstrained
	log(CH4)	unconstrained	>−2.4
	log(H2O)	N/A	−3.1
Self-consistent	log(CO)	N/A	−9.5
(0.1 bar)	log(CH4)	N/A	−3.3

Note. There are no uncertainties for the self-consistent runs, which have been performed with Pyrat Bay only.

Download table as:  ASCIITypeset image

Given the relatively similar results both codes obtained in the cloud-free retrieval case (Section 6.1), differences seen here between the retrieved parameters' values within Pyrat Bay and TauREx can be mostly attributed to different cloud parameterization schemes between the two codes: higher complexity of one cloud model compared to the other, different treatment of thermal scattering, larger number of free parameters, more freedom in the shape of the log-normal particle distribution, flexibility in the cloud base location and the shape of the cloud, and possible differences in the resolutions of the output models.

6.2.2. MgSiO₃ Clouds

For the MgSiO₃ clouds, we performed the same "free" retrieval runs as for the MnS clouds. However, in this case, "free" retrievals did not manage to converge to a physically realistic solution for both TauREx and Pyrat Bay. The codes appear to misinterpret the bump at around 8 μm as an emission feature, cut the contribution functions at the very top of the atmosphere, and fit the spectrum with a high CH₄ abundance and an inverted temperature profile (see Figure 15).

Zoom In Zoom Out Reset image size

To overcome the encountered issue, we guided the temperature–pressure model to explore only noninverted solutions in both codes. TauREx was unable to converge to a physical solution, while Pyrat Bay achieved a match between the retrieved temperature and cloud parameters and the input model. In Figure 16, we show the results of this retrieval. To explore only the noninverted temperature profiles, in Pyrat Bay we fixed the γ₂ and α parameters to zero, the γ₁ to values smaller than zero, and let κ and β parameters to be free. The retrieved condensate particle size is 1 μm as the input value, ${\mathrm{log}}_{10}({r}_{\mathrm{eff}})=-4.069\pm 0.509\,\mathrm{cm}$, and the cloud mole fraction upper limit is around 10⁻². H₂O abundance again closely matches the Gaussian prior, log₁₀(H₂O) = −3.58 ± 0.305, and the contribution functions reveal the location of the cloud at 10⁻³ bar, at the same level as the input model. The retrieved best-fit spectrum, median temperature profile, and the condensate particle size (${\mathrm{log}}_{10}({r}_{\mathrm{eff}})=-4.039\pm 0.128$ cm) for the Pyrat Bay "self-consistent" retrieval are almost identical to the constrained "free" retrieval case, with the cloud mole fraction having similar upper limit, but with the contribution functions located at lower pressures (10⁻⁵ bar). The chemical abundances values for all MgSiO₃ cases are gathered in Table 4. Apart from the prior value retrieved for the water, neither Pyrat Bay nor TauREx are able to provide useful constraints on the chemical abundances. Particularly, in the free retrieval case, Pyrat Bay provides a biased measurement of the CH₄ abundance. Despite these shortcomings, Pyrat Bay provides a surprisingly good match of the particle size and a realistic upper limit on the cloud height that cannot be obtained with TauREx.

Zoom In Zoom Out Reset image size

Table 4.  Comparison of the Retrieved Abundances Obtained by TauREx and Pyrat Bay for the MgSiO₃ Clouds Cases

Case	Molecule	TauREx	Pyrat Bay
	log(H2O)	not converged	−3.78 ± 0.29
Free	log(CO)	not converged	<−1.36
	log(CH4)	not converged	>−3.4
	log(H2O)	not converged	−3.58 ± 0.31
Constrained	log(CO)	not converged	unconstrained
T-P	log(CH4)	not converged	unconstrained
	log(H2O)	N/A	−3.0
Self-consistent	log(CO)	N/A	−9.0
(0.1 bar)	log(CH4)	N/A	−3.3

Download table as:  ASCIITypeset image

In conclusion, based on the cloudy retrieval analysis, distinguishing between cloud-free and cloudy atmospheres in JWST/MIRI data could present a challenge without a careful approach. Even for the 1 μm synthetic models with MgSiO₃ clouds, which show a noticeable silicate feature around 10 μm, we saw a degenerate solution with high CH₄ abundance, inverted temperature profile, and no traces of clouds. The challenge will become even higher for particle sizes larger than 10 μm, as the silicate feature becomes even less pronounced (Wakeford & Sing 2015).

To overcome this issue, one must closely examine the physical background of the retrieved model (e.g., temperature profile, spectral features seen in emission rather than in absorption, abundances of the retrieved species) and discard nonphysical solutions. In addition, the complexity of the underlying parameterized cloud models plays a crucial role. More advanced and realistic cloud models that, for example, include thermal scattering, more freedom in the cloud particle distribution, and the location and shape of the cloud, have higher chances to discard false models and allow Bayesian sampling to land in a realistic phase space. Greene et al. (2016) also points out the importance of the synergy between different instruments on board JWST to obtain wider wavelength coverage and avoid abundances or degeneracies in the cloud properties. Combining NIRCam, NIRISS, and MIRI instruments, covering the 0.6–11 μm range, would be our best way forward (Schlawin et al. 2018). In particular, Mai & Line (2019) show that NIRCam wavelength range (λ = 2.5–5 μm) is critical in inferring atmospheric properties (precise compositional constraints are possible due to the presence of CO and CO₂ features at these wavelengths), while NIRISS + MIRI (λ = 0.6–2.5, 5–11 μm) are necessary for constraints on cloud parameters (NIRISS is required to constrain the scattering slope at shorter wavelengths, while MIRI is important due to the existence of mid-IR resonance features). Finally, a 2.5D approach such as the one of Irwin et al. (2020) would further improve the ability to retrieve the cloud properties.

Overall, our retrieval models can distinguish between a cloudy and a cloudless nightside, but if the clouds are too thick or have a single-scattering albedo in the infrared that is too large, such as in the MgSiO₃ case, our state-of-the art retrieval methods have difficulty analyzing the complexity of the planned JWST/MIRI observations. If not carefully guided, they can produce biased detection of a high CH₄ abundance, regardless of a well-known specific absorption feature at 10 μm. More work is needed to leverage the planned observations and quantify the nightside abundances in the presence of nightside clouds.