© K. L. Chubb and M. Min 2022

1 Introduction

The spectra of a star that hosts a transiting planet can offer a wealth of information on that planet’s atmosphere, both via transmission spectroscopy during primary transit, and by comparing emission spectra before and during secondary transit. Information about different regions of an exoplanet’s atmosphere can be deduced by analysing these complementary spectra. Various ‘1D’ assumptions have been typically employed to model and retrieve parameters from such spectral observations of exo-planet atmospheres. 1D transmission retrievals, for example, typically assume an isothermal pressure-temperature structure, which is the same at both morning and evening terminator regions of the planet. 1D emission retrievals will often make assumptions such as the pressure-temperature structure being modelled the same across the entire dayside of the planet, and that molecular abundances are constant throughout the atmosphere. These can often be reasonable simplifications, particularly when analysing transmission spectra from a narrow wavelength region that comes from a layer of the atmosphere where the pressure-temperature structure is close to isothermal. The retrieval of emission spectra, which are more sensitive to the pressure-temperature profile than transmission spectra, requires a non-isothermal structure, but this will usually be parametrised (see, for example, Hansen 2008; Madhusudhan & Seager 2009; Guillot 2010; Waldmann et al. 2015a; Barstow & Heng 2020). One reason for employing such simplifying assumptions is that the computation of many thousands of forward models in a retrieval necessitates each forward model to be computationally efficient. This is particularly useful when analysing multiple datasets in population studies such as Sing et al. (2016); Tsiaras et al. (2018); Fisher & Heng (2018); Pinhas et al. (2018); Welbanks et al. (2019).

There are a number of works that have explored the biases introduced into atmospheric models and retrievals that assume 1D rather than 3D (or closer to 3D) atmospheres, for example: Rocchetto et al. (2016); Blecic et al. (2017); Feng et al. (2016); Taylor et al. (2020); Pluriel et al. (2020); MacDonald et al. (2020). These works make clear the requirement to develop modelling tools that take higher dimension effects into account (for example, allowing molecular abundances, pressure-temperature profiles, and cloud coverage to vary as a function of longitude, latitude, and altitude). The development of various such tools, which range from 1.5D to 3D effects, has been the focus of studies such as Irwin et al. (2020); Caldas et al. (2019); Changeat & Al-Refaie (2020); Feng et al. (2020); Beltz et al. (2020); Cubillos et al. (2021); Lee et al. (2022); Falco et al. (2022); MacDonald & Lewis (2022). These studies build on early examples of work into developing 3D transmission spectra models, such as Fortney et al. (2010); Burrows et al. (2010); Dobbs-Dixon et al. (2012).

WASP-43b is a very well-studied “hot Jupiter” exo-planet (e.g. Louden & Kreidberg 2018; Gandhi & Madhusudhan 2017; Komacek et al. 2017; Keating & Cowan 2017; Gandhi & Jermyn 2020; Feng et al. 2020; Lustig-Yaeger et al. 2021; Garai et al. 2021). It has a rich variety of available observational data, including an emission phase curve (Stevenson et al. 2014, 2017) and a transmission spectrum that has been found to contain strong evidence of molecular signatures (e.g. Kreidberg et al. 2014; Fisher & Heng 2018; Chubb et al. 2020a). The close-in giant exoplanet is assumed to be tidally locked, so information on atmospheric variability across the planet’s surface can be inferred by utilising the phase curve data to generate spectra at different phases during the planet’s orbit. This makes WASP-43b an ideal candidate for detailed studies of atmospheric circulation models, and for testing retrieval models that extend beyond the traditional 1D formalism, such as the one presented in this work. It has already been the target for such retrievals that go beyond the traditional 1D formalism for retrieving properties of exo-planet atmospheres, for example: Irwin et al. (2020); Feng et al. (2020); Changeat & Al-Refaie (2020); Cubillos et al. (2021). The presence of strong equatorial jets has been suggested by previous studies, such as Kataria et al. (2015), largely due to the strong day-night temperature contrast (Stevenson et al. 2014; Kataria et al. 2015; Gandhi & Madhusudhan 2017; Irwin et al. 2020). Another recent study from Carone et al. (2020) has also investigated the potential presence of deep atmospheric wind jets in WASP-43b, which are linked to retrograde (westwards) airflow to explain an apparent small hotspot shift and a large day-night temperature gradient.

The publicly available transmission and emission phase data for WASP-43b are primarily a result of observations by the Wide Field Camera 3 (WFC3) instrument on the Hubble Space Telescope (HST) and the Infrared Array Camera (IRAC) instrument on board the Spitzer Space Telescope (Spitzer). HST/WFC3 was used to observe three full-orbit phase curves, three primary transits, and two secondary eclipses (proposal ID 13467, PI: Jacob Bean, Bean 2013). Spitzer/IRAC was used to observe three broadband photometric phase curves of WASP-43b (proposal IDs 10169 and 11001, PI: Kevin Stevenson, Stevenson et al. 2017). The light curve fitting of the transmission data was first carried out by Kreidberg et al. (2014) and later, independently, by Tsiaras et al. (2018). We consider the former in this work. The Spitzer IRAC emission phase curve data measured at 3.6 μm and 4.5 μm by Stevenson et al. (2017) were later reduced independently by Mendonça et al. (2018a); Morello et al. (2019); May & Stevenson (2020); Bell et al. (2021). Previous analyses of the WFC3 data by Kreidberg et al. (2014) and Stevenson et al. (2017) have found evidence of H₂O in both transmission and emission, with Stevenson et al. (2017) deducing the presence of CO and/or CO₂ based on the Spitzer data points. Weaver et al. (2019) also find evidence of H₂O in transmission. CH₄ was found to vary in emission with phase by Stevenson et al. (2017), with some caution on the derived abundances demonstrated by Feng et al. (2016). Chubb et al. (2020a) recently found evidence for AlO in the transmission spectra of WASP-43b, based on a systematic search for evidence of all species that have available opacity data with spectral features in the HST/WFC3 wavelength region.

The main aim of the present study is to create a retrieval framework that encapsulates the 3D nature of an exoplanet’s atmosphere. We have developed and implemented such a framework into the atmospheric modelling and Bayesian retrieval package ARCiS (ARtful modelling Code for exoplanet Science) (Min et al. 2020; Ormel & Min 2019). We utilise available observed emission phase curve data and transmission spectra of WASP-43b. We perform a series of retrievals on these spectra in order to test various input assumptions on retrieval outcomes.

This paper is structured as follows. Section 2 outlines the theory behind our 3D retrieval model. This model is applied to WASP-43b using retrieval code ARCiS, as outlined in Sect. 3. A total of eight different retrieval setups are detailed here, the results of which are outlined in Sect. 4 (and corresponding figures in Appendix A and the supplementary information document associated with this work¹). Discussions on various aspects of these results are also included in Sect. 4. Finally, our conclusions and future outlook are presented in Sect. 5.

2 3D Model for Use in Retrievals: Theory

ARCiS is an atmospheric modelling and Bayesian retrieval package. General details can be found in Min et al. (2020) and Ormel & Min (2019). ARCiS can handle simultaneous retrieval of transmission and emission spectra, including phase-dependent spectra. A benchmark of the emission spectra, based on that of Baudino et al. (2017), can be found in Appendix B. We have recently updated ARCiS to allow for a 3D characterisation of exoplanet atmospheres during this simultaneous retrieval, as outlined in this work. There is no limit to the amount of observational data that can be included in such a retrieval, although the computational time does scale with the number of observed data points. Although the 3D setup implemented in ARCiS is used in this work, we also provide a general routine for the 3D setup used². This is with the intention that it can be implemented into other general retrieval codes.

In the 3D retrieval setup, we have to model how the incoming stellar flux on the dayside of the planet is spread throughout the atmosphere. Here we model this using diffusion and horizontal winds. In general the normalised incoming stellar flux at the top of the atmosphere as a function of longitude Λ and latitude Φ is $${\beta }_{\star }=\cos \Lambda cos\Phi$$(1)

on the dayside, and $${\beta }_{\star }=0$$(2)

on the nightside. For a static atmosphere (with no diffusion or winds), this allows us to compute the temperature structure at each location. We know this is unphysical. The atmosphere is not static and part of the energy incoming on the dayside is transported to the nightside. We therefore define a new parameter β, which specifies how the energy incident from the host star, β_⋆, is spread throughout the atmosphere. Again, β is a function of longitude Λ and latitude Φ. We model this spread of energy by a diffusion equation on β: $$\nabla \cdot v\beta -K_{\Lambda \Phi }{\nabla }^2\beta =S$$(3)

