Testing 2D temperature models in Bayesian retrievals of atmospheric properties from hot Jupiter phase curves

Spectroscopic phase curves of transiting hot Jupiters are spectral measurements at multiple orbital phases, giving a set of disc-averaged spectra that probe multiple hemispheres. By fitting model phase curves to observations, we can constrain the atmospheric properties of hot Jupiters such as molecular abundance, aerosol distribution and thermal structure, which offer insights into their dynamics, chemistry, and formation. In this work, we propose a novel 2D temperature scheme consisting of a dayside and a nightside to retrieve information from near-infrared phase curves, and apply the scheme to phase curves of WASP-43b observed by HST/WFC3 and Spitzer/IRAC. In our scheme, temperature is constant on isobars on the nightside and varies with cos$^n$(longitude/$\epsilon$) on isobars on the dayside, where $n$ and $\epsilon$ are free parameters. We fit all orbital phases simultaneously using the radiative transfer package NEMESISPY coupled to a Bayesian inference code. We first validate the performance of our retrieval scheme with synthetic phase curves generated from a GCM, and find our 2D scheme can accurately retrieve the latitudinally-averaged thermal structure and constrain the abundance of H$_2$O and CH$_4$. We then apply our 2D scheme to the observed phase curves of WASP-43b and find: (1) the dayside temperature-pressure profiles do not vary strongly with longitude and are non-inverted; (2) the retrieved nightside temperatures are extremely low, suggesting significant nightside cloud coverage; (3) the H$_2$O volume mixing ratio is constrained to $5.6\times10^{-5}$--$4.0\times10^{-4}$, and we retrieve an upper bound for CH$_4$ at $\sim$10$^{-6}$.


INTRODUCTION
Exoplanet surveys suggest planets are common around stars in our galaxy (Winn & Fabrycky 2015).The diversity in their characteristics from system architecture to bulk properties poses challenging questions in the theory of planetary formation (Mordasini et al. 2009).Gaseous giant planets with close-in orbits (period < 10 days), dubbed hot Jupiters, are a key piece of the puzzle for two reasons: (1) they likely undergo significant orbital migration and play an important role in shaping planetary system architecture (Dawson & Johnson 2018), and (2) they are the easiest targets for spectroscopic characterisation, and the constraints on their atmospheric properties give valuable insights into planetary formation (Madhusudhan et al. 2017;Mordasini et al. 2016).
The spectral appearance of hot Jupiters is determined by the opacity structure and the thermal structure of their atmospheres.Conversely, by fitting spectra generated from atmospheric models ★ E-mail: jingxuan.yang@hertford.ox.ac.uk to observations, we could constrain the atmospheric properties of these planets in a process known as atmospheric retrievals (Irwin et al. 2008;Madhusudhan & Seager 2009;Line et al. 2013;Changeat & Al-Refaie 2020;Cubillos & Blecic 2021;MacDonald & Batalha 2023).Two observing methods are widely used: (1) transmission spectroscopy, which measures the stellar light filtered through the planetary limb during primary transits (Barstow 2017;Sing et al. 2016); (2) eclipse spectroscopy, which extracts the dayside emission spectra by monitoring the combined stellar and planetary flux during secondary eclipses (Lee et al. 2012;Mansfield et al. 2021).Such observations have been done at both high resolution (e.g., Brogi & Line 2019) from ground-based facilities and at low-resolution using space telescopes (e.g., Wakeford et al. 2017), resulting in a myriad of atomic and molecular detections (e.g., Fe : Hoeĳmakers et al. 2018;Na: Snellen et al. 2008; H 2 O: Evans et al. 2016; CO 2 : Team et al. 2022).
A major challenge in the analysis of low-resolution emission spectra is the degeneracy between thermal structure and molecular abundance; in particular, modelling highly inhomogeneous thermal structure with a single temperature-pressure (TP) profile can lead to significant biases in retrieved molecular abundance (Feng et al. 2016;Blecic et al. 2017;Taylor et al. 2020).This degeneracy can be partially broken by measuring the emission spectra of transiting hot Jupiters at multiple orbital phases.As hot Jupiters are likely tidally locked to their stars due to the short orbital separations, it is straightforward to relate the orbital phases to the central longitudes of the visible hemisphere.Note that we assume the planets have edge-on orbits, and we define the equator to be in the orbital plane.Such observations, called 'phase curves,' allow us to constrain the thermal structure of the atmospheres better, resulting in better constraints on the chemical abundance as well.
In order to retrieve information from a set of phase curves using Bayesian inference, we first need to construct appropriate parametric atmospheric models.If we want to analyse all orbital phases simultaneously, then the models must describe the entire observable atmosphere to generate disc-averaged emission spectra at multiple orbital phases.Crucially, we need multidimensional temperature models that can capture the longitudinal variation of thermal structure in the pressure ranges probed by emission spectroscopy.The model should also contain as few parameters as possible to ensure Bayesian parameter estimation can be done in reasonable time, and to avoid overfitting.There are two main approaches to this problem.The first approach is to split the atmosphere into disjoint regions, where the thermal structure in each region is modelled with a onedimensional TP model.For example, Feng et al. (2020) split the atmosphere into a dayside and a nightside, whereas Changeat et al. (2021) further model an additional hot spot within the dayside, and Irwin et al. (2020) divide the atmosphere into meridian bands with linearly interpolated temperature maps between the bands.The second approach is to construct highly-simplified three-dimensional analytical atmospheric models.For example, Dobbs-Dixon & Blecic (2022) propose a 3D model by separating radiative and convective components from Global Circulation Model (GCM) outputs, whereas Chubb & Min (2022) prescribe a 3D model by restricting heat transfer to diffusion and zonal winds.
The studies summarised above offer unique and valuable perspectives on the analysis of hot Jupiter phase curves.However, it is difficult to compare the different retrieval schemes for several reasons.Firstly, the ways in which the retrieval schemes are validated differ significantly across studies.The most direct way to assess the performance of a retrieval scheme is first to create synthetic data from a model atmosphere, then test how well the retrieval scheme can recover the input atmospheric properties.Out of the studies that include such validation tests, the synthetic data are often generated from toy models resembling the temperature parameterisations of the retrieval schemes, so it needs to be clarified how well the retrieval models can perform on data generated from more realistic atmospheric models.Secondly, there are multiple modelling steps within each retrieval study, for example, the modelling of thermal structure, the modelling of chemical abundance, and the modelling of radiative transfer, which can all vary across studies.Thirdly, the studies that perform analysis of real observations often do not retrieve on exactly the same data set, which hinders the comparison of the retrieved constraints.
In this work, we propose a novel 2D retrieval scheme (model 4 in section 2.4) that can be used to retrieve chemical abundance and thermal structure from hot Jupiter phase curves.In this model, the temperature is a function of pressure and longitude.The model is split into a dayside and a nightside: on the dayside, the temperature varies with cos  (longitude/) on isobars, where  and  are free parameters, and on the nightside the temperature is con-stant on isobars.We use this scheme, together with several other simpler 2D retrieval schemes for comparison, to retrieve molecular abundance and latitudinally averaged thermal structure from phase curves of WASP-43b observed by Hubble Space Telescope/Wide Field Camera 3 (HST/WFC3) and Spitzer/Infra-Red Array Camera (Spitzer/IRAC).We first validate the performance of the 2D schemes by retrieving atmospheric properties from synthetic phase curves generated from a GCM-based model of WASP-43b, where the 'ground truth' is known, so we can assess the accuracy of the retrieved properties.We then apply the 2D schemes to the observed HST/WFC3 and Spitzer/IRAC phase curves of WASP-43b.We also compare the 2D approach to the phase-by-phase approach, where the spectrum at each orbital phase is analysed separately, in appendix B.
This paper is structured as follows.In section 2, we describe our routine for simulating spectroscopic phase curves, our 2D temperature models, and our retrieval set-up.Section 3 demonstrates that simplified atmospheric models can reproduce the synthetic data generated from a GCM.Section 4 presents the retrieval results on synthetic GCM phase curves using our 2D temperature models, followed by application to the observed phase curves in section 5. We discuss the implications of our results and compare them with past studies in section 6 and end with a conclusion in section 7.

