Two-dimensional Eclipse Mapping of the Hot-Jupiter WASP-43b with JWST MIRI/LRS

We present eclipse maps of the two-dimensional thermal emission from the dayside of the hot-Jupiter WASP-43b, derived from an observation of a phase curve with the JWST MIRI/LRS instrument. The observed eclipse shapes deviate significantly from those expected for a planet emitting uniformly over its surface. We fit a map to this deviation, constructed from spherical harmonics up to order ℓmax=2 , alongside the planetary, orbital, stellar, and systematic parameters. This yields a map with a meridionally averaged eastward hot-spot shift of (7.75 ± 0.36)°, with no significant degeneracy between the map and the additional parameters. We show the latitudinal and longitudinal contributions of the dayside emission structure to the eclipse shape, finding a latitudinal signal of ∼200 ppm and a longitudinal signal of ∼250 ppm. To investigate the sensitivity of the map to the method, we fix the parameters not used for mapping and derive an “eigenmap” fitted with an optimized number of orthogonal phase curves, which yields a similar map to the ℓmax=2 map. We also fit a map up to ℓmax=3 , which shows a smaller hot-spot shift, with a larger uncertainty. These maps are similar to those produced by atmospheric simulations. We conclude that there is a significant mapping signal which constrains the spherical harmonic components of our model up to ℓmax=2 . Alternative mapping models may derive different structures with smaller-scale features; we suggest that further observations of WASP-43b and other planets will drive the development of more robust methods and more accurate maps.


INTRODUCTION
Eclipse mapping measures the two-dimensional spatial features of an object when it is eclipsed by another object.The eclipsed object is covered along one axis during the eclipse ingress, and is revealed along another axis during the eclipse egress.Combining the information from these two axes reveals two-dimensional information about the surface of the eclipsed object (de Wit et al. 2012;Majeau et al. 2012).
Eclipse maps have been derived for objects like Pluto during its eclipse by Charon (Stern 1992), or the white dwarf BD +16 • 516B during its eclipse in a binary system (Warner et al. 1971).Williams et al. (2006) proposed applying this technique to exoplanets, and Rauscher et al. (2007) showed that the James Webb Space Telescope (JWST) would have sufficient precision to accurately map exoplanets.Eclipse mapping is currently the only method by which 2D information can be measured for exoplanets, as out-of-eclipse "phase curves" provide lowresolution information as a function of longitude only.2D spatial information is crucial for understanding atmospheric circulation (Lewis & Hammond 2022), chemical composition (Taylor et al. 2020;Yang et al. 2023), and cloud structure (Parmentier et al. 2016).
de Wit et al. (2012) and Majeau et al. (2012) derived eclipse maps of the hot Jupiter HD 189733b using observations with the Spitzer Space Telescope.Their analyses were restricted to fitting large-scale shapes with a high degree of uncertainty due to the limited precision of the measurements, uncertainty in and degeneracy with orbital parameters, and the dependence of the result on the mapping model.de Wit et al. (2012) described the uncertainty in the mapping signal due to the impact parameter, stellar density, eccentricity, and argument of periastron in particular.Rauscher et al. (2018), Mansfield et al. (2020), and Challener & Rauscher (2022) developed more advanced methods to fit eclipse maps of 2D and 3D structure, based on fitting orthogonal phase curves rather than orthogonal surface maps, and applied these methods to the previous observations of HD 189733b.Coulombe et al. (2023) presented the first eclipse map measured with JWST, deriving a map of thermal emission from observations of the hot Jupiter WASP-18b from 0.85 to 2.85 µm with the NIRISS/SOSS instrument.They found no longitudinal shift in the hotspot, sharp gradients in brightness towards the terminator, and no clear measurement of latitudinal structure.
This study presents eclipse maps derived from broadband JWST MIRI/LRS (Kendrew et al. 2015) observations from 5 to 10.5 µm of a full phase curve of the hot Jupiter WASP-43b containing two eclipses and one transit (Bell et al. 2024).WASP-43b is a "hot Jupiter" exoplanet with strong thermal emission from its permanent dayside, and a short, tidally-locked orbit around a K7 main sequence star that exhibits low variability (Hellier et al. 2011;Scandariato et al. 2022).These properties make it ideal for precise time-series observations of thermal emission.Bell et al. (2024) analysed this MIRI/LRS phase curve, finding a large difference in dayside and nightside brightness temperatures and evidence for water absorption.We fit its inclination to be 82.11+0.050  −0.052 , and the ratio of its semi-major axis to stellar radius to be 4.859 +0.013  −0.012 , which corresponds to an impact parameter of 0.667 +0.006  −0.006 .This means that the edge of the star crosses the planet at an angle of ∼42 • (the "stellar edge angle" defined in Boone et al. (2023)).The longitudinal and latitudinal features of the day-side therefore affect the eclipse mapping signal almost equally, with a small bias towards longitudinal features.This geometry makes WASP-43b especially well suited for eclipse mapping.
We fit an eclipse map to this dataset, finding the clearest eclipse mapping signal to date and separating the latitudinal and longitudinal parts of this signal for the first time.In Section 2 we describe the methods used to fit the eclipse maps.Section 3 then shows the different maps we fit to the observations, and compares their structures and statistical evidence.These maps are then compared to numerical simulations in Section 4. In Section 5 we conclude that the spherical harmonic components up to ℓ max = 2 are constrained well, but that further observations are needed to precisely constrain higher-order mapping structures.

DATA AND METHODS
A full orbit of WASP-43b with the JWST MIRI/LRS slitless (SL) mode was observed as part of the JWST-ERS-1366 program, performing target acquisition with the F1500W filter and using the SLITLESSPRISM subarray for the science observation.The science observation lasted 26.5 hours, consisting of 9316 integrations lasting 10.34 s each.The two eclipses show a downwards systematic trend with time, similar to that identified in other MIRI/LRS time-series observations (Bouwman et al. 2023).We discard the data beyond 10.5 µm due to the "shadowed region effect" described in Bell et al. (2024) alongside a more detailed analysis of the data acquisition.
Figure 1 shows the raw MIRI/LRS dataset (black points) in units of planetary flux divided by stellar flux.We use the Eureka!pipeline (Bell et al. 2022) to reduce the observed data, starting from the Eureka!Stage 4 output from the "Eureka v1" reduction in Bell et al. (2024).We trim the initial 780 integrations (not plotted) where instrumental systematic effects are the strongest.The only subsequent methodological difference to Bell et al. (2024) is the method used to model the planetary emission in the phase curve -we used a starry model of a planet with a 2D emission map over its surface (Luger et al. 2019), instead of a Fourier series model of the phase curve.

Astrophysical Model
We fit a 2D map of thermal emission to this phase curve, simultaneously with the orbital, planetary, stellar, and instrumental systematic parameters.Given the significant effect of the planetary emission map on the observed eclipse shape, this should derive more accurate parameters than those derived using a Fourier series model in Bell et al. (2024).This simultaneous fit also tests if there are any degeneracies between the system parameters and the eclipse map, such as between the eclipse timing and the longitudinal emission offset (Williams et al. 2006).
We fit the planet-to-star radius ratio, the linear ephemeris, the inclination, the ratio of the semi-major axis to the stellar radius, and the two parameterised stellar limb-darkening parameters in Kipping (2013).
We set the orbital period constant at 0.813474 days as it is known with sufficient precision already (Kokori et al. 2023).We set the obliquity of the planet to be 0 as we assume a tidally locked orbit.We set the stellar radius to be constant at 0.665 R ⊙ (Bell et al. 2024), although this is an arbitrary value that serves only to give our model a dimensional form, as the ratio F p /F S is only sensitive to the ratios R p /R * and a/R * which we fit separately.The stellar mass does affect the signal (de Wit et al. 2012) but is entirely determined in our model by the fitted value of a/R * and the fixed value of the orbital period.
de Wit et al. (2012) showed how orbital eccentricity can be degenerate with an eclipse mapping signal.Gillon et al. (2012) constrained the eccentricity of WASP-43b to be 0.0035 +0.0060 −0.0025 , which has generally been used to assume zero eccentricity for this planet (Bell et al. 2024).To check this assumption, we separately fitted the model including the eccentricity and periastron, finding an eccentricity of 0.001078 0.00185 −0.00025 , and no meaningful constraint on periastron, as discussed in Appendix B. We take this to verify the assumption of zero eccentricity, and proceed to fit the models with eccentricity set to zero.This is consistent with the circularization timescale of 3 Myr that we estimate using Adams & Laughlin (2006) (assuming Q P = 10 6 ), which is much shorter than the age of the system estimated by Hellier et al. (2011) to be 400 +200 −100 Myr.For the instrumental systematic effects, we fit a uniform baseline and linear trend in time, the magnitude and timescale of an exponential ramp, and trends of the spatial position and width of the data on the detector.Table 1 shows the new parameters fitted for the ℓ max = 2 eclipse map model, compared to the parameters fitted with an n = 2 Fourier series model (like in Bell et al. (2024)).
The difference in the time it takes light to travel from the planet and the star has an important effect on this dataset, so we include it in our model.In eclipse, the light from the planet takes about 8.5 seconds longer to reach the observer than the light from the star.For the average gradient in eclipse ingress and egress of roughly 8 ppm per second (see Figure 1), this corresponds to a maximum effect of around 70 ppm on the eclipse shape.Figure 2 shows that this corresponds to about 15% of the maximum deviation in eclipse shape due to the eclipse mapping signal itself.