Here, β now specifies how the energy incident from the central star is spread through the atmosphere, and is found as a function of longitude Λ and latitude Φ. For simplicity we convert to a unit sphere and dimensionless velocity and diffusion coefficients. The diffusion is computed over this unit sphere assuming the velocity vector υ is in the longitudinal direction. It is known from global circulation models (GCMs) that these longitudinal winds are strong around the equator and much weaker at the poles (see e.g. Kataria et al. 2015; Carone et al. 2020). Therefore, we also use a coefficient n_ϕ to specify the latitudinal extent of the winds: $v=V_{\Lambda }{\cos }^{n_{\varphi }}\Phi$. In spherical coordinates, Eq. (3) is then given by: $$\frac{V_{\Lambda }{\cos }^{n_{\varphi }}\Phi }{\cos \Phi }\frac{\partial \beta }{\partial \Lambda }-\frac{K_{\Lambda \Phi }}{\cos \Phi }\frac{\partial }{\partial \Phi }\left(\cos \Phi \frac{\partial \beta }{\partial \Phi }\right)-\frac{K_{\Lambda \Phi }}{{\cos }^2\Phi }\frac{{\partial }^2\beta }{\partial {\Lambda }^2}=S,$$(4)

where the source term is given by $$S=f_{\text{day}}{\beta }_{\star }$$(5)

and where f_day is the efficiency with which the planet absorbs the stellar radiation (i.e. in principle f_day = 1 − A_b, where A_b is the planetary bond albedo).

To aid the interpretation of the variables in our equations, we substitute a few parameters with those that are better tailored to the interpretation of observations. First we introduce a new parameter $f_v=\frac{V_{\Lambda }}{K_{\Lambda \Phi }}$ as the ratio of the equatorial wind speed over the horizontal diffusion parameter. Equation (4) then becomes: $$f_v\frac{{\cos }^{n_{\varphi }}\Phi }{\cos \Phi }\frac{\partial \beta }{\partial \Lambda }-\frac{1}{\cos \Phi }\frac{\partial }{\partial \Phi }\left(\cos \Phi \frac{\partial \beta }{\partial \Phi }\right)-\frac{1}{{\cos }^2\Phi }\frac{{\partial }^2\beta }{\partial {\Lambda }^2}=\frac{S}{K_{\Lambda \Phi }}.$$(6)

Finally, we define a parameter f_red, which is the integrated value of β on the nightside over the integrated value of β on the dayside. This gives a measure of the heat distribution across the atmosphere. For a static atmosphere, we have f_red = 0, while for a fully mixed atmosphere we have f_red = 1. We iteratively adjust the value of K_ΛΦ to reach the given value of f_red. Numerically this is rather easy since with all other parameters fixed f_red is always changing in the same direction as K_ΛΦ.

To summarise, our base 3D setup (before including other free parameters related to molecular abundance or temperature profile) consists of only four parameters: (f_red, f_v, f_day, and n_Φ). For these parameters we iteratively solve for K_ΛΦ using Eq. (6) and solve for β at each location on the planet’s globe. Next, at each Λ, Φ coordinate, we compute a ID temperature structure. We allow some of the parameters in our model (or their logarithm) to vary linearly with β, while other parameters are kept fixed over the entire planet. This has the advantage that for each value of β, we only have to compute a single ID structure and position it at different latitudinal and longitudinal locations. Example structures of β as a function of longitude and latitude are shown in Fig. 1.

In the retrievals presented in this paper, we use ray tracing to compute the outgoing flux from the parametrised pressure-temperature structure. Ray tracing is done on a 3D structure, with no assumptions on the geometry. This means that effects such as limb darkening are automatically taken into account. In addition, effects of the planet becoming aspherical because of the larger vertical extend of hotter regions, are also automatically computed. For both scattering and temperature computations, where non-Monte Carlo, ray-tracing techniques are used, we achieve good agreement with literature benchmarks, as can be seen in Appendix B.

Fig. 1 Example structures of β as a function of latitude and longitude. Top row. two examples for a pure diffusive case (no significant longitudinal wind), while the lower panels show two examples for significant longitudinal winds. The red vertical lines indicate the position of the longitude 0, while the blue vertical lines indicate the location of the terminator. For all four cases, we take fday = 1, nΦ = 3, where fday and nΦ are the efficiency with which the planet absorbs incoming stellar radiation, and the coefficient specifying the latitudinal extent of the longitudinal winds, respectively.

2.1 Pressure-Temperature Profiles

In principle, knowing β (where β specifies how the energy incident from the central star is spread through the atmosphere; see Eq. (3)) at each location allows the pressure temperature structure to be computed from radiative transfer. However, for a retrieval model this is very computationally demanding. Therefore, here we use the pressure-temperature profile parametrisa-tion of Guillot (2010) to model the atmospheric structure at a given longitude and latitude. It is possible to have a completely free temperature profile with our setup; however, this introduces many extra free parameters into the retrieval. We instead opt here for the Guillot profile, which requires fewer free parameters but is still relatively flexible. In the future, when more data is available from missions such as the James Webb Space Telescope (JWST), we may be able to constrain the free pressure-temperature profile also. For the Guillot parametrisation, the temperature T at a given atmospheric layer, defined by the optical depth τ, is given by: $$T^4=\frac{3T_{\mathrm{int}}^4}{4}\left[\frac{2}{3}+\tau \right]+\frac{3T_{\text{irr}}^4}{4}\beta \left[\frac{2}{3}+\frac{1}{\gamma \sqrt{3}}+\left(\frac{\gamma }{\sqrt{3}}-\frac{1}{\gamma \sqrt{3}}\right){\text{e}}^{-\gamma \tau \sqrt{3}}\right]$$(7)

Here, $\gamma =\frac{{\kappa }_{\text{V}}}{{\kappa }_{\text{IR}}}$ (i.e. the ratio of mean optical to infra-red opacities). τ is the optical depth at IR wavelengths, related to pressure P by gravity g and the mean infra-red opacity κ_ir: $$\tau =\frac{{\kappa }_{\text{IR}}P}{g}$$(8)

These equations were derived for a static atmosphere where the stellar flux received by the planet at the top of the atmosphere is given by $\frac{\sigma T_{\text{irr}}^4{\beta }_{\star }}{4\pi }$ The irradiation from the host star is computed using $T_{\text{irr}}=T_{\star }\sqrt{R_{\star }/d_p}$, where T_⋆ is the stellar temperature, R_⋆ the stellar radius, and d_p the distance between the planet and host star. The internal heat flux which is emitted by the planet is given by $\sigma T_{\mathrm{int}}^4$, where σ is the Stefan-Boltzmann constant. We apply this equation also for our non-static atmosphere by using β instead of β_⋆ at each longitude and latitude (see the text surrounding Eq. (4)). We therefore compute a pressure-temperature profile for each of these grid points using our Λ– and Φ– dependent value of β together with Eq. (7). Using this setup, the full 3D temperature structure is determined by seven parameters: (f_red, f_v, f_day, n_Φ, T_int, κ_ir, and, γ).

2.2 Retrieval Setup

We are using the Bayesian retrieval package ARCiS (Min et al. 2020; Ormel & Min 2019), which utilises the MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009, 2013) Monte Carlo nested sampling algorithm, to sample the specified parameter space for the region of maximum likelihood. We use 500 live points, an evidence tolerance of 0.5, and a sampling efficiency of 0.4 in our retrievals. The parameters and associated priors relevant for the present study are given in Table 1. A benchmark of the transmission spectra computed by ARCiS can be found in Ormel & Min (2019), and a benchmark for the emission spectra, based on that of Baudino et al. (2017), can be found in Appendix B.

Above, we described the 3D setup of the forward model. Some of the parameters in Table 1 are global, while some others (or their logarithmic values) are varied linearly with β (the parameter which determines how incoming stellar radiation is spread throughout the atmosphere; see Eq. (4)). For these parameters which are found as a function of β, we have a night value, occurring at the lowest value of β, and a day value, corresponding to the highest value of β. The values in between these two extremes are varied linearly. For example, in the free molecular abundance retrievals, two free parameters exist (day and night abundances) and the longitudinal abundances are computed by interpolating with respect to β. In the equilibrium case, the metallicity, [Z], and carbon-to-oxygen ratio, C/O, are conserved globally, but the abundances change longitudinally with the local temperature. The cloud parameters (we use a simple cloud layer and allow the pressure-top level, cloud opacity, and cloud albedo to vary; see Sect. 3.1), are also retrieved using two free parameters. Similar to the free molecular abundances, these are related to the minimum and maximum values of β (night and day), with the values in between computed by interpolating linearly with respect to β. All other retrieval parameters (see the section titled Parameters included in all setups in Table 1) are retrieved globally (i.e. there is only one free parameter). The pressure-temperature profiles, computed using Eq. (7), depend on the defined parameters used to describe the Guillot profile, and on β, which leads to different temperature structures as a function of longitude and latitude.