METHODOLOGY
This work aims to assess the performance of a novel 2D parametric temperature model (model 4) in retrieving atmospheric properties from low-resolution spectroscopic phase curves.We also test three simpler models (model 1, 2, and 3) for comparison.The data we model are the phase curves of WASP-43b, as presented in Stevenson et al. (2017), and synthetic phase curves of the same resolution simulated from a GCM-based model.The planet WASP-43b is a hot Jupiter around a K7 star discovered by Hellier et al. (2011), with planetary parameters of 2.034±0.052Jupiter mass and 1.036±0.019Jupiter radii as given by Gillon et al. (2012).Due to its short 19.5hour orbit and large planet-to-star flux ratio, WASP-43b is a prime target for phase-resolved spectroscopic observations.The observation contains 15 phase curves from the HST/WFC3 instrument, which are binned in equally spaced bins of width 0.035 μm spanning the wavelength range 1.1425-1.6325μm, and two phase curves from Spitzer/IRAC broad channels centred at 3.6 and 4.5 μm.
In this section, we describe the GCM used for simulating synthetic phase curves in 2.1, our procedure for the radiative transfer calculation in 2.2, and our method for calculating disc-averaged spectra in 2.3.We describe our 2D atmospheric temperature models in 2.4, our retrieval set-up in 2.5, and our Bayesian inference scheme in 2.6.

GCM data
We use a GCM of WASP-43b to simulate synthetic phase curves to validate our retrieval schemes.We can directly assess the performance of our retrieval schemes by comparing the atmospheric properties retrieved from the synthetic data with the input GCM.The GCM is a cloud-free model calculated using SPARC/MITgcm (Showman et al. 2009) based on the set-up of Parmentier et al. (2016), and used for validating the 2.5D retrieval scheme by Irwin et al. (2020).We plot temperature as a function of longitude and latitude at three pressure levels of the GCM in Fig. 1 .Temperature (Kelvin) as a function of longitude and latitude at three pressure levels in the WASP-43b GCM.The super-rotating equatorial jet is clearly visible and shifts the 'hot spot' eastward of the substellar point (where the star would be perceived to be directly overhead).Note that the substellar point is at 0 degree longitude.Such jet-like features would cause the phase curve amplitudes to peak before secondary eclipses.Note that the latitudinal distance is weighted by cos(latitude) to mimic the effect that polar latitudes would appear foreshortened to us because we observe WASP-43b from above the equator.
of pressure and wavelength channel number in Fig. 2. The model is H 2 /He dominated and contains four spectrally active gas species: H 2 O, CO, CO 2 and CH 4 , which are expected to be the dominant opacity sources in the atmosphere of WASP-43b in the observed wavelengths.The chemical abundance in the GCM is initially set according to chemical equilibrium, resulting in significant variation in CH 4 abundance from dayside to nightside in the photospheric pressures.However, disequilibrium chemistry processes such as horizontal quenching are expected to smooth out such inhomogeneity (Cooper & Showman 2006;Agúndez et al. 2014).Hence, we reset the abundance of all molecules to be the latitudinally averaged abundances in the 0.1 to 1-bar pressure region (using cos(latitude) as the weight) at the sub-stellar meridian, following Irwin et al. (2020).We then use this model to simulate the synthetic data.By resetting the chemical abundance to uniform values, we isolate the effect of temperature parameterisation in retrievals.We discuss the motivation and implication of using uniform abundance in section   1.The synthetic phase curves are thus generated from the thermal structure of the GCM and the uniform VMRs, using the disc-averaging scheme described in section 2.3.The chemical abundance and thermal structure of this GCM-based model are seen as the 'ground-truths' for our retrieval tests in section 4.

Radiative transfer calculation
We use the correlated-k method (Lacis & Oinas 1991) to accurately and efficiently implement our radiative transfer calculations, following Irwin et al. (2008).Consider the mean transmission for a homogeneous path of absorber amount  in the wavelength bin [,  + Δ], where  () is the absorption cross-section. is the key quantity that links opacity structure and thermal structure in thermal emission calculations.The cross-section  is a rapidly varying function of , so it is computationally expensive to numerically calculate equation (1).However, since the ordering of  in the wavelength bin [,  + Δ] does not affect the value of equation ( 1), we sort  in ascending order within each wavelength bin, which gives a monotonic distribution of  that is easier to handle in quadrature schemes.Mathematically, let the cumulative frequency distribution of  be (), then the inverse of (), which we denote as  (), is well-defined and monotonic.The function  () is called the -distribution, and can be tabulated on a grid of pressures and temperatures for each spectrally active molecule before calculations.
During radiative transfer calculations, the -distributions of multiple gases are combined with the random-overlapping-line approximation (Lacis & Oinas 1991).Such approximation gives residuals insignificant compared with measurement error, as found by Irwin et al. (2020) and Mollière et al. (2015).To calculate  through an inhomogenous path, we first split the path into multiple subpaths (Irwin et al. 2008) that are sufficiently homogeneous, then model each sub-path with the absorber-amount weighted averaged sub-path properties.In the monochromatic case, the transmission of each sub-path can be multiplied together to give the transmission of the total path.However, to use the -distribution technique where we have reordered , we need to additionally assume that the wavelength regions of high opacity are correlated for all sub-paths, which is the correlated-k approximation.This is a good approximation for our set-up, as we are assuming constant vertical distribution of chemical abundance in our atmospheric model.Using the -distribution technique, the mean transmission of a sub-path as defined in equation ( 1) can be well-approximated with a Gaussian quadrature scheme: where   is the th quadrature point, Δ  the corresponding weight, and   the number of quadrature points.The total transmission of an inhomogenous path is then where we have multiplied the transmission of all  layer layers together.We use -distribution look-up tables ('-tables') with   = 20 generated from the line data summarised in Table 1.
To speed up calculations, we use channel-averaged -tables for the 15 channels of HST/WFC3 and the 2 channels of Spitzer/IRAC.Irwin et al. (2020) find that such an approach produces an excellent approximation, resulting in residuals much less than typical measurement uncertainties of observations using these facilities.Apart from the above molecular opacity, we additionally include collisioninduced absorption of H 2 -H 2 pairs and H 2 -He pairs using the coefficients of Borysow et al. (1989) and Borysow & Frommhold (1989), as well as Rayleigh scattering for a H 2 /He dominated atmosphere using data from Allen (1976).
. Illustration of our disc averaging scheme at 45 degree orbital phase (0 degree being the primary transit and 180 degree the secondary eclipse).
The visible region of the illuminated dayside at this orbital phase is shaded pink, whereas the visible nightside region is shaded in grey.The crosses mark the quadrature points for disc integration, and the dashed circles mark the positions of the zenith angle quadratures.

