Exoplanet atmosphere retrievals in 3D using phase curve data with ARCiS: application to WASP-43b

Aims. To create a retrieval framework which encapsulates the 3D nature of exoplanet atmospheres, and to apply it to observed emission phase curve and transmission spectra of hot Jupiter exoplanet WASP-43b. Methods. We present our 3D framework, which is freely available as a stand-alone module from GitHub. We use the atmospheric modelling and Bayesian retrieval package ARCiS (ARtful modelling Code for exoplanet Science) to perform 8 3D retrievals on simultaneous transmission (HST/WFC3) and phase-dependent emission (HST/WFC3, Spitzer/IRAC) observations of WASP-43b as a case study. We assess how input assumptions affect our retrieval outcomes. In particular we look at constraining equilibrium chemistry vs a free molecular retrieval, the case of no clouds vs parametrised clouds, and using Spitzer phase data that have been reduced from two different literature sources. For the free chemistry retrievals, we retrieve abundances of H2O, CH4, CO, CO2, AlO, and NH3 as a function of phase, with many more species considered for the equilibrium chemistry retrievals. Results. We find consistent super-solar C/O (0.6-0.9) and super-solar metallicities (1.7-2.9 dex) for all retrieval setups that assume equilibrium chemistry. Atmospheric heat distribution, hotspot shift (15.6 vs 4.5 for different Spitzer datasets), and temperature structure are influenced by the Spitzer emission phase data. Comparisons are made with other studies of WASP-43b, including available GCM simulations. Conclusions. The parametrised 3D setup we have developed provides a valuable tool to analyse extensive observational datasets such as spectroscopic phase curves. Near-future observations with missions such as the James Webb Space Telescope (JWST) and Ariel will greatly improve our understanding of the atmospheres of exoplanets such as WASP-43b.


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 exoplanet 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 pressuretemperature structure is close to isothermal. The retrieval of emission spectra, which are more sensitive to the pressuretemperature profile than transmission spectra, requires a nonisothermal structure, but this will usually be parametrised (see, for example, Hansen 2008;Madhusudhan & Seager 2009;Guil-klc20@st-andrews.ac.uk m.min@sron.nl lot 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.5 D to 3 D effects, has been the focus of studies such as Irwin et al. (2020) 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(Stevenson et al. , 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 exoplanet 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 daynight 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 2 O in both transmission and emission, with Stevenson et al. (2017) deducing the presence of CO and/or CO 2 based on the Spitzer data points. Weaver et al. (2019) also find evidence of H 2 O in transmission. CH 4 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 Section 3. A total of eight different retrieval setups are detailed here, the results of which are outlined in Section 4 (and corresponding figures in Appendix A and the supplementary information document associated with this work 1 ). Discussions on various aspects of these results are also included in Section 4. Finally, our conclusions and future outlook are presented in Section 5.

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 2 . 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 on the dayside, and 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 β: 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 #» v 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: In spherical coordinates, Eq. 3 is then given by: where the source term is given by and where f day is the efficiency with which the planet absorbs the stellar radiation (i.e. in principle 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 = V Λ K ΛΦ as the ratio of the equatorial wind speed over the horizontal diffusion parameter. Eq. 4 then becomes: 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 1D 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 1D 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 pressuretemperature 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.

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 parametrisation 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: Here, γ = κ v κ 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 : 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 σT 4 irr β 4π . The irradiation from the host star is computed using T irr = T R /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 σT 4 int , 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 Equation 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, γ).

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. 2009Feroz et al. , 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 Equation 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 Section 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 Section 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 Sections 3.3 and 3.4.

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. (2014Stevenson et al. ( , 2017 can be used to essentially map the atmosphere across the planetary surface. AR-CiS 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 Section 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. 2009Feroz et al. , 2013 algorithm to sample the specified parameter space for the Table 1. Description of the parameters retrieved by ARCiS using the different 3D setups in this study, as described in the text.

Name Description Priors
Parameters included in all setups R p Planet radius 0.94 to 1.13 R J f day Proportion of incoming stellar flux absorbed by the planet's dayside atmosphere. f day = 1 − A b a 0 (no flux absorbed) to 1 (all flux absorbed) f red Heat redistribution factor. Integrated energy on the nightside divided by integrated energy on the dayside 0 (no heat redistributed from day to night) to 1 (all heat redistributed from day to night) Base-10 logarithm of the ratio of equatorial wind speed V Λ over horizontal diffusion parameter K ΛΦ b -1 for negligible winds to 8 for strong winds n φ Exponent determining the latitudinal extent of equatorial winds, given by cos n φ Φ 3 for far-reaching winds to 10 for a narrow band of wind T int Temperature at an optical depth τ = 2 3 as caused by internal heat from the planet 10 to 3000 K log(γ) Base-10 logarithm of the ratio between optical and IR opacities i.e. γ = κ v κ IR -2 to 2 log(κ IR ) Base-10 logarithm of the mean IR opacity κ IR varies from 10 −4 to 10 4 cm 2 /g log(g p ) Base-10 logarithm of the planet surface gravity c 3.56 to 3.83 (with g given in cgs units) Parameters only included in setups with clouds log(P cloud ) Cloud top pressure (function of phase) P cloud = 10 −3 bar for high-altitude clouds to 10 3 bar for deep-atmosphere clouds log(κ cloud ) A measure of cloud opacity (function of phase) κ cloud = 10 −7 cm 2 /g for transparent clouds to 10 2 cm 2 /g for very opaque clouds A cloud A measure of cloud albedo (function of phase) 0 to 1 Parameters only included in free molecular retrieval setups log(VMR(H 2 O)) Base-10 logarithm of VMR of H 2 O as a function of phase 10 −12 for negligible to 1 for very abundant log(VMR(CO)) Base-10 logarithm of VMR of CO as a function of phase 10 −12 for negligible to 1 for very abundant log(VMR(CO 2 )) Base-10 logarithm of VMR of CO 2 as a function of phase 10 −12 for negligible to 1 for very abundant log(VMR(CH 4 )) Base-10 logarithm of VMR of CH 4 as a function of phase 10 −12 for negligible to 1 for very abundant log(VMR(NH 3 )) Base-10 logarithm of VMR of NH 3 as a function of phase 10 −12 for negligible to 1 for very abundant log(VMR(AlO)) Base-10 logarithm of VMR of AlO as a function of phase 10 −12 for negligible to 1 for very abundant Parameters only included in equilibrium chemistry retrieval setups C/O Atmospheric carbon-to-oxygen ratio (global) 0.1 (sub-solar) to 1.5 (super-solar) [Z] Atmospheric metallicity (global) -3 (sub-solar) to 3 (super-solar) dex a : A b is the bond albedo.
b : The value of K ΛΦ is adjusted until the given value of f red is reached. c : We use a Gaussian prior on the mass based on radial velocity measurements of WASP-43b (Gillon et al. 2012). We compute the mass based on R p and log(g p ) and then place the prior on this derived parameter.
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 Article number, page 5 of 31 A&A proofs: manuscript no. WASP43b_3D_phases_after_lang_corr_twocol Table 2. Observed data used in the WASP-43b retrieval setups described in this work. All phase data from HST are from Stevenson et al. (2014). The phase data observed using Spitzer are from Stevenson et al. (2017) and were independently reduced by two different sources: Stevenson et al. (2017) or Mendonça et al. (2018a). We use slightly different phases in this work for these two cases, as detailed in the footnote of this 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.

Clouds vs no clouds
Clouds are predicted to be present on WASP-43b, particularly on the nightside (Helling et al. , 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 (cloudtop pressure, ranging from 10 −3 bar for high-altitude clouds or hazes to 10 3 bar for deep-atmosphere clouds), κ cloud (a measure of cloud opacity from 10 −7 cm 2 /g for transparent clouds to 10 2 cm 2 /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 Section 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 Section 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(Feroz et al. , 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 700,000). 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 Section 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.

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 ) website 3 . Opacity formats for other retrieval codes are also available from ExoMolOP: Tau , 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]: Here, X i is the number density of element i.

Free molecular chemistry
For the free molecular retrievals, we include the base set of molecules (H 2 O, CO, CO 2 , CH 4 , and NH 3 ) 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 AR-CiS 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 2 O, CH 4 , CO, CO 2 , AlO, and NH 3 . All line lists used are the same as those listed in Section 3.3. 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 Section 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 5, 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 Section 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. Table 3. Description of the different retrieval setups presented in this work. We are considering setups: (i) with and without clouds, (ii) free chemistry vs equilibrium chemistry, and (iii) using different sets of Spitzer points. The retrievals with clouds included were run in the constant efficiency mode of Multinest, and those without clouds in the full efficiency mode of Multinest, with the exception of setup I. See text of Section 3.2 for more information.  −0.12 in the scenario where they allow a completely free temperature profile. This temperature profile setup is different to our parametrised Guillot profile (see Section 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 4 , 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, Table 5. Retrieved parameters found using different retrieval setups. See Table 1 for definitions of the parameters. The different retrieval setups are explained in the footnote to this table. The error bounds on R p are all in the region of 1 × 10 −4 (in units of Jupiter radii). could also be an indication of cloud formation, which typically de-oxidises the gas phase (Helling 2019). 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 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.

Retrieved spectra at different phases
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. 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 Section 4.2 and the transmission spectra shown in Figure 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 metallicities 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 con-trolled 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 Section 3.2), and so the posterior distributions given in the corresponding plots are more constrained than they should be.

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 Section 2). This value of β can be used in Equation 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 Figures A.5 and 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 Figure 2 (d) and (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 Figure 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 Figure A.5). Figure 3 of Kataria et al. (2015) gives the longitude-latitude temperature cut at various pressure levels for 1 × solar metallicity and for 5 × solar in their Figure 7. The temperature cuts at a pressure layer of 10.63 mbar (≈ 0.01 bar ) can be compared to our Figure 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 1 × 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 2 (Ib) of Carone et al. (2020) gives the longitudelatitude 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 Figure 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 exoplanets (see e.g. Cho et al. 2003Cho et al. , 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: 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.005 −0.001 by integrating their retrieved model spectra. This corresponds to a value of f red = 0.002 +0.010 −0.002 , 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 Section 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 recirculation) 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

Influence of Spitzer phase curve data reduction
As mentioned in Section 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 Section 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 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.

Hotspot shift
All panels of Figures A.4, A.5, and 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 Figure 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.

Variation of molecular abundances across the atmosphere
Figures A.7 -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 Figures A.7 -A.12. The same data used as presented for these six species are available for each species included in the equilibrium chemistry setup (see Section 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 4 , CO, CO 2 , H 2 O, and NH 3 ). The plots of Figures A.8 -A.11 (which correspond to CH 4 , CO, CO 2 , and H 2 O, respectively) can be compared to Figure 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.

AlO
As can be seen in Figure 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.

CH 4
We find a very consistent trend of CH 4 abundance as a function of phase, particularly in the retrievals cases which use the Stevenson et al. (2017) Spitzer data (see Figure A.8). We find CH 4 to vary between ∼10 −4 on the dayside and ∼10 −3 on the nightside. This can be compared to Figure 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.

CO and CO 2
We retrieve relatively high abundances of both CO and CO 2 (see Figures A.9 and 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 −4 . 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 2 also has some weaker features in the HST/WFC3 region.

H 2 O
We find the same trend with all of our retrieval setups (with one exception; see Figure A.11): a higher abundance of H 2 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 2 O abundance for the free retrievals that include clouds are artificial, due to the constant efficiency mode of Multinest (see Section 3.2).

NH 3
Although NH 3 is predicted to be relatively abundant by the equilibrium chemistry retrievals (see Figure 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.  May & Stevenson (2020) and Bell et al. (2021) 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 Figure 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 Figure 6). We can also compare to GCM outputs, for example see the first panel of Figure 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 Figure 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 Section 3.2), but likely also because of the shape of the parametrised profile. We also include the results of re-trieval setup I as part of Figure A.13, which is the same as setup D but using the full efficiency as opposed to the constant efficiency mode of Multinest (see Section 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 τ = 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 10 (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 Section 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.

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 Section 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 Section 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. 2020Helling et al. , 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 AR-CiS (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.

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 1f 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

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 pressuretemperature profile on a grid of longitude-latitude points. The routine is provided as a stand-alone code, which can be freely downloaded via GitHub 4 . 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 phasedependent 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 4 , where is it found to be higher on the nightside and decreasing on the dayside in all retrievals. Generally H 2 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 3 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 2 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 Section 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 T int 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. Aluminiumbearing 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.

Supplementary material
The  Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.  Stevenson et al. (2017) and Mendonça et al. (2018a), respectively.         Variation of retrieved cloud parameters as a function of longitude ( • ). The first column gives P cloud (cloud-top pressure), the second column κ cloud (cloud opacity), and the third column A cloud (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 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.

P[bar]
PETITcode ARCiS GJ436b WASP-12b . 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.