We compute the phase-averaged temperature-profiles for each specified phase, which are better linked to the observed data points and can be directly compared to the results of other studies and the outputs of GCMs. The phase-averaged outputs are centred about the latitudinal line on the atmosphere with the same value in degrees. Outputs given as a function of longitude give the exact value at that point on the atmosphere, whereas outputs given as a function of phase give the average value integrated over the visible disk of the planet, centred on the given latitude (see Sect. 4 for all outputs). The parameters retrieved as free parameters in the 3D ARCiS retrieval are given in Table 1, with a short description of each. Further details on the molecules included in the equilibrium chemistry and free molecular retrieval setups are given in Sects. 3.3 and 3.4.

3 Case Study: 3D Retrievals of WASP-43b

WASP-43b was discovered by Hellier et al. (2011) around an active K7V star, with deduced planetary parameters from radial velocity measurements and transit observations of 2.034 ± 0.052 M_j and 1.036 ± 0.019 R_j (Gillon et al. 2012). It is assumed to be tidally locked. Therefore, the emission phase observations from Stevenson et al. (2014, 2017) can be used to essentially map the atmosphere across the planetary surface. ARCiS allows a combination of emission, phase, and transmission spectra to be included simultaneously in the retrieval. The data we included in this case study of WASP-43b is summarised in Table 2. We follow the convention that the morning terminator corresponds to 270° or 0.75 phase, and the evening terminator to 90° or 0.25 phase. We choose not to use Spitzer transmission data points in this study, as there may be potential issues with combining data from different instruments (see, for example, Yip et al. 2020). Nevertheless, as will be shown by our results, and as pointed out in previous studies, the Spitzer emission data are very important in order to distinguish between different interpretations of the atmospheric composition and dynamics. We thus include the Spitzer data from the emission phase curves along with the HST/WFC3 emission phase curve data. Having observations with more extensive wavelength coverage taken from a single instrument, such as JWST, will greatly aid future retrievals.

As mentioned in Sect. 2.2, we are using Bayesian retrieval package ARCiS (Min et al. 2020; Ormel & Min 2019), which uses the MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009, 2013) algorithm to sample the specified parameter space for the region of maximum likelihood. In order to fully explore the effects of including different physical parameters in a retrieval (in particular, including parametrised clouds, or allowing free vs. equilibrium chemistry), and also the effects of using different analyses of the Spitzer phase curve emission data, we performed a series of various setups of retrievals, as described in Table 3. The physical parameters of the WASP-43 system that we used in all retrieval setups are given in Table 4. Below, we explain briefly the differences between these eight retrieval setups. We included isotropic scattering from stellar and planetary flux in all retrieval setups.

3.1 Clouds VS. No Clouds

Clouds are predicted to be present on WASP-43b, particularly on the nightside (Helling et al. 2020, 2021). We therefore performed a set of retrievals with clouds included in a parametrised way. We used a simple cloud setup, for which we introduced 3 parameters to describe a cloud or haze layer; P_cloud (cloud-top pressure, ranging from 10⁻³ bar for high-altitude clouds or hazes to 10³ bar for deep-atmosphere clouds), κ_cloud (a measure of cloud opacity from 10⁻⁷ cm² g⁻¹ for transparent clouds to 10² cm² g⁻¹ for very opaque clouds), and A_cloud (cloud albedo, ranging from 0 for completely unreflective clouds to 1 for fully reflective clouds). As outlined in Sect. 2, this cloud setup introduces six more free parameters in the retrieval setups that include clouds than in those that don’t. This is because we introduced values of each of the three parameters that correspond to both minimum and maximum values of β (i.e. the hottest dayside and the coolest nightside parts of the planet), and then allowed the parameters to vary linearly in between. These six extra free parameters inevitably increased the computation time for our retrievals. One way to make these computations more feasible was to explore the settings and modes available in the Multinest algorithm used in ARCiS to sample parameter space. We thus used the constant efficiency mode of Multinest, which noticeably improved the efficiency. A short description of this is given below in Sect. 3.2.

3.2 Computational Demands: Constant Efficiency Mode for Cloud Setups

The constant efficiency mode of Multinest can vastly speed up retrievals computations (Feroz et al. 2009, 2013). We tested this for our 3D retrieval setup and found significant speed improvements. To put into context the increased computational effort of the 3D retrieval setup in comparison to a standard 1D retrieval using ARCiS: a standard 1D retrieval (of combined emission and transmission spectra) computes approximately 5 models per second (on a MacBook, 2 GHz Quad-Core Intel Core i5, 16 GB memory), whereas the 3D model computes approximately 1 model per second. A retrieval of just transmission spectra using a 1d setup, with no scattering included, can compute around 33 models per second. Feroz et al. (2009) note that although Multinest’s constant efficiency mode could potentially miss sampling some regions of parameter space, it may nevertheless produce reasonably accurate posterior distributions for parameter estimation purposes. The evidence estimates in this mode, however, should not be relied upon. We therefore do not use these values as a way to judge between our different retrieval setups, such as in Trotta (2008). In order to test whether MULTINEST was missing degenerate solutions in this mode, we ran a test by setting up the same retrieval as setup D (with clouds, equilibrium chemistry, and Stevenson et al. 2017 Spitzer analysis; see Table 3), but using full efficiency instead of constant efficiency. We refer to this as retrieval setup I. The number of model calls is over ten times larger for the full efficiency mode in comparison to the equivalent constant efficiency mode (comparing retrieval setup D at 60 000 and retrieval setup I at 700000). This had a big impact on computational time. We did find some differences between some of the retrieved parameters of the two retrievals, but many were within σ error bounds (see the posterior plots contained in the supplementary information document of this work. This document also contains a comparison of results from setups D and I). The main difference that we note is that the parameters retrieved in the constant efficiency case appear to be better constrained than those in the full-efficiency case. This has a noticeable effect on the pressure-temperature profiles (see Sect. 4.4). We conclude this is artificial, and it should be kept in mind when viewing our results that we do not expect them to be as well constrained as they often appear in the cases where clouds are included in the setups.

3.3 Equilibrium Chemistry

In the equilibrium chemistry retrievals, the molecular abundances are computed assuming equilibrium chemistry using the GGchem code (Woitke et al. 2018), which is fully integrated into ARCiS. Two parameters are used to parametrise the elemental abundances: carbon-to-oxygen ratio, C/O, and metallicity, [Z]. Opacities were all computed in ARCiS format as part of the Exo-MolOP database (Chubb et al. 2021), and are publicaly available on the ExoMol (Tennyson et al. 2020) website³. Opacity formats for other retrieval codes are also available from ExoMolOP: TauREx3 (Waldmann et al. 2015b,a; Al-Refaie et al. 2021), NEMESIS (Irwin et al. 2008), and petitRADTRANS (Mollière et al. 2019). We include the following set of species, given along with the source of the line-list used for the corresponding opacities: H₂O (Polyansky et al. 2018), CO₂ (Yurchenko et al. 2020), CH₄ (Yurchenko et al. 2017), CO (Li et al. 2015), AlO (Patrascu et al. 2015), NH₃ (Coles et al. 2019), H⁻ (John 1988), SiO (Barton et al. 2013), FeH (Wende et al. 2010), TiO (McKemmish et al. 2019), VO (McKemmish et al. 2016), SH (Gorman et al. 2019), H₂S (Azzam et al. 2016), C₂H₂ (Chubb et al. 2020b), CN (Brooke et al. 2014), CP (Ram et al. 2014), H₂CO (Al-Refaie et al. 2015), H₂O₂ (Al-Refaie et al. 2016), K (Allard et al. 2016, 2019; Kramida et al. 2013), Na (Allard et al. 2019; Kramida et al. 2013), OH (Brooke et al. 2016; Yousefi et al. 2018; Bernath 2020), PH₃ (Sousa-Silva et al. 2015), and ScH (Lodi et al. 2015). There are other species available from ExoMolOP, such as Yorke et al. (2014); Yurchenko et al. (2018); Mant et al. (2018), but we limit the species included in our chemistry retrievals to those listed above. This is because combining the k-tables of multiple species has a significant effect on computing time for emission spectra in particular. The number of molecular species included is less of an issue for the constrained chemistry case than the free molecular retrieval, where in the latter the number of free parameters also has a significant effect on computational time. We therefore made a selection of which species to include based on those that are expected to be there, and those with notable absorption features in the wavelength region able to be constrained by our data. We would consider including additional species in the future if we found that they could affect our resulting C/O.