Disc-average scheme
The disc-averaged spectral radiance (W m −2 sr −1 μm −1 ) of an inhomogeneous atmosphere for a distant observer is where  = cos() is the cosine of the zenith angle 1 , and  is the azimuth angle.To carry out sampling-based Bayesian parameter estimation, we need to evaluate R() for many different model atmospheres to approximate the posterior distributions, so it is important to have a numerical integration scheme for equation ( 4) that is both accurate and computationally inexpensive.We use the method of Irwin et al. (2020): the zenith integration with respect to  is done with a Gauss-Lobatto quadrature scheme with   quadrature points, while the azimuthal integration with respect to  is done with a Trapezium rule quadrature scheme with   quadrature points.For the trapezium rule integration, the quadrature points are placed on the circles corresponding to each zenith quadrature point such that the arc-length between neighbouring points is ∼  plt /  (Fig. 3).Overall, our numerical integration scheme for equation ( 4) is given by: R() = 2 where   are the Lobatto quadrature points for the zenith integration, Δ  are the corresponding quadrature weights,    are the Trapezium rule quadrature points for the azimuth integration (which are different for each zenith angle), and    are the Trapezium rule quadrature weights for the th azimuthal angle and the th zenith angle.We assume that the orbit of the planet is exactly edge-on, and that the atmosphere is symmetric about the equatorial plane, which is defined to be in the orbital plane.As a result, we only need to evaluate the integration over half of the visible disc and multiply the result by a factor of two.The quadrature points for   = 2, 3, 4, 5 are shown in Fig. 3.For atmospheric models that partition the atmosphere in longitude into several regions each modelled with a single TP and abundance profile, our disc average routine can be further simplified.For example, if the planet is divided into a uniform dayside and a uniform nightside, as illustrated in Fig. 3, then all the quadrature points on the same zenith angle ring (blue dashed circles) in the same region have the same radiance.The azimuth integration is then a matter of calculating what fraction of the zenith angle rings are in each region.This greatly speeds up the disc averaging routine for the simple models in 2.4.

2D atmospheric temperature models
We describe four parametric temperature models for the atmospheres of hot Jupiters.Model 4 is our proposed model, and the other three simpler models are included for comparison.While all of the models use the one-dimensional analytical TP profile of Guillot (2010) to describe temperature as a function of pressure, they can also be easily interfaced with other parametric TP profiles.The Guillot profile is given by equation ( 29) of Guillot ( 2010) where  is the infrared optical depth defined by The TP profile has four free parameters:  th is the mean infrared opacity,  is the ratio between the mean visible and mean infrared opacities, T int is the internal heat flux, and  is a catch-all parameter of order unity that models the effects of albedo (on the dayside) and the redistribution of stellar flux due to atmospheric circulation (on both the dayside and the nightside).We assume the change in gravity  is negligible in the pressure range probed by emission spectroscopy, so that  is linear in .Finally,  irr is the irradiation temperature defined by where  is the orbital semi-major axis, and  star and  star are the host star radius and temperature.In section 3.1, we show that the Guillot profile is able to approximate the temperature-pressure profiles found in the WASP-43b GCM well enough to reproduce the synthetic data.For all models, we place the sub-stellar point at the origin in our longitude-latitude coordinate system.We denote lon-gitude by Λ.All models describe the thermal structure with two TP profiles: a representative dayside profile and a representative nightside profile.Note that temperature is set to be uniform as a function of latitude on isobars, and we treat the retrieved temperature profiles as latitudinally averaged temperature profiles.In our models, the centres of 'dayside' and the 'nightside' can shift away from the substellar point and the anti-stellar point, respectively.

Model 1
Model 1 (Fig. 4, top panel) our simplest model in which the dayside and nightside are both set to span 180 • in longitude.The centre of the dayside region ( ′ ) is allowed to shift eastward or westward by some longitude  relative to the substellar point, representing the effect of atmospheric dynamics, in particular equatorial jets, on redistributing heat around the atmosphere.The dayside is thus bound by the meridians Λ =  − 90 • and Λ =  + 90 • .Within each region, temperature is only a function of pressure.Note that in all our models, the 'dayside' denotes the region of the atmosphere modelled by the dayside TP profile and does not necessarily coincide with the physically illuminated dayside.The model contains 9 parameters: 4 parameters for each of the two TP profiles and 1 parameter  for the longitudinal offset.Note that model 1 is equivalent to the '2TP-Crescent' approach of Feng et al. (2020) when applied to all phases simultaneously.In summary, model 1 is given by: (9)

Model 2
Model 2 (Fig. 4, lower panel) is similar to model 1, with the only difference being that the 'width' (longitudinal extent) of the dayside is now a free parameter.We introduce a scaling parameter , so that the dayside is now bounded by the meridians Λ =  − 90 • ×  and Λ =  + 90 • ×  and spans 180 • ×  in longitude.The model contains 10 parameters: 4 parameters for each TP profile, 1 parameter  for the longitudinal offset, and 1 parameter  for the dayside width.Parameterising the dayside area fraction has been shown to be effective in analysing disc-averaged emission spectrum of tidally-locked hot Jupiters by Taylor et al. (2020), and this approach has been applied to phase curve analysis by Feng et al. (2020) in their '2TP-Free' model, albeit only in the phase-by-phase approach.Model 2 is our way of implementing the dayside area fraction parameterisation self-consistently when fitting all phases of phase curves simultaneously.In summary, model 2 is given by: (10)

Model 3
Model 3 (Fig. 5, upper panel) is an extension of model 2, where the temperature is now a continuous function of longitude across the dayside boundary.The dayside is bound by the meridians Λ =  − 90 • ×  and Λ =  + 90 • × .Within the dayside, temperatures at each pressure level vary with the cosine of longitude.Outside the dayside, the TP is set to be a single nightside profile.In summary, model 3 is given by: (11)

Model 4
Model 4 is a generalisation of model 3, where we allow the exponent of the cosine term in equation ( 11) to be a variable, so that we parameterise how strongly temperatures vary with longitude on isobars: The parameters for the temperature models, together with the other parameters of our atmospheric models, are summarised in Table 2.

Retrieval set-up
We run retrievals on two sets of data: (1) synthetic HST/WFC3 and Spitzer/IRAC phase curves simulated from the GCM-based model of WASP-43b described in section 2.1, and (2) observed HST/WFC3 and Spitzer/IRAC phase curves of WASP-43b as presented in Stevenson et al. (2017).For each set of data, we run four retrievals using each of the atmospheric models described in section 2.4.We fit the spectra at all phases simultaneously using spectra generated from the parametric atmospheric models.We use Nested Sampling (Feroz & Hobson 2008) to calculate the posterior distribution of the atmospheric model parameters and the Bayesian evidence of the model, described in section 2.6.In appendix B, we also compare our retrieval results to the phase-by-phase retrieval approach, where the spectrum at each orbital phase is analysed independently.For all of our retrievals, the atmospheric model is defined from 20 to 10 −3 bar, on 20 points equally spaced in log pressure.The atmospheric models have two components: a temperature model and a chemical abundance model.In each of our retrieval schemes, we test a different temperature model described in section 2.4.On the other hand, all of our retrieval schemes share the same chemical abundance model, which assumes a H 2 /He dominated atmosphere and contains four spectrally active gases: H 2 O, CO 2 , CO, CH 4 .The abundance model is parameterised by the volume mixing ratios of the spectrally active gases, assumed to be constant with respect to pressure, longitude and latitude.Furthermore, we do not include clouds/hazes in any of our models, and we assume aerosols with no significant spectral features to be degenerate with the other components of our atmospheric models.The limitations of these assumptions are discussed in section 6.The model parameters and their prior ranges are listed in Table 2, and we prescribe uniform priors for all of our model parameters.