Eclipse Mapping
The eclipse mapping method works by constructing a map from a number of "basis maps", which each have an associated basis phase curve.We match the observed light curve with a sum F (t) constructed from the the basis phase curves f i (t) weighted by coefficients c i : (1) Figure 2. The observed data in the ingress and egress of the eclipses, subtracted by a phase curve fitted using an n = 2 Fourier series model.We fitted this Fourier series model using the orbital parameters derived using an eclipse map in Table 1, so that the residual signal depends entirely on the different models.This leaves a residual "eclipse mapping signal" showing the effect of partial stellar coverage of the non-uniform planetary emission.Green points show the first eclipse, purple points show the second eclipse, and black points show their mean.The blue shaded region labelled "Latitude-Longitude Map Fit" shows the range of fitted phase curves from the ℓmax = 2 eclipse map in Figure 3, with the two shaded regions showing the first and second quantiles, containing 68.27% and 95.45% of the posterior distribution.The red shaded region labelled "Longitude-Only Map Fit" shows the fitted eclipse of the ℓmax = 2 eclipse map with flat latitudinal structure shown in Figure 4.The poor fit of this model shows the presence of latitudinal information in the dataset, and the need to fit the latitudinal structure of the map.This "longitude-only" residual is around 250 ppm at its largest, while the "latitude-only" signal (estimated from the difference between the 1D and 2D map fits) is around 200 ppm.The total eclipse mapping signal is around 450 ppm at its largest.
We derive the actual eclipse map Z(θ, ϕ) from the fitted coefficients c i , as a sum of the spherical harmonic basis maps z i (θ, ϕ): (2) Luger et al. (2019) shows the phase curves f i (t) of each spherical harmonic z i (θ, ϕ) (where θ and ϕ are longitude and latitude).We model the astrophysical signal using starry1 (Luger et al. 2019) and fit the coefficients c i using PyMC3 (Salvatier et al. 2016).We use "pixel sampling" to enforce a positive emission map globally (e.g. as used in Gorski et al. (2005)).This method samples the brightness of pixels distributed uniformly over the mapped surface and transforms them to spherical harmonic coefficients c i to compute the actual phase curve.10, for the Fourier series fit described in Section 3 and the ℓmax = 2 eclipse map fit plotted in Figure 3.The limb darkening parameters q1 and q2 are as described in Kipping (2013).The polynomial parameters C0 and C1 describe the constant baseline and linear trend with time, as used in Bell et al. (2022).The ramp parameters describe the magnitude and timescale of a linearly decaying exponential r0e (−r 1 t) (where t = 0 at the start of the observation) as used in Bell et al. (2022).The spatial position and PSF width are as described in Bell et al. (2022).The uncertainty scaling describes a multiplicative parameter to inflate the expected errors to be consistent with the residual between the data and the fitted model.
To fit the ℓ max = 2 eclipse map alongside the orbital and systematic parameters, we sample 16 pixels evenly spaced on a Mollweide projection to represent the spherical harmonic space with 4ℓ 2 pixels (McEwen & Wiaux 2011).We use a (natural) log-normal prior for the pixels with a mean magnitude of 6000 ppm and a standard deviation of 3000 ppm (transformed to log-space).
After fitting the pixels representing the ℓ max = 2 eclipse map simultaneously with the orbital, planetary, stellar, and systematic parameters, we fix these additional parameters to their derived values, and re-fit the map with a variety of methods.We fit maps with: 1. ℓ max = 3 spherical harmonics 2.An eigenmap as described in Rauscher et al. (2018) 3. ℓ max = 2 spherical harmonics to the first eclipse alone, the second eclipse alone, and both eclipses The next section presents the results from each of these re-fitted maps, and describes how the details of each fitting procedure affect the derived map.It also presents the statistical evidence for an eclipse mapping signal, and compares the different mapping models.