For a C/O above solar, we adjust the oxygen abundance and keep carbon at solar abundances, while for a C/O below solar, we adjust the carbon abundance and keep oxygen at solar abundances. This simulates the removal of these elements during the formation or evolution of the planet. After the relative abundances of all elements heavier than He are computed using this C/O, the H and He abundances are then scaled using the metallicity parameter [Z]: $$\left[Z\right]={\log }_{10}\left[\left(\frac{{\displaystyle {\sum }_{i\ne \text{H,He}}{X}_i}}{{\displaystyle {\sum }_{i=\text{H,He}}{X}_i}}\right)\cdot {\left(\frac{{\displaystyle {\sum }_{i\ne \text{H,He}}{X}_i}}{{\displaystyle {\sum }_{i=\text{H,He}}{X}_i}}\right)}_{\text{Solar}}\right].$$(9)

Here, X_i is the number density of element i.

3.4 Free Molecular Chemistry

For the free molecular retrievals, we include the base set of molecules (H₂O, CO, CO₂, CH₄, and NH₃) which are typically included in previous studies of WASP-43b (for example, Stevenson et al. 2017; Irwin et al. 2020; Kreidberg et al. 2014). We also include AlO due to the studies of Chubb et al. (2020a), who found some potential evidence for AlO in the Kreidberg et al. (2014) transmission spectra of WASP-43b. As with the equilibrium chemistry setups, all opacity data were formatted for ARCiS as part of the publicly available ExoMolOP database (Chubb et al. 2021). The species used in the free molecular retrieval setup are as follows: H₂O, CH₄, CO, CO₂, AlO, and NH₃. All line lists used are the same as those listed in Sect. 3.3.

3.5 Using Different Analyses of the Spitzer Phase Curve

Although we do not consider them all here, there have been multiple reductions of the same observed Spitzer phase curves of WASP-43b; for both 3.6 and 4.5 μm (programmes 10169 and 11001, PI: Kevin Stevenson). These include the two which are used in this study: Stevenson et al. (2017) and Mendonça et al. (2018a). Subsequent reductions have been published by Morello et al. (2019), May & Stevenson (2020), and Bell et al. (2021).

There are some disagreements between these different analyses. For example, Morello et al. (2019) found higher nightside temperatures and smaller hotspot offsets than Stevenson et al. (2017). The Spitzer points from the latter are used in a number of retrieval studies of WASP-43b, including Irwin et al. (2020) and Changeat et al. (2021). We therefore chose to use them in this work in order to aid comparisons. We also use those from Mendonça et al. (2018a) that represent spectra where the apparent day-night temperature difference is less pronounced, in order to observe the effect on our retrieval results. We do not consider the other reductions of the Spitzer phase curves that are mentioned in the present study, but we note that this choice can have an effect on the retrieval results, as will be illustrated in Sect. 4.

4 Results and Discussion of 3D Retrievals of WASP-43b

As previously mentioned, we performed eight different retrievals using slightly different setups, as detailed in Table 3. This allowed us to explore the parameter space and the effects that different prior assumptions have on retrieval results. Table 5 gives the values of a selection of retrieved parameters found from the different retrieval setups, with parameters all described in Table 1. Further results can be found in Appendix A, as detailed below. We also include all corner plots in the supplementary information of this work. A reminder of the eight retrieval setups is given in the footnote of Table 5. We also include a ninth setup, labelled I, which is the same as setup D but using the full efficiency mode of Multinest (see Sect. 3.2).

Here, we focus on some key findings. In particular, we highlight the retrieved parameters that appear consistent between different setups and could therefore be considered more robust. We also discuss degeneracies and covariance between different parameters. Various figures are presented for these eight different outputs, as detailed in the text of this section. These can be found in the appendix to this paper, with corner plots in a separate supplementary information document.

4.1 C/O and Metallicity

Knowledge of atmospheric C/O and [Z] can have implications for inferring information about planet formation history (see, for example, Miller & Fortney 2011; Mordasini et al. 2016; Öberg & Bergin 2016; Madhusudhan et al. 2017; Shibata et al. 2020, 2022; Khorshid et al. 2021; Cridland et al. 2019; Turrini et al. 2021). It can be seen from Table 5 that both C/O and [Z] are consistently retrieved to be super-solar for all the equilibrium chemistry retrieval setups. Super-solar values are consistent with those found in the full retrieval by Changeat et al. (2021), who find $\text{Z=1}{.81}_{-0.17}^{+0.19}$ and $\text{C}/\text{O}\text{=0}{.68}_{-0.12}^{+0.11}$ in the scenario where they allow a completely free temperature profile. This temperature profile setup is different to our parametrised Guillot profile (see Sect. 2.1). Similar to the present study, Changeat et al. (2021) test various retrieval setups when performing simultaneous phase curve retrievals of WASP-43b observations. While they retrieve super-solar metallicity for all runs, they find the C/O to be more dependent on the exact setup and, in particular, the thermal profile used. Kataria et al. (2015) find that an atmospheric metallicity that is 5 × solar (without TiO/VO) is the best match to the WASP-43b observations. Irwin et al. (2020) use optimal estimation for their 2.5D retrieval of WASP-43b, but recommend a more extensive Bayesian sampling of parameter spaces, such as nested sampling, as is used in the present work. Based on their high retrieved abundances of CO and CH₄, they infer a C/O of ~0.91, which is similar to the values found in our equilibrium chemistry retrievals. They use the Spitzer emission data from Stevenson et al. (2017). Irwin et al. (2020) also point out that any errors in the offset between HST and Spitzer data could lead to erroneous molecular abundances, which could then have an effect on both the C/O and metallicity. In the context of planet formation theory, a high metallicity typically indicates that the planet was formed further out from the star and migrated inwards. However, this is also generally assumed to be linked to a low C/O because of the assumed oxygen-rich nature of the accreted solids that leads to the high metallicity (Khorshid et al. 2021; Turrini et al. 2021). A combination of high metallicity and a high C/O is therefore potentially unexpected, based on predictions from planet formation theory. A high C/O, however, could also be an indication of cloud formation, which typically de-oxidises the gas phase (Helling 2019).

4.2 Retrieved Spectra at Different Phases

Figure 2 shows the observed emission data for different phases plotted alongside the retrieved spectra from ARCiS for retrieval setup A. The shading represents the 1, 2, and 3 σ bounds. On the left of each panel is the corresponding heat distribution for each phase, with the cooler nightside in purple and the dayside hotspot in red. Figures A.1 and Figures A.2 show retrieved best-fit spectra at three different phases (90° (evening), 180° (dayside), 270° (morning)) for each of the eight different retrieval setups, as detailed in each panel.

It can be seen that there is some obvious asymmetry in the terminators from the observed data, with one terminator apparently hotter, and therefore emitting more flux than the other. This is particularly apparent for the setups A-D which use the Spitzer data from Stevenson et al. (2017). Setups E-H which use the Spitzer data from Mendonça et al. (2018a) do not have the same extent of asymmetry. It is apparent that all retrieved spectra fit the observed data reasonably well, with no drastic outliers.

4.3 Retrieved Transmission Spectra

Figure A.3 shows the retrieved and observed transmission spectra for all eight of our retrieval setups, as detailed in each panel. Retrieval setup I (which is the same as setup D, but in the full efficiency mode of Multinest) is also included. Our 3D model, which we use for all retrievals simultaneously, retrieves both the emission phase spectra presented in Sect. 4.2 and the transmission spectra shown in Fig. A.3, with equal weight assigned to all observed data points. The fit is generally reasonable, although it does not show the absorption features associated with AlO as was found in our previous retrieval of the transmission spectra only (Chubb et al. 2020a). The flatter transmission spectra retrieved in setups F and H (equilibrium chemistry, using Spitzer data from Mendonça et al. (2018a), without and with clouds included, respectively) is likely due to the high metallici-ties retrieved for both these cases, as given in Table 5. Such high metallicities are associated with heavier atmospheres, which thus reduce the scale height of the atmosphere and flatten absorption features. The absorption features in the emission spectra are not affected in the same way as they are more controlled by the temperature structure of the atmosphere. It should also be noted that setup H (including clouds) uses the constant efficiency mode of Multinest (see Sect. 3.2), and so the posterior distributions given in the corresponding plots are more constrained than they should be.

4.4 Atmospheric Heat Distribution