Bayesian parameter estimation
We extract information from phase curves using Bayesian inference.Consider a set of phase curve data  that we wish to analyse.Suppose we have a parametric atmospheric model  with parameter space Θ, so that for each point  ∈ Θ we can calculate model phase curves  ( ()), where  is our 'forward model' that encapsulates all the modelling steps required to generate model phase curves from an atmospheric model.The probability distribution of the parameters of  given  is where () ≡ Pr(,  |) is the posterior distribution, L () ≡ Pr(|, ) is the likelihood, () ≡ Pr(|) is the prior, and  ≡ Pr(|) is the evidence.To proceed, we define the log likelihood function as where  obs is the total number of data points in the observed phase curves,   and  ( ())  are the th points of the observed and model phase curves, respectively, and   is the associated measurement uncertainty.The Bayesian evidence is then given by where the integral is over the -dimensional parameter space Θ.
We can approximately calculate the posterior distribution and the evidence by sampling the parameter space Θ, which is a computationally expensive task for phase curve retrievals because Θ is high-dimensional and  ( ()) is expensive to calculate.The Nested Sampling algorithm (Feroz & Hobson 2008) is an efficient way to carry out these tasks, which starts by rewriting the evidence in terms of the prior volume , defined as The evidence is then given by where L () is a monotonic function and can be evaluated with simple quadrature schemes, and the likelihood contours L (  ) are approximated by sampling the parameter space within nested ellipsoids.Numerically, the evidence is given by where   ∈ [0, 1] are a decreasing sequence of quadrature points (starting from 1) and   the corresponding weights, and  iter is determined by some convergence criterion.The nested sampling algorithm discards the point with the lowest-likelihood at each iteration, so the posterior can be generated by assigning weights to those points by For all of our retrievals, we use the Python interface PyMultiNest (Buchner et al. 2014) to implement nested sampling with 1000 sampling live points.

PRELIMINARY TESTS
Before testing our retrieval schemes on synthetic HST/WFC3 and Spitzer/IRAC phase curves, we investigate if simplified atmospheric models can reproduce the phase curves generated from the GCM described in section 2.1.We show that we can reproduce these GCM phase curves to well within realistic measurement uncertainties: (1) if we replace the TP profiles in the GCM at all locations with best-fit 1D Guillot profiles; (2) if we replace the TP profiles in the GCM on the same meridian with the latitudinally averaged TP profile of that meridian; and (3) if we replace the volume mixing ratios of all gases with a uniform profile.These results justify our use of 2D models coupled with the Guillot profile in our retrievals, and in section 4, we show that such 2D models can adequately model synthetic HST/WFC3 and Spitzer/IRAC phase curves.

Replace GCM TP profiles with 1D model fits
We use the Guillot TP profile (Guillot 2010) as the basis of our 2D temperature models described in section 2.4.We demonstrate here that this profile is flexible enough to approximate the range of TP profiles found in our WASP-43b GCM with the following procedure.First, we directly fit the 1D TP profile to the TP profiles of the GCM on all longitude-latitude grid points in the pressure range 20-10 −3 bar, which covers the support of the transmission weighting function.We find that extending the pressure range has negligible effects on the spectra.We then generate phase curves from the total collection of best-fit 1D profiles, and compare them to those generated directly from the GCM.Both sets of phase curves are simulated using the volume mixing ratios of the original GCM, which are set via chemical equilibrium.In Fig. 6, we compare the phase curves simulated from the best-fit 1D profiles (blue curves) with the phase curves simulated directly from the GCM (black curves).We overplot the measurement errors of Stevenson et al. (2017) on the phase curves simulated directly from the GCM.We see that the phase curves simulated from the 1D best-fit profiles can match the GCM phase curves to within error at almost all phases.1, which are the synthetic data we retrieve on in section 4. The orange curves are simulated with the abundance listed in Table 1, but with the methane abundance multiplied by 40, which illustrates the fact that models with uniform gas abundance can match the phase curves produced from a chemical equilibrium model.Note that the error bars on the GCM phase curves are the estimated observational uncertainties of Stevenson et al. (2017).

Replace GCM thermal structure with latitudinally-averaged thermal structure
In the 2D models described in section 2.4, the atmospheric temperature varies with longitude and pressure only, and we prescribe that temperature is constant with respect to latitudinal variation.The key point is that we interpret the retrieved thermal structure as a latitudinally averaged thermal structure, as we have very limited sensitivity to latitudinal variation of atmospheric properties in the data.We now demonstrate that the latitudinally-averaged thermal structure of the GCM can reproduce the synthetic data, which justifies our choice of 2D temperature models.To this end, we replace all the TP profiles on the same meridian with some latitudinally averaged TP structure.We find that if we pick the TP profile averaged from 0 degree latitude to 45 degree latitude using cos(latitude) as the weight, the resulting 2D temperature model could produce phase curves (green curves, Fig. 6) that agree with the phase curves simulated directly from the GCM (black curves, Fig. 6) to measurement uncertainties quoted in Stevenson et al. (2017).Note that both sets of phase curves are simulated using the chemical equilibrium VMRs of the original GCM.In conclusion, we find that a 2D temperature model can reproduce HST/WFC3 and Spitzer/IRAC quality phase curves simulated from our WASP-43b GCM.Furthermore, we expect the retrieved thermal structure using a 2D model to resemble the latitudinally averaged thermal structure of the atmosphere.In Fig. 7, we plot the GCM temperature as a function of pressure and longitude at the equator on the top, compared to the latitudinally averaged temperature (from 0 degree latitude to 45 degree latitude using cos(latitude) as the weight) on the bottom.In section 4, we compare the retrieved thermal structure from the synthetic data with this latitudinally averaged GCM thermal structure.