Eclipse map fitted with orbit and systematics
The red line in Figure 1 shows the phase curve resulting from the ℓ max = 2 eclipse map fitted alongside the orbital parameters, stellar parameters, and systematic model.The blue line in Figure 1 shows the model of instrumental systematics fitted at the same time.We normalised the F p /F S values (the ratio of planetary to stellar flux) in Figure 1 so that the systematic model has zero mean over the course of the fitted observations.
Appendix B shows the posterior distribution for the fitted orbital, planetary, stellar, and systematic parameters.Some of the parameters are consistent with those derived using a Fourier series fit like that used in Bell et al. (2024), but some are not consistent.There are small but statistically significant differences in the planetary radius, eclipse timing, and inclination, all of which we expect to be more accurately fitted by the more realistic eclipse shape in the eclipse mapping model.The posterior distribution of the stellar limb darkening parameters are different but we found that, when combined, these resulted in consistent limb darkening profiles.The stellar limb darkening is relatively weak at these long wavelengths, so is poorly constrained but is consistent between the two models.The systematic parameters are different to those produced by a Fourier series fit (Bell et al. 2024) but result in an almost identical model of systematics, as a short exponential ramp plus a linear trend is almost identical to a long exponential ramp.
Figure 1 shows that this ℓ max = 2 eclipse map model fits the observed phase curve well.The information about the 2D eclipse map structure is contained within the shapes of the ingress and egress of each eclipse.Fitted max = 2 map, orbit, and systematics Figure 3.An ℓmax = 2 eclipse map fitted to the data in Figure 1.First panel: an eclipse map constructed with spherical harmonics up to ℓmax = 2, fitted simultaneously with the orbital, planetary, stellar, and systematic parameters, via "pixel sampling" as described in Section 3.1.Note that the nightside of the map contains no latitudinal information, with longitudinal information from the phase curve only (out of eclipse).The plotted 2D map uses the median of the posterior distribution of each fitted spherical harmonic coefficient.Second panel: the posterior distribution of the longitudinal structure along the equator of the map, showing the first and second quantiles.Third panel: the posterior distribution of the latitudinal structure through the substellar point.Fitted max = 2 map with uniform latitude structure, with fixed orbit and systematics Figure 4.An eclipse map fitted to the dataset in Figure 1, using the parameters listed in Table 1 derived with the ℓmax = 2 eclipse map model, and constructed from ℓmax = 2 spherical harmonics averaged in latitude.This is distinct from the "Fourier series" phase curve model used in Table 1, which models the phase curve directly and does not represent the partial coverage of the map in eclipse.The fit to the ingress and egress data is shown in red in Figure 2; as described in Section 3.2, this map cannot accurately fit the observed eclipse shape, showing the presence of latitudinal information in the observations.purple points), and a "control" model representing the phase curve as a Fourier series up to order n = 2.This Fourier series was fitted with orbital and systematic parameters fixed to the values derived when fitting the ℓ max = 2 eclipse map, listed in Table 1, so that the residual is only due to the different emission models.
The green and purple points in Figure 2 show the residual of the observed data in the first and second eclipses, binned every 8 points, and the black points show the average of the two eclipses.The ℓ max = 2 eclipse mapping model and the n = 2 Fourier series model fit the out-of-eclipse phase curve as well as each other, so the residual outside the range plotted in Figure 2 is determined by the uncertainty in the data.This is not surprising, as both models have access to almost identical fitting functions for the out-of-eclipse phase curve -the ℓ max = 2 harmonics produce the same phase curves as the n = 2 Fourier series, apart from a small modification from the orbital obliquity.
The size of the residual between the observed data and the n = 2 control model shows the size of the eclipse mapping signal, which is the effect of partial stellar coverage of non-uniform emission from the planet (de Wit et al. 2012).For example, when the star covers an area that is emitting more than the average of the planetary disk, the observer measures less flux than would be expected for a uniform disk, so the residual signal is negative.
The blue shaded area shows the range of phase curves fitted using the ℓ max = 2 eclipse map model, with the Fourier series control model subtracted.The two shaded regions show the first and second quantiles, containing 68.27% and 95.45% of the posterior distribution respectively.The eclipse map model matches the residual signal in ingress and egress, with a uncertainty comparable to the error bars of the observed data.This implies a good fit to a robust eclipse mapping signal, which we quantify in the next section.
Figure 3 shows the ℓ max = 2 eclipse map itself.The plotted 2D map is constructed from the median values of the posterior distribution of each of the spherical harmonic coefficients.This median map has a hotspot near the substellar point, with a small meridionally averaged shift of (7.75 +0.36  −0.36 ) • eastwards and a shift of (10.72 +4.14  −4.68 ) • southwards2 .We present the maps as a ratio of planetary flux to stellar flux, as a conversion to temperature is more complex and requires modelling or assumptions about the temperature structure.We make a simple estimate of the brightness temperature assuming planetary blackbody emission B T,λ in the bandpass λ = 5 to 10.5 µm, and a PHOENIX (Allard & Hauschildt 1995;Hauschildt et al. 1999;Husser et al. 2013) stellar model spectrum with effective temperature T eff = 4300 K and surface gravity log g = 4.50 as used in Bell et al. (2024).The stellar radius that we assumed to be 0.665 R ⊙ does affect the value of the brightness temperature.We do not include effect of the uncertainty of this and other stellar parameters on the derived brightness temperature, in order to highlight the uncertainty due to mapping alone.This determines that this ℓ max = 2 map corresponds to a brightness temperature of (1790.0+23.0 −29.0 ) K at the substellar point, (1293.0+27.0 −34.0 ) K at the equatorial east terminator, (1114.0+30.0 −36.0 ) K at the equatorial west terminator, (1100.0+128.0 −109.0 ) K at the south pole, and (1011.0+111.0 −108.0 ) K at the north pole (where the distinction between the poles is our arbitrary choice).
The posterior distribution of the fitted maps is shown east-west along the equator and north-south through the substellar point.This range of fitted maps includes any degeneracies with the orbital or systematic parameters as these are fitted simultaneously.We confirmed that the "median map" plotted in 2D corresponds closely to the center of the east-west and north-south posterior distributions (not shown).The maximum likelihood map is very similar to the median map, except near the poles where its latitudinal structure deviates slightly from the median but remains within the first quantile.
Figure 4 shows an ℓ max = 2 eclipse map fitted to the same dataset, but with all of its fitting basis maps averaged in latitude, removing the ability of the map to fit latitudinal structure.This makes some of the basis maps uniform everywhere, so we remove them from the fitting process (reducing the number of parameters).The resulting map is very tightly constrained by its need to fit the out-of-eclipse phase curve.Its residual signal in ingress and egress is shown by the red region in Figure 2, which fails to match the observed residual points.This shows the presence of a latitudinal mapping signal of around 200 ppm (the difference between this longitudeonly fit, and the observed residual), compared to the longitudinal mapping signal of around 250 ppm (the magnitude of the residual of this longitude-only fit).

Eclipse Mapping Signal
We can quantify the significance of the signal in Figure 2 by comparing the Bayesian evidence of the ℓ max = 2 eclipse map model (the blue shaded region) and the n = 2 Fourier series control model (the zero line) in Figure 2. As applied in Placek et al. (2017) to a similar analysis of exoplanet phase curves, the log-likelihood function is: for data D i at N times t i fitted by a model F (t) with variance σ 2 .The log-odds ratio ln O = ln L 1 −ln L 2 then represents the relative performance of models 1 and 2, with a value of ln O above 1 corresponding to a better fit to model 1.The log-odds ratio for the ℓ max = 2 map model compared to the Fourier series model (with the eclipse shape of a uniformly emitting disk) is 52.6, showing that the eclipse model is overwhelmingly preferred according to the criterion in Placek et al. (2017).
Both models fit the out-of-eclipse data very similarly as they have access to almost identical out-of-eclipse fitting functions, so this difference in likelihood is entirely due to the improved fit in the ingress and egress of the eclipses shown in Figure 2. Later, we use the Bayesian Information Criterion (BIC) to compare the fit quality of the eclipse map models while taking the number of fitted parameters into account.We cannot use the BIC to compare a map model to the Fourier series model as they are not nested.Notably, the log-odds ratio of an n = 2 Fourier series fit with free orbital, planetary, stellar, and systematic parameters (used to derive the parameters in Table 1) is 20.7 compared to the ℓ max = 2 map model.This shows how the non-mapping parameters are modified by the Fitted max = 2 eigenmap with 6 components, with fixed orbit and systematics Figure 5. Two alternative mapping methods to Figure 3. First row: an eclipse map using spherical harmonics up to order ℓmax = 3, fitted using constant systematic and orbital parameters derived by the map fitted in Figure 3 (listed in Table 1) The median map, longitudinal structure, and latitudinal structure are laid out as in Figure 3. Bottom row: an eclipse map fitted with eigenmaps up to order ℓmax = 3 as described in Section 3.4, for six basis maps selected by optimising the BIC, using constant systematic and orbital parameters derived by the fit in Figure 3.
Fourier series model to slightly better fit the ingress and egress shape.However, this makes a minor difference overall and the eclipse map model is still very strongly preferred.
Comparing the eclipse map and Fourier series models confirms that there is a strong eclipse mapping signal in this dataset.This does not immediately imply the presence of 2D information in the dataset, as Coulombe et al. (2023) only found evidence for longitudinal information in their eclipse map of WASP-18b.We can show that latitudinal information is present by comparing the 2D ℓ max = 2 map to the "1D" ℓ max = 2 map where all of the basis maps are averaged in latitude, shown in Figure 4.The red region in Figure 2 shows the residual eclipse mapping signal for this 1D map, which does not fit the real dataset as well as the 2D map.
The log-odds ratio of the 2D map to the 1D map fits in Figure 2 is 9.8, showing that the 2D map with variable latitudinal structure is greatly preferred.Notably, the uncertainty on the residual of the 1D map in Figure 2 is very small, as the longitudinal structure is well constrained by the additional information from the out-of-eclipse phase curve.This implies that most of the uncertainty on the 2D model in Figure 2 is related to uncertainty about the latitudinal structure of the map.Coulombe et al. (2023) conducted the same comparison between 1D and 2D eclipse maps of WASP-18b and found no evidence that the 2D map was preferred over the 1D map.They therefore concluded that the eclipse mapping signal revealed longitudinal information only.In contrast, the eclipse mapping signal in this study of WASP-43b is sensitive to both longitudinal and latitudinal information, as would be expected from the higher impact parameter of the orbit.This is therefore the first time that latitudinal information has been shown to be detected on an exoplanet.