Figure A.4 gives a 2D projection (as a function of longitude and latitude) of heat distribution, with contours representing different values of β (where β is a function of longitude and latitude, and gives the heat due to stellar radiation at the top of the atmosphere after diffusion has been taken into account; see Sect. 2). This value of β can be used in Eq. (7), along with the retrieved Guillot parameters, to find the pressure-temperature profile at a given longitude and latitude point of the atmosphere, and at a given altitude (as determined by atmospheric pressure). This is illustrated in Figs. A.5 and Figs. A.6, which give the temperature as a function of longitude and latitude for a cut of the atmosphere at a pressure layer of 0.01bar and 0.12 bar, respectively. Some comparison can be made to outputs from GCMs of WASP-43b such as Kataria et al. (2015), Mendonça et al. (2018b), and Carone et al. (2020). See for example Figs. 2d,h of Mendonça et al. (2018b), which show horizontal temperature maps at a pressure-layer cut of 0.01bar for atmospheres which do include and do not include a thermal inversion in the upper atmosphere, respectively. It should be noted that the longitude axis is from 0 to 360^°, whereas ours is from –180 to +180°. 0.01 bar is the pressure layer where their jet reaches maximum strength. A comparison can be made to our Fig. A.5, which gives our retrieved temperature maps at a pressure layer of 0.01 bar. For the case where a thermal inversion is included in the GCM of Mendonça et al. (2018b) (labelled G = 2.0), their model appears qualitatively to match best to our retrieval setup F (equilibrium chemistry, no clouds, and Spitzer data from Mendonça et al. 2018b). For the GCM of Mendonça et al. (2018b) without a thermal inversion (labelled G = 0.5), the temperature map appears most similar to the results of our retrieval setup G (free chemistry, with clouds, and Spitzer data from Mendonça et al. 2018b; see Fig. A.5).

Figure 3 of Kataria et al. (2015) gives the longitude-latitude temperature cut at various pressure levels for l× solar metallic-ity and for 5× solar in their Fig. 7. The temperature cuts at a pressure layer of 10.63 mbar (≈0.01 bar) can be compared to our Fig. A.5, which gives temperature maps at a pressure cut of 0.01 bar. The results from our retrieval setups B and D (i.e. equilibrium chemistry setups with Spitzer data from Stevenson et al. 2017, without and with clouds, respectively) appear most similar to both the l× solar and 5× solar GCM models of Kataria et al. (2015). Kataria et al. (2015) assume local equilibrium chemistry (accounting for condensate rainout, the same as in the present work) in their models, and ignore clouds and hazes. Their models could only be compared to HST/WFC3 data from Stevenson et al. (2014) at the time they were published as the Spitzer data of Stevenson et al. (2017) had not been observed. Our results indicate that their GCM models offer a good match to our retrieval results using Stevenson et al. (2017) Spitzer data and assuming equilibrium chemistry.

Figure 2Ib of Carone et al. (2020) gives the longitude-latitude temperature map as a result of their GCM at a cut of 0.12 bar, which can be compared to our retrieved temperature maps cut at the same pressure in Fig. A.6. The results of our retrieval setup C (i.e. free chemistry, with clouds, and using Spitzer data from Stevenson et al. 2017), match the GCM of Carone et al. (2020) most closely at this pressure layer. Beneath this pressure layer, the GCM models of Carone et al. (2020) predict retrograde equatorial wind jets. The deep circulation framework of all these GCMs are far more complex than the wind parametrisation used in this work, and the models of Carone et al. (2020) show potential retrograde equatorial wind flow due to deep atmospheric wind jets in WASP-43b. Other variabilities and dynamical processes such as storms are expected to have impacts on the atmosphere of hot Jupiter exo-planets (see e.g. Cho et al. 2003, Cho et al. 2021; Skinner & Cho 2021). Our retrieved spectra generally fit very well to the observed data, demonstrating that the retrieval model outlined in this work is well placed to deal with observations of WASP-43b using HST/WFC3 and Spitzer/IRAC. We expect it to also be well placed to deal with observations from the recently launched JWST, which will become available in the near future. This will be explored in future work. This, as well as the fact that we show very reasonable (and in some cases very good) matches to the GCM model outputs of Kataria et al. (2015), Mendonça et al. (2018b), and Carone et al. (2020), is promising and gives weight to the ability of the GCMs to predict atmospheric conditions on exoplanets such as WASP-43b.

We note that the coefficient n_ϕ which we use to specify the latitudinal extent of the winds is hard to constrain, because the observed data we are retrieving from is averaged across all latitudes.

One of the parameters we retrieve, f_red, gives an indication of how heat is distributed from the day to the nightside of the planet (see Table 1; f_red can vary between 0 (no heat redistributed from day to night) to 1 (all heat redistributed from day to night)). Stevenson et al. (2017) derive a related parameter, F, to represent this heat redistribution: $$F=\frac{1}{2}(1+f_{red}),$$(10)

with F varying between 0.5 (low heat redistribution) and 1 (high heat redistribution). Stevenson et al. (2017) find a value of $F={0.501}_{-0.001}^{+0.005}$ by integrating their retrieved model spectra. This corresponds to a value of $f_{\text{red}}={0.002}_{-0.002}^{+0.010}$, which can be compared to the values of f_red found in our eight different retrieval setups, as can be found in Table 5. In general our retrieved values of F agree very well, particularly for the setups which use the Spitzer data from Stevenson et al. (2017), as may be expected. See Sect. 4.5 for further discussion on the heat distribution parameter f_red and potential degeneracies with other parameters.

Keating et al. (2019) find a heat recirculation factor (which varies from 0 for no recirculation to 1 for perfect heat recir-culation) of 0.27 ± 0.07 for WASP-43b. They infer a dayside temperature of 1664 ± 69 K and a relatively hot nightside temperature of 984 ± 67 K; this hotter nightside is more consistent with the results we get from the retrieval setups which use the Mendonça et al. (2018a) Spitzer data. Keating et al. (2019) also use the Mendonça et al. (2018a) analysis of the Spitzer phase curve in their study.

Fig. 2 Emission spectra retrieved at different phases plotted with the observed data for each phase (as given in the bottom left of each panel) for retrieval setup A (no clouds, free chemistry, and Stevenson et al. 2017 Spitzer analysis). The left of each panel shows the heat distribution over the surface of the planet, as computed by ARCiS during the retrieval.

4.5 Influence of Spitzer Phase Curve Data Reduction

As mentioned in Sect. 3.5, there are various different reductions of the Spitzer phase curves available in the literature. We find the choice to influence our interpretations of the data through our 3D retrievals. This can be seen in Table 5 and as described in Sect. 4.4.

In general, f_red, the flux redistribution factor, is very low or zero for the retrievals using the Spitzer data from Stevenson et al. (2017). It is typically slightly higher for the retrievals that instead use the Mendonça et al. (2018a) Spitzer data. There is a relation here with f_v, with a higher value of f_v representing a scenario with a higher ratio of wind speed to diffusion (i.e. a greater hotspot shift). It appears that using the Stevenson et al. (2017) Spitzer data results in a greater hotspot shift (which represents a greater asymmetry between the East and West of the planet), but overall very little heat distributed across from the dayside to the nightside. Using the Mendonça et al. (2018a) Spitzer data, however, generally results in a very small hotspot shift, but a larger amount of heat distributed across from the dayside to the nightside. This is evident when looking at the spectra (e.g. see Fig. A.1), where the difference in flux is not so extreme between night and day when using the Mendonça et al. (2018a) Spitzer data in comparison to the Stevenson et al. (2017) Spitzer data. The outlier here is the retrieval setup where Spitzer data is from Mendonça et al. (2018a), equilibrium chemistry is enforced, and no clouds are included (setup F). It is unclear exactly why this is, but this hints that there are degeneracies in the interpretation of the observed data used in this study with our model. We note that due to data availability, we do not include Spitzer data at 90° (i.e centred on the evening terminator) for the Mendonça retrievals, and only include data from HST/WFC3 at this phase (see Table 2). However, due to the way the 3D retrieval is setup, and the fact that we use Spitzer data for phases either side of 90^° (67.5° and 112.5°), we do not expect this to influence our findings significantly.

4.6 Hotspot Shift

All panels of Figs. A.4-Figs. A.6 show a hotspot shifted slightly away from the middle of the dayside. The value of this hotspot shift for the different retrieval setups is given in Table 6, which can be compared to values found in previous studies, which are also given in Table 6. The caveat is that these hotspot shifts were all derived based on the 4.5 μm Spitzer data only, which could give limited information. The values determined in the present work are based on our full 3D retrievals, which include effects due to the geometry of the planet. Nevertheless, a comparison is still of interest. In general there is a greater shift in hotspot for those retrievals that use Spitzer data from Stevenson et al. (2017) as opposed to Spitzer data from Mendonça et al. (2018a). This could be linked to the higher asymmetry between the morning and evening terminator regions in the spectra of Stevenson et al. (2017, see Fig. A.1). In general, our values broadly agree with those found in the other studies noted in Table 6. However, we note that those studies do only consider the 4.5 μm Spitzer data to infer the hotspot shifts given in Table 6.

4.7 Variation of Molecular Abundances Across the Atmosphere