Replace chemical equilibrium VMRs with constant VMRs
As mentioned in section 2.1, the chemical equilibrium abundance of the original WASP-43b GCM is expected to be homogenised by horizontal quenching (Cooper & Showman 2006;Agúndez et al. 2014), and Irwin et al. (2020) reset the GCM gas abundances at all altitudes and locations to be the latitudinally averaged abundances in the 0.1-1-bar pressure region (using cos(latitude) as the weight) at the sub-stellar meridian.We now compare the phase curves simulated using uniform abundance as those used by Irwin et al. (2020) (yellow curves) with those simulated using the chemical equilibrium abundance (black curves) in Fig. 6.We see the main difference is that the different distributions of CH 4 result in significantly different phase curves at 3.6 μm.This effect has been investigated by Steinrueck et al. (2019), who explore if disequilibrium effects such as the quenching of CH 4 can explain why GCMs systematically overestimate phase curve amplitudes compared to observations.We echo the finding of Steinrueck et al. (2019) that phase curves observed in the wavelength range covered by the 3.6 μm Spitzer channel are an effective diagnostic for disequilibrium methane chemistry on hot Jupiters.
The synthetic phase curves we use to validate our retrieval schemes are simulated from uniform gas VMRs as listed in Table 1 and are the same synthetic phase curves in Irwin et al. (2020).To further justify our use of uniform gas abundance in our model, we show that if we multiply the CH 4 abundance of Irwin et al. (2020) by a factor of 40, the resultant phase curves (orange curves) agree well with the phase curves simulated from the equilibrium chemistry VMRs (black curves).This suggests that while using uniform VMRs can adequately fit HST/WFC3 and Spitzer/IRAC quality phase curves, this approach can lead to significantly biased CH 4 abundance.The limitation of the constant chemistry assumption is discussed in 6.3.3.

APPLICATION TO SYNTHETIC PHASE CURVES
We validate our retrieval schemes with synthetic data simulated from the GCM-based model of WASP-43b described in section 2.1.Each retrieval scheme is identified with one of the temperature models described in section 2.4, while all other aspects of the retrieval schemes are identical, so we refer to each retrieval scheme by the temperature model used.The synthetic phase curves are simulated at the same wavelengths as those presented in Stevenson et al. (2017), and we use their measurement uncertainties to set the uncertainties of the synthetic data.We do not add random noise to the synthetic phase curves.We assess the retrieval schemes on three criteria: (1) the goodness of fit to the synthetic phase curves; (2) the accuracy of the retrieved chemical abundance; (3) the goodness of fit of the retrieved thermal structure to the latitudinally averaged GCM thermal structure in Fig. 7. Overall, apart from model 1, all other models can fit the synthetic phase curves within measurement uncertainties at almost all wavelengths and all orbital phases, and can accurately constrain the abundance of H 2 O and CH 4 , as well as retrieve the latitudinally averaged thermal structure of the GCM.
We plot the phase curves generated from the medians of the posterior distributions of the model parameters in Fig. 8.It is clear that model 1, where the dayside 2 width is fixed at 180 • , gives the worst spectral fit to the synthetic data.The phase curve amplitudes 2 We reiterate that the 'dayside' in our models denotes the region of the at- Here we plot the best-fit model phase curves calculated from the posterior medians and compared them to the synthetic data.The synthetic data is shown with the measurement uncertainties of Stevenson et al. (2017).The models are described in section 2.4.retrieved by model 1 are too small at most wavelengths, whereas the other models mostly retrieve the correct phase curve amplitudes.It can be seen in Fig. 10 that model 1 is not flexible enough to approximate the thermal structure of the GCM, which has a hot region significantly narrower than 180 degree in longitude.This also leads model 1 to retrieve a biased high H 2 O abundance, as the model tries to match the amplitudes of the synthetic phase curves by pushing the photosphere higher, where the day/night flux contrast is larger.We thus demonstrate that the dayside area fraction is an important parameter in phase curve retrievals, and the exclusion of its implementation can lead to significant biases in retrieved molecular abundance.
We plot the posterior distributions of gas VMRs using different retrieval schemes in Fig. 9.As mentioned previously, model 1 retrieves biased H 2 O abundance due to the inflexibility of the temperature parameterisation, namely the fixed dayside fraction.The other models produce precise and accurate constraints on H 2 O, as well as accurate upper bounds on CH 4 .None of the models can constrain CO and CO 2 from the data, as CO and CO 2 have weak mosphere modelled by the dayside profile and does not necessarily coincide with the permanently illuminated 'physical dayside'.opacities in the HST/WFC3 wavelengths and their retrieved abundance are mainly driven by the two Spitzer wavelengths.Hence, their abundance are more susceptible to degeneracy with the thermal structure and are therefore poorly constrained.
We plot the retrieved thermal structures, calculated with the median parameters of the posterior distributions, in Fig. 10.We compare the retrieved thermal structures to the appropriate latitudinally-averaged TP structure of the GCM, since we have shown that the appropriately averaged GCM TP structure can reproduce the synthetic phase curves generated directly from the GCM.As described in 3.2, the average is between 0 and 45 degree latitude and with cos(latitude) as the weight.It is now clear that model 1 performs badly because the width of the hot region in the GCM is significantly narrower than the dayside width prescribed by model 1.Model 2 can accurately approximate the typical dayside and nightside TP profiles; however, since model 2 contains a discontinuity at the dayside/nightside boundary, there are large jumps in temperature on isobars around the day/night boundary.Hence, the fits at those regions deviate significantly from the GCM.Model 3 and 4, by virtue of being continuous models, avoids this problem and can approximate the latitudinally averaged GCM structure to well within ± 300 K at most pressures and longitudes.However, both models perform relatively poorly in the deep atmosphere as the retrieved deep atmosphere temperatures are too low (right column, Fig. 10).This leads to incorrectly retrieved phase curve amplitudes at the wavelengths with the deepest photospheres.By computing the transmission weighting function in Fig. 2, we find the three channels at 1.2475 μm, 1.2825μm, 1.3175μm are sensitive mainly to the deep atmosphere at close to 10 bar, whereas most other channels are sensitive to lower pressure levels.In Fig. 8, we can see that the phase curve fits at these three channels by model 3 and model 4 have lower flux than the synthetic phase curves.The retrieved deep atmosphere temperature is biased because there are more data points constraining the TP profile at lower pressure levels, and the TP profile used cannot satisfy all constraints equally well.We expect that the biased deep atmosphere temperature can be resolved by using a more sophisticated TP profile.However, this issue does not lead to significantly biased retrieved abundance and our precision is still in line with past studies.

APPLICATION TO REAL PHASE CURVES
In section 4, we test the performance of our 2D retrieval schemes against synthetic data.We find that we can accurately constrain the abundance of H 2 O and CH 4 from synthetic data simulated from a GCM, provided that the temperature model is flexible enough to approximate the thermal structure of the GCM atmosphere.We now apply the retrieval schemes to the observed HST/WFC3 and Spitzer/IRAC phase curves of Stevenson et al. (2017).We again look at the spectral fits and the retrieved chemical abundance and thermal structure.
We plot the phase curves generated from the medians of the posterior distributions of model parameters in Fig. 11.The model fits to the real data are worse than the model fits to the synthetic data, partly because we do not add random noise to our synthetic data.Furthermore, there are two interesting points of comparison with the fits to synthetic data.First, in the case of real data, model 1 can fit the phase curves almost as well as model 2, whereas in the case of synthetic data, model 1 fits the phase curves markedly worse than model 2. The reason that model 1 cannot fit the synthetic data is because the GCM used to generate the data has a hot region that is significantly narrower in longitudinal extent than the width of dayside region in model 1, which is fixed at 180 • .However, the thermal structure retrieved from the real data by both models is consistent with a hot dayside region that spans approximately 180 • in longitude.This is the reason why model 1 could perform as well as model 2 on the real data, while failing to do so on the synthetic data.Second, model 3 produces the poorest fits to the real data, and often produces phase curve maxima that are too large while simultaneously under-produce the amplitudes of the intermediate phases.Since model 3 prescribes that temperatures on isobars vary sinusoidally with longitude on the dayside, the misfits suggest that temperature must vary less strongly with longitude on isobars.This is confirmed by the retrieved thermal structure of model 4.
We plot the posterior distributions of gas VMRs using different retrieval schemes in Fig. 12.The retrieved H 2 O abundances are consistent across models.We take model 4 as our fiducial model, which gives a constraint of 5.6 × 10 −5 -4.0 × 10 −4 at 1.As with the synthetic data, we cannot constrain the abundance of CO and CO 2 , but we can place an upper bound on the VMR of CH 4 at ∼10 −6 .We plot the retrieved thermal structures, calculated with the median parameters of the posterior distributions, in Fig. 13.The retrieved dayside temperatures of model 4 suggest that the dayside thermal structure of WASP-43b is relatively homogeneous, meaning that temperature does not vary strongly as a function of longitude on isobars.On the other hand, the retrieved temperatures on the nightside are extremely cold, which is likely due to thick cloud coverage that lifts the photosphere to lower pressure levels.We discuss the results of our retrievals in the next section.