Eclipse map fitted with ℓ max = 3 spherical harmonics
The spherical harmonic order of the ℓ max = 2 map fit to the dataset was limited by the number of parameters that could be fitted simultaneously with the orbital, stellar, and systematic parameters.Figure 10 in Appendix B shows the posterior distribution for these parameters, demonstrating that there are no significant degeneracies.We therefore fix these parameters to the values (listed in Table 1) derived using the ℓ max = 2 eclipse map model and re-fit the map with a variety of methods.Bell et al. (2024) found that an n = 1 Fourier series fitted the out-of-eclipse phase curve very poorly, so we do not consider any ℓ max = 1 spherical harmonic fits.The top row of Figure 5 shows the eclipse map re-fitted using spherical harmonics up to order ℓ max = 3.We fix the orbital, planetary, stellar, and systematic parameters to the median values derived from the fit in Figure 3.We found that the ℓ max = 3 maps required too many pixels to sample effectively with PyMC3, so in this case we sampled the spherical harmonic coefficients directly with mc3 (Cubillos et al. 2017).This method allows us to impose positivity by excluding negative maps from the fitting process, instead of sampling pixels with positive priors.We set the spherical harmonic coefficients to have Gaussian priors with a mean of 0 and a standard deviation equal to the mean of the observed planetary flux.
The resulting map in the top row of Figure 5 shows a similar median structure to the map in Figure 3, with more uncertainty due to the increased degrees of freedom.There is a longitudinal offset of (0.50 +14.79  −8.04 ) • , which is consistent with the ℓ max = 2 fit value of (7.75 +0.36  −0.36 ) • .Its median value is almost zero, which may seem unusual as the overall phase curve has an offset of (7.3 +0.4  −0.4 ) • .This apparent discrepancy is due to the more complex shapes allowed by the ℓ max = 3 spherical harmonics, which produce a dayside hot-spot which peaks at the substellar point but extends further east than west (as shown in Figure 5) resulting in an overall shift in the phase curve.We discuss this distinction between hot-spot position and phase curve offset in Section 4, showing how the phase curve offsets and hot-spot positions can differ in numerical simulations.
The latitudinal structure of the ℓ max = 3 map is somewhat different to the ℓ max = 2 map, being flatter near the equator with no significant latitudinal hot-spot shift.We will also discuss in Section 4 how latitudinal structure could be degraded by using low-order basis maps to fit the data.For all these reasons, the spatial resolution of the map used to fit the data is a key question.In Section 3.6, we will therefore discuss how many parameters are justified to fit the dataset.

Eclipse map fitted with eigenmapping
The bottom row of Figure 5 shows a map fitted using "eigenmapping" (Rauscher et al. 2018) with ThERESA (Challener & Rauscher 2022), with the orbital and sys-tematic parameters again fixed to the median values derived from the ℓ max = 2 fit in Figures 1 and 3.
Briefly, eigenmapping starts with spherical harmonic phase curves up to some degree ℓ max , orthogonalizes them with principal component analysis to create "eigencurves", and fits the observed phase curve as a sum of N E highest-variance eigencurves.ℓ max and N E are selected to minimize the BIC.Each eigencurve has a corresponding eigenmap, and the fitted map is the corresponding sum of these eigenmaps.We again enforce a positive-flux constraint on the fitted map globally.Eigenmapping has been used to map observations of HD 189733b with the Spitzer Space Telescope (Rauscher et al. 2018;Challener & Rauscher 2022), and observations of WASP-18b with JWST NIRISS/SOSS (Coulombe et al. 2023).
We find the best BIC for ℓ max = 2 and N E = 6, a larger number of eigencurves than used for WASP-18b in Coulombe et al. (2023), showing the improved precision of this WASP-43b dataset.We find a slight eastward hot-spot offset of (7.5 +0.5 −0.5 ) • .This eigenmap achieves the best BIC score in Table 2; its fitted map is almost the same as the ℓ max = 2 map in Figure 3, suggesting that it achieves this improved BIC score by discarding the map components that contribute the least to the observed phase curve.This means that it produces the same fit quality using fewer parameters (see Section 4 for a related discussion of the mapping "null space").

Eclipse maps fitted with single eclipses
The maps in Figures 3 and 5 are fitted to the entire dataset, including 2D information from both eclipses as well as 1D information from the rest of the phase curve.
To isolate the 2D mapping information in each eclipse, Figure 6 shows the result of re-fitting the ℓ max = 2 eclipse map using the first and second eclipses only.As in Figure 5, we fix the orbital and systematic parameters to the median values derived with the model in Figures 3, listed in Table 1.The exceptions are the magnitude of the exponential ramp and the magnitude of the linear trend, which we re-fit as these require different values given the different starting points in time The first row in Figure 6 shows an ℓ max = 2 eclipse map fitted using a ∼ 2.2 hour section of data containing the first eclipse only with a small section of phase curve either side.We only show the map on the dayside of the planet, as there is no information about the nightside contained in this limited dataset.The map is similar to Figure 3  Fitted max = 2 map, both eclipses only, with fixed orbit and systematics Figure 6.Eclipse maps fitted using each eclipse separately.First row: an eclipse map fitted using spherical harmonics up to ℓmax = 2, using a ∼ 2.2 hour section of the data centered on the first eclipse.The orbital and systematic parameters are fixed to those derived using the map fitted in Figure 3 (except the time-dependent systematics, which are re-fitted), listed in Table 1.Only the dayside of the map is shown as there is no information about the nightside in this limited dataset.The map has a similar longitudinal structure to Figure 3, but finds more latitudinal asymmetry.Second row: an eclipse map fitted using a ∼ 2.2 hour section of the data centered on the second eclipse.This finds a similar map to that in Figure 3.These two maps show that both eclipses are mostly consistent, although the first eclipse implies more latitudinal asymmetry.Third row: an eclipse map fitted using both the first and second eclipses as defined above, excluding the rest of the phase curve.This combined fit is similar to the individual fits, especially in its longitudinal structure.
be insufficient to derive an accurate eclipse map.There is a greater range of fitted maps at −90 • latitude than +90 • latitude because the inclination of the orbit angles this pole away from the observer.
The second row in Figure 6 shows an eclipse map fitted to the second eclipse only, with a smaller latitudinal hot-spot shift than the fit to the first eclipse only.The third row in Figure 6 shows an eclipse map fitted to both eclipses together.We suggest that the original fit to the entire dataset in Figure 3 which found a latitudinal shift of (−10.72 +4.14  −4.68 ) • is better constrained in latitude than this two-eclipse fit (despite containing no additional latitudinal information) because the additional measurement of the rest of the phase curve provides independent longitudinal information.This breaks degeneracies between the longitudinal and latitudinal information in the eclipse mapping signal (Boone et al. 2023).
The maps in Figure 6 show that the eclipse mapping signals in both eclipses are mostly consistent with each other, as would be expected from the similar residual signals in Figure 2.However, they have different degrees of latitudinal asymmetry, which also manifests in the overall fit in Figure 3. Figure 6 shows that the first eclipse implies more latitudinal asymmetry.This may be a real effect, but it may also be due to the increased effect of instrumental systematics on this eclipse.The blue line in Figure 1 shows that the systematic signal is still relatively large during the first eclipse, but that it is almost negligible during the second eclipse.This shows the utility of observing a full phase curve where the periodicity of the astrophysical signal isolates it from the systematic signal.

Model Selection
To compare our fitted eclipse maps, we calculate the Bayesian Information Criterion (BIC) (Schwarz 1978) and the Akaike Information Criterion (AIC) (Akaike 1981).The BIC is: where k is the number of model parameters, N is the number of data points, and L is the model likelihood: for a sum over all N data points D i with uncertainty σ i , fitted by a model M i .This uncertainty is derived from the residuals of the originally fitted phase curve, by scaling the expected error by a factor referred to as the "Uncertainty Scaling Factor" in Figure 10 (see Bell et al. (2024)).
A smaller BIC implies a better model, with a better fit or fewer fitting parameters.The difference in BIC is the relevant quantity for model comparison, so Table 2 shows the ∆BIC for each model compared to Model M5, the eigenmap in Figure 5.This criterion allows comparison of nested models where it is assumed that the true model is inside the set of tested models (Burnham & Anderson 2004).
We also consider the AIC of each model, which is: Due to the large number of points N in this dataset, the BIC applies a stronger penalty for each parameter than the AIC.We also present the "weights" of each model, w BIC and w AIC, which are the relative probabilities of each compared model in a set, based on their ∆BIC or ∆AIC.The w BIC for a particular model is: where ∆BIC is compared to the best-performing model, and there are R total models being compared with a score of ∆BIC r each.w AIC depends on ∆AIC in the same way.Table 2 summarises our compared models, showing their χ 2 values, their ∆BIC and ∆AIC relative to the model with the best BIC and AIC, and the implied probability of each model.
Model M1 is the Fourier series model used in Bell et al. (2024).It achieves a significantly worse χ 2 value than all the eclipse mapping models, due to a worse fit to the eclipse shape.Model M2 is the eclipse map model shown in Figure 4, where the latitudinal structure is fixed to be flat.This achieves a better χ 2 value than the Fourier series model, as it includes the effect of longitudinal map structure on the eclipse shape.However, its χ 2 value is worse than the fully two-dimensional eclipse map fits, as shown by its poor fit in Figure 2. Interestingly, due to its smaller number of parameters, it achieves the secondbest BIC score of all the eclipse mapping models due to the strong penalty applied to the number of parameters by the BIC.However, it achieves the worst AIC score of the eclipse map models, where the number of parameters is penalised less strongly.We can confidently reject models M1 and M2 compared to the 2D eclipse mapping models.
Model M3 is the ℓ max = 2 eclipse mapping model in Figure 3.It achieves a better χ 2 value than models M1 and M2, as shown by its good fit in Figure 2. It has a worse BIC score than the eigenmap (model M5), as they have similar χ 2 values, but model M3 uses two more parameters.It has a more similar AIC score to the eigen-map, where these additional parameters are penalized less heavily.
Model M4 is the ℓ max = 3 eclipse mapping model in Figure 5.It achieves the best χ 2 value due to its increased degrees of mapping freedom, but has poor BIC and AIC scores due to its increased number of parameters.These metrics imply that the data quality does not justify this number of degrees of freedom, which is also implied by the large uncertainty on the ℓ max = 3 map in Figure 5.
Model M5 is the ℓ max = 2 eigenmap using 6 basis maps.It achieves a similar χ 2 to the other eclipse mapping models, and has the best BIC and AIC scores due to its reduced number of parameters.It achieves these by discarding mapping structures that contribute weakly to the observed phase curve, so is able to match the observations with fewer degrees of freedom.
It is not clear whether the BIC or AIC is the better metric for model comparison, or if another metric would be more appropriate.The BIC assumes that the set of fitted models includes the "true" physical system, while the AIC assumes that all of the tested models are inexact representations of this system (Burnham & Anderson 2004).It could be argued that eclipse maps fitted to real planets do include the "true" map, as spherical harmonics form a complete basis set on the sphere.On the other hand, information about the real map is inevitably lost in the mapping process due to the need to truncate fits to low-order harmonics, and due to the presence of a "null space" (see Section 4 for discussion of both issues).
To summarise our model comparison, 1.The Fourier series model (M1), and the eclipse map with flat latitudinal structure (M2) perform badly on all statistical metrics, so we confidently reject them.
2. The latitude-longitude ℓ max = 2 eclipse map (M3) is favoured over the longitude-only eclipse map (M2), showing the presence of latitudinal information in the data.
3. The ℓ max = 3 eclipse map achieves a better χ 2 value than model M3, but at the cost of many more parameters, so has worse BIC and AIC scores.
4. The ℓ max = 2 eigenmap achieves the best BIC and AIC scores, producing a fit of comparable quality to model M3 with two fewer parameters Therefore, the eigenmap (model M5) is the best performing model on these metrics.However, it requires model M3 to fit the orbital, stellar, and systematic parameters simultaneously first, and then both models produce a similar result in the end.The eigenmap achieves a better AIC and BIC by discarding structures that do not contribute strongly to the map, producing a very similar map to model M3 overall.We suggest that more work can be done on the process of fitting eclipse maps and comparing levels of model complexity.Cross-validation may provide a more practical metric for eclipse mapping model comparison, such as the leave-one-out cross-validation applied by Challener et al. (2023), which is asymptotically equivalent to the AIC but provides advantages such as an absolute measurement of model predictive power (Stone 1977).