Figures A.7-Figures A.12 give the variation of molecular abundances as a function of phase for each of the eight different retrieval setups and six different molecules that are included in the free molecular retrievals. These are the volume mixing ratio (VMR) abundances which are averaged over the visible part of the disk, centred at each of the given phases, and taken at a pressure-layer cut of 0.18 bar at the equator. These phase-averaged plots allow us to compare the molecular abundances found in the free molecular retrieval case with those found in the corresponding equilibrium chemistry case; see Figs. A.7-Figs. A.12. The same data used as presented for these six species are available for each species included in the equilibrium chemistry setup (see Sect. 3.3), but we only present here those which can be compared to the six species which are included in the free molecular retrieval setups (i.e. AlO, CH₄, CO, CO₂, H₂O, and NH₃). The plots of Figs. A.8-Figs. A.11 (which correspond to CH₄, CO, CO₂, and H₂O, respectively) can be compared to Fig. 7 of Stevenson et al. (2017). We generally find a similar pattern for these four species that we can compare to Stevenson et al. (2017), as described below. We also discuss briefly the similarities between our eight different retrievals, and highlight where the same trend is seen in all or most retrieval setups.

4.7.1 AlO

As can be seen in Fig. A.7, AlO is predicted by all the equilibrium chemistry retrievals to be most abundant on the dayside, where the atmosphere is hottest. Although we find this same trend in two of the free chemistry retrievals, we also find it to be unconstrained or not very abundant in the other two free chemistry retrievals. A potential indication for AlO was found in WASP-43b’s transmission spectra by Chubb et al. (2020a). We note that if there really is AlO present at one or more of the terminators, then we would expect to also see it on the dayside of the planet. It is possible that the emission data is not precise enough to allow for a strong determination of its presence. In any case, observing a wider wavelength region would really help to confirm whether this species is present in the atmosphere of WASP-43b, and if it is present, to place tighter constraints on its abundance.

4.7.2 CH₄

We find a very consistent trend of CH₄ abundance as a function of phase, particularly in the retrievals cases which use the Stevenson et al. (2017) Spitzer data (see Fig. A.8). We find CH₄ to vary between ~10^-4 on the dayside and ~10⁻³ on the nightside. This can be compared to Fig. 7 of Stevenson et al. (2017), where a higher abundance is also found on the nightside than on the dayside, but with different overall abundances. The abundance trend as a function of phase for the four retrievals that use the Stevenson et al. (2017) Spitzer data are particularly consistent.

4.7.3 CO and CO₂

We retrieve relatively high abundances of both CO and CO₂ (see Figs. A.9 and Figs. A.10). Stevenson et al. (2017) combine these two species together. They also find a relatively consistent abundance across day and night, at a VMR of around 10⁻⁴. This is not too dissimilar from our findings. In particular, our equilibrium chemistry retrievals find an extremely high abundance of CO across all phases. As pointed out by Irwin et al. (2020), any errors in the offset between HST and Spitzer data could lead to erroneous molecular abundances; in particular of CO. Although there is likely to be some degeneracy between these species, which makes it potentially hard to distinguish between them with our given dataset, we choose to keep them separate when performing our retrievals. Both species have a very strong absorption feature at 4.5 μm, but CO₂ also has some weaker features in the HST/WFC3 region.

4.7.4 H₂O

We find the same trend with all of our retrieval setups (with one exception; see Fig. A.11): a higher abundance of H₂O on the dayside than on the nightside. This appears for both the equilibrium and free chemistry retrievals. This trend is also found by the 1D set of retrievals from Stevenson et al. (2017). We note that the apparent tight constraints on H₂O abundance for the free retrievals that include clouds are artificial, due to the constant efficiency mode of Multinest (see Sect. 3.2).

4.7.5 NH₃

Although NH₃ is predicted to be relatively abundant by the equilibrium chemistry retrievals (see Fig. A.12), it remains unconstrained by our free chemistry retrievals. A wider wavelength region is likely needed in order to properly constrain the abundances of this species, if it is present in the observable atmospheric region.

4.8 Temperature-pressure Profiles

As previously mentioned, we also retrieve pressure-temperature profiles for each phase, averaged over the visible disk of the planet. Figure A.13 gives these for the dayside (180°), nightside (315^°), morning (270^°) and evening (90^°) terminators, with each panel giving the temperature profiles retrieved as a result of each of our different retrieval setups.

We use our 3D retrievals to find the temperature structure at each of the points specified on the grid of 36 longitude and 18 latitude points. We then average for each specified phase in order to facilitate comparisons with other works. In Fig. A.13 we present four phases only: 90^° (evening), 180^° (dayside), 270^° (morning), and 315° (nightside). The shading of the plots represents the 1, 2, and 3 σ regions. These plots can also be compared to other works such as those from the 1D retrievals of Stevenson et al. (2017, e.g. see their Fig. 6). We can also compare to GCM outputs, for example see the first panel of Fig. 2 of Helling et al. (2020), which presents the results of 3D GCM models of WASP-43b (Parmentier et al. 2018). The horizontal lines added to the plots given in Fig. A.13 represent the approximate region outside which the temperature cannot be constrained from the observed data. It is important to note that the temperature may appear to be artificially well constrained, in particular for the setups which include clouds, largely due to the mode of Multinest used to make some of the computations feasible (see Sect. 3.2), but likely also because of the shape of the parametrised profile. We also include the results of retrieval setup I as part of Fig. A.13, which is the same as setup D but using the full efficiency as opposed to the constant efficiency mode of Multinest (see Sect. 3.2). It is very apparent, by comparing these two panels, that the retrieval using the constant efficiency mode (setup D) is too artificially well constrained compared to using the full efficiency mode (setup I).

There is a link here to the retrieved parameter T_int (see Table 1), which in the absence of stellar radiation (i.e. due to internal planet heat) gives the atmospheric temperature at an optical depth of $\tau =\frac{2}{3}$.It can be seen from the posterior distributions contained in the supplementary information that this parameter is well constrained as having a relatively high value of log₁₀ (T_int) = 2.69 ± 0.02, whereas the same parameter in retrieval setup I remains unconstrained. T_int appears in Eq. (7), which determines the temperature at a given atmospheric layer in the Guillot (2010) temperature-pressure parametrisation that we use. A similar trend in terms of how constrained a parameter is is seen for other parameters, including log(γ) and log(κ_IR), which also determine the pressure-temperature structure.

Changeat et al. (2021) point out that the Guillot profile (which we use in this work; see Sect. 2.1) may not represent the nightside atmospheric structure of the planet so well, and so they test for allowing the use of a free profile in their retrievals of WASP-43b. This has also been explored in other studies such as Blecic et al. (2017); Madhusudhan & Seager (2009). Our retrieval setup allows for the use of different pressure-temperature profiles, such as a free temperature profile, but this would come with an increase in computational time. Allowing for different pressure-temperature profiles while performing our 3D analysis is something we will investigate further in the future, particularly for the nightside of the planet. When using a Guillot profile, Changeat et al. (2021) retrieve a non-inverted temperature profile and a chemistry structure similar to that of Feng et al. (2020). However, when allowing for a free temperature profile, as opposed to the parametrised Guillot profile, which we also use in the present study, they find that the temperature structure (at least on the dayside and at hotspot regions), is consistent with a thermal inversion.

4.9 Clouds

Figure A.14 gives the retrieved parameters that represent our cloud setup for each of the four retrievals for which clouds are included, as detailed in each panel (see Sect. 3.1 for details). We retrieve P_cloud (cloud-top pressure), κ_cloud (cloud opacity), and A_cloud (cloud albedo). For each we retrieve the values corresponding to the minimum and maximum values of β, and allow them to vary linearly in between (see Sect. 3.1). In general, most of the cloud parameters in all of the retrievals that include clouds are not well constrained, particularly the albedo. This is perhaps unsurprising, considering the narrow wavelength coverage of the observations we are retrieving, and particularly because we do not include any data in the optical region. We find decent fits to the data without clouds included, which could indicate their inclusion is not necessary to fit the observed data. Nevertheless, clouds are expected throughout the atmosphere (see, for example, Helling et al. 2020, 2021). In particular, the detailed cloud and haze models of WASP-43b computed by Helling et al. (2020) predict relatively transparent haze at very high altitudes on the dayside of WASP-43b, and more opaque clouds on the nightside. It could therefore be that the current observed data does not contain enough information to retrieve information on these clouds predicted to be present. In order to properly explore such a trend in a retrieval, we would not only need data that covers a wider wavelength region, but we would also ideally use more detailed cloud formation models than those used in the present work. There is a cloud formation model in ARCiS (Ormel & Min 2019; Min et al. 2020) that we did not use for the current study in order to make the retrievals more computationally feasible. This model allows the C/O to vary locally throughout the atmosphere, with an explicit link to cloud formation throughout the atmosphere. This more sophisticated approach to treating the presence of clouds in the atmosphere, including allowing for patchy clouds, has obvious benefits, with the familiar caveat of increased computational time.