DISCUSSIONS
We compare our results to previous retrieval studies of WASP-43b, discuss the effects of nightside clouds, and detail the limitations of our retrieval model in this section.

Comparison with previous retrievals
We present a summary of past retrieval studies of WASP-43b in Table 3.We focus on the studies which analyse the HST transmission, secondary eclipse and phase curve data (GO Program 13467, PI: Jacob Bean, Kreidberg et al. 2014), Spitzer secondary eclipse data (Program ID 70084, Blecic et al. 2017), and Spitzer phase curve data (Programs 10169 and 11001, PI: Kevin Stevenson, Stevenson et al. 2017).We present the constraints on the abundance of H 2 O at 1 for the studies that publish such a result.While all of the H 2 O abundance constraints overlap, our constraint is on the lower end compared to past studies.Additionally, we find two points of discussion.
Firstly, we find that multiple studies support the hypothesis that WASP-43b has no optically-thick clouds on the dayside.Kreidberg et al. (2014) analyse the transmission spectrum from HST/WFC3, and find that the day-night terminator of WASP-43b contains no significant clouds at the pressure levels probed by transmission spectroscopy.Since the dayside is hotter than the terminator region and likely has comparable chemical inventory, this suggests that the dayside may be free from significant cloud coverage as well.More directly, Fraine et al. ( 2021) find a very low dayside geometric albedo (<0.06) using HST WFC3/UVIS secondary eclipse data in the optical wavelengths, and report a non-detection of clouds on the dayside at P > 1 bar.Furthermore, Stevenson et al. (2014) estimate the Bond albedo of WASP-43b to be 0.18 +0.07 −0.12 by computing the day-and night-side bolometric fluxes from the model spectra retrieved from the HST/WFC3 phase curves.These findings provide justifications for our cloud-free retrieval scheme on the dayside, in accordance with the prediction from modelling work that a cloud-free hot spot should dominate the dayside of hot Jupiters (Parmentier et al. 2016).
Secondly, we find that while some studies find variable H 2 O abundance as a function of longitude, this might be the result of the 1D phase-by-phase approach used by these studies.For example, Stevenson et al. (2017) analyse the same HST/WFC3 and Spitzer phase curves as in this work, and apply a 1D phase-by-phase retrieval  .The emission data analysed by Kreidberg et al. (2014) refer to the HST/WFC3 secondary eclipse data as presented in Kreidberg et al. (2014) and the Spitzer secondary eclipse data from Blecic et al. (2014).The transmission data refer to the HST/WFC3 primary transit data as presented in Kreidberg et al. (2014).The phase curves refer to the HST/WFC3 and Spitzer phase curves as presented in Stevenson et al. (2017).We note the constraints on the abundance of H 2 O at 1 for the studies that publish such a result.For the studies that analysed phase curves, we note the retrieval methods.approach, where the spectrum at each orbital phase is retrieved independently and assuming uniform abundance and TP structure for each phase.They find the H 2 O abundance varies between the dayside and the nightside phases, and give a constraint of 2.5×10 −5 -1.1 × 10 −4 at 1 on H 2 O for the nightside and a constraint of 1.4 × 10 −4 -6.1 × 10 −4 at 1 for the dayside.However, the 1D retrieval approach, where uniform atmospheric condition is assumed for the visible atmosphere under observation, is now known to give biased results in the analysis of disc-averaged hot Jupiter spectra (Blecic et al. 2017;Taylor et al. 2020).Furthermore, the phase-by-phase approach is not geometrically self-consistent and under-utilises the constraints from the fact that hemispheres observed at neighbouring phases overlap (Irwin et al. 2020).In Fig. B1 in the appendix, we plot our phase-by-phase retrieval results of the synthetic data, and show that the retrieved H 2 O abundance varies as a function of orbital phase even though the true abundance is uniform across the planet in the GCM used to simulate the synthetic data.In Fig. B2, we plot our phase-by-phase retrieval results of the observed HST/WFC3 and Spitzer phase curves, and find a similar result that the retrieved H 2 O abundance is higher on the dayside than on the nightside as Stevenson et al. (2017).

The influence of nightside clouds
By comparing the synthetic phase curves simulated from the cloudfree WASP-43b GCM with the observed phase curves, we see that the GCM phase curves under-predict the phase curve amplitudes and over-predict the phase curve maximum offsets.The mismatch in phase curve offsets suggest that the strength of heat circulation is weaker on WASP-43b than predicted by the GCM, and the low nightside brightness temperatures further suggest significant nightside cloud coverage.According to Parmentier et al. (2020), when nightside clouds are present, the day-to-night heat transport becomes extremely inefficient, and the nightside photosphere is lifted to higher altitude.This could explain the low nightside temperatures and small phase curve offsets observed on WASP-43b.

Limitations and future work
We have introduced a new 2D retrieval scheme (model 4), where the atmospheric temperature is parameterised by equation ( 12).We now discuss the limitations of our retrieval scheme and directions for future work.

Aerosol model
We do not explicitly model the effects of clouds in our retrieval scheme.Past studies (e.g., Burningham et al. 2017;Mollière et al. 2020) have shown that flexible TP profiles can mimic the spectral contribution of clouds in low-resolution spectroscopy.We have assumed clouds with uniform-with-wavelength spectral features in the observed wavelengths are fully degenerate with thermal structure and chemical abundance.Disentangling this degeneracy is beyond the scope of this work.We recognise that the lack of cloud parameterisation is likely the most significant source of error in our retrieved atmospheric properties, though we expect the retrieved dayside properties are reliable as both transmission spectroscopy and broadband emission observation of WASP-43b find no evidence of clouds on WASP-43b (Kreidberg et al. 2014;Fraine et al. 2021).Recent studies have shown that clouds play an important role in shaping the phase curves of hot Jupiters when they are present on the nightside (Parmentier et al. 2020;Roman et al. 2021), and we plan to include aerosols in our retrieval scheme and validate such scheme with cloudy GCMs in future work.