COMPARISON TO GENERAL CIRCULATION MODELS
In this section, we interpret the eclipse maps by comparing them to three-dimensional General Circulation Models (GCMs).We use some of the GCMs presented in Bell et al. (2024), selecting those that matched the out-of-eclipse phase curve well.We include a simulation from four models: THOR (Mendonça et al. 2016(Mendonça et al. , 2018a,b,b The THOR and RM-GCM simulations use semi-grey radiative transfer, while the PCM and expeRT/MITgcm simulations use multi-band correlated-k schemes.The THOR, RM-GCM, and PCM simulations feature clouds (on the nightside only for THOR), whereas the ex-peRT/MITgcm simulation is free of clouds.Output from each model was post-processed using several multiband radiative transfer codes to calculate spectrally resolved emission from each column, which was then integrated from 5 to 10.5 µm, weighted by the MIRI/LRS throughput.Appendix A and Bell et al. (2024) give more detail on each model.

Observable Features of GCMs
Figure 7 shows thermal emission maps post-processed from each GCM simulation.The top row shows the thermal emission integrated from 5 to 10.5 µm.These maps contain small-scale structure that is not present in the observed eclipse maps in Section 3. A chevronlike structure is present in all the GCM maps, which is commonly attributed to the temperature structure associated with equatorial Kelvin and Rossby waves (Matsuno 1966 Not all spatial brightness patterns produce a signal in phase curve space (Luger et al. 2021).For instance, rotational phase curves of spatially unresolved objects (i.e., exoplanet phase curves) are insensitive to latitudinal structures.If the object is eclipsed, as is the case for WASP-43b, then the shape of eclipse ingress and egress breaks many of these degeneracies, but still leaves a "null space" of unobservable patterns (Challener & Rauscher 2023).This null space means that retrieved eclipse maps will not always match thermal emission output from GCM simulations.This effect may account for some of the discrepancy between the GCM maps shown in the top row of Figure 7 and the observed eclipse maps in Section 3.
To illustrate this effect, we separate our GCM maps into their observable and null components following Challener & Rauscher (2023).The second row of Figure 7 shows the observable component for each GCM.The most significant difference between the observable GCM maps and the original GCM maps is the lack of spatial variation on the nightside of the observable maps.In this region, we only have information from the phase curve, so our mapping capabilities are limited to largescale longitudinal variation.This means that the fine structure associated with, for example, nightside equatorial Rossby waves, cannot be constrained by observations.By contrast, the dayside is scanned by the eclipse, which means that some features remain observable, such as the broad shape and location of the hot-spot.While the observable maps appear to have some new structures compared to the original GCM maps, this is simply the result of removing the high-order structures of the null space.
There are a number of differences between the observable parts of the GCMs, and the observed eclipse maps.Each observable GCM map has a large hot-spot that is offset eastwards from the substellar point, while the eclipse maps 3 have smaller or negligible hot-spot sifts.In addition, the dayside emission in the GCMs varies less with latitude near the equator than the emission in any of the ℓ max = 2 eclipse maps (Figures 3 and 6), or the eclipse map derived using the eigenmapping method (bottom row of Figure 5).
None of the GCM maps in the top row of Figure 7 show a significant latitudinal hot-spot offset, which is present in the observed ℓ max = 2 full phase curve map, the ℓ max = 2 map derived from the first eclipse only, and the eigenmap.It is unsurprising that the GCMs do not show a latitudinal offset, given that the forcing and boundary conditions for each model are hemispherically symmetric.However, the observable maps for each do have small latitudinal offsets due to asymmetries introduced by the inclined viewing angle.This implies that asymmetries could be introduced by mapping at high orders.
For the purposes of comparison to the ℓ max = 2 and ℓ max = 3 maps we fit to the observations, it is important to note that ℓ max = 2 and ℓ max = 3 maps have no null space for this observation.The non-zero null spaces identified in Figure 7 are due to the high-order spherical harmonic bases used to represent the small-scale structure of the GCM results.Therefore, while the presence of a null space places a theoretical upper limit on the accuracy of eclipse mapping it should not have a direct effect on our fitted ℓ max = 2 and ℓ max = 3 maps, and should not introduce a latitudinal asymmetry at these low orders.The eigenmap basis may have some latitudinal asymmetry introduced by its calculation of the structures which contribute most strongly to the light curve.

THOR
. A comparison of four GCM simulations of WASP-43b and how their spatial distributions could appear in an eclipse map with restricted spherical harmonic order.First row: the modelled outgoing longwave radiation (OLR) in the 5 to 10.5 µm MIRI/LRS bandpass for each of the four GCM simulations, expressed as a ratio of planetary to stellar flux.Second row: observable emission from the GCM simulations, with the null space removed as described by Luger et al. (2021) and Challener & Rauscher (2023).Third row: the OLR represented using spherical harmonics up to ℓmax = 3, showing how structure is lost; this is the best map that could be achieved by a fit like the ℓmax = 3 fit in Figure 5. Fourth row: the OLR represented using a basis of the eigenmaps used in Figure 5. Fifth row: the OLR represented using spherical harmonics up to ℓmax = 2; this is the best map that could be achieved by a fit like that in Figure 3.All of the simulations match the dayside and nightside amplitudes of the observed phase curve and eclipse map fairly well, as well as the small eastward phase curve shift.The succession of plots here shows how using low-order spherical harmonics limits the structures that can be fitted with an eclipse map.

Effect of Mapping Order
The eclipse maps in Section 3 are fitted with low-order spherical harmonics up to ℓ max = 2 and ℓ max = 3, as including higher-order harmonics could lead to overfitting.This means that some realistic structures cannot be fitted accurately -for example, a sharp brightness temperature gradient caused by the onset of cloud formation needs high-order harmonics to fit it accurately (Parmentier et al. 2016).To illustrate this effect, the third and fourth rows of Figure 7 show emission maps from each GCM, truncated to use information from harmonics up to ℓ max = 2 and ℓ max = 3 only (computed using starry).
Latitudinal and longitudinal cross-sections of the truncated maps (taken at the sub-stellar longitude, and along the equator, respectively) are shown in Figure 8, where they are compared to the relevant observed eclipse maps.The eigenmaps fitted in Figure 5 are also compared to each GCM expressed using a basis of the eigenmaps from this particular fit.
Truncation of the GCM emission to ℓ max = 2 forces the dayside emission to be more strongly peaked on the equator than in the higher-order representations.This difference is consistent with the different latitudinal structure suggested by the ℓ max = 2 and ℓ max = 3 eclipse maps from Section 3 (comparing, e.g., Figures 3 and the top row of Figure 5).The ℓ max = 2 fit is therefore consistent with the real WASP-43b either having i) a peaked latitudinal structure in reality, or ii) a flatter latitudinal structure in reality (like the GCMs), which is being masked by the spherical harmonic truncation.While the flatter fits in the ℓ max = 3 map in Figure 5 are in better agreement with the GCM simulations, the ℓ = 3 modes are not constrained well enough to conclude that the flatter structure is a better representation of the "real" emission from WASP-43b.
Turning to longitudinal structure, Figure 7 shows that restricting the emission to ℓ max = 2 broadens the longitudinal structure of the hot-spot, most notably for the THOR and expeRT/MITgcm simulations.It also removes sharp gradients in emission at the terminators of the THOR and RM-GCM simulations, which are associated with the formation of clouds on the nightside.As with the latitudinal structure, these differences between ℓ max = 2 and ℓ max = 3 are consistent with those present in the eclipse maps presented in Section 3 and Figure 8.

Longitudinal Offsets
The eastwards hot-spot shift of the temperature structure of tidally locked planets is one of their key observable features, traditionally derived from the phase offset of the maximum of their phase curve.This phase curve offset is not, however, identical to the actual shift of the 2D emission structure.Instead, it represents the integration of the emission over the observed hemisphere, weighted by viewing angle.Figure 9 shows how the simulated phase curves for each GCM have different phase curve offsets than the offsets of their emission map, which we define as the maximum longitude θ max of +π/2 −π/2 F p (θ, ϕ) cos ϕ dϕ.We refer to this averaging as "meridional averaging" from now on, and always consider this to be the "longitudinal offset" of a 2D map, for comparison to the "peak offset" of a phase curve.
By this measure, the THOR GCM has a phase curve offset of 11.1 • but an emission map offset of 17.5 • (meridionally averaged).Figure 9 shows how this difference is especially pronounced for this simulation, due to its sharp gradients shown in the top row of Figure 7.These gradients are poorly captured by the ℓ max = 2 representation, so the phase curve has a very different offset, as it is dominated by low-order modes (Cowan & Agol 2008).This also results in different emission map offsets in its ℓ max = 2 and ℓ max = 3 representations, which are 13.2 • and 26.5 • respectively.The ℓ max = 2 emission map offset is similar to the phase curve offset, as expected.This shows how low-order representations of emission maps can bias measurements of the largescale offsets, as well as discard small-scale information.
Figure 9 shows that the meridionally averaged emission map offset for the ℓ max = 2 map in Figure 3 is (7.75 +0.36  −0.36 ) • .This is consistent with the emission map offset for the first, second, and combined eclipses in Figure 6, which are (7.03+7.57−5.41 ) • , (6.67 +5.77 −4.32 ) • , and (4.86 +3.6  −2.88 ) • respectively.The fit to the full phase curve is more precise due to the increased longitudinal information present in the out-of-eclipse phase curve, which constrains the low-order longitudinal structure very well and break degeneracies with the latitudinal structure (Boone et al. 2023).The ℓ max = 3 map fit to the whole dataset shown in Figure 5 has a map offset of (0.50 +14.79  −8.04 ) • that is consistent with the ℓ max = 3 eclipseonly fits, as well as with the ℓ max = 2 map fits.As discussed in Section 3.3, this map has a hot-spot offset with a median position on the substellar point, but a non-zero phase offset of around 7 degrees, consistent with the overall phase curve, showing the difference between these two metrics.The ℓ max = 3 fit also has much more uncertainty due to its increased degrees of mapping freedom, resulting in the large range of fitted maps in Figure 5.
The longitudinal offsets of the n = 2 Fourier series fit, and the ℓ max = 2 map fit to the full dataset, appear very precise compared to the ℓ max = 3 map fit to the full dataset.This reflects the choice of fitting   2024)), compared to the phase curve offset simulated from each GCM.Right: The (meridionally longitudinal offset in each 2D map.The ℓmax = 2 and ℓmax = 3 offsets (using all the dataset) are from the map fits in Figures 3 and 5.The ℓmax = 2 map fits using the eclipses only correspond to the fits in Figure 6 (the black point shows the fit to both eclipses, and the grey points show the fits to the first and second eclipses).The eclipse-only eigenmap fits and the eclipse-only ℓmax = 3 map fits are not plotted in full elsewhere.The GCM points show the cos ϕ-weighted offset for the GCMs truncated to ℓmax = 2 and ℓmax = 3 representations shown in Figure 7, as well as for the maps with no truncation.
functions rather than a truly increased precision.The out-of-eclipse phase curve in Figure 1 is very well constrained, so the n = 2 Fourier series model has a very precise peak as shown in Figure 9. Similarly, most of the ℓ max = 2 modes in the map in Figure 3 are tightly constrained by the out-of-eclipse phase curve, resulting in an apparently more precise value in Figure 9.This precision is due to the limited fitting functions rather than actual statistical certainty.This is shown by the greater uncertainty on the ℓ max = 3 map in Figure 5, which has more degrees of mapping freedom.It can be consistent with the out-of-eclipse phase curve with many different maps, resulting in the larger uncertainty on its meridionally averaged longitudinal shift in Figure 9.In summary, the derived emission map offset varies between around 0 and 20 degrees based on the mapping method used.The map offsets in the GCMs used for comparison are larger, varying between around 10 and 30 degrees depending on the mapping basis.We conclude that while the ℓ max = 2 map and eigenmap are the best constrained maps, they may be discarding important structures from the real map.The ℓ max = 3 map has access to more structures, but has too large an uncertainty to provide precise conclusions.Observations of more eclipses would provide the information needed to better constrain the ℓ = 3 modes.

CONCLUSIONS
In this study, we have presented eclipse maps of the thermal emission of WASP-43b derived from a JWST MIRI/LRS phase curve.We fitted a map using ℓ max = 2 spherical harmonics simultaneously with the parameters of the system and the systematic parameters of the in-strument.There is a clear residual mapping signal in the ingress and egress of the eclipses, which a 2D map model fits much better than a Fourier series model.This residual signal is at most 450 ppm, composed of a signal ∼250 ppm due to longitudinal structure and ∼200 ppm due to latitudinal structure.This is the first time that the magnitude of the signal of the latitudinal structure has been statistically identified in an eclipse map of an exoplanet.
Fitting this dataset with an eclipse map model derived statistically significantly different parameters to those derived with a Fourier series model as in Bell et al. (2024).We suggest these updated parameters are more accurate due to a more accurate fit to the eclipse shapes, but note that they have little effect on the eclipse map itself, as shown in Figure 11.There were no significant degeneracies between the fitted parameters.
Figure 3 shows that the ℓ max = 2 map found a small eastward hot-spot shift of (7.75 +0.36  −0.36 ) • (defined using a meridional average), compared to an eastward phase shift of (7.3 +0.4  −0.4 ) • derived from fitting a Fourier series to the phase curve.This ℓ max = 2 fit finds a sharply peaked latitudinal structure, with a small latitudinal offset of (−10.72 +4.14  −4.68 ) • .Figure 5 shows an ℓ max = 3 eclipse map fitted with the non-mapping parameters fixed to those derived from the ℓ max = 2 map fit.This found a smaller but more uncertain hot-spot shift of (0.5 +14.79 −8.04 ) • degrees east, and a flatter latitudinal structure near the equator than the ℓ max = 2 map.It achieved a better χ 2 than the ℓ max = 2 model, but a worse BIC and AIC due to its higher number of parameters.
In addition, we fitted an "eigenmap" (with the orbital and systematic parameters fixed), which fitted the data with six eigenmaps with orthogonal phase curves from an ℓ max = 2 basis.Figure 5 shows that this derived a map very similar to the original ℓ max = 2 eclipse map in Figure 3, but with fewer parameters as it discarded the mapping structures that contribute least to the phase curve.Table 2 shows how this model achieved the best BIC score, with a similar fit quality but fewer parameters.
As described in Section 3.1, a simple estimate of brightness temperature from the observed broadband flux corresponds to a temperature at the substellar point of (1790.0+23.0 −29.0 ) K for the ℓ max = 2 map.This is consistent with the brightness temperature of the eigenmap at the substellar point of (1783.0+51.0 −70.0 ) K, as well as with the brightness temperature of the ℓ max = 3 map which is (1822.0+75.0 −92.0 ) K at the substellar point.We tested the mapping information in each eclipse separately by re-fitting a map to each eclipse individually.Figure 6 shows that the second eclipse produced a dayside map similar to that derived from the full phase curve, but the first eclipse produced one with a larger latitudinal offset.This may be due to increased systematic effects on the first eclipse; it could also be a real difference in thermal emission structure but we suggest that such large-scale variability over one orbit is unlikely.We recommend that those looking to use MIRI/LRS for eclipse mapping observations allow for significant settling time at the start of their observations to reduce the impact of the large instrumental systematics.
Figures 7 and 8 compare the fitted eclipse maps to four GCM simulations from Bell et al. (2024).In general, the GCMs match the eclipse maps well given the effect of truncating the spherical harmonic order of the maps.The ℓ max = 2 map is largely consistent with the ℓ max = 2 GCM representations.The ℓ max = 3 map has almost no hot-spot shift, different to the GCMs.It has a flatter latitudinal structure near the equator, which appears more consistent with the structure of the GCMs, but the higher uncertainty on the ℓ max = 3 map makes it difficult to compare exactly.
We conclude that there is a strong eclipse mapping signal in this MIRI/LRS observation of WASP-43b, which can strongly constrain the ℓ max = 2 spherical harmonic components of the planetary emission.Our fiducial map is the ℓ max = 2 eclipse map in Figure 3, fitted simultaneously with the orbital, stellar, and systematic parameters.This fiducial map finds a (meridionally averaged) hot-spot shift of (7.75 +0.36  −0.36 ) • eastward (shown in Figure 9), compared to a shift in the peak of the phase curve of (7.3 +0.4  −0.4 ) • eastward.
We also highlight our statistically preferred map (according to the BIC and AIC), which is the eigenmap fitted using ℓ max = 3 and N E = 6 basis maps; this produces a very similar map to the fiducial ℓ max = 2 map using fewer fitting parameters.However, the limited spatial order of the fitting functions means we may be missing important structure as shown in Figure 7.The limitation of our maps in general to ℓ max = 2 structures is a very model-dependent result; there may be information about smaller-scale structures in the data that the spherical harmonics are not well suited to capture.Future studies could investigate mapping methods with more degrees of freedom to relax the dependence of our conclusions on the form of the low-order spherical harmonics (Horne 1985).Measuring further eclipses would allow better constraints on the smaller-scale structures.In general, we suggest that observing multiple eclipses, and ideally observing a full phase curve, is necessary to obtain reliable eclipse maps for even the best targets like WASP-43b.
In summary, the dayside eclipse maps fitted to this dataset are generally symmetric about the equator, have a small hot-spot shift eastward from the substellar point, and vary smoothly away from the substellar point.This structure implies a weak eastward heat transport by atmospheric dynamics (Hammond & Lewis 2021), no strong latitudinal asymmetries driven by magnetic fields or atmospheric dynamics (Rogers & Komacek 2014;Skinner & Cho 2022), and no dayside homogeneity driven by magnetic fields (Beltz et al. 2021).They are generally consistent with our four numerical simulations of WASP-43b, which have small eastwards hemispheric hot-spot shifts, and are symmetric about the equator.However, the eastward hot-spot shift measured by the eclipse maps is generally smaller than that predicted by the simulations.This dataset contains spectroscopic information that we averaged out in this broadband analysis; future studies could derive maps in particular wavelength bands or regions of molecular absorption, in order to measure three-dimensional temperature structure or the distribution of particular chemical species.However, the accompanying decrease in precision may weaken the mapping signal and the accuracy of the resulting maps.More precise mapping may also reveal subtler spatial features that are not present in this analysis, and mapping of planets with different properties may reveal different spatial structures on their daysides.

GCMS AND POST-PROCESSING
We compare four GCM simulations from Bell et al. (2024) to our derived eclipse maps.The first is THOR (Mendonça et al. 2016) (Simulation 31 in Bell et al. 2024), with the same model configuration as was used to simulate WASP-43b previously (Mendonça et al. 2018a,b).The simulation presented here uses semi-grey radiative transfer, and a simple parameterized cloud scheme on the nightside of the planet (Mendonça et al. 2018a).It was run for roughly 9400 orbits, and the output data was averaged over the last 500 days.The emission from each column was calculated by post-processing the results with a multiwavelength radiative transfer model (Mendonça et al. 2015), assuming 1× solar metallicity and equilibrium chemical species concentration calculated with the FastChem model (Stock et al. 2018).
The second is expeRT/MITgcm (Simulation 20 in Bell et al. 2024).This model uses the dynamical core of the MITgcm (Adcroft et al. 2004) on a C32 cubed-sphere grid, coupled to a non-grey radiative transfer scheme based on petitRADTRANS (Mollière et al. 2019).It follows a setup used in Carone et al. (2020) and Schneider et al. (2022) to investigate the deep dynamics of hot Jupiters, with radiative transfer as described in Schneider et al. (2022).The only difference between the configuration used here and that in Schneider et al. (2022) is the omission of TiO and VO.Eleven frequency bins are used, with 16 k-coefficients for each.The simulation shown here has 47 vertical levels and 10× solar metallicity, and was run for 1500 days, with the results shown here averaged over the last 100 days.The spectrally resolved emission was post-processed using petitRADTRANS (Mollière et al. 2019) and prt phasecurve (Schneider et al. 2022) using a spectral resolution of R = 100.
The third is the Generic Planetary Climate Model (Generic PCM) (Simulation 7 in Bell et al. 2024), which has been used to model exoplanets (Charnay et al. 2015;Turbet et al. 2016;Teinturier et al. 2023) and the planets of the Solar System (Turbet et al. 2021;Spiga et al. 2020).The simulation presented here uses a horizontal resolution of 64×48 with 40 vertical levels.The Generic PCM treats clouds as radiatively active tracers of fixed radii.The correlated-k  1 lists the numerical values.Some of the parameters have statistically significant differences, although Figure 11 shows that these differences do not produce large changes in the eclipse map.We expect that the blue posteriors derived by the ℓmax = 2 eclipse map model are more accurate, due to its more accurate model of the eclipse shape.1.The map is almost identical to that in Figure 3, showing that the small (but statistically significant) differences between the two fits in Figure 10 do not meaningfully affect the resulting ℓmax = 2 eclipse map.
radiative transfer uses 27 shortwave and 26 longwave bins, each with 16 k-coefficients.The model is run for 2000 orbits, and the data presented is an average of the model output over the last 100 days.The simulation presented here has 1× solar metallicity and models Mg 2 SiO 4 clouds with radius 1µm.The simulation was post-processed with the Pytmosph3R code (Falco et al. 2022) to calculate the spectrally resolved emission from each column, which was then integrated from 5 to 10.5 µm.
The fourth is the RM-GCM (Simulation 26 in Bell et al. 2024), which was adapted from the GCM of Hoskins & Simmons (1975) by Menou & Rauscher (2009), Rauscher &Menou (2010), andRauscher &Menou (2012), and has been applied to investigations of exoplanets using semi-grey radiative transfer and aerosol scattering (Roman & Rauscher 2017;Roman et al. 2021).The simulation presented here has 1x solar metallicity and condensate clouds represented by aerosols as described in Bell et al. (2024), with 50 vertical levels, and was run for over 3500 orbits.It was post-processed following Zhang et al. (2017) and Malsky et al. (2021) to produce the spectrally resolved emission.

B. ORBITAL AND SYSTEMATIC PARAMETERS
Figure 10 shows the orbital, planetary, stellar, and systematic parameters fitted with the n = 2 Fourier series model described in Section 3, compared to the parameters fitted with the ℓ max = 2 eclipse map model shown in Figure 3.We find no significant degeneracies between these parameters, or between these parameters and the Fourier series coefficient or map pixels used to fit the light curve shape as described in Section 2 (not shown).There are some degeneracies between the parameters of the systematic model, due to inherent degeneracies between a linear slope and a long-timescale exponential ramp.
As described in Section 3, the eclipse map model produces a better fit to the data (the smaller χ 2 value listed in Table 2).This is as expected, as it should provide a more realistic model of the eclipse shape.We therefore assume that the parameters derived by the eclipse map model (coloured blue in Figure 10) are more accurate.The eclipse timing is very well constrained despite the possibility for degeneracy between eclipse timing and eclipse mapping (Williams et al. 2006).We suggest this is due to the high precision and cadence of the observations, the independent measurement of longitudinal structure from the phase curve, and the presence of a transit in the dataset.
The Fourier series model finds a statistically significantly smaller planetary radius, higher orbital inclination, and larger semi-major axis, than the eclipse map model.The Fourier series model also finds different stellar limb darkening parameters to the eclipse mapping model, but the resulting limb darkening profile posteriors (not shown) are almost the same.This reflects degeneracies in the limb darkening model, and the limited limb darkening information available at these longer wavelengths (Morello et al. 2017;Morello 2018).The two fits in Figure 10 also find different systematic parameters, but these result in almost exactly the same systematic model.The Fourier series model fits a weak exponential ramp with a long timescale, plus a strong decreasing linear trend.The eclipse map model fits a strong exponential ramp with a short timescale, plus a weak decreasing linear trend.These two systematic models are almost identical within observational uncertainty when the linear and exponential trends are combined.
These statistically significant differences in the parameters fitted by the two models do not produce significant differences in the resulting eclipse map. Figure 11 shows an ℓ max = 2 eclipse map fitted using orbital, planetary, stellar, and systematic parameters fixed to those derived using the n = 2 Fourier series model (listed in Table 1).The resulting map is very similar to that fitted simultaneously with the additional parameters in Figure 3, showing that the differences in parameters in Figure 10 do not have a significant effect on the eclipse map.We suggest that it is still best to fit the orbital and systematic parameters simultaneously with an eclipse map model in general, to derive more accurate parameters and to search for degeneracies with the fitted map.

Figure 1 .
Figure 1.The JWST MIRI/LRS dataset, showing the observed phase curve (black points), the fitted ℓmax = 2 eclipse map model (red line) shown in Figure 3, and the systematics model (blue line) fitted alongside the eclipse map.
Figure  2shows the residual difference in ingress and egress between the observations (shown as black, green, and and 5, but has a larger latitudinal hot-spot shift than the map in Figure3.The latitudinal peak of the fit to the first eclipse only is weakly constrained, showing how an observation of only a single eclipse is likely to Figure 8.A comparison of the GCM simulation results in Figure 7 to the eclipse maps fitted to the data in Figure 1.First row: the GCM maps in the first row of Figure 7. Second row: the ℓmax = 2 map in Figure 3 compared to the ℓmax = 2 GCM maps in the fourth row of Figure 7. Third row: the eigenmap in Figure 5 compared to the eigenmap representations of the GCMs.Fourth row: the ℓmax = 3 map in Figure 3 compared to the ℓmax = 3 GCM maps in the third row of Figure 7.

Figure 9 .
Figure9.Left: the phase curve offset derived from the n = 2 fit using a Fourier series model (giving the same result asBell et al. (2024)), compared to the phase curve offset simulated from each GCM.Right: The (meridionally longitudinal offset in each 2D map.The ℓmax = 2 and ℓmax = 3 offsets (using all the dataset) are from the map fits in Figures3 and 5.The ℓmax = 2 map fits using the eclipses only correspond to the fits in Figure6(the black point shows the fit to both eclipses, and the grey points show the fits to the first and second eclipses).The eclipse-only eigenmap fits and the eclipse-only ℓmax = 3 map fits are not plotted in full elsewhere.The GCM points show the cos ϕ-weighted offset for the GCMs truncated to ℓmax = 2 and ℓmax = 3 representations shown in Figure7, as well as for the maps with no truncation.

Figure 10 .
Figure10.A corner plot of the posterior distributions of the orbital, planetary, stellar, and systematic parameters fitted to the observations in Figure1, produced with Foreman-Mackey (2016).The red posteriors correspond to the n = 2 Fourier series model described in Section 3, and the blue posteriors correspond to the ℓmax = 2 eclipse map model in Figure3.Table1lists the numerical values.Some of the parameters have statistically significant differences, although Figure11shows that these differences do not produce large changes in the eclipse map.We expect that the blue posteriors derived by the ℓmax = 2 eclipse map model are more accurate, due to its more accurate model of the eclipse shape.
Figure11.An ℓmax = 2 eclipse map fitted using the orbital and systematic parameters derived by fitting the Fourier series model, shown in Figure10and listed in Table1.The map is almost identical to that in Figure3, showing that the small (but statistically significant) differences between the two fits in Figure10do not meaningfully affect the resulting ℓmax = 2 eclipse map.

Table 2 .
(Challener & Rauscher 2022)wis & Ham- A comparisonof the models used to fit the full dataset.The offset value for the Fourier Series model is the offset of the phase curve, and the offset value for the eclipse maps is the meridionally averaged longitudinal offset as defined in Section 4.2.The ∆BIC and ∆AIC scores are calculated relative to the eigenmap model.When calculating the χ 2 , BIC, and AIC values, we exclude the transit because ThERESA(Challener & Rauscher 2022)does not model this explicitly.We do not calculate a comparative BIC or AIC between the Fourier series model and the eclipse map models as they are not nested.
This work is based on observations made with the NASA/ESA/CSA JWST.The data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127.The specific observations analyzed are associated with program JWST-ERS-01366 and can be accessed via DOI: 10.17909/kj9a-8d81 Support for this program was provided by NASA through a grant from the Space Telescope Science Institute.The results reported herein benefited during the design phase from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate.M.H. acknowledges support from Christ Church, University of Oxford.L.T. acknowledges access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.T.J.B. acknowledges funding support from the NASA Next Generation Space Telescope Flight Investigations program (now JWST) via WBS 411672.07.05.05.03.02. APPENDIXA.