Fraine et al. (2021) analyse reflected light observations in the optical, taken using the HST WFC3/UVIS instruments. They do not detect clouds on the dayside of the planet within the pressure layers probed (P > 1 bar). They rule out the presence of a high-altitude, bright, uniform cloud layer, but they cannot rule out the presence of patchy clouds. Irwin et al. (2020) retrieve a colder nightside of WASP-43b using their 2.5D-setup than predicted from GCMs. They assume from their retrievals of WASP-43b that there is a thick cloud on the nightside with a cloud-top pressure <0.2 bar.

4.10 Albedo

In theory we should be able to compare our retrieved bond albedo values, indicating the proportion of energy reflected as opposed to absorbed, with those deduced from previous studies by 1 – f_day. However, we do not compute the temperature structure of the atmosphere self-consistently, which would take into account the absorption properties of the atmospheric species and their retrieved abundances. ARCiS does have the capability to do this, but including it in these 3D retrievals is currently computationally unfeasible.

The values of bond albedo of WASP-43b found in the literature are generally much lower than the values which we can deduce from our retrievals using 1 – f_day (see Table 5). For example, Keating et al. (2019) find a bond albedo for WASP-43b of ${0.22}_{-0.12}^{+0.13}$, and Stevenson et al. (2014) derive a bond albedo of 0.078−0.262 based on computed day and night fluxes.

5 Conclusions and Future Outlook

In this work we have outlined the general theory behind a 3D retrieval setup for use in characterising exoplanet atmospheres. This setup models incoming radiation from the host star and how it spreads throughout the planetary atmosphere via diffusion and longitudinal winds, which is used to find the vertical pressure-temperature profile on a grid of longitude-latitude points. The routine is provided as a stand-alone code, which can be freely downloaded via GitHub⁴. Here, we have implemented it into the Bayesian retrieval package ARCiS and run a series of 3D retrievals on spectral observations of WASP-43b as a demonstration. In running eight different retrieval setups, we are able to compare a number of different input assumptions; in particular, constraining equilibrium chemistry vs. a free molecular retrieval, testing the case of no clouds vs. parametrised clouds, and using Spitzer phase data which has been reduced from two different sources (Stevenson et al. 2017 and Mendonça et al. 2018a). We give some examples of how results from our 3D retrievals can be compared to GCM simulations, with indications of some good agreements. We hope to expand on this in the future, in a more quantitative way.

In summary, we find a consistent super-solar metallicity and super-solar C/O for all retrieval setups which assume equilibrium chemistry (see Table 5 for a summary of retrieved parameters). We find that all eight retrieval setups fit the observed spectra reasonably well, with no drastic outliers. This indicates that there are degeneracies in how the 3D setup can fit the observed data. We find that the retrievals that use observed phase-dependent emission Spitzer data from Stevenson et al. (2017) have a low heat redistribution factor, while the value is retrieved to be higher when using Spitzer data from Mendonça et al. (2018a). There is a hotspot shift found in all retrievals, with a higher shift for those using the Stevenson et al. (2017) Spitzer data than for those using the Mendonça et al. (2018a) Spitzer data. The Mendonça et al. (2018a) Spitzer dataset apparently requires a more diffusion-based solution rather than strong equatorial winds, which are more favoured by the Stevenson et al. (2017) Spitzer dataset. The trend of molecular abundances varying across the atmosphere is most consistent in the case of CH₄, where is it found to be higher on the nightside and decreasing on the dayside in all retrievals. Generally H₂O is found to consistently be lower on the nightside and increasing on the dayside, with one exception; AlO is generally expected to be most abundant on the dayside, as found in all equilibrium chemistry retrievals and in two of the free molecular abundance retrievals. The AlO abundance in the other two free retrievals remains unconstrained. NH₃ is also predicted to be relatively abundant if assuming equilibrium chemistry, but is unconstrained in all the free molecular retrievals, most likely due to the observed data wavelength coverage. CO and CO₂ are generally found to be present in high abundances for all retrievals, although again we find they are less constrained for the free molecular abundance retrievals. We find the temperature-pressure profiles to be particularly influenced by the mode of Multinest used (see Sect. 3.2). Artificially well-constrained parameters lead to seemingly very well-constrained profiles. T_int is in general hotter for all retrievals where equilibrium chemistry is assumed. Although such values of Tint are feasible based on the theoretical predictions of Thorngren et al. (2019) for an exoplanet such as WASP-43b, we note that setup I has a lower value than setup D (both setups are identical, including clouds, and assume equilibrium chemistry, but they use different modes of MULTINEST). It can be seen from the respective posterior plots in the supplementary information document of this work that this is due to T_int being unconstrained in setup I, but artificially constrained in setup D. We therefore advise caution when interpreting such results. Our retrieved cloud parameters remain unconstrained for all retrievals where clouds were included. We therefore conclude that clouds cannot be inferred from the currently available observed data using our retrieval setup with a relatively simple cloud setup (we allow for clouds to vary linearly from day to night, but do not allow for patchy clouds, for example). This can be improved by more data and a more sophisticated cloud setup. We do not infer that the atmosphere is cloud-free from our results, but perhaps our cloud-free retrievals can be considered more reliable in the current study due to the mode of MULTINEST used.

It was mentioned a few times in this work that knowledge of the atmosphere of WASP-43b would be significantly improved by more data observed over a wider wavelength range. We find that the choice of phase-dependent Spitzer emission data has a big impact on the retrieved parameters, particularly those relating to heat distribution. More observations across this wavelength region will be vital to improve our understanding of the atmosphere of WASP-43b. Fortunately, observations of WASP-43b are expected to be obtained as a result of several JWST Cycle 1 guaranteed time observations (GTO) and Early Release Science (ERS) programmes. Amongst these, JWST will observe a phase curve of WASP-43b in the 5–12 μm region using the MIRI LRS instrument, as detailed in Venot et al. (2020) and also in Bean et al. (2018). This will also be complemented by another phase curve observation by JWST, using the NIRSpec instrument to cover the region around 3-5 μm (GTO programme 1224, PI: S. Birkmann). The increased wavelength coverage that JWST will observe is expected to really help in breaking degeneracies and constraining our retrieved parameters. Ideally these proposed observations will be complemented by observations of WASP-43b at shorter wavelengths, including the visible, which will particularly help with constraining cloud parameters and temperature-pressure profiles through the ratio of mean optical to infra-red opacity. Observations over a large wavelength region using one instrument will also deal with any potential calibration issues between different instruments (in particular between Spitzer and HST in this work). This will greatly help in interpreting the physical properties of the atmosphere. We demonstrate that when such observations become available in the near future, we have a powerful tool to help interpret the observations and compare them to current theoretical predictions and properties inferred from currently available observations.

Some studies predict that there are disequilibrium chemistry effects affecting the atmosphere, and therefore affecting the observed spectra of WASP-43b (see, for example, Mendonça et al. 2018b). Although we ran a set of free chemistry retrievals that allowed for disequilibrium chemistry effects, it would be of interest in the future to treat these effects in a parametrised way by including the disequilibrium chemistry parametrisation that has recently been included in ARCiS (Kawashima & Min 2021). However, there are currently only a limited number of species included in this disequilibrium chemistry scheme. Aluminium-bearing species such as AlO are currently not included, largely due to the availability of suitable data. This is something that could be explored further in the future.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement 776403, ExoplANETS-A. K.L.C. acknowledges funding from STFC, under project number ST/V000861/1.

Appendix A Figures of Retrieval Results

Fig. A.1 Retrieved spectra for the 90° (evening), 180° (dayside) and 270° (morning) phases for eight different retrieval setups, as detailed in each panel. The shading represents the 1, 2, and 3 σ bounds. Observed data (HST and Spitzer) points are plotted with error bars. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Fig. A.2 Retrieved spectra showing the HST/WFC3 data region only for the 90° (evening), 180° (dayside) and 270° (morning) phases, for eight different retrieval setups, as detailed in each panel. The shading represents the 1, 2, and 3 σ bounds. Observed data (HST) points are plotted with error bars. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Fig. A.3 Retrieved and observed HST/WFC3 transmission spectra for our eight different retrieval setups (plus additional setup I), as detailed in each panel. The shading represents the 1, 2, and 3 σ bounds. Observed data (HST/WFC3) points are plotted with error bars. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Fig. A.4 Contour plots giving retrieved β (where β specifies how the energy incident from the central star is spread through the atmosphere; see Eq. 3) as a function of longitude and latitude for the eight different retrieval setups, as detailed above each panel, β can be used in Eq 7 to compute the temperature at a given pressure layer as a function of longitude and latitude; see Figures A.5 and A.6. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Fig. A.5 Contour plots giving retrieved temperature at an altitude of 0.01 bar (see Eq. 7) as a function of longitude and latitude for the eight different retrieval setups, as detailed above each panel. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Fig. A.6 Same as Fig A.5, but for an altitude level of 0.12 bar