Temperature model
Our atmospheric temperature model is strongly parameterised to keep the retrieval timescale tractable.We discuss here the limitations of the parameterisation.Firstly, our model is 'two dimensional', meaning that temperature varies with pressure and longitude, but not with latitude.We then interpret the retrieved thermal structure as a latitudinally averaged thermal structure weighted towards the low latitude regions, as described in section 3.2.Irwin et al. (2020) show that the HST + Spitzer phase curves of WASP-43b do not allow the retrieval of the latitudinal variation of atmospheric properties, so the thermal structure of WASP-43b remains poorly constrained.However, JWST phase curves may allow latitudinal variation to be probed, especially when analysed in conjunction with the eclipse mapping technique (e.g., Rauscher et al. 2018).The joint analysis of eclipse mapping data and phase curves will provide the most detailed constraints on the 3D structure of hot Jupiter atmospheres.We plan to upgrade our current 2D model to include latitudinal variation in order to analyse phase curves and eclipse maps jointly in our future work, in particular the JWST/MIRI data of WASP-43b.
Secondly, we assume both north-south symmetry about the equator and east-west symmetry about the dayside central meridian in our temperature model.The assumption of north-south symmetry is based on the GCM of Parmentier et al. (2016) described in section 2.1, which exhibits negligible differences between the northern and southern hemispheres.The symmetry between the northern and southern hemispheres has been seen in other hot Jupiter GCM studies as well (Amundsen et al. 2016;Roman et al. 2021;Mendonça 2020).However, the simulations of Cho et al. (2021) suggest that the atmospheres of hot Jupiter may be highly turbulent, and the atmospheric thermal structures may exhibit significant time variability that breaks the north-south symmetry.Repeated observations of hot Jupiters can detect such time-variation.Since we have not validated our retrieval scheme on such turbulent atmospheres, our temperature model would require further testing when time-variation is present.The east-west symmetry, on the other hand, limits the flexibility of our model to capture certain thermal structures seen in GCMs, for example, temperature decreasing faster with longitude in the westward direction than in the eastward direction, or pressure-dependent hot spot shift (when the hottest hemispheres in each pressure level do not align).In the case of retrieving synthetic phase curves, these two effects only mildly limit our spectral fits, as shown in Fig. 8. Furthermore, as seen in Fig. 9 and Fig. 10, the chemical abundance and thermal structure retrieved are not significantly biased, and are in line with literature results in terms of precision (e.g., Feng et al. 2020;Irwin et al. 2020).Nevertheless, we expect such a approach would not be appropriate for data of higher quality than the ones we have analysed.As we plan to apply our retrieval schemes to JWST-quality data, particularly the MIRI observation of WASP-43b, we plan to make modifications to our schemes to include such secondary structures.

Chemistry model
We assume that the chemical abundance of gas species are constant with location and constant with pressure in our atmospheric model as we focus on the parameterisation of atmospheric temperature in this work.The validity of this assumption depends on the gas species, the characteristics of the planetary atmosphere and the pressure range we are interested in.The chemical and dynamical modelling work of Cooper & Showman (2006) and Agúndez et al. (2014) suggest that for H 2 O and CO, which are predicted to be the most abundant spectrally active molecules on hot Jupiters in the pressure ranges probed by low-resolution spectroscopy (∼10-10 −3 bar), the constant abundance assumption is valid as atmospheric circulation effectively homogenises their abundance.The assumption holds less well for CO 2 , though we do not expect this to be a significant source of error, as the variation in CO 2 abundance in the models of Agúndez et al. (2014) is less than one order of magnitude.The case of CH 4 is most problematic: although its abundance shows negligible variation with longitude, the vertical variation is significant (Agúndez et al. 2014).More recently, Baeyens et al. (2021) calculated a grid of pseudo-2D chemistry models for hot Jupiter atmospheres, and their results enable us to directly assess the constant chemistry assumption for a WASP-43b-like atmosphere (Figure 18, Baeyens et al. 2021).Their model suggests that the variations of H 2 O, CO, and CO 2 abundance are well within one-order of magnitude in the pressure range 100-10 −4 bar, both with respect to longitude and with respect to pressure.This pressure range should more than cover the pressure range probed by low-resolution infrared emission spectroscopy of WASP-43b.Since typically even the best molecular abundance constraints from HST + Spitzer data have uncertainties around one order of magnitude, we conclude that for H 2 O, CO, and CO 2 it is valid to assume constant chemistry for HST + Spitzer data based on the cited modelling work.As for CH 4 , the model of Baeyens et al. (2021) suggests that for a WASP-43b like atmosphere while the CH 4 abundance is significant in the deep atmosphere (at pressure level greater than about 1 bar), its abundance rapidly decreases with decreasing pressure.If this is the case, then our retrieved upper bound of CH 4 abundance at around 10 −6 would not reflect the true CH 4 abundance of the atmosphere, which would be highly pressure-dependent.Apart from the interplay between circulation and equilibrium chemistry, photochemistry and molecular dissociation can also affect the spatial distribution of molecular abundance.However, we only expect dissociation to be significant in much more strongly irradiated planets than WASP-43b, and we expect photochemical products to be insignificant in the pressure region probed by emission spectroscopy.For JWST data, we expect that the constant abundance approximation could still be valid for H 2 O, CO, and CO 2 when analysing emission spectroscopy of hot Jupiters, based on the current modelling work.However, more sophisticated parameterisation of abundance variation is necessary for CH 4 , and for joint analysis of transmission and emission data where a large pressure range is probed.

CONCLUSIONS
We propose a novel 2D temperature parameterisation for the retrievals of hot Jupiter phase curves, which is described by equations ( 12).The temperature model is a function of pressure and longitude, and can be used to retrieve the chemical composition and latitudinally averaged thermal structure of hot Jupiters atmospheres from phase curves.The model is built on two TP profiles, signifying the representative profiles for the dayside and the nightside.In our model, the temperature is uniform on isobars on the nightside, and varies with cos  (longitude/) on isobars on the dayside, where  and  are free parameters.Both the dayside central longitude and dayside fraction (longitudinal extent) are free parameters of our model.We first apply our proposed retrieval scheme, together with several other 2D models for comparison, to synthetic phase curves simulated from a cloud-free GCM of WASP-43b, representing a more realistic atmospheric model than the typically simple models used for validating retrieval schemes in the literature.We find that the models that allow variable dayside longitudinal extent can fit synthetic HST/WFC3 and Spitzer/IRAC phase curves to within typical measurement uncertainties, as well as accurately and precisely constraining the water and methane abundance.The retrieved ther-mal structures using these models are good approximations to the latitudinal-average of the GCM thermal structure weighted towards the low latitude regions.We then apply our retrieval schemes to retrieve information from the observed phase curves of WASP-43b presented in Stevenson et al. (2017).We constrain the abundance of water to be 5.6 × 10 −5 -4.0 × 10 −4 at 1 using model 4, as well as an upper bound on CH 4 at ∼ 1 × 10 −6 .We find that the latitudinally averaged dayside TP structure of WASP-43b is likely to be homogeneous (meaning that temperature does not vary strongly as a function of longitude on isobars) and non-inverted.We expect the nightside of WASP-43b to be covered by thick clouds due to the extremely low retrieved nightside temperature, in agreement with previous studies.We have thus demonstrated the efficacy of our retrieval scheme, which simultaneously fits all orbital phases of a set of phase curves at a modest computation cost.
as new developments, including the implementation of several 2D retrieval schemes for analysing exoplanet phase curves.The most computationally expensive routines are compiled to machine code at runtime using a just-in-time (JIT) compiler 3 , so that the speed of our code is on par with compiled languages.

APPENDIX B: COMPARISON WITH PHASE-BY-PHASE RETRIEVALS
We show the retrieval results using the 1D phase-by-phase approach, where the spectrum at each orbital phase is independently analysed with a uniform TP and abundance profile.The TP profile used for the 1D retrievals is the same Guillot profile that we use for the 2D retrievals.In Fig. B1, we present the retrieved molecular abundance from the synthetic phase curves, where the truths are marked by horizontal black lines.We see that the retrieved H 2 O abundance using the 1D phase-by-phase approach varies significantly with orbital phases, and is biased at several orbital phases.We echo the findings of past studies (Blecic et al. 2017;Taylor et al. 2020) that the 1D phase-by-phase approach could lead to significantly biased molecular abundance.In Fig. B2, we present the retrieved molecular abundance from the real phase curves presented in Stevenson et al. (2017).We find a similar trend as Stevenson et al. (2017), that H 2 O abundance is higher for the dayside phases than for the nightside phases.We plot the abundance constraints from model 4 on the figures for comparison.We can see from Fig. B1 that model 4 performs markedly better than the 1D approach on synthetic phase curves both in terms of accuracy and precision in retrieved molecular abundance.

APPENDIX C: FULL POSTERIOR PLOTS
We include the full posterior distributions of all models for the retrievals of the synthetic data.We additionally include the full posterior distribution of model 4 for the retrievals of the observed data.2.     2.
Figure1.Temperature (Kelvin) as a function of longitude and latitude at three pressure levels in the WASP-43b GCM.The super-rotating equatorial jet is clearly visible and shifts the 'hot spot' eastward of the substellar point (where the star would be perceived to be directly overhead).Note that the substellar point is at 0 degree longitude.Such jet-like features would cause the phase curve amplitudes to peak before secondary eclipses.Note that the latitudinal distance is weighted by cos(latitude) to mimic the effect that polar latitudes would appear foreshortened to us because we observe WASP-43b from above the equator.

Figure 2 .
Figure 2. Transmission weighting function of the WASP-43b GCM atmosphere at the substellar point as a function of pressure and wavelength channel.The first 15 HST/WFC3 channels are in equally spaced bins of width 0.035 μm spanning the wavelength range 1.1425-1.6325μm, and the last two Spitzer/IRAC broad channels are centred at 3.6 and 4.5 μm, respectively.

Figure 4 .
Figure 4. Schematics of model 1 and model 2. Model 1 (top panel) is defined by equation (9) and divides the atmosphere into a dayside and a nightside.Each region is then modelled with a single representative TP profile.The dayside central longitude  is allowed to vary, and the dayside width (longitudinal extent) is fixed to be 180 • .Model 2 (bottom panel) is defined by equation (10) and generalises model 1 by allowing the dayside width to vary, which now spans 180 • ×  in longitude.

Figure 5 .
Figure 5. Temperature as a function of pressure and longitude in two examples of model 3 and model 4. Model 3 (top panel) is defined by equation (11) .Model 4 (bottom panel) is defined by equation (12); in this example  is set to be 1.75.Note that model 4 is equivalent to model 3 if  is being set to 1.

Figure 6 .
Figure 6.Comparison of the phase curves simulated from simplified models to the phase curves simulated from the original WASP-43b GCM (black), where abundance is set by chemical equilibrium.The blue curves are simulated from the best-fit Guillot profiles.The green curves are simulated from latitudinally-averaged TP profiles from 0 • to 45 • , using cos(latitude) as the weight.The yellow curves are simulated with uniform abundance listed in Table1, which are the synthetic data we retrieve on in section 4. The orange curves are simulated with the abundance listed in Table1, but with the methane abundance multiplied by 40, which illustrates the fact that models with uniform gas abundance can match the phase curves produced from a chemical equilibrium model.Note that the error bars on the GCM phase curves are the estimated observational uncertainties ofStevenson et al. (2017).

Figure 7 .
Figure7.Top panel: temperature as a function of longitude and pressure in the WASP-43b GCM at the equator.Lower panel: latitudinally averaged TP profiles from 0 to 45 degree latitude, using cos(latitude) as the weight.Since the retrieved thermal structure resemble the latitudinally averaged profiles, we would underpredict the hot spot offset at the equator for typical hot Jupiter GCMs from synthetic data.

Figure 8 .
Figure 8. Results from the retrievals of synthetic phase curves generated from the WASP-43b GCM.Here we plot the best-fit model phase curves calculated from the posterior medians and compared them to the synthetic data.The synthetic data is shown with the measurement uncertainties ofStevenson et al. (2017).The models are described in section 2.4.

Figure 9 .
Figure9.Results from the retrievals of synthetic phase curves generated from the WASP-43b GCM.Here we plot the posterior distributions of the retrieved gas VMRs using different retrieval schemes, and we mark the abundance used to simulated the synthetic data with black lines ('truths').

Figure 10 .
Figure 10.Results from the retrievals of synthetic phase curves generated from the WASP-43b GCM.The rows from top to bottom correspond to model 1, model 2, model 3, and model 4, respectively.Left column: retrieved temperature structures, which are calculated using the median parameters of the posterior distributions.Middle column: latitudinally averaged TP profile of the GCM from 0 to 45 degree latitude, using cos(latitude) as the weight.Right column: difference between right and middle columns.

Figure 11 .
Figure 11.Retrieval results of the real phase curves.The retrieved phase are calculated from the posterior medians and compared to the real data.The models are described in 2.4.

Figure 12 .
Figure 12.Retrieval results of the real phase curves.Posterior distributions of the retrieved gas VMRs.

Figure 13 .
Figure 13.Retrieval results of the real phase curves.The rows from top to bottom correspond to model 1, model 2, model 3, and model 4, respectively.Left column: retrieved temperature structures, which are calculated using the median parameters of the posterior distribution.Middle column: latitudinally averaged TP profile of the GCM from 0 to 45 degree latitude, using cos(latitude) as the weight.Right column: difference between right and middle columns.

Figure B1 .Figure B2 .
Figure B1.Retrieved molecular abundance from the synthetic phase curves using the 1D phase-by-phase approach (blue), compared to the retrieved abundance using model 4 (red).The vertical lines mark the 1  confidence intervals.The true abundances used to generate the data are marked by the black horizontal lines.

Figure C1 .
Figure C1.Full posterior distribution of the model parameters of model 1 for the retrievals of the synthetic phase curves.The parameters are summarised in Table2.

Figure C2 .
Figure C2.Full posterior distribution of the model parameters of model 2 for the retrievals of the synthetic phase curves.The parameters are summarised in Table2.

Figure C3 .
Figure C3.Full posterior distribution of the model parameters of model 3 for the retrievals of the synthetic phase curves.The parameters are summarised in Table2.

Figure C4 .
Figure C4.Full posterior distribution of the model parameters of model 4 for the retrievals of the synthetic phase curves.The parameters are summarised in Table2.

Figure C5 .
Figure C5.Full posterior distribution of the model parameters of model 4 for the retrieval of the observed phase curves.The parameters are summarised in Table2.

Table 2 .
Parameters of our atmospheric models.All models have 4 parameters for gas VMRs, with the other parameters specifying the thermal structure.Model 1 has 13 parameters in total, whereas model 2 and model 3 have 14 parameters in total.Model 4 has 15 parameters in total.

Table 3 .
Summary of previous retrieval studies of