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

Mark Hammond, Taylor J. Bell, Ryan C. Challener, Neil T. Lewis, Megan Weiner Mansfield, Isaac Malsky, Emily Rauscher, Jacob L. Bean, Ludmila Carone, João M. Mendonça, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mark Hammond, Taylor J. Bell, Ryan C. Challener, Neil T. Lewis, Megan Weiner Mansfield, et al. Two-dimensional Eclipse Mapping of the Hot-Jupiter WASP-43b with JWST MIRI/LRS. The Astronomical Journal, 2024, 168, 10.3847/1538-3881/ad434d. insu-04729881

HAL Id: insu-04729881

https://insu.hal.science/insu-04729881v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

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

Mark Hammond¹, Taylor J. Bell²,³, Ryan C. Challener⁴,⁵, Neil T. Lewis⁶, Megan Weiner Mansfield⁷,<fnref n="32" />, Isaac Malsky⁵, Emily Rauscher⁵, Jacob L. Bean⁸, Ludmila Carone⁹, João M. Mendonça¹⁰, Lucas Teinturier¹¹,¹², Xianyu Tan¹³,¹⁴, Nicolas Crouzet¹⁵, Laura Kreidberg¹⁶, Giuseppe Morello¹⁷,¹⁸,¹⁹, Vivien Parmentier²⁰, Jasmina Blecic²¹,²², Jean-Michel Désert²³, Christiane Helling²⁴,²⁵, Pierre-Olivier Lagage²⁶, Karan Molaverdikhani²⁷,²⁸, Matthew C. Nixon²⁹, Benjamin V. Rackham³⁰,³¹, and Jingxuan Yang¹

¹ Atmospheric, Oceanic, and Planetary Physics, Department of Physics, University of Oxford, Oxford, UK

² BAER Institute, NASA Ames Research Center, Moffett Field, CA, USA

³ Space Science and Astrobiology Division, NASA Ames Research Center, Moffett Field, CA, USA

⁴ Department of Astronomy, Cornell University, Ithaca, NY, USA

⁵ Department of Astronomy, University of Michigan, Ann Arbor, MI, USA

⁶ Department of Mathematics and Statistics, University of Exeter, Exeter, Devon, UK

⁷ Steward Observatory, University of Arizona, Tucson, AZ, USA

⁸ Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL, USA

⁹ Space Research Institute, Austrian Academy of Sciences, Graz, Austria

¹⁰ DTU Space, Technical University of Denmark, Kgs. Lyngby, Denmark

¹¹ LESIA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CNRS, Meudon, France

¹² Laboratoire de Météorologie Dynamique, IPSL, CNRS, Sorbonne Université, Ecole Normale Supérieure, Université PSL, Ecole Polytechnique, Institut Polytechnique de Paris, Paris, France

¹³ Tsung-Dao Lee Institute, Shanghai Jiao Tong University, 520 Shengrong Road, Shanghai, People's Republic of China

¹⁴ School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai, People's Republic of China

¹⁵ Leiden Observatory, Leiden University, Leiden, The Netherlands

¹⁶ Max Planck Institute for Astronomy, Heidelberg, Germany

¹⁷ Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden

¹⁸ Instituto de Astrofísica de Canarias (IAC), La Laguna, Tenerife, Spain

¹⁹ Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain

²⁰ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

²¹ Department of Physics, New York University Abu Dhabi, Abu Dhabi, UAE

²² Center for Astro, Particle and Planetary Physics (CAP3), New York University Abu Dhabi, Abu Dhabi, UAE

²³ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, The Netherlands

²⁴ Space Research Institute, Austrian Academy of Science, Schmiedlstrasse 6, 8042 Graz, Austria

²⁵ Institute for Theoretical Physics and Computational Physics, Graz University of Technology, Petersgasse 16/II, 8010 Graz, Austria

²⁶ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, Gif-sur-Yvette, France

²⁷ Universitäts-Sternwarte, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, D-81679 München, Germany

²⁸ Exzellenzcluster Origins, Boltzmannstraße 2, 85748 Garching, Germany

²⁹ Department of Astronomy, University of Maryland, College Park, MD, USA

³⁰ Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA

³¹ Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA, USA

Received 2024 January 24; revised 2024 April 22; accepted 2024 April 23; published 2024 June 4

Abstract

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 $l_{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 \pm 0.36)^\circ$, 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 $l_{max} = 2$ map. We also fit a map up to $l_{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 $l_{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.

<figure><img src="image_4.png" /></figure>

<footnote n="32">³² NHFP Sagan Fellow.</footnote>
*Unified Astronomy Thesaurus concepts: Exoplanet atmospheres (487); Exoplanet atmospheric structure (2310); Exoplanet atmospheric dynamics (2307); Extrasolar gaseous giant planets (509); Hot Jupiters (753)*

*Supporting material: figure set*

# 1. Introduction

Eclipse mapping measures the two-dimensional (2D) 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 2D 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 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 low-resolution 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 three-dimensional (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 hot spot, sharp gradients in brightness toward 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) analyzed 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.052}^{+0.050}$, and the ratio of its semimajor axis to stellar radius to be 4.859$_{-0.012}^{+0.013}$, 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. 2024). The longitudinal and latitudinal features of the dayside therefore affect the eclipse mapping signal almost equally, with a small bias toward the longitudinal features. This geometry makes WASP-43b especially well suited for eclipse mapping.

We fit an eclipse map to this data set, 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 $l_{max} = 2$ are constrained well, but that further observations are needed to precisely constrain higher-order mapping structures.

# 2. Data and Methods

A full orbit of WASP-43b with the JWST MIRI/LRS slitless 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 hr, 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 data set (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.

## 2.1. 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 semimajor axis to the stellar radius, and the two parameterized stellar limb-darkening parameters in Kipping (2013).
<figure><img src="image_5.png" /><figcaption>**Figure 1.** The JWST MIRI/LRS data set, showing the observed phase curve (black points), the fitted $\ell_{max} = 2$ eclipse map model (red line) shown in Figure 3, and the systematics model (blue line) fitted alongside the eclipse map.</figcaption></figure>

<table><caption>**Table 1**
The Orbital and Systematic Parameters Shown in Figure 10, for the Fourier Series Fit Described in Section 3 and the $\ell_{max} = 2$ Eclipse Map Fit Plotted in Figure 3</caption><thead><tr><th>Parameter</th><th>Fourier Series Fit</th><th>Eclipse Map Fit</th></tr></thead><tbody><tr><td>Planet-star radius ratio R<sub>p</sub> (R<sub>*</sub>)</td><td>0.15734<sup>+0.00017</sup><sub>-0.00015</sub></td><td>0.15839<sup>+0.00025</sup><sub>-0.00040</sub></td></tr><tr><td>Transit time (BMJD)</td><td>55934.292296<sup>+1.2×10<sup>-5</sup></sup><sub>-1.2×10<sup>-5</sup></sub></td><td>55934.292283<sup>+1.1×10<sup>-5</sup></sup><sub>-1.1×10<sup>-5</sup></sub></td></tr><tr><td>Inclination (deg)</td><td>82.277<sup>+0.050</sup><sub>-0.050</sub></td><td>82.106<sup>+0.050</sup><sub>-0.052</sub></td></tr><tr><td>Semimajor axis a (R<sub>*</sub>)</td><td>4.881<sup>+0.012</sup><sub>-0.012</sub></td><td>4.859<sup>+0.013</sup><sub>-0.012</sub></td></tr><tr><td>Limb-darkening parameter q<sub>1</sub></td><td>0.0565<sup>+0.0079</sup><sub>-0.0077</sub></td><td>0.0182<sup>+0.0081</sup><sub>-0.0045</sub></td></tr><tr><td>Limb-darkening parameter q<sub>2</sub></td><td>0.043<sup>+0.068</sup><sub>-0.032</sub></td><td>0.059<sup>+0.280</sup><sub>-0.350</sub></td></tr><tr><td>Constant baseline C<sub>0</sub> (ppm)</td><td>-2577.0<sup>+37</sup><sub>-50</sub></td><td>-2881.0<sup>+30</sup><sub>-30</sub></td></tr><tr><td>Linear trend C<sub>1</sub> (ppm/day)</td><td>-910<sup>+100</sup><sub>-80</sub></td><td>-240<sup>+60</sup><sub>-60</sub></td></tr><tr><td>Ramp magnitude r<sub>0</sub> (ppm)</td><td>758.0<sup>+80</sup><sub>-63</sub></td><td>1319.0<sup>+65</sup><sub>-67</sub></td></tr><tr><td>Ramp time constant r<sub>1</sub> (1/day)</td><td>9.7<sup>+2.0</sup><sub>-1.9</sub></td><td>3.7<sup>+0.3</sup><sub>-0.3</sub></td></tr><tr><td>Spatial position trend</td><td>0.0121<sup>+0.0012</sup><sub>-0.0012</sub></td><td>0.0122<sup>+0.0012</sup><sub>-0.0013</sub></td></tr><tr><td>Spatial point-spread function width trend</td><td>-0.0362<sup>+0.0070</sup><sub>-0.0069</sub></td><td>-0.0385<sup>+0.0071</sup><sub>-0.0072</sub></td></tr><tr><td>Uncertainty scaling factor</td><td>1.2252<sup>+0.0098</sup><sub>-0.0093</sub></td><td>1.2225<sup>+0.0092</sup><sub>-0.0092</sub></td></tr></tbody></table>

**Notes.** The limb-darkening parameters $q_1$ and $q_2$ are as described in Kipping (2013). The polynomial parameters $C_0$ and $C_1$ 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 $r_0e^{(-r_1t)}$ (where $t = 0$ at the start of the observation) as used in Bell et al. (2022). The spatial position and point-spread function 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.

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 zero 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<sup>+0.0060</sup><sub>-0.0025</sub>, 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<sup>+0.00185</sup><sub>-0.00025</sub> 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<sup>+200</sup><sub>-100</sub> 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 $\ell_{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 data set, so we include it in our model. In eclipse, the light from the planet takes about 8.5 s longer to reach the observer than the light from the star. For the average gradient in eclipse ingress and egress of roughly 8 ppm s<sup>-1</sup> (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.
<figure><img src="image_6.png" /><figcaption>**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 nonuniform planetary emission. Green points show the first eclipse, purple points show the second eclipse, and black points show their mean. The blue shaded region labeled “Latitude-Longitude Map Fit” shows the range of fitted phase curves from the *l*<sub>max</sub> = 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 labeled “Longitude-Only Map Fit” shows the fitted eclipse of the *l*<sub>max</sub> = 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 data set, 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.</figcaption></figure>

## 2.2. 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$:

$$F(t) = \sum_{i} c_i f_i(t). \qquad (1)$$

We derive the actual eclipse map $Z(\theta, \phi)$ from the fitted coefficients $c_i$ as a sum of the spherical harmonic basis maps $z_i(\theta, \phi)$:

$$Z(\theta, \phi) = \sum_{i} c_i z_i(\theta, \phi). \qquad (2)$$

Luger et al. (2019) shows the phase curves $f_i(t)$ of each spherical harmonic $z_i(\theta, \phi)$ (where $\theta$ and $\phi$ are longitude and latitude). We model the astrophysical signal using starry (Luger et al. 2019)<fnref n="33" /> 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. To fit the $l_{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 $4l^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 $l_{max} = 2$ eclipse map simultaneously with the orbital, planetary, stellar, and systematic parameters, we fix these additional parameters to their derived values, and refit the map with a variety of methods. We fit maps with the following:

1. $l_{max} = 3$ spherical harmonics.

2. An eigenmap as described in Rauscher et al. (2018).

3. $l_{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 refitted 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.

# 3. Results

## 3.1. Eclipse Map Fitted with Orbit and Systematics

The red line in Figure 1 shows the phase curve resulting from the $l_{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 normalized 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, though 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.

<footnote n="33">³³ starry.readthedocs.io</footnote>
Fitted $ℓ_{max} = 2$ map, orbit, and systematics

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

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. Figure 2 shows the residual difference in ingress and egress between the observations (shown as black, green, and 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 eight 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 nonuniform 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 hot spot near the substellar point, with a

small meridionally averaged shift of $(7.75_{-0.36}^{+0.36})^{\circ}$ eastwards and a shift of $(10.72_{-4.68}^{+4.14})^{\circ}$ southwards.34

We present the maps as a ratio of planetary flux to stellar flux, as a conversion to temperature is more complex and requires modeling or assumptions about the temperature structure. We make a simple estimate of the brightness temperature assuming planetary blackbody emission $B_{T,\lambda}$ in the bandpass $\lambda = 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_{\text{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$_{\odot}$ does affect the value of the brightness temperature. We do not include the 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_{-29.0}^{+23.0})$ K at the substellar point, $(1293.0_{-34.0}^{+27.0})$ K at the equatorial east terminator, $(1114.0_{-36.0}^{+30.0})$ K at the equatorial west terminator, $(1100.0_{-109.0}^{+128.0})$ K at the south pole, and $(1011.0_{-108.0}^{+111.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 data set, 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

<footnote n="34">34 The observations cannot distinguish between an orbit about a vector pointed “up” or “down” in the sky, so our choice of north and south poles is made arbitrarily and kept consistent for plotting. The eastward direction corresponds to the direction of rotation of the planet.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** An eclipse map fitted to the data set in Figure 1, using the parameters listed in Table 1 derived with the $\ell_{max} = 2$ eclipse map model, and constructed from $\ell_{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.</figcaption></figure>

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 longitude-only 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).

## 3.2. Eclipse Mapping Signal

We can quantify the significance of the signal in Figure 2 by comparing the Bayesian evidence of the $\ell_{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

$$ \log L = -\frac{1}{2\sigma^2}\sum_{i=1}^{N}(F(t_i) - D_i)^2 - \frac{N}{2}\log 2\pi\sigma^2, \quad (3) $$

for data $D_i$ at $N$ times $t_i$ fitted by a model $F(t)$ with variance $\sigma^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 $\ell_{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 $\ell_{max} = 2$ map model. This shows how the nonmapping parameters are modified by the 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 data set. This does not immediately imply the presence of 2D information in the data set, 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 $\ell_{max} = 2$ map to the “one-dimensional” (1D) $\ell_{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 data set 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.

## 3.3. Eclipse Map Fitted with $\ell_{max} = 3$ Spherical Harmonics

The spherical harmonic order of the $\ell_{max} = 2$ map fit to the data set 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 $\ell_{max} = 2$ eclipse map model and refit the map with a variety of methods.
Fitted $ℓ_{max}$ = 3 map, with fixed orbit and systematics

<figure><img src="image_9.png" /></figure>

Fitted $ℓ_{max}$ = 2 eigenmap with 6 components, with fixed orbit and systematics

<figure><img src="image_10.png" /><figcaption>**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 optimizing the BIC, using constant systematic and orbital parameters derived by the fit in Figure 3.</figcaption></figure>

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 refitted 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 zero 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_{-0.04}^{+14.79})^\circ$, which is consistent with the $ℓ_{max} = 2$ fit value of $(7.75_{-0.36}^{+0.36})^\circ$. 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})^\circ$. 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 data set.

### 3.4. 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 systematic 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).
<table><caption>**Table 2**
A Comparison of the Models Used to fit the Full Data Set</caption><thead><tr><th>Model</th><th>Order</th><th>Longitude Offset</th><th>Parameters</th><th>χ<sup>2</sup></th><th>ΔBIC</th><th>ΔAIC</th><th>wBIC</th><th>wAIC</th></tr></thead><tbody><tr><td>M1: Fourier series</td><td>n = 2</td><td>(7.3<sup>+0.4</sup><sub>-0.4</sub>)°</td><td>5</td><td>8121.0</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>M2: Uniform latitude</td><td>ℓ<sub>max</sub> = 2</td><td>(7.3<sup>+0.4</sup><sub>-0.4</sub>)°</td><td>5</td><td>8039.9</td><td>6.4</td><td>20.4</td><td>3.9%</td><td>0.0%</td></tr><tr><td>M3: Spherical harmonics</td><td>ℓ<sub>max</sub> = 2</td><td>(7.75<sup>+0.36</sup><sub>-0.36</sub>)°</td><td>10</td><td>8016.5</td><td>19.0</td><td>5.0</td><td>0.0%</td><td>7.5%</td></tr><tr><td>M4: Spherical harmonics</td><td>ℓ<sub>max</sub> = 3</td><td>(0.50<sup>+14.79</sup><sub>-8.04</sub>)°</td><td>17</td><td>8014.9</td><td>80.2</td><td>17.4</td><td>0.0%</td><td>0.0%</td></tr><tr><td>M5: Eigenmap</td><td>ℓ<sub>max</sub> = 2, N<sub>E</sub> = 6</td><td>(7.5<sup>+0.5</sup><sub>-0.5</sub>)°</td><td>8</td><td>8015.5</td><td>0</td><td>0</td><td>96.1%</td><td>92.5%</td></tr></tbody></table>

**Notes.** 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 χ², 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.

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 data set. We find a slight eastward hot-spot offset of $(7.5^{+0.5}_{-0.5})^\circ$. 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”).

### 3.5. Eclipse Maps Fitted with Single Eclipses

The maps in Figures 3 and 5 are fitted to the entire data set, 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 refitting 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 refit 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 hr 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 data set. The map is similar to Figures 3 and 5, but has a larger latitudinal hot-spot shift than the map in Figure 3. 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 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 data set 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. 2024).

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.

### 3.6. Model Selection

To compare our fitted eclipse maps, we calculate the BIC (Schwarz 1978) and the Akaike information criterion (AIC; Akaike 1981). The BIC is

$$ \text{BIC} = k \ln(N) - 2 \ln(\hat{L}), \quad (4) $$

where k is the number of model parameters, N is the number of data points, and $\hat{L}$ is the model likelihood:

$$ \ln(\hat{L}) = -\frac{1}{2} \sum_{i}^{N} \left( \frac{M_i - D_i}{\sigma_i} \right)^2, \quad (5) $$

for a sum over all N data points $D_i$ with uncertainty $\sigma_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

$$ \text{AIC} = 2k - 2 \ln(\hat{L}). \quad (6) $$

Due to the large number of points N in this data set, the BIC applies a stronger penalty for each parameter than the AIC. We also present the “weights” of each model, wBIC and wAIC, which are the relative probabilities of each compared model in
<figure><img src="image_11.png" /><figcaption>Fitted $l_{max}$ = 2 map, second eclipse only, with fixed orbit and systematics</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fitted $l_{max}$ = 2 map, both eclipses only, with fixed orbit and systematics</figcaption></figure>

<figure><img src="image_13.png" /></figure>

**Figure 6.** Eclipse maps fitted using each eclipse separately. First row: an eclipse map fitted using spherical harmonics up to $l_{max}$ = 2, using a ~2.2 hr 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 refitted), listed in Table 1. Only the dayside of the map is shown as there is no information about the nightside in this limited data set. 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 hr 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.

a set, based on their ΔBIC or ΔAIC. The wBIC for a particular model is

$$wBIC = \frac{\exp(-\Delta BIC/2)}{\sum_{r=1}^{R} \exp(-\Delta BIC_r/2)}, \quad (7)$$

where ΔBIC is compared to the best-performing model, and there are R total models being compared with a score of ΔBIC_r each. wAIC depends on ΔAIC in the same way. Table 2 summarizes our compared models, showing their $\chi^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 $\chi^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 $\chi^2$ value than the Fourier series model, as it includes the effect of longitudinal map structure on the eclipse shape. However, its $\chi^2$ value is worse than the fully 2D eclipse map fits, as shown by its poor fit in Figure 2. Interestingly, due to its smaller number of parameters, it achieves the second-best 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 penalized less strongly. We can confidently reject models M1 and M2 compared to the 2D eclipse mapping models.

Model M3 is the $l_{max} = 2$ eclipse mapping model in Figure 3. It achieves a better $\chi^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 $\chi^2$ values, but model M3 uses two more parameters. It has a more similar AIC score to the eigenmap, where these additional parameters are penalized less heavily.

Model M4 is the $l_{max} = 3$ eclipse mapping model in Figure 5. It achieves the best $\chi^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 $l_{max} = 3$ map in Figure 5.

Model M5 is the $l_{max} = 2$ eigenmap using six basis maps. It achieves a similar $\chi^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 summarize 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 $l_{max} = 2$ eclipse map (M3) is favored over the longitude-only eclipse map (M2), showing the presence of latitudinal information in the data.

3. The $l_{max} = 3$ eclipse map achieves a better $\chi^2$ value than model M3, but at the cost of many more parameters, so has worse BIC and AIC scores.

4. The $l_{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).

# 4. Comparison to General Circulation Models

In this section, we interpret the eclipse maps by comparing them to 3D 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, 2018a, 2018b; Simulation 31 in Bell et al. 2024), expeRT/MITgcm (Carone et al. 2020; Schneider et al. 2022; Simulation 20 in Bell et al. 2024), the Generic Planetary Climate Model (Generic PCM; Teinturier et al. 2024; Simulation 7 in Bell et al. 2024), and the RM-GCM (Rauscher & Menou 2012; Roman & Rauscher 2017; Roman et al. 2021; Simulation 24 in Bell et al. 2024).

The THOR and RM-GCM simulations use semi-gray radiative transfer, while the PCM and expeRT/MITgcm simulations use multiband correlated-k schemes. The THOR, RM-GCM, and PCM simulations feature clouds (on the nightside only for THOR), whereas the expeRT/MITgcm simulation is free of clouds. The output from each model was postprocessed 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.

## 4.1. Observable Features of GCMs

Figure 7 shows thermal emission maps postprocessed 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 chevron-like 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; Showman & Polvani 2011; Lewis & Hammond 2022). THOR, expeRT/MITgcm, and the RM-GCM each show cold “lobes” on the nightside, which are likely associated with stationary equatorial Rossby waves (Matsuno 1966; Showman & Polvani 2011; Lewis & Hammond 2022). A Rossby wave–like structure also appears to dominate the dayside emission in the RM-GCM simulation.

### 4.1.1. Observable and Null Maps

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
<figure><img src="image_14.png" /><figcaption>**Figure 7.** 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 modeled outgoing longwave radiation (OLR) in the 5–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.</figcaption></figure>

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 large-scale 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 in Section 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 the $ℓ_{max} = 2$ and $ℓ_{max} = 3$ maps have no null space for this observation. The nonzero 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.

### 4.1.2. 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 substellar 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, for example, Figure 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.

## 4.2. 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 $\int_{-π/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 large-scale 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})^\circ$. 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})^\circ$, $(6.67^{+5.77}_{-4.32})^\circ$, and $(4.86^{+3.6}_{-2.88})^\circ$, 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 breaks degeneracies with the latitudinal structure (Boone et al. 2024). The $ℓ_{max} = 3$ map fit to the whole data set shown in Figure 5 has a map offset of $(0.50^{+14.79}_{-8.04})^\circ$ that is consistent with the $ℓ_{max} = 3$ eclipse-only fits, as well as with the $ℓ_{max} = 2$ map fits. As discussed in Section 3.3, this map has a hot-spot offset
GCM maps along equator and through substellar point

<figure><img src="image_15.png" /></figure>

Compare $l_{max} = 2$ map fit to GCM $l_{max} = 2$ maps

<figure><img src="image_16.png" /></figure>

Compare $l_{max} = 2$ eigenmap fit to GCM $l_{max} = 2$ eigenmaps, both with 6 components

<figure><img src="image_17.png" /></figure>

Compare $l_{max} = 3$ map fit to GCM $l_{max} = 3$ maps

<figure><img src="image_18.png" /><figcaption>**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 $l_{max} = 2$ map in Figure 3 compared to the $l_{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 $l_{max} = 3$ map in Figure 3 compared to the $l_{max} = 3$ GCM maps in the third row of Figure 7.</figcaption></figure>

with a median position on the substellar point, but a nonzero phase offset of around 7°, consistent with the overall phase curve, showing the difference between these two metrics. The

$l_{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.
<figure><img src="image_19.png" /><figcaption>**Figure 9.** Left: the phase curve offset derived from the *n* = 2 fit using a Fourier series model (giving the same result as Bell et al. 2024), compared to the phase curve offset simulated from each GCM. Right: the (meridionally averaged) longitudinal offset in each 2D map. The *ℓ*<sub>max</sub> = 2 and *ℓ*<sub>max</sub> = 3 offsets (using all the data set) are from the map fits in Figures 3 and 5. The *ℓ*<sub>max</sub> = 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 gray points show the fits to the first and second eclipses). The eclipse-only eigenmap fits and the eclipse-only *ℓ*<sub>max</sub> = 3 map fits are not plotted in full elsewhere. The GCM points show the cos φ-weighted offset for the GCMs truncated to the *ℓ*<sub>max</sub> = 2 and *ℓ*<sub>max</sub> = 3 representations shown in Figure 7, as well as for the maps with no truncation.</figcaption></figure>

The longitudinal offsets of the *n* = 2 Fourier series fit and the
*ℓ*<sub>max</sub> = 2 map fit to the full data set appear very precise
compared to the *ℓ*<sub>max</sub> = 3 map fit to the full data set. This
reflects the choice of fitting 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
*ℓ*<sub>max</sub> = 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 *ℓ*<sub>max</sub> = 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° based on the mapping method used. The
map offsets in the GCMs used for comparison are larger,
varying between around 10° and 30° depending on the
mapping basis. We conclude that while the $\ell_{\text{max}}$ = 2 map
and eigenmap are the best constrained maps, they may be
discarding important structures from the real map. The
$\ell_{\text{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 $\ell$ = 3 modes.

5. 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 $\ell_{\text{max}}$ = 2 spherical
harmonics simultaneously with the parameters of the system
and the systematic parameters of the instrument. 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 data set 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 $\ell_{\text{max}}$ = 2 map found a small eastward hot-spot shift of $(7.75_{-0.36}^{+0.36})^{\circ}$ (defined using a meridional average), compared to an eastward phase shift of $(7.3_{-0.4}^{+0.4})^{\circ}$ derived from fitting a Fourier series to the phase curve. This $\ell_{\text{max}}$ = 2 fit finds a sharply peaked latitudinal structure, with a small latitudinal offset of $(-10.72_{-4.68}^{+4.14})^{\circ}$.

Figure 5 shows an $\ell_{\text{max}}$ = 3 eclipse map fitted with the nonmapping parameters fixed to those derived from the $\ell_{\text{max}}$ = 2 map fit. This found a smaller but more uncertain hot-spot shift of $(0.5_{-8.04}^{+14.79})^{\circ}$ degrees east, and a flatter latitudinal structure near the equator than the $\ell_{\text{max}}$ = 2 map. It achieved a better $\chi^2$ than the $\ell_{\text{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 $\ell_{\text{max}}$ = 2
basis. Figure 5 shows that this derived a map very similar to the
original $\ell_{\text{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<sup>+23.0</sup><sub>-29.0</sub>) K for the
ℓ<sub>max</sub> = 2 map. This is consistent with the brightness temper-
ature of the eigenmap at the substellar point of (1783.0<sup>+51.0</sup><sub>-70.0</sub>) K,
as well as with the brightness temperature of the ℓ<sub>max</sub> = 3 map,
which is (1822.0<sup>+75.0</sup><sub>-92.0</sub>) K at the substellar point.

We tested the mapping information in each eclipse
separately by refitting 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 $l_{max} = 2$ map is largely consistent with the $l_{max} = 2$ GCM representations. The $l_{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 $l_{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 $l_{max} = 2$ spherical harmonic components of the planetary emission. Our fiducial map is the $l_{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})^\circ$ eastward (shown in Figure 9), compared to a shift in the peak of the phase curve of $(7.3_{-0.4}^{+0.4})^\circ$ eastward.

We also highlight our statistically preferred map (according to the BIC and AIC), which is the eigenmap fitted using $l_{max} = 3$ and $N_E = 6$ basis maps; this produces a very similar map to the fiducial $l_{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 $l_{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 data set 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 & Komaček 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 data set 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 3D 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.

## Acknowledgments

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 NASAs 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 Région Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the program 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. This research was funded in whole, or in part, by a UKRI STFC studentship.

## Appendix A

### General Circulation Models and Postprocessing

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, 2018b). The simulation presented here uses semi-gray 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 were averaged over the last 500 days. The emission from each column was calculated by postprocessing the results with a multiwavelength radiative transfer model (Mendonça et al. 2015), assuming $1 \times$ 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 nongray 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 \times$ solar metallicity, and was run for
1500 days, with the results shown here averaged over the last 100 days. The spectrally resolved emission was postprocessed using petitRADTRANS (Mollière et al. 2019) and prt_pha-securve (Schneider et al. 2022) using a spectral resolution of R = 100.

The third is the *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. 2024) and the planets of the solar system (Spiga et al. 2020; Turbet et al. 2021). 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 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₂SiO₄ clouds with radius 1 µm. The simulation was postprocessed 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), and Rauscher & Menou (2012), and has been applied to investigations of exoplanets using semi-gray radiative transfer and aerosol scattering (Roman & Rauscher 2017; Roman et al. 2021). The simulation presented here has 1× 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 postprocessed following Zhang et al. (2017) and Malsky et al. (2021) to produce the spectrally resolved emission.

## Appendix 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 $l_{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 $\chi^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 (colored 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 data set.

The Fourier series model finds a statistically significantly smaller planetary radius, higher orbital inclination, and larger semimajor 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 $l_{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><img src="image_20.png" /><figcaption>**Figure 10.** A corner plot of the posterior distributions of the orbital, planetary, stellar, and systematic parameters fitted to the observations in Figure 1, 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 *ℓ*<sub>max</sub> = 2 eclipse map model in Figure 3. Table 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 *ℓ*<sub>max</sub> = 2 eclipse map model are more accurate, due to its more accurate model of the eclipse shape. This figure is part of a figure set available in the online journal, which shows alternative fits that also include the eccentricity, argument of periastron, and stellar radius as free parameters, to explain why we fix these parameters to constant values in our fiducial model.</figcaption></figure>

(The complete figure set (3 images) is available.)
Fitted $l_{max} = 2$ map, with fixed orbit and systematics from Fourier series fit

<figure><img src="image_21.png" /><figcaption>**Figure 11.** An $l_{max} = 2$ eclipse map fitted using the orbital and systematic parameters derived by fitting the Fourier series model, shown in Figure 10 and listed in Table 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 $l_{max} = 2$ eclipse map.</figcaption></figure>

## ORCID iDs

Mark Hammond https://orcid.org/0000-0002-6893-522X

Taylor J. Bell https://orcid.org/0000-0003-4177-2149

Ryan C. Challener https://orcid.org/0000-0002-8211-6538

Megan Weiner Mansfield https://orcid.org/0000-0003-4241-7413

Isaac Malsky https://orcid.org/0000-0003-0217-3880

Emily Rauscher https://orcid.org/0000-0003-3963-9672

Jacob L. Bean https://orcid.org/0000-0003-4733-6532

Ludmila Carone https://orcid.org/0000-0001-9355-3752

João M. Mendonça https://orcid.org/0000-0002-6907-4476

Lucas Teinturier https://orcid.org/0000-0002-0797-5746

Nicolas Crouzet https://orcid.org/0000-0001-7866-8738

Laura Kreidberg https://orcid.org/0000-0003-0514-1147

Giuseppe Morello https://orcid.org/0000-0002-4262-5661

Jasmina Blecic https://orcid.org/0000-0002-0769-9614

Jean-Michel Désert https://orcid.org/0000-0002-0875-8401

Christiane Helling https://orcid.org/0000-0002-8275-1371

Karan Molaverdikhani https://orcid.org/0000-0002-0502-0428

Matthew C. Nixon https://orcid.org/0000-0001-8236-5553

Benjamin V. Rackham https://orcid.org/0000-0002-3627-1676

Jingxuan Yang https://orcid.org/0009-0006-2395-6197

## References

Adams, F. C., & Laughlin, G. 2006, ApJ, 649, 1004

Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. 2004, MWRv, 132, 2845

Akaike, H. 1981, *J. Econometrics*, 16, 3

Allard, F., & Hauschildt, P. H. 1995, ApJ, 445, 433

Bell, T., Ahrer, E.-M., Brande, J., et al. 2022, *JOSS*, 7, 4503

Bell, T. J., Crouzet, N., Cubillos, P. E., et al. 2024, *NatAs*, in press

Beltz, H., Rauscher, E., Roman, M. T., & Guilliat, A. 2021, AJ, 163, 35

Boone, S., Grant, D., & Hammond, M. 2024, *MNRAS*, 528, 596

Bouwman, J., Kendrew, S., Greene, T. P., et al. 2023, *PASP*, 135, 038002

Burnham, K. P., & Anderson, D. R. 2004, *Sociol. Methods Res.*, 33, 261

Carone, L., Baeyens, R., Mollière, P., et al. 2020, *MNRAS*, 496, 3582

Challener, R. C., & Rauscher, E. 2022, AJ, 163, 117

Challener, R. C., & Rauscher, E. 2023, AJ, 166, 176

Challener, R. C., Welbanks, L., & McGill, P. 2023, AJ, 166, 251

Charnay, B., Meadows, V., Misra, A., Leconte, J., & Arney, G. 2015, ApJL, 813, L1

Coulombe, L.-P., Benneke, B., Challener, R., et al. 2023, *Natur*, 620, 292

Cowan, N. B., & Agol, E. 2008, ApJL, 678, L129

Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3

de Wit, J., Gillon, M., Demory, B. O., & Seager, S. 2012, A&A, 548, A128

Falco, A., Zingales, T., Pluriel, W., & Leconte, J. 2022, A&A, 658, A41

Foreman-Mackey, D. 2016, *JOSS*, 1, 24

Gillon, M., Triaud, A., Fortney, J., et al. 2012, A&A, 542, A4

Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hammond, M., & Lewis, N. T. 2021, *PNAS*, 118, e2022705118

Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377

Hellier, C., Anderson, D., Cameron, A. C., et al. 2011, arXiv:1104.2823

Horne, K. 1985, *MNRAS*, 213, 129

Hoskins, B. J., & Simmons, A. J. 1975, *QJRMS*, 101, 637

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Kendrew, S., Scheithauer, S., Bouchet, P., et al. 2015, *PASP*, 127, 623

Kipping, D. M. 2013, *MNRAS*, 435, 2152

Kokori, A., Tsiaras, A., Edwards, B., et al. 2023, *ApJS*, 265, 4

Lewis, N. T., & Hammond, M. 2022, ApJ, 941, 171

Luger, R., Agol, E., Foreman-Mackey, D., et al. 2019, AJ, 157, 64

Luger, R., Foreman-Mackey, D., Hedges, C., & Hogg, D. W. 2021, AJ, 162, 123

Majeau, C., Agol, E., & Cowan, N. B. 2012, ApJL, 747, L20

Malsky, I., Rauscher, E., Kempton, E. M. R., et al. 2021, ApJ, 923, 62

Mansfield, M., Schlawin, E., Lustig-Yaeger, J., et al. 2020, *MNRAS*, 499, 5151

Matsuno, T. 1966, *JMeSJ*, 44, 25

McEwen, J. D., & Wiaux, Y. 2011, *ITSP*, 59, 5876

Mendonça, J. M., Grimm, S. L., Grosheintz, L., & Heng, K. 2016, ApJ, 829, 115

Mendonça, J. M., Malik, M., Demory, B.-O., & Heng, K. 2018a, AJ, 155, 150

Mendonça, J. M., Read, P. L., Wilson, C. F., & Lee, C. 2015, *P&SS*, 105, 80

Mendonça, J. M., Tsai, S.-m., Malik, M., Grimm, S. L., & Heng, K. 2018b, ApJ, 869, 107

Menou, K., & Rauscher, E. 2009, ApJ, 700, 887

Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67

Morello, G. 2018, AJ, 156, 175

Morello, G., Tsiaras, A., Howarth, I. D., & Homeier, D. 2017, AJ, 154, 111

Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

Placek, B., Angerhausen, D., & Knuth, K. H. 2017, AJ, 154, 154

Rauscher, E., & Menou, K. 2010, ApJ, 714, 1334

Rauscher, E., & Menou, K. 2012, ApJ, 750, 96

Rauscher, E., Menou, K., Seager, S., et al. 2007, ApJ, 664, 1199

Rauscher, E., Suri, V., & Cowan, N. B. 2018, AJ, 156, 235

Rogers, T. M., & Komacek, T. D. 2014, ApJ, 794, 132

Roman, M., & Rauscher, E. 2017, ApJ, 850, 17

Roman, M. T., Kempton, E. M. R., Rauscher, E., et al. 2021, ApJ, 908, 101

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, *PeerJ Comp. Sci.*, 2, e55

Scandariato, G., Singh, V., Kitzmann, D., et al. 2022, A&A, 668, A17

Schneider, A. D., Carone, L., Decin, L., et al. 2022, A&A, 664, A56

Schwarz, G. 1978, *AnSta*, 6, 461

Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71

Skinner, J., & Cho, J. Y. 2022, *MNRAS*, 511, 3584

Spiga, A., Guerlet, S., Millour, E., et al. 2020, *Icar*, 335, 113377

Stern, S. A. 1992, *ARA&A*, 30, 185

Stock, J. W., Kitzmann, D., Patzer, A. B. C., & Sedlmayr, E. 2018, *MNRAS*, 479, 865
Stone, M. 1977, J. R. Stat., 39, 44

Taylor, J., Parmentier, V., Irwin, P. G., et al. 2020, MNRAS, 493, 4342

Teinturier, L., Charnay, B., Spiga, A., et al. 2024, A&A, 683, A231

Turbet, M., Bolmont, E., Chaverot, G. et al. 2021, Natur, 598, 276

Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112

Warner, B., Robinson, E. L., & Nather, R. E. 1971, MNRAS, 154, 455

Williams, P. K., Charbonneau, D., Cooper, C. S., Showman, A. P., & Fortney, J. J. 2006, ApJ, 649, 1020

Yang, J., Irwin, P. G., & Barstow, J. K. 2023, MNRAS, 525, 5146

Zhang, J., Kempton, E. M. R., & Rauscher, E. 2017, ApJ, 851, 84