Fig. A.7 Variation of retrieved AlO abundance (VMR) as a function of phase (°) for the eight different retrieval setups, as detailed in each panel. The VMRs are averaged over the visible part of the disk, centred at each of the given phases, and taken at a pressure-layer cut of 0.18 bar at the equator. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Fig. A.8 Same as Fig A.7, but for CH4.

Fig. A.9 Same as Fig A.7, but for CO

Fig. A.10 Same as Fig A.7, but for CO2

Fig. A.11 Same as Fig A.7, but for H2O

Fig. A.12 Same as Fig A.7, but for NH3

Fig. A.13 Retrieved pressure-temperature profiles for the 90° (evening), 180° (dayside), 270° (morning), and 315° (nightside) phases for eight different retrieval setups, as labelled in each panel. We note that the pressure-temperature profiles are artificially too well constrained in many of these figures, particularly in the lower part of the atmosphere.

Fig. A.14 Variation of retrieved cloud parameters as a function of longitude (°). The first column gives Pcloud (cloud-top pressure), the second column κcloud (cloud opacity), and the third column Acloud (cloud albedo). Each row is for a different retrieval setup that include clouds, as detailed in each panel. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

Appendix B Emission Benchmark

We have tested the emission spectra as computed by ARCiS against the benchmark computations presented in Baudino et al. (2017). Here we present both the results for the computations with fixed temperature structure, as well as those where the temperature structure is computed self-consistently. The computations presented in Baudino et al. (2017) are for homogeneous planets (so ID structures only). We therefore run ARCiS in the mode where f_red = 1 and v_Λ = 0. This results in a homogeneous planet. We adjust the other parameters to reproduce the benchmark planets presented in that study.

Appendix B.1 Fixed PT Profiles

Five cases are presented with a fixed pressure temperature profile. Figure B.l shows the resulting emission spectra compared to the petitCODE (Mollière et al. 2015; Mollière et al. 2017) spectra presented in Baudino et al. (2017). We show the spectra computed with ARCiS in the mode where the chemistry is computed using GGchem and in the mode where the chemical abundances are read from the ones provided by Baudino et al. (2017). We see that there are small differences in the spectra. Some of them can be assigned to differences in the computed chemical abundances, others cannot. It is most likely that these small remaining differences are caused by different molecular opacities; it can be seen that the line lists used by petitCODE, given in Table 1 of Baudino et al. (2017), are all different to those used in the present study (see Section 3.3). In particular, the Na and K lines in the optical show significant differences. It is important to note that in all five of these cases, the parameters for the central star were not given in the paper, and so we assume that scattering from the central star was also not taken into account in the spectrum from the planet.

Appendix B.2 Self Consistent PT Profiles

In Baudino et al. (2017) cases are also presented where the pressure temperature profile is computed self-consistently. We also implemented this mode in the ARCiS code. In Fig. B.2 we show the resulting PT profiles for two irradiated planets from Baudino et al. (2017). We see that the structure largely agrees. For WASP-12b we see significant differences. These were also found by Malik et al. (2019) for the same object. We conclude that likely the chemistry in the upper atmosphere is in a very delicate balance and is disturbed by the slightest difference in the opacities or the way the chemistry is computed. Nevertheless, the overall structure is still similar. In Fig. B.3 we plot the resulting spectra for these two cases. It is also important to note the computed scattered light part for these two. This is significant, as can be seen by the difference when we switch the scattering part off, and agrees well with the scattered light computed by PETITcode.

We conclude that, overall, ARCiS successfully reproduces the benchmark results from Baudino et al. (2017). Any differences are likely caused by differences in chemistry and molecular opacities.

Fig. B.1 Comparison of the emission spectra of the benchmark planets with fixed pressure temperature structures and different effective temperatures. Shown are the spectra from PETITcode (black), ARCiS reading the precomputed molecular opacities (green), and ARCiS when molecular opacities are computed using GGchem (red).

Fig. B.2 Self-consistent pressure temperature structures as computed using PETITcode Baudino et al. (2017) and ARCiS.

Fig. B.3 Spectra from the self-consistent models using PETITcode Baudino et al. (2017) and ARCiS. In green we also show the spectra as computed using ARCiS but without the scattered light component.

References

All Tables

All Figures

Fig. 1 Example structures of β as a function of latitude and longitude. Top row. two examples for a pure diffusive case (no significant longitudinal wind), while the lower panels show two examples for significant longitudinal winds. The red vertical lines indicate the position of the longitude 0, while the blue vertical lines indicate the location of the terminator. For all four cases, we take fday = 1, nΦ = 3, where fday and nΦ are the efficiency with which the planet absorbs incoming stellar radiation, and the coefficient specifying the latitudinal extent of the longitudinal winds, respectively.

In the text

	Fig. 2 Emission spectra retrieved at different phases plotted with the observed data for each phase (as given in the bottom left of each panel) for retrieval setup A (no clouds, free chemistry, and Stevenson et al. 2017 Spitzer analysis). The left of each panel shows the heat distribution over the surface of the planet, as computed by ARCiS during the retrieval.
In the text

Fig. A.1 Retrieved spectra for the 90° (evening), 180° (dayside) and 270° (morning) phases for eight different retrieval setups, as detailed in each panel. The shading represents the 1, 2, and 3 σ bounds. Observed data (HST and Spitzer) points are plotted with error bars. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

In the text

Fig. A.2 Retrieved spectra showing the HST/WFC3 data region only for the 90° (evening), 180° (dayside) and 270° (morning) phases, for eight different retrieval setups, as detailed in each panel. The shading represents the 1, 2, and 3 σ bounds. Observed data (HST) points are plotted with error bars. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

In the text

Fig. A.3 Retrieved and observed HST/WFC3 transmission spectra for our eight different retrieval setups (plus additional setup I), as detailed in each panel. The shading represents the 1, 2, and 3 σ bounds. Observed data (HST/WFC3) points are plotted with error bars. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

In the text

Fig. A.4 Contour plots giving retrieved β (where β specifies how the energy incident from the central star is spread through the atmosphere; see Eq. 3) as a function of longitude and latitude for the eight different retrieval setups, as detailed above each panel, β can be used in Eq 7 to compute the temperature at a given pressure layer as a function of longitude and latitude; see Figures A.5 and A.6. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

In the text

	Fig. A.5 Contour plots giving retrieved temperature at an altitude of 0.01 bar (see Eq. 7) as a function of longitude and latitude for the eight different retrieval setups, as detailed above each panel. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.
In the text

	Fig. A.6 Same as Fig A.5, but for an altitude level of 0.12 bar
In the text

Fig. A.7 Variation of retrieved AlO abundance (VMR) as a function of phase (°) for the eight different retrieval setups, as detailed in each panel. The VMRs are averaged over the visible part of the disk, centred at each of the given phases, and taken at a pressure-layer cut of 0.18 bar at the equator. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

In the text

	Fig. A.8 Same as Fig A.7, but for CH4.
In the text

	Fig. A.9 Same as Fig A.7, but for CO
In the text

	Fig. A.10 Same as Fig A.7, but for CO2
In the text

	Fig. A.11 Same as Fig A.7, but for H2O
In the text

	Fig. A.12 Same as Fig A.7, but for NH3
In the text

	Fig. A.13 Retrieved pressure-temperature profiles for the 90° (evening), 180° (dayside), 270° (morning), and 315° (nightside) phases for eight different retrieval setups, as labelled in each panel. We note that the pressure-temperature profiles are artificially too well constrained in many of these figures, particularly in the lower part of the atmosphere.
In the text

Fig. A.14 Variation of retrieved cloud parameters as a function of longitude (°). The first column gives Pcloud (cloud-top pressure), the second column κcloud (cloud opacity), and the third column Acloud (cloud albedo). Each row is for a different retrieval setup that include clouds, as detailed in each panel. Stevenson and Mendonça refer to the use of Spitzer data analysed by Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.

In the text

	Fig. B.1 Comparison of the emission spectra of the benchmark planets with fixed pressure temperature structures and different effective temperatures. Shown are the spectra from PETITcode (black), ARCiS reading the precomputed molecular opacities (green), and ARCiS when molecular opacities are computed using GGchem (red).
In the text

	Fig. B.2 Self-consistent pressure temperature structures as computed using PETITcode Baudino et al. (2017) and ARCiS.
In the text

	Fig. B.3 Spectra from the self-consistent models using PETITcode Baudino et al. (2017) and ARCiS. In green we also show the spectra as computed using ARCiS but without the scattered light component.
In the text