# Identifying and Fitting Eclipse Maps of Exoplanets with Cross-Validation

Mark Hammond,1 Neil T. Lewis,2 Sasha Boone,1 Xueqing Chen,3 João M. Mendonça,4
Vivien Parmentier,6 Jake Taylor,5 Taylor Bell,7,8 Leonardo dos Santos,9
Nicolas Crouzet,10 Laura Kreidberg,11 Michael Radica,12 and Michael Zhang13
1Atmospheric, Oceanic and Planetary Physics, Department of Physics, University
of Oxford, Oxford, UK
2Department of Physics and Astronomy, University of Exeter, Exeter, UK
3Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh, EH9 3HJ, UK
4National Space Institute, Technical University of Denmark, Elektrovej, 2800
Kgs. Lyngby, Denmark
5Astrophysics, Department of Physics, University of Oxford, Oxford, UK
6Université Côte d’Azur, Av. Valrose, 06000 Nice, France
7BAER Institute, NASA Ames Research Center, Moffett Field, CA, USA
8Space Science and Astrobiology Division, NASA Ames Research Center, Moffett
Field, CA, USA
9Space Telescope Science Institute, Baltimore, MD, USA
10Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The
Netherlands
11Max Planck Institute for Astronomy, Heidelberg, Germany
12Trottier Institute for Research on Exoplanets, Université de Montréal,
Montréal, Québec, Canada
13Department of Astronomy and Astrophysics, The University of Chicago,
Chicago, IL 60637, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Eclipse mapping uses the shape of the eclipse of an exoplanet to measure its
two-dimensional structure. Light curves are mostly composed of longitudinal
information, with the latitudinal information only contained in the brief
ingress and egress of the eclipse. This imbalance can lead to a spuriously
confident map, where the longitudinal structure is constrained by out-of-
eclipse data and the latitudinal structure is wrongly determined by the priors
on the map. We present a new method to address this issue. The method tests
for the presence of an eclipse mapping signal by using k-fold cross-validation
to compare the performance of a simple mapping model to the null hypothesis of
a uniform disk. If a signal is found, the method fits a map with more degrees
of freedom, optimising its information content. The information content is
varied by penalising the model likelihood by a factor proportional to the
spatial entropy of the map, optimised by cross-validation. We demonstrate this
method for simulated datasets then apply it to three observational datasets.
The method identifies an eclipse mapping signal for JWST MIRI/LRS observations
of WASP-43b but does not identify a signal for JWST NIRISS/SOSS observations
of WASP-18b or Spitzer Space Telescope observations of HD 189733b. It is
possible to fit eclipse maps to these datasets, but we suggest that these maps
are overfitting the eclipse shape. We fit a new map with more spatial freedom
to the WASP-43b dataset and show a flatter east-west structure than previously
derived.

###### keywords:

methods: observational — planets and satellites: atmospheres

††pubyear: 2024††pagerange: Identifying and Fitting Eclipse Maps of Exoplanets
with Cross-Validation–A

##  1 Introduction

Planetary atmospheres are controlled by multi-dimensional processes.
Observations of transiting exoplanets are generally limited to measurements of
point sources, where the flux from the entire stellar and planetary system is
observed at once. For such a system, eclipse mapping is currently the only
method by which both two-dimensional longitudinal (east-west) and latitudinal
(north-south) information can be measured over the surface of an exoplanet.
This information is derived from the shape of the eclipse of an exoplanet by
its star. The eclipse shape is a function of the spatial distribution of flux
from the planet and of the geometry of the coverage of the planet by its star.
If we know the geometry of the eclipse coverage, or can fit a model of it, we
can derive a map of the flux from the planet.

Eclipse mapping has been used in other contexts to derive images of the
surface brightness of accretion disks (Horne, 1985) and maps of the single-
scattering albedo of Pluto and Charon (Buie et al., 1992). Rauscher et al.
(2007) showed the possibility of mapping exoplanets with sufficiently high-
precision observations, highlighting the suitability of JWST in particular.
Eclipse mapping has previously been applied to observations of three “hot
Jupiter” exoplanets: HD 189733b with the Spitzer Space Telescope (de Wit et
al., 2012; Majeau et al., 2012), WASP-18b with JWST NIRISS/SOSS (Coulombe et
al., 2023), and WASP-43b with JWST MIRI/LRS (Hammond et al., 2024). Majeau et
al. (2012) used spherical harmonics up to first order to map HD 189733b. They
found that the brightest part of the planet was shifted slightly eastwards,
and that its latitudinal position was consistent with the equator. de Wit et
al. (2012) used more flexible mapping functions to derive a more localised
hot-spot from this same dataset and highlighted degeneracies between this map
and the orbital parameters of the planet. (Coulombe et al., 2023) derived a
map of WASP-18b with a notably uniform longitudinal structure near the
substellar point, but they could not find evidence for latitudinal information
in the dataset. Hammond et al. (2024) detected latitudinal information in the
eclipse shape of WASP-43b and found an eastward longitudinal hot-spot shift.

Deriving reliable eclipse maps from observational datasets has two intrinsic
issues. First, the problem is ill-posed – one spatial pattern can produce the
same flux time-series as another, or even no flux at all (Cowan & Agol, 2008;
Challener & Rauscher, 2023). Alongside the uncertainty inherent in
observational data, this means that an observed light curve implies a range of
possible maps instead of a single unique solution. Second, the data that
contain two-dimensional information are normally a very small part of an
overall time-series observation – the ingress and egress of an eclipse being
typically at least an order of magnitude shorter than the whole orbit of a
planet. These issues lead to two statistical challenges.

The first challenge is to identify when it is justifiable to fit an eclipse
map to a particular observation. In general, an eclipse map model has more
degrees of freedom than a model of a planet with uniform emission. It will,
therefore, normally produce a better fit to the data in an eclipse, even when
the deviations from the shape of a uniform eclipse are driven by random noise,
or when unresolved systematic errors produce an eclipse shape that does not
depend on the true emission map. This means that it is easy to fit spurious
eclipse maps to noisy or inaccurate datasets. We therefore present a
statistical test using k-fold cross-validation to determine when an “eclipse
mapping signal” is present. This statistical test avoids both underfitting
(where the model cannot match the observations), and overfitting (where the
model derives spurious information from noise).

The second statistical challenge is how much complexity to allow in the fitted
map, if an eclipse mapping signal is indeed present. Fitting with a high
degree of spatial freedom results in highly uncertain maps with spurious
small-scale features due to uncertainty and degeneracy in the observational
data. Section 2 therefore presents a method to derive the appropriate spatial
complexity of the fitted map, again using k-fold cross-validation. This method
is based on the pixel sampling and entropy maximisation methods presented in
Horne (1985) for eclipse mapping of accretion disks.

We optimise the information content of our fitted maps like Chen et al.
(2024), by penalising the likelihood of the fitted model by a regularisation
parameter proportional to the entropy of the map. The entropy is a function of
the brightness of the pixels composing the map and corresponds to the
information content of the map relative to a uniformly bright image. Welbanks
et al. (2023) shows a similar method for fitting spectroscopic models to
exoplanet observations using leave-one-out cross validation. Challener &
Rauscher (2023) showed the application of leave-one-out cross-validation to
eclipse mapping, but we chose not to use the leave-one-out method as we seek
to test the predictive power of a model for a longer duration of omitted data.

We demonstrate these methods with simulated data in Section 3, showing that
they correctly identify when fitting a simple low-order eclipse map is
justified. We show that alternative statistical metrics are not as reliable
for identifying the presence of a robust signal. We then show that, if a
mapping signal is present, the method can then identify the appropriate
information content for a map composed of higher-order spherical harmonics,
depending on the spatial scale of the true map and the precision of the data.

We then demonstrate this method for real observational data in Section 4,
applying it to observations of WASP-43b with JWST MIRI/LRS, WASP-18b with JWST
NIRISS/SOSS, and HD 189733b with the Spitzer Space Telescope. The method
identifies a clear eclipse mapping signal for observations of WASP-43b with
JWST MIRI/LRS, consistent with the findings of Hammond et al. (2024). We then
derive a new eclipse map for WASP-43b, composed of higher-order spherical
harmonics than the map presented in Hammond et al. (2024). Our method does not
identify a significant eclipse mapping signal for observations of WASP-18b
with NIRISS SOSS or the Spitzer Space Telescope observations of HD 189733b;
this differs from the conclusions of Coulombe et al. (2023), Majeau et al.
(2012), and de Wit et al. (2012).

We conclude that our proposed method can reliably determine the presence of an
eclipse mapping signal in observational data, and then determine the
appropriate spatial scale of the map to fit to this data, addressing the two
statistical issues presented above.

Figure 1: First panel: the emission map from the THOR GCM used to simulate the
time-series observations in Figures 2 and 3. Second panel: the true map
represented by $\ell_{\mathrm{max}}=2$ spherical harmonics, showing how these
are not sufficient to capture its spatial structure. Third panel: the true map
represented by $\ell_{\mathrm{max}}=4$ spherical harmonics, showing how this
captures most of the true structure. The brightest point of the true map is at
$40^{\circ}$ east of the substellar point. The brightest point of the
$\ell_{\mathrm{max}}=2$ representation is incorrect, at $24^{\circ}$ east. The
$\ell_{\mathrm{max}}=4$ representation is more accurate, with its brightest
point at $43^{\circ}$ east.

##  2 Methods

This section describes the methods we use to simulate and retrieve eclipse
mapping data. It lists the models used to simulate observations of thermal
emission light curves, our sources of real observational light curves, the
models used to fit these light curves with eclipse maps, and the cross-
validation metric we propose to compare the quality of different models.

###  2.1 Simulated time-series thermal emission

We simulate time-series thermal emission to represent the typical precision
and cadence of observations suitable for eclipse mapping with JWST. We set up
our simulation to emulate observations of the hot Jupiter WASP-43b (Hellier et
al., 2011) with the MIRI/LRS instrument on JWST (Bell et al., 2024). Our
conclusions do not depend on this particular choice, and we vary our model
choices later. WASP-43b orbits a K7, 0.717 M$\odot$ star with an orbital
period of 19.2 hours, a mass of 2.034 $M_{J}$, and a semi-major axis of 0.0153
AU (Gillon et al., 2012). It has an estimated equilibrium temperature of 1440
K (Blecic et al., 2014).

We use a simulation from the General Circulation Model (GCM) THOR (Mendonça et
al., 2016; Deitrick et al., 2020) of a cloudless atmosphere on WASP-43b to
model observations. This simulation was previously used to study the
atmospheric temperature structure, cloud cover, and chemistry of WASP-43b
(Mendonça et al., 2018a, b). THOR is based on a dynamical core that solves the
non-hydrostatic compressible Euler equations on an icosahedral grid. The
simulation uses a simple two-band formulation of radiative transfer calibrated
to reproduce the results from more complex non-grey models of WASP-43b
(Kataria et al., 2015; Malik et al., 2017). More details about the model
configuration can be found in Mendonça et al. (2018a). We are not concerned
about the time-varying behaviour of the atmosphere (which is small in any
case) as we are using this as an exemplar case, so we post-process a single
snapshot of the GCM temperature field.

To simulate the observed thermal emission in the MIRI/LRS bandpass, the GCM
simulation was post-processed using the HELIOS code (Malik et al., 2017) as
described in Bell et al. (2024) to produce a map of the thermal emission
averaged over the entire bandpass. We do not represent limb darkening for this
latitude-longitude map, as we fit a two-dimensional map that varies only due
to the rotation and occultation of the planet.

This simulation results in the map of thermal emission shown in the left-hand
panel of Figure 1. The figure also shows representations of this “true map”
using spherical harmonics up to a given order $\ell_{\mathrm{max}}$,
calculated using starry (Luger et al., 2019). This shows how an
$\ell_{\mathrm{max}}=2$ representation fails to capture the spatial structure
of the map; its brightest point is on the equator $24^{\circ}$ east of the
substellar point, compared to the position of the brightest point in the
original map at $40^{\circ}$ east. This offset is due to the low-order
representation omitting many of the higher-order structures in the true map
(Hammond et al., 2024).

The $\ell_{\mathrm{max}}=4$ representation in Figure 1 captures most of the
spatial structure, with its brightest point at $43^{\circ}$ east, more closely
matching the brightest point of the original map. This implies that fitting an
eclipse map composed of $\ell_{\mathrm{max}}=2$ harmonics would not produce an
accurate fit to the real structure, while the $\ell_{\mathrm{max}}=4$
representation would in theory be able to match the real structure closely.
For real observations where we do not know the spatial structure in advance,
it is desirable to use a high degree of spatial freedom, constrained only by
observational precision.

To simulate the light curve with starry (Luger et al., 2019), we use the full
map represented precisely by harmonics up to $\ell_{\mathrm{max}}=10$.
WASP-43b is known to be tidally locked to its star (Hammond et al., 2024) so
the fitted eclipse map corresponds to its permanent dayside. The simulated
time series runs from orbital phase -0.6 to 0.4, covering a single orbit with
an eclipse and a transit. It has 7000 data points, giving a cadence of 10.04
seconds per point, similar to the cadence of the JWST observations of WASP-43b
in Bell et al. (2024).

We add Gaussian noise with a standard deviation $\sigma$, testing data with
$\sigma=$ 150 ppm, 250 ppm, and 2000 ppm. We chose the value of 250 ppm to be
comparable to the precision achieved by averaging the two eclipses of WASP-43b
observed with JWST MIRI/LRS in Bell et al. (2024) and mapped in Hammond et al.
(2024). The value of 150 ppm represents a best-case scenario that could be
achieved by averaging five such eclipses. The value of 2000 ppm represents
low-precision data where there is no discernible mapping signal in the eclipse
shape, but the rest of the phase curve can be resolved.

###  2.2 Observed time-series thermal emission

We also apply our method to three different observational datasets, plotted
later in Section 4. They are all time-series observations of thermal emission
from hot Jupiters. The first is a JWST MIRI/LRS observation of WASP-43b,
averaged from 5 to 10.5 $\mu$m (Bell et al., 2024). The second is a JWST
NIRISS/SOSS observation of WASP-18b, averaged from 0.85 to 2.85 $\mu$m
(Coulombe et al., 2023). The final is a Spitzer Space Telescope observation of
HD 189733b at 8 $\mu$m, compiled from seven separate eclipse observations
(Majeau et al., 2012). Although the JWST datasets contain spectroscopic
information, we average over all available wavelengths to achieve the highest
possible precision.

A systematic model of instrumental effects and an astrophysical model of the
orbital parameters was fitted to each dataset in each of the original studies
to produce the final light curves. In this study, we do not explore the
possibility of degeneracies between the systematic model, the astrophysical
model, and the map model that we fit. Instead, we assume that the first two
models are correctly known in advance, and only fit for the map (or a Fourier
series model). In reality, there may be degeneracies between these models
(Hammond et al., 2024). The JWST MIRI/LRS light curve we use for WASP-43b had
its systematic and astrophysical model fitted alongside an
$\ell_{\mathrm{max}}=2$ eclipse map, which should at least partially mitigate
such degeneracies (Hammond et al., 2024).

###  2.3 Fitting thermal emission with a Fourier series

The simplest model of a light curve of a planetary system uses a Fourier
series for the variation in flux with time, multiplied by the eclipse and
transit shapes of a planet emitting uniformly from its surface. We refer to
this as a “Fourier series model” with flux $F(t)$:

$F(t)=\biggl{(}A_{0}+\sum_{n=1}^{n_{\mathrm{max}}}\biggl{[}A_{n}\sin(2\pi
n(t-t_{0})/P)\\\ +B_{n}\cos(2\pi
n(t-t_{0})/P)\biggr{]}\bigg{)}F_{p,0}(t)+F_{s}(t),$ (1)

where $t_{0}$ is the transit time, $P$ is the orbital period, $A_{n}$ and
$B_{n}$ are the fitted parameters up to order $n_{\mathrm{max}}$, $F_{p,0}$ is
the time-series flux from a planet with uniform emission, and $F_{s}$ is the
stellar flux including the transit model. We refer to a Fourier series model
by its maximum order, so a model with maximum order $n_{\mathrm{max}}=2$ is
$\mathcal{M}(n_{\mathrm{max}}=2)$.

This fits the out-of-eclipse light curve as well as any eclipse mapping model.
However, it does not represent the effect of non-uniform planetary flux on the
eclipse ingress and egress shapes. We will use this as our null hypothesis
when fitting datasets. To justify fitting an eclipse map, the map must fit the
dataset better than this null hypothesis. Crucially, the map must not achieve
this better fit by overfitting to noise in the eclipse shape and must actually
be a better model of the true eclipse shape. We use an
$\mathcal{M}(n_{\mathrm{max}}=2)$ Fourier series model as our null hypothesis
as this is sufficient to capture almost all the out-of-eclipse variation of a
phase curve, due to the rapid decreasing contribution to a light curve from
higher-order modes (Cowan & Agol, 2008).

###  2.4 Fitting eclipse maps

For data of sufficient precision, we can fit an eclipse map model of the two-
dimensional flux from the dayside of the planet using the information encoded
in the shape of the eclipse. We refer to an “eclipse map model” by its maximum
spherical harmonic degree $\ell_{\mathrm{max}}$, so an eclipse map model with
maximum degree $\ell_{\mathrm{max}}=2$ is referred to as
$\mathcal{M}(\ell_{\mathrm{max}}=2)$.

We define the map in terms of pixels evenly spaced in planetary surface area,
known as a Mollweide projection (Snyder, 1987). These pixels uniquely define a
spherical harmonic representation which continuously covers the mapped
surface. starry uses this spherical harmonic representation to calculate the
time-series emission from the map. Two-dimensional information is only
accessible in the part of the surface eclipsed in both ingress and egress;
this approximately corresponds to the dayside, corrected slightly for
planetary rotation during the eclipse. We simplify this by plotting the
posterior distribution of the fitted map on the dayside only, and plot the
median fitted map globally for completeness.

The pixel representation and spherical harmonic representation of the map are
interchangeable. Sampling the pixels makes it easier to require positive
emission from the planet, as we simply place positive log-normal priors on the
brightness of each pixel. These have a mean value of the peak brightness of
the lightcurve divided by $\pi$, and a scale factor of 0.1 to set a very wide
prior spanning multiple orders in – see the documentation of PyMC3 for the
precise form of this prior (Salvatier et al., 2016) . This enforces positivity
at the pixel locations only, so sometimes leads to the fitting of small areas
of negative emission between pixels.

If we sampled the spherical harmonics directly instead, we could not impose
positivity as a prior, and would need to measure the positivity of the map at
each step while it is sampled, which we found to be computationally intensive.
We also found that sampling the pixel representation detected small-scale
features more effectively than the spherical harmonic representation. The
dense mass matrix approach available in PyMC3 (Salvatier et al., 2016) (first
implemented in Foreman-Mackey et al. (2020)) allows us to efficiently sample
the degeneracies in the posterior distribution of the brightness of the
pixels.

The number of mapped pixels is defined by the maximum spherical harmonic
order, increased by an oversampling factor that adds additional pixels to
better cover the planetary surface. Increasing the spherical harmonic order
increases the complexity of the features represented by the modelled light
curve, while increasing the oversampling factor improves the accuracy of the
measurement of map entropy in Section 2.7. Increasing both of these factors is
desirable but requires more computation time.

For our high-order maps, designed to have high spatial freedom, we found that
a spherical harmonic basis of $\ell_{\mathrm{max}}=4$ and an oversampling
factor of 3 (resulting in 62 pixels to sample) produced a good compromise
between spatial freedom and computational time. For our low-order maps,
designed to compare directly to the null hypothesis of an $n_{\mathrm{max}}=2$
Fourier series, we use a basis of $\ell_{\mathrm{max}}=2$ and an oversampling
factor of 3 (corresponding to 16 pixels). We define the pixels for a given map
as a vector of brightness values $\mathbf{p}$. The sampled pixel values
$\mathbf{p}$ define the vector $\mathbf{c}$ of the coefficients $c_{i}$ of the
spherical harmonic representation via a matrix $\mathbf{M}$:

$\mathbf{c}=\mathbf{M}\cdot\mathbf{p}.$ (2)

The matrix $\mathbf{M}$ is pre-calculated for our method by starry as it is
unique to each combination of spherical harmonic order $\ell_{\mathrm{max}}$
and degree of oversampling. The 2D eclipse map $Z(\theta,\phi)$ is then
constructed as a sum of the spherical harmonic basis maps $z_{i}(\theta,\phi)$
(where $\theta$ and $\phi$ are longitude and latitude), weighted by the
coefficients $c_{i}$ of the vector $\mathbf{c}$:

$Z(\theta,\phi)=\sum_{i}c_{i}z_{i}(\theta,\phi).$ (3)

The time-series emission $F(t)$ from the map defined by the coefficients
$c_{i}$ is then a sum of the time-series emission $f_{i}(t)$ from each
individual harmonic $z_{i}(\theta,\phi)$ (see Luger et al. (2019) for their
forms) weighted by its corresponding coefficient:

$F(t)=\sum_{i}c_{i}f_{i}(t).$ (4)

The modelled light curve $F(t)$ is therefore ultimately a function of the
vector of the pixel brightness values $\mathbf{p}$, which are the parameters
that we sample.

We fit a model $\mathcal{M}$ to the data points $D$ in the light curve using
PyMC3 (Salvatier et al., 2016). We sample the posterior distribution of the
brightness values of the pixels $\mathbf{p}$ with 250 samples in 8 chains,
thinned by a factor of 4. This ensured that the Gelman-Rubin statistic was
below 1.1 for our fitted maps. Ideally we would have used more samples, but
compromised with this smaller number due to the need to re-run the entire
sampling process on the order of 100 times to calculate the cross-validation
score. It would be possible to fit the orbital and instrumental parameters at
the same time as this map (Hammond et al., 2024), but in this study we assume
these parameters are known in advance.

###  2.5 Model comparison with AIC and BIC

Eclipse maps with more parameters (more surface pixels or a higher spherical
harmonic order) can fit smaller-scale features. This lets them fit an observed
time series more closely, but they may use this freedom to overfit to the
observational noise, mapping spurious small-scale features that do not really
exist and increasing the uncertainty of the map. On the other hand, maps with
too few parameters may not be able to fit the observed eclipse shape.
Comparing models by their likelihood will always prefer more complex models
that may overfit the data. We therefore need statistical metrics that penalise
using too many parameters to overfit the data.

The Bayesian Information Criterion (BIC) (Schwarz, 1978) and the Akaike
Information Criterion (AIC) (Akaike, 1981) are metrics that address this
issue, as they penalise the model likelihood by a function of the number of
model parameters. Previous studies have used the BIC to address issues with
eclipse map models overfitting small-scale features to eclipse shapes
(Rauscher et al., 2018; Mansfield et al., 2020; Challener & Rauscher, 2022).
The BIC is:

$\mathrm{BIC}=k\ln(N)-2\ln p,$ (5)

where $k$ is the number of model parameters, $N$ is the number of data points,
and $p$ is the model likelihood:

$\ln
p=-\frac{1}{2}\sum_{i}^{N}\left(\frac{\mathcal{M}_{i}-D_{i}}{\sigma}\right)^{2}-N\ln\sigma-\frac{N}{2}\ln(2\pi),$
(6)

for a sum over all $N$ data points $D_{i}$ with uncertainty $\sigma$, fitted
by a model $\mathcal{M}$ with values $\mathcal{M}_{i}$. A smaller BIC implies
a better model, either through a good fit to the data or a small number of
parameters. The BIC 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:

$\mathrm{AIC}=2k-2\ln p.$ (7)

Unlike the BIC, the AIC assumes that the true model is only approximated by
the set of tested models (Burnham & Anderson, 2004). This is essentially a
philosophical difference to the BIC, and there is not a clear answer about
which metric is more appropriate for our case. This is part of our motivation
for using the data-driven cross-validation approach instead.

###  2.6 Model comparison with cross-validation

We propose the use of k-fold cross-validation to compare models of time-series
thermal emission from exoplanet eclipses, as it penalises both underfitting
and overfitting. Cross-validation is a metric of the ability of a model to
predict data outside the dataset it is fitted to. In our case, the method
removes a section of observed data with size $k$, re-fits the model, and tests
how well the re-fitted model predicts the omitted data (Hastie et al., 2009).
We found that k-fold cross-validation was better suited to tuning a
regularisation parameter (details of which are below), compared to leave-one-
out cross-validation (Challener & Rauscher, 2023). We suggest this is because
omitting a k-fold (with a longer duration than a single point) produces a
stronger test of the predictive ability of a model due to the local
autocorrelation in the time-series data.

There is no objectively correct way to choose the size of a k-fold, although
$10\%$ of the size of a dataset is often used (Hastie et al., 2009). We chose
a k-fold size of $20\%$ of the size of the ingress of an eclipse (i.e. $10\%$
of the total dataset containing two-dimensional information) which we found to
be a good test of the predictive ability of a model. For the typical cadence
and eclipse duration of our datasets, this results in k-folds containing
roughly 10 data points each, in line with Arlot & Lerasle (2016).

This choice of the width of the k-folds will limit the size of the spatial
features that can be resolved with our mapping method. Each fold omits a
stripe approximately $36^{\circ}$ wide from either ingress or egress, and
tests the ability of the model to predict the brightness of this stripe. This
width is not a hard limit in itself, as the stripe is still measured across a
different angle by the other ingress or egress, and the model still has access
to the magnitude and gradient of the brightness either side of the omitted
stripe. However, it must places some upper limit on the precision of the map
and so this choice could be varied in future work.

Measuring the k-fold cross-validation score is a time-consuming process as the
eclipse map must be re-fitted every time a new k-fold is removed. We greatly
speed up this process by only testing k-folds in the ingress and egress of the
eclipses, as well as sections of the light curve on either side of the eclipse
with durations equal to the eclipse ingress and egress. This omits most of the
out-of-eclipse phase curve from the test, but we found that this makes no
practical difference. For typical datasets, the model only fails to fit the
phase cure if it has far too few degrees of freedom, which our small test set
easily identifies. For datasets with multiple eclipses, we apply these k-folds
periodically, removing the same section of each eclipse (or either side of
it). This is the same as stacking the eclipses (averaging them together) and
removing the k-fold from the averaged eclipse. We already effectively stack
the eclipses by fitting a periodic model that cannot vary from one eclipse to
another.

We split the tested part of the dataset (the eclipse ingress and egress, and
sections of the phase curve of equal duration either side of the eclipse) into
$K=20$ sets, each of length $N$ points. When leaving out the $k^{\mathrm{th}}$
fold, we label the dataset as $D_{-k}$, and the model as $\mathcal{M}^{-k}$.
The cross validation score for the $k^{\mathrm{th}}$ fold is:

$CV^{k}=\frac{1}{N}\sum^{N}_{i=0}\log\left(\frac{1}{S}\sum_{s=1}^{S}p(D_{i}|\mathcal{M}^{-k})\right),$
(8)

where $p(D_{i}|\mathcal{M}^{-k})$ is the likelihood of each point $D_{i}$ in
the k-fold, given the model $\mathcal{M}^{-k}$ fitted to the dataset $D_{-k}$.
For each data point, we sum over the $S$ samples drawn from the posterior
distribution as in Vehtari et al. (2017), and then sum over each of the $N$
data points in the $k^{\mathrm{th}}$ fold. $CV^{k}$ therefore tests the
ability of the model to predict the points in the $k^{\mathrm{th}}$ fold when
it is fitted without access to these points. In the context of eclipse
mapping, it tests the ability of a model to predict the brightness of a small
slice across the planetary disk when the contribution of that slice to the
eclipse shape is removed. When we use this method to optimise a fitted map, we
always plot the map fitted to the full dataset and never explicitly plot maps
fitted to datasets with data removed from k-folds.

We find the overall cross-validation score $CV_{\mathcal{M}}$ for a model
$\mathcal{M}$ by averaging over all $K$ folds:

$CV_{\mathcal{M}}=\frac{1}{K}\sum_{k=1}^{K}CV^{k}.$ (9)

The cross-validation score penalises underfitting (when the model cannot fit
the data) as the predictive likelihood $p(D_{i}|\mathcal{M}^{-k})$ is low for
a poorly fitting model. Importantly, it also penalises overfitting because a
model $\mathcal{M}^{-k}$ with too many degrees of freedom will produce a
posterior distribution in the $k^{\mathrm{th}}$ fold that is much wider than
the true spread of the omitted data $D_{k}$. This means that the average
predictive likelihood $p(D_{i}|\mathcal{M}^{-k})$ will be low when summed over
all $S$ samples, resulting in a low cross-validation score.

Conversely, a well-fitting model will predict the omitted datapoints
accurately and precisely as it is fitting the underlying physical system
(instead of overfitting to noise), so the data points in one particular k-fold
are implied by those in the others. This will result in a model
$\mathcal{M}^{-k}$ with a posterior distribution with a similar mean and
standard deviation to the real omitted data $D_{k}$, producing a high
predictive likelihood $p(D_{i}|\mathcal{M}^{-k})$ and a high cross-validation
score.

To compare the cross-validation score for two competing models
$\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, we take the difference between them:

$\Delta
CV_{\mathcal{M}_{1}\mathcal{M}_{2}}=CV_{\mathcal{M}_{1}}-CV_{\mathcal{M}_{2}}.$
(10)

The standard deviation of this estimate is calculated from the variance $V$ of
the difference in the scores in each k-fold (Chen & Yang, 2021; Welbanks et
al., 2023):

$\mathrm{SD}(\Delta
CV_{\mathcal{M}_{1}\mathcal{M}_{2}})=\sqrt{V_{k=1}^{K}(CV^{k}_{\mathcal{M}_{1}}-CV^{k}_{\mathcal{M}_{2}})}$
(11)

where $CV^{k}_{\mathcal{M}_{1}}$ and $CV^{k}_{\mathcal{M}_{1}}$ are calculated
from Equation 8 for each individual k-fold. The standard error is then
calculated from the standard deviation (Chen & Yang, 2021):

$\mathrm{SE}(\Delta
CV_{\mathcal{M}_{1}\mathcal{M}_{2}})=\frac{\mathrm{SD}(\Delta
CV_{\mathcal{M}_{1}\mathcal{M}_{2}})}{\sqrt{K}}$ (12)

Due to the assumptions made in deriving this standard error estimate, we
caution against using it as an “multiple-sigma” preference for one model over
another, and restrict its use to identifying a simple preference one way or
the other. The standard error of the difference in cross-validation score
between two models tells us if the difference is significant or not. If two
models have consistent cross-validation scores within this standard error
estimate, we assess that they are as good as each other (and therefore we
prefer the simpler model).

The exception to this is when we later tune a regularisation parameter to find
the best cross-validation score; in that case, we select the model with the
highest cross-validation score as the best-performing, even though other
nearby values of the regularisation parameter will necessarily be within one
standard error. We chose not to use the the “one standard error” rule (Chen &
Yang, 2021) in this case, as we have no a priori tendency to prefer a higher
or lower regularisation parameter.

###  2.7 Optimising the spatial scale of fitted maps

The eclipse maps fitted with the method in Section 2.4 have a spatial scale
determined by their number of pixels (corrected for oversampling). Maps with a
small number of pixels (or low-order spherical harmonics) can only fit very
large-scale features with similarities to the shapes of the low-order
spherical harmonics. Figure 1 shows how a low-order spherical harmonic
representation may not be able to fit the true shape of a map. Maps with large
numbers of pixels (or high-order spherical harmonics) can fit spatial features
more accurately, as shown in Figure 1. However, they may also produce spurious
small-scale features due to the noise and degeneracy in the observed light
curve. To address this issue, we use a large number of pixels to allow freedom
in fitting small-scale features, but apply the method of Horne (1985) (applied
to brown dwarfs in Chen et al. (2024)) to select maps with an appropriate
spatial scale given the precision of the data.

Figure 2: Two simulated datasets with clear eclipse mapping signals. First
row: a simulated observation (with precision $\sigma=150$ ppm) of time-series
thermal emission from the simulation of WASP-43b in the THOR GCM. The left-
hand panel shows the dataset fitted by an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$. The right-hand panel shows the residual
eclipse mapping signal, which is the data (binned every 10 points) subtracted
by a fit with a Fourier series model $\mathcal{M}(n_{\mathrm{max}}=2)$ that
assumes a uniform planetary disk. There is a clear residual signal (dashed red
line) that the data (black points) resolves, which the eclipse map model (blue
region) fits better than the Fourier series model (dashed black line). The
blue region showing the map model contains two shaded regions showing the
first and second quantiles, containing 68.27% and 95.45% of the posterior
distribution. Second row: a simulated observation (with precision $\sigma=250$
ppm) of time-series thermal emission from the simulation of WASP-43b in the
THOR GCM. As in the first row, there is a clear eclipse mapping signal which
the eclipse mapping model fits better than the Fourier series model, although
with slightly less certainty due to the lower-precision data.

We measure the information content of a map from the spatial entropy of its
pixels (Gull & Daniell, 1978). We seek to match the information content of the
map to the information content of the dataset. A map with more information
than the dataset will overfit, and a map with less information than the
dataset will underfit. We follow Horne (1985) and consider the entropy $S$ of
the pixels defining a map $I$ composed of pixels denoted by $j$:

$S=-\sum_{j}I(j)\ln\left[\frac{I(j)}{D(j)}\right]\text{. }$ (13)

where $D$ is the “default image” defined by Horne (1985). This default image
depends on a weighting function $w(k,j)$ which defines the relation of all
pixels $j$ to an individual pixel $k$:

$D(k)=\frac{\sum_{j}w(k,j)\ I(j)}{\sum_{j}w(k,j)}.$ (14)

We use a weighting function $w(k,j)=1$ so that the entropy tracks deviations
from a uniform map. This matches our null hypothesis, which is the eclipse
shape of a uniformly bright planetary disk as defined in Section 2.3.

We include the information content of the map in the fitting process by
penalising the model likelihood $p$ by a factor $2\alpha S$, following Chen et
al. (2024):

$\hat{p}=p-2\alpha S,$ (15)

where $\alpha$ is a regularisation parameter. A large value of $\alpha$
promotes maps with low entropy and large-scale features, which may underfit
the data. A small value of $\alpha$ allows maps with high entropy and small-
scale features, which may overfit the data. We aim to find an optimal value of
$\alpha$, which achieves the best k-fold cross-validation score. We therefore
vary $\alpha$ until we find the value that results in the best cross-
validation score. We expect that this will produce a map with an appropriate
information content given the precision of the data. This method is analogous
to fitting a spline to a one- or two-dimensional dataset, where there is a
single regularisation parameter that gives an optimal cross-validation score
(Wahba, 1990).

##  3 Simulated Results

In this section, we apply the methods from Section 2 to simulated data. We
first demonstrate the process of identifying eclipse mapping signals in
simulated data using cross-validation. We then show how fitting an eclipse map
with spherical harmonics alone can lead to under- or over-fitting. To remedy
this, we show how we can optimise the information content (smoothness) of a
high-order spherical harmonic map using cross-validation, and derive an
optimal map for this simulated dataset. Note that in Sections 3.1 and 3.2 we
do not apply a penalty to the likelihood based on the entropy of the fitted
map (effectively, we keep $\alpha=0$ until Section 3.3).

###  3.1 Identifying eclipse mapping signals in simulated data

Figure 3: Two simulated datasets that are too inaccurate or too imprecise to
fit an eclipse map. First row: a simulated observation (with precision
$\sigma=150$ ppm) of time-series thermal emission from the simulation of
WASP-43b in the THOR GCM, where the map is fitted with an incorrect value for
the transit time of $t_{0}^{\prime}=t_{0}+10\mathrm{s}$. The layout of the
panels is the same as in Figure 2. The incorrect astrophysical model means
that the longitudinal structure imposed by the phase curve produces an
incorrect fit to the data in the eclipse, shown by the mismatch between the
black data points and the blue fitted model. Second row: a simulated
observation (with precision $\sigma=2000$ ppm) of time-series thermal emission
from the simulation of WASP-43b in the THOR GCM. The low precision of the data
means that it does not resolve the eclipse mapping signal, so the
$\ell_{\mathrm{max}}=2$ model does not fit this signal any better than the
$n_{\mathrm{max}}=2$ Fourier series model (the dashed black zero line). The
red dashed lines showing the true eclipse mapping signals are identical to
those in Figure 2.

As described in Section 2.6, it is possible to fit spurious eclipse maps to
datasets that are not accurate enough to support them. In this section, we
therefore show how the cross-validation method described in Section 2.6 can
quantitatively identify the presence or absence of an eclipse mapping signal
in a dataset.

Figure 4: An eclipse map model $\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the 150 ppm simulated data in Figure 2. The first panel shows a map defined by the median values of the sampled distribution at every location. The second and third panels show the fitted map as a function of latitude and longitude through the substellar point, with two shaded regions showing the first and second quantiles, containing 68.27% and 95.45% of the posterior distribution. The limited freedom of the $\ell_{\mathrm{max}}=2$ map means that it cannot exactly fit the true map from the GCM simulation (the dashed black lines, shown in full in Figure 1) which has more complex structures composed of higher-order harmonics. This results in the brightest point on the equator being measured at $(24^{+0}_{-0})^{\circ}$ east, compared to the true value of $40^{\circ}$ east shown in Figure 1. The lack of uncertainty on this measurement is due to the strong constraint from the full phase curve on the longitudinal structure of the fitted low-order spherical harmonics. Dataset | $\Delta\chi^{2}$ |  $\Delta$BIC |  $\Delta$AIC |  $\Delta$CV
---|---|---|---|---
WASP-43b GCM Simulation ($\sigma=150$) | +126.04 | +81.77 | +116.04 | $+0.148\pm 0.016$
WASP-43b GCM Simulation ($\sigma=250$) | +73.37 | +29.1 | +63.37 | $+0.076\pm 0.007$
WASP-43b GCM Simulation ($\sigma=150$, $t_{0}^{\prime}=t_{0}+10\mathrm{s}$) | +54.79 | +10.52 | +44.79 | $-0.002\pm 0.013$
WASP-43b GCM Simulation ($\sigma=2000$) | -0.31 | -44.58 | -10.31 | $-0.004\pm 0.001$
Table 1: Statistical metrics for the eclipse map models
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the datasets in Figures 2 and
3, relative to the metrics for the Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$ (the null hypothesis). The $\chi^{2}$ value
is the lowest $\chi^{2}$ value in the posterior of the fitted model, which is
also used to derive the BIC and the AIC. A positive value for a metric means
that the eclipse mapping model is preferred over the Fourier series model.

Figure 5: Eclipse map models $\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to
the simulated datasets in Figure 3, with low accuracy and low precision
respectively. First row: the map fitted to the inaccurate simulated data in
the first row of Figure 3 with an incorrect transit time
$t_{0}^{\prime}=t_{0}+10\mathrm{s}$. The layout of the panels is the same as
Figure 4. The systematic error on the eclipse timing means that the fitted map
does not match the true map (the dashed black lines, shown in full in Figure
1) but has spuriously high confidence. Second row: the map fitted to the
imprecise data in the second row of Figure 3. The fitted map has very high
uncertainty but does include the true map inside its posterior distribution as
it has no systematic errors like the error in eclipse time present in the map
in the first row. Neither dataset is suitable for fitting an eclipse map,
shown by the cross-validation score in Table 1 and Figure 6.

Figures 2 and 3 show light curves produced from the GCM simulation of WASP-43b
as described in Section 2.1. Their left-hand panels show the time-series
thermal emission. Their right-hand panels show the same data, binned every 10
points, and subtracted by a model of a planet with a uniform dayside emission.
This leaves a residual signal which corresponds to the deviation of the
eclipse shape from the eclipse shape of a uniform disk — the “eclipse mapping
signal”, as discussed in Hammond et al. (2024). A consistent deviation in the
data points from the zero line implies an eclipse mapping signal. The dashed
red lines show the true residual signal in eclipse before the Gaussian noise
is added. A positive residual signal implies that the strip across the
planetary disk covered by the stellar edge at that moment is brighter than a
uniform disk producing the same total eclipse depth and vice versa.

The first row of Figure 2 shows the simulated data for $\sigma=$ 150 ppm. This
represents a best-case scenario for observations with JWST, corresponding to
an average of five eclipses of WASP-43b with MIRI/LRS (Bell et al., 2024).
There is a clear eclipse mapping signal, shown by the large residual in the
data in the right-hand panel, fitted well by the $\ell_{\mathrm{max}}=2$
eclipse map model.

Table 1 shows that the eclipse mapping model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ achieves a significantly better
$\chi^{2}$ value for the $\sigma=150$ ppm data than the Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$, due to its better fit to the eclipse shape
in the first row of Figure 2. This does not necessarily demonstrate the
presence of two-dimensional information in the eclipse shape, as the eclipse
mapping model has more degrees of freedom so may achieve a better $\chi^{2}$
score by overfitting to noise in the eclipse. Both models fit the out-of-
eclipse phase curve as well as each other, as the $\ell_{\mathrm{max}}=2$
spherical harmonics produce the same out-of-eclipse phase curve components as
the $n_{\mathrm{max}}=2$ Fourier components, apart from a small correction for
inclination.

Figure 4 shows the eclipse map model $\mathcal{M}(\ell_{\mathrm{max}}=2)$
fitted to the $\sigma=$ 150 ppm dataset in Figure 2. The first panel shows the
median of the posterior distribution of fitted maps (note that the two-
dimensional structure is only constrained on the dayside). We plot the median
at every point because we sample the pixels composing the map; this plot
therefore shows the median value of the sampled posterior distribution. The
next two panels of Figure 4 show slices of the posterior distribution of the
fitted maps on the dayside through the substellar point, as functions of
longitude and latitude. The true map is composed of many high-order spherical
harmonic components beyond $\ell=2$, so this fitted $\ell_{\mathrm{max}}=2$
map cannot match the shape of the true map very closely. It does, at least,
fit the residual signal accurately, so we suggest that this map is
statistically supported, if not especially accurate.

The second row of Figure 2 shows the simulated data for $\sigma=$ 250 ppm.
This represents a typical scenario for observations with JWST, corresponding
to an average of two eclipses of WASP-43b with MIRI/LRS (Bell et al., 2024).
The large residual shows an eclipse mapping signal, although resolved less
clearly than in the first row of Figure 2. Table 1 again shows that the
eclipse mapping model $\mathcal{M}(\ell_{\mathrm{max}}=2)$ achieves a much
better $\chi^{2}$ value than the Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$ for this dataset.

The first row of Figure 3 shows the simulated data for $\sigma=$ 150 ppm, with
the Fourier series and eclipse map models fitted using an incorrect value for
the transit time $t_{0}^{\prime}=t_{0}+10\mathrm{s}$. As we assume a circular
orbit, this produces the same offset of $+10\mathrm{s}$ in the assumed eclipse
time. This simple change produces a systematic error in the astrophysical
model, which means that the eclipse mapping signal and the eclipse map are
fitted poorly regardless of the precision of the data. An error in assumed
eclipse time could also arise in other ways, such as a slightly eccentric
orbit producing an offset in eclipse timing relative to a measured transit
time. Fitting the model with this delayed eclipse timing results in an eclipse
ingress that is too high, and an eclipse egress that is too low. The model
wrongly adjusts the eclipse map to compensate for this error.

This example is intended to represent any source of systematic error in the
fit of the model to the data, such as instrumental systematic effects or other
incorrect orbital parameters. Both of these can produce significant deviations
in the modelled eclipse shape (de Wit et al., 2012; Hammond et al., 2024). The
modelled eclipse map residual in the first row of Figure 3 fails to match the
eclipse data accurately due to the systematic error in the astrophysical model
caused by the incorrect eclipse timing. This is a distinctive signature of a
map fitted with a systematic error — the longitudinal structure imposed by the
phase curve produces a residual shape that is not consistent with the observed
residual; the latitudinal structure of the map is then wrongly adjusted to try
to compensate for this error. The first row of Figure 5 shows the eclipse map
fitted to this dataset, which fails to match the true map, especially in its
latitudinal structure.

Table 1 shows that the eclipse mapping model still achieves a better
$\chi^{2}$ value than the Fourier series model for this dataset. This is
despite the highly inaccurate map and despite the fact that its posterior
distribution is highly statistically different from the true residual signal
(compare the blue posterior distribution to the red dashed line in the first
row of Figure 3). The better $\chi^{2}$ value is due to the fact that the
eclipse map model has increased freedom to fit the eclipse shape, even though
this does not correctly correspond to a mapping signal.

Figure 6: Left panel: The relative cross-validation scores (Equation 10) for
eclipse map models $\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the
simulated datasets in Figures 2 and 3, compared to the scores for Fourier
series models $\mathcal{M}(n_{\mathrm{max}}=2)$. The scores are normalised by
their standard errors (Equation 12). The eclipse map models
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the higher precision data in
Figure 2 have significantly better cross-validation scores than the Fourier
series models, indicating the presence of eclipse mapping signals. The eclipse
map models $\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the inaccurate and
imprecise data in Figure 3 achieve cross-validation scores that are lower than
(or consistent with) the Fourier series model, indicating the absence of
eclipse mapping signals. Right panel: the difference in cross-validation score
for eclipse map models $\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$ fitted to
the simulated datasets in Figure 2, with their likelihood penalised by the
regularisation parameter $\alpha$ as described in Section 2.7. This identifies
an optimal information content (spatial scale) for the map in each case. The
cross-validation scores are plotted relative to the optimal cross-validation
score for each dataset, so the maximum value of each curve is zero.

The second row of Figure 3 shows the simulated data for $\sigma=$ 2000 ppm.
This represents a light curve that is too imprecise for eclipse mapping. The
posterior distribution of the eclipse mapping model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ does not match the true residual shape
(the dashed red line). The second row of Figure 5 shows the eclipse map fitted
to this dataset, which is very uncertain due to the low precision of the data.
Table 1 shows that the mapping model achieves almost the same $\chi^{2}$ value
as the $\mathcal{M}(n_{\mathrm{max}}=2)$ model (actually a slightly worse
value). Figure 3 shows how the two models fit the data essentially as well as
each other as the precision of the data is too low to reveal the effect of the
two-dimensional structure on the eclipse shape.

These simulated datasets show different pitfalls that could be encountered
when fitting eclipse maps. We need statistical metrics that will strongly
penalise mapping models that use too many parameters to fit the eclipse shape
of inaccurate or imprecise data.

Figure 7: A high-order eclipse map model $\mathcal{M}(\ell_{\mathrm{max}}=4)$
fitted to the 150 ppm simulated data in Figure 2. The layout of the panels is
the same as Figure 4. The higher degree of spatial freedom means that it is
possible in theory for the fitted map to closely match the true map. However,
in practice the increased spatial freedom allows too large a range of
solutions, making the uncertainty on the map very large. This makes the
measurement of the brightest point on the equator very imprecise at
$(56^{+65}_{-58})^{\circ}$ east, compared to the true value of $40^{\circ}$
east shown in Figure 1.

Figure 8: Three eclipse map models $\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$
fitted to the simulated $\sigma=150$ ppm data in Figure 2 with variable values
of $\alpha$. First row: an eclipse map model
$\mathcal{M}(\alpha=2154,\ell_{\mathrm{max}}=4)$. This low value of $\alpha$
weakly penalises the likelihood of maps with a high information content
(Equation 15), producing an uncertain map with a low cross-validation score in
Figure 6. The fitted equatorial brightness maximum is
$(49^{+11}_{-36})^{\circ}$ east, compared to the true value of $40^{\circ}$
shown in Figure 1. Second row: a model
$\mathcal{M}(\alpha=21544,\ell_{\mathrm{max}}=4)$ fitted to the same data.
This allows an appropriate information content in the map given the precision
of the input data, giving the optimal cross-validation score in Figure 6. This
results in an accurate and precise measurement of the equatorial brightness
maximum at $(45^{+7}_{-7})^{\circ}$ east, compared to the true value of
$40^{\circ}$. Third row: a model
$\mathcal{M}(\alpha=215443,\ell_{\mathrm{max}}=4)$ fitted to the same data.
This high value of $\alpha$ strongly penalises the likelihood of maps with
high information content, which allows too little information in the map and
achieves a poor cross-validation score in Figure 6. It still achieves an
accurate but overly precise measurement of the position of maximum brightness
on the equator, of $(38^{+7}_{-4})^{\circ}$ east, compared to the true value
of $40^{\circ}$.

###  3.2 Model comparison in simulated data

We need a model comparison metric to identify the presence of an eclipse
mapping signal. We define that the metric must quantitatively prefer the
simplest eclipse mapping model $\mathcal{M}(\ell_{\mathrm{max}}=2)$ over the
simplest Fourier series model $\mathcal{M}(n_{\mathrm{max}}=2)$. We want the
metric to identify an eclipse mapping signal in the simulated datasets in
Figure 2, but not in the simulated datasets in Figure 3. Section 2.5 describes
the BIC, which has been used previously for eclipse mapping model comparison
(Challener & Rauscher, 2022), and the AIC, which is an alternative metric
(Burnham & Anderson, 2004). These metrics penalise the fit likelihood by a
function of the number of model parameters to penalise both overfitting and
underfitting.

Table 1 lists the difference in AIC and BIC between an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ and a Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$ fitted to the four simulated datasets. A
positive $\Delta$AIC or $\Delta$BIC means that the metric prefers the eclipse
mapping model. This table shows that both the AIC and the BIC correctly
identify an eclipse mapping signal in the 150 ppm and 250 ppm simulated data
in Figure 2. This is because the advantage in $\chi^{2}$ for the eclipse
mapping model is so high for each dataset that it outweighs the penalty for
its increased number of parameters in each case.

Both the AIC and BIC also correctly identify the lack of an eclipse mapping
signal in the imprecise 2000 ppm data in Figure 3, which is fitted just as
well by the Fourier series model. In this case, the $\chi^{2}$ value is
essentially the same for the eclipse mapping model and the Fourier series
model (the best value is actually slightly better for the Fourier series
model, but the difference is not significant). The penalty for the increased
number of parameters in the eclipse mapping model then gives it much worse BIC
and AIC scores.

The picture is more complex for the 150 ppm data fitted with an incorrect
value of $t_{0}$ in Figure 3. The eclipse map model does fit the observed data
more closely than the Fourier series model, because there are elements of the
eclipse shape that it models correctly. However, it also incorrectly adjusts
its latitudinal structure to compensate for the error in the eclipse timing,
producing a highly incorrect eclipse map with spuriously high confidence. We
suggest that a model comparison metric should not prefer the eclipse mapping
model in this case, but both the BIC and the AIC strongly prefer the eclipse
mapping model. This is because they cannot distinguish when a model achieves a
good fit by an incorrect route. The penalty for the increased number of
parameters is then not sufficient to outweigh the spuriously improved
$\chi^{2}$ score. We suggest that this is a strength of the cross-validation
score (see below), which identifies the fact that the eclipse map model is not
fitting a robust and consistent signal.

The BIC and AIC also produce different answers for the relative weighting of
models, and there is no clear answer about which metric is more appropriate a
priori as discussed in Section 2.5. Despite these issues, the AIC and BIC are
useful improvements on the $\chi^{2}$ metric as they do attempt to penalise
for overfitting. However, crucially for our purposes, the AIC and BIC cannot
be used for the regularisation approach described in Section 2.7. We therefore
turn to model comparison by cross-validation. Cross-validation lets us
calculate the predictive power of the model, which applies a penalty for
overfitting without needing to specify a number of model parameters. We
calculate the k-fold cross-validation score as described in Section 2.6 for
the simulated datasets in Figures 2 and 3, for the simplest eclipse mapping
model $\mathcal{M}(\ell_{\mathrm{max}}=2)$ compared to the simplest Fourier
series model $\mathcal{M}(n_{\mathrm{max}}=2)$. For this test, we do not
penalise the model likelihood by the map entropy as described in Section 2.7
(i.e. $\alpha=0$ here).

The left-hand panel of Figure 6 compares the cross-validation score difference
between the eclipse mapping model and the Fourier series model for these three
datasets. We normalise each relative cross-validation score by its standard
error, so that a difference greater than 1 is significant. This shows that the
eclipse map models achieve significantly better cross-validation scores for
the datasets with precisions of $\sigma=150$ ppm and $\sigma=250$ ppm (Figure
2). As expected, this shows that the fitted eclipse map models are much better
at predicting out-of-sample data than the Fourier series models, so they are
not achieving better $\chi^{2}$ values by overfitting to noise.

Figure 6 hows that an eclipse map model fitted to the dataset in Figure 3 with
an incorrect value of $t_{0}^{\prime}=t_{0}+10\mathrm{s}$ achieves a lower
cross-validation score than a Fourier series model. The standard error of the
difference between the two scores shows that their scores are statistically
consistent with each other, so the eclipse mapping model is no better or worse
at predicting out-of-sample data than the Fourier series model. This means
that there is no statistically robust eclipse mapping signal. This is an
improvement on the BIC and the AIC, which preferred the eclipse mapping model
in this case in Table 1.

The eclipse mapping model also achieves a significantly worse cross-validation
score than the Fourier series model for the dataset with precision
$\sigma=2000$ ppm in the second row of Figure 3. This is because the mapping
model predicts a wide range of implausible solutions that can be very
different to the omitted data, while the Fourier series model is restricted to
a small range of solutions matching the shape of a uniform eclipse.

We conclude that, for these simulated datasets, the cross-validation score
correctly identifies if the fitted model is tracking the true eclipse mapping
signal, instead of overfitting to any imprecise or inaccurate data points. The
BIC and AIC perform less well at identifying the presence or absence of
eclipse mapping signals and – unlike the cross-validation score – cannot be
used with a regularisation parameter to optimise information content.

###  3.3 Optimising the information content of fitted maps

Having used a simple eclipse mapping model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ to establish the presence of an eclipse
mapping signal in the dataset with $\sigma=$ 150 ppm in Figure 2, we can now
fit a map with more complex structure. Figure 7 shows an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=4)$ fitted to this dataset using spherical
harmonics up to order $\ell_{\mathrm{max}}=4$. This high-order map is able to
fit the true map in theory, but in practice it overfits a wide range of
solutions, resulting in an impractically large uncertainty.

So, while there is a clear eclipse mapping signal, both the low-order map
(Figure 4) and the high-order map (Figure 7) fail to match the true map. This
issue is even worse for real data as we do not know the true map in advance,
so do not know how many spherical harmonic orders are enough to fit
accurately. This is the motivation behind the method to fit a high-order map
with optimised information content described in Section 2.7. The right-hand
panel of Figure 6 shows how the cross-validation score for an eclipse map
model $\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$ varies as the regularisation
parameter $\alpha$ in Equation 15 is varied. For both datasets, there is a
significant difference between the cross-validation scores of the lowest value
of $\alpha$ (which tends towards the score for $\alpha=0$) and the optimal
value of $\alpha$.

When $\alpha$ is very small, the fitted map is not strongly penalised for a
high entropy (high spatial variability), so the very imprecise map in the
first row of Figure 8 produces the best likelihood (see Equation 15), and the
data is overfitted so the cross-validation score is low. When $\alpha$ is very
large, spatial variation is strongly penalised so the fit in the third row of
Figure 8 produces the best likelihood, and the data is underfitted so the
cross-validation score is low. For the optimal value of $\alpha$, the second
row of Figure 8 produces the best likelihood, achieving the best cross-
validation score as it has an appropriate information content. The fitted map
matches the data and also allows an appropriate amount of variation as well as
appropriate average gradients and feature sizes for the true map. Therefore,
cross-validation gives us a data-driven method to find the optimal value of
$\alpha$, and the corresponding optimal information content for the map, with
no prior knowledge of the true map. We suggest that as a shortcut, choosing
and applying a “reasonable” value of $\alpha$ is a quick way to produce a
better result than a $\mathcal{M}(n_{\mathrm{max}}=4)$ model, especially if
the reasonable value is derived from a previous similar observation. However,
deriving the optimal value of $\alpha$ manually using cross-validation for
each dataset will be best, although time-consuming.

Figure 9: The thermal emission from 5 to 10.5 $\mu$m from the hot Jupiter
WASP-43b observed by JWST MIRI/LRS (Bell et al., 2024). The layout of panels
is the same as Figure 2. The black points show the first eclipse and the grey
points show the second eclipse, binned every 10 points. The eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ fits the eclipse shape much better than
the Fourier series model $\mathcal{M}(n_{\mathrm{max}}=2)$. This is confirmed
by the higher cross-validation score for the eclipse map model for this
dataset shown in Figure 12 and listed in Table 2. The data points in the
right-hand panel appear slightly different to the data points in Hammond et
al. (2024) because here we bin them with a slightly longer time cadence.

Figure 10: The thermal emission from 0.8 to 2.8 $\mu$m from the hot Jupiter
WASP-18b observed by JWST NIRISS/SOSS (Coulombe et al., 2023). The data in the
right-hand panel is binned every 10 points. There is some residual deviation
in the data from the Fourier series model $\mathcal{M}(n_{\mathrm{max}}=2)$,
but it is not fitted well by the eclipse mapping model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$. Figure 12 suggests there is no
significant eclipse mapping signal for this dataset, as the eclipse mapping
model is no better than the Fourier series model at predicting out-of-sample
sections of the eclipse.

Figure 11: The observed thermal emission at 8 $\mu$m from the hot Jupiter HD 189733b observed by the Spitzer Space Telescope (Majeau et al., 2012). The data in the right-hand panel is binned every 5 points. There is some residual deviation from the Fourier series model $\mathcal{M}(n_{\mathrm{max}}=2)$, but it is not fitted well by the eclipse mapping model $\mathcal{M}(\ell_{\mathrm{max}}=2)$. Figure 12 suggests there is no significant eclipse mapping signal for this dataset, as the eclipse mapping model is no better than the Fourier series model at predicting out-of-sample sections of the eclipse. Dataset | $\Delta\chi^{2}$ |  $\Delta$BIC |  $\Delta$AIC |  $\Delta$CV
---|---|---|---|---
WASP-43b (JWST MIRI/LRS) | +103.4 | +58.2 | +93.4 | $+0.060\pm 0.008$
WASP-18b (JWST NIRISS SOSS) | -18.2 | -57.69 | -28.2 | $-0.054\pm 0.006$
HD 189733b (Spitzer) | +0.92 | -32.98 | -9.08 | $-0.067\pm 0.008$
Table 2: Statistical metrics comparing eclipse mapping models
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ to Fourier series models
$\mathcal{M}(n_{\mathrm{max}}=2)$, for the observational datasets in Figures
9, 11, and 11. We propose that $\Delta$CV is the best metric for identifying
an eclipse mapping signal, as discussed in Section 3. Only the JWST MIRI/LRS
WASP-43b dataset in Figure 9 has statistical evidence for an eclipse mapping
signal by this metric. Note that the BIC and AIC values for the WASP-43b
dataset are slightly different to the equivalent models fitted in Hammond et
al. (2024) because here we do not discard the data in the transit, and do not
include an additional parameter to fit the baseline of the eclipses. The
WASP-18b dataset prefers the Fourier series model (the null hypothesis, with
no mapping information) by all metrics. The HD 189733b dataset prefers the
Fourier series model by all metrics except $\Delta\chi^{2}$, but we suggest
this is due to overfitting. Figure 12: Left panel: the difference in cross-
validation score (Equation 10) for each dataset for the eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ compared to the Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$. Each score is normalised by its standard
error (Equation 12). Of the three observational datasets, only the JWST
MIRI/LRS observations of WASP-43b in Figure 9 achieves a better cross-
validation score with the eclipse map model. Right panel: the cross-validation
score of the eclipse map model $\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$
fitted to the JWST MIRI/LRS observations of WASP-43b in Figure 9, with its
likelihood penalised by the entropy scaled by $\alpha$ as described in Section
2.7. The best cross-validation score is achieved by setting $\alpha=4641.6$,
which produces the map in Figure 14.

##  4 Observational Results

In this section, we take the method demonstrated on simulated data in Section
3, and apply it to three observational datasets. These datasets are the JWST
MIRI/LRS observations of WASP-43b in Figure 9 (Bell et al., 2024), the JWST
NIRISS/SOSS observations of WASP-18b in Figure 11 (Coulombe et al., 2023), and
the Spitzer Space Telescope observations of HD 189733b in Figure 11 (Majeau et
al., 2012). The left-hand panel of these plots shows the raw observations,
fitted by an eclipse map model $\mathcal{M}(\ell_{\mathrm{max}}=2)$. The
right-hand panel of these plots shows the residual signal in ingress and
egress.

###  4.1 Identifying eclipse mapping signals in observed data

Figure 13: An eclipse map fitted to the data in Figure 9, using a model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$. Table 2 and Figure 12 show how this
model is strongly preferred over the Fourier series model (the null
hypothesis) by the cross-validation metric. Hammond et al. (2024) discussed
how the structure of this $\ell_{\mathrm{max}}=2$ map is limited by its
limited spatial freedom. Figure 14: An eclipse map fitted to the data in
Figure 9 using a model $\mathcal{M}(\alpha=4641.6,\ell_{\mathrm{max}}=4)$.
This value of $\alpha$ gives the optimal cross-validation score in Figure 12.
We suggest that this value imposes an information content which reflects a
combination of the spatial scale of the true map and the level of mapping
precision achievable given the precision of the data itself.

We apply our statistical method defined in Section 2 to determine if there is
quantitative evidence for an eclipse map in each dataset. The left-hand panel
of Figure 12 shows the difference in cross-validation score for an eclipse map
model $\mathcal{M}(\ell_{\mathrm{max}}=2)$ compared to a Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$, for each of the observational datasets in
Figures 9, 11, and 11. We do not apply a penalty to the likelihood of a fitted
map based on its entropy for these $\mathcal{M}(\ell_{\mathrm{max}}=2)$
models; we only later consider a non-zero $\alpha$ once we have established
the presence of a mapping signal with this initial test. The WASP-43b dataset
has a positive $\Delta$CV value, showing significant evidence for an eclipse
mapping signal, matching the clear fit to the large residual in the data in
Figure 9. Figure 13 shows an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to this dataset. The posterior
distribution of this fitted map is slightly different to the
$\ell_{\mathrm{max}}=2$ map in Hammond et al. (2024) because the prior placed
on the brightness of each pixel is wider in Figure 13.

On the other hand, the WASP-18b dataset has a negative $\Delta$CV value in
Figure 12, suggesting that there is no statistical evidence for an eclipse
mapping signal by this metric. The WASP-18b dataset has a comparable signal-
to-noise ratio to the WASP-43b dataset, so we would expect that it would be
possible to resolve a mapping signal. The mismatch between the fitted and
observed residual signal in Figure 11 implies that the longitudinal structure
imposed by the out-of-eclipse phase curve is inconsistent with the structure
implied by the eclipse shape. This may indicate an unresolved systematic error
in the astrophysical model or the instrumental model, like the inaccurate
simulated data in the top row of Figure 3. It is also possible that the day-
side has a complex structure that produces a smoothly varying phase curve but
a largely flat structure that requires higher-order harmonics than
$\ell_{\mathrm{max}}=2$ to resolve. The eclipse mapping model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ (shown in Supplementary Figure 18)
actually achieves a worse $\chi^{2}$ score than the Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$ in Table 2 for this dataset.

The HD 189733b dataset also has a negative $\Delta$CV values in Figure 12,
again suggesting that there is no statistical evidence for an eclipse mapping
signal by this metric. This conclusion matches the visual evidence in Figures
11 and 11, where there is no clear fit to a residual signal. We suggest that
the precision of the HD 189733b dataset may be too low to identify an eclipse
mapping signal, like the imprecise simulated dataset in Figure 3. Table 2
shows that the eclipse mapping model $\mathcal{M}(\ell_{\mathrm{max}}=2)$
(shown in Supplementary Figure 19) achieves a slightly better $\chi^{2}$ score
for the HD 189733b data than the Fourier series model
$\mathcal{M}(n_{\mathrm{max}}=2)$. We suggest this is due to overfitting to
noise in the eclipse shape, rather than matching an underlying signal. This
overfitting is penalised by the cross-validation score, which prefers the
simpler Fourier series model. This suggests that the two-dimensional structure
implied by this map is due to a combination of the out-of-eclipse phase curve
providing longitudinal structure and limits on the latitudinal structure due
to the priors (including the need for positive brightness) placed on the
sampled parameters.

The signs of the relative cross-validation scores in Table 2 match the signs
of the BIC and AIC scores and are consistent with the presence or absence of
visual signals in Figures 9, 11, and 11. We conclude that the $\Delta$CV score
provides strong statistical evidence for an eclipse mapping signal in the
WASP-43b dataset only. We now proceed to the second part of our method for
this dataset, fitting an $\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$ model
with its information content optimised by the method described in Section 2.7.

###  4.2 An optimised map fitted to observations of WASP-43b

As discussed in Section 3, spherical harmonics up to $\ell_{\mathrm{max}}=2$
have limited spatial freedom and may not be able to fit the true map for the
observed WASP-43b dataset. We therefore fit a new model
$\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$ with spherical harmonics up to
$\ell_{\mathrm{max}}=4$, with its information content determined by a
regularisation parameter $\alpha$ (Equation 15). The right-hand panel of
Figure 12 shows how the cross-validation score for this eclipse map varies as
a function of $\alpha$. There is a peak at $\alpha=4641.6$, corresponding to
the optimal amount of map information that allows an accurate fit to the data
without overfitting to the noise.

Figure 14 shows the resulting $\ell_{\mathrm{max}}=4$ map with
$\alpha=4641.6$. The map has a flatter emission structure near the substellar
point in longitude than the $\ell_{\mathrm{max}}=2$ map. Its location of
maximum brightness is at $(15^{+27}_{-38})^{\circ}$ longitude and
$(-3^{+13}_{-13})^{\circ}$ latitude, relative to the substellar point. There
appear to be two peaks in the dayside brightness distribution in the median
map in Figure 14, but we suggest that the posterior distribution shows this is
not a significant feature. The large uncertainty of the longitudinal position
of the brightness maxima reflects the flat emission structure along the
equator near the substellar point, rather than a range of distinct peaks, so
the assumed Gaussian distribution that leads to the large uncertainty may not
be a good assumption here. The latitudinal structure contributes significantly
to the eclipse shape, as identified in Hammond et al. (2024), and is
consistent with a brightness maximum on the equator. There is a small
asymmetry in the latitudinal structure which is not statistically significant,
given the width of the posterior distribution. The latitudinal structure
varies more smoothly than the longitudinal structure, which transitions from a
low gradient near the substellar point to sharp gradients near the
terminators.

Hammond et al. (2024) suggested that a uniform longitudinal structure near the
substellar point could be caused by dayside stationary Rossby waves. These
waves appear in some of the GCM simulations in Hammond et al. (2024) and
produce a less sharply peaked longitudinal structure (Lewis & Hammond, 2022).
However, they also produce a two-peaked latitudinal structure with maxima
north and south of the equator (Matsuno, 1966), but the data are not
sufficiently precise to resolve or rule out a two-peaked structure in Figure
14. More observations might constrain the dayside structure strongly enough to
detect such structures.

###  4.3 Comparison to theoretical expectations

To put these results in context, we compare them to the precision predicted by
the analytic “Eclipse Mapping Metric” in Boone et al. (2024). To calculate the
predicted $\mathrm{EMM}$ for the most precise simulated WASP-43b data in
Figure 2, we set the error per point to $\sigma=150$ ppm (note that Boone et
al. (2024) denotes the error per point as $\sigma_{1}$), the number of points
to $N_{\mathrm{points}}=7000$, the peak flux as $F_{0}=7000$ ppm, the
inclination to $i=82.106^{\circ}$, the semi-major axis to $a=4.859R_{*}$, and
the planetary radius to $R_{p}=0.15839R_{*}$. This results in an eclipse
mapping metric of $\mathrm{EMM}=23^{\circ}$, which predicts the smallest scale
resolvable by an eclipse map. This is consistent with Figure 14, where the map
can resolve variations in latitude and longitude on the scale of tens of
degrees.

We use the same set of parameters to calculate the $\mathrm{EMM}$ for the real
WASP-43b data in Figure 9, except for the precision per point which we change
to be $\sigma=250$ ppm, representing the accuracy achieved by averaging two
eclipses (the precision of a single point is $\sigma\approx 370$ ppm (Bell et
al., 2024)). This gives $\mathrm{EMM}=27^{\circ}$, implying that the map in
Figure 14 should also be able to resolve variations on the scale of tens of
degrees, although with less precision than the best-case simulated map. To
achieve the same eclipse mapping metric, we would need to observe five such
eclipses of WASP-43b.

To estimate the $\mathrm{EMM}$ for the WASP-18b dataset in Figure 11, we set
the error per point to $\sigma=121$ ppm, the number of points for a full phase
curve to $N_{\mathrm{points}}=9184$, the peak flux as $F_{0}=1000$ ppm, the
inclination to $i=84.39^{\circ}$, the semi-major axis to $a=3.483R_{*}$, and
the planetary radius to $R_{p}=0.0943R_{*}$. This results in
$\mathrm{EMM}=62^{\circ}$, which is significantly hampered by the low impact
parameter of the eclipse. Note that this is a less precise $\mathrm{EMM}$ than
the predicted value of $\mathrm{EMM}=47^{\circ}$ for WASP-18b in Boone et al.
(2024), which used an estimated value for the precision of NIRISS observations
of WASP-18b that was more precise than the actual observations. The eclipse
mapping metric for the longitudinal direction on WASP-18b is
$\mathrm{EMM}_{x}=21^{\circ}$, predicting that the effect of the longitudinal
structure on the eclipse shape should be resolvable in this case. The mismatch
between the longitudinal structure implied by the phase curve, and the
observed eclipse shape, implies that there may be unresolved systematic errors
affecting the eclipse shape, or that the dayside may have a structure that is
very uniform and is poorly represented by $\ell_{\mathrm{max}}=2$ harmonics.
Observing a full phase curve including a transit would help in fitting
astrophysical and instrumental models for an observation of WASP-18b, as
separating time-correlated noise from physical phase variations is difficult
without a full phase curve constrained by periodicity.

We also estimate the $\mathrm{EMM}$ for the HD 189733b data in Figure 11. We
set the error per point to $\sigma=252$ ppm, the number of points for a full
phase curve to $N_{\mathrm{points}}=7160$ (calculated from the higher-cadence
data which covers the eclipse), the peak flux to $F_{0}=3500$ppm, the
inclination to $i=85.68^{\circ}$, the semi-major axis to $a=8.863R_{*}$, and
the planetary radius to $R_{p}=0.153R_{*}$. This results in
$\mathrm{EMM}=36^{\circ}$, which is lower than the WASP-43b dataset but still
implies that mapping may be possible. The lack of a discernible residual shape
in Figure 11 and the related lack of a detection of a mapping signal via
cross-validation implies that the precision of the data is too low or there
are unresolved systematic errors.

Based on these results, we suggest that the eclipse mapping metric defined by
Boone et al. (2024) should be at least $\mathrm{EMM}=30^{\circ}$ to resolve an
eclipse map. Observing three additional eclipses of WASP-43b would make its
potential resolution comparable to our highest-precision simulated data.

##  5 Conclusions

In this study, we have presented a new method for fitting eclipse maps to
observations of exoplanets. This method provides two improvements on previous
work. Firstly, it provides a data-driven method to assess when an eclipse
mapping signal is resolvable in a dataset. Secondly, it provides a data-driven
method to assess the appropriate amount of spatial information to include when
fitting a map to this signal. The method has two stages:

  1. 1.

Test for the presence of an eclipse mapping signal by fitting a low-order
Fourier series model $\mathcal{M}(n_{\mathrm{max}}=2)$ and a low-order eclipse
map model $\mathcal{M}(\ell_{\mathrm{max}}=2)$, and comparing their cross-
validation scores. There is an eclipse mapping signal if the eclipse map model
has a better cross-validation score. The residual signal in ingress and egress
should confirm this. We suggest an $\ell_{\mathrm{max}}=2$ eclipse map model
and an $n_{\mathrm{max}}=2$ Fourier series model because these are simple
enough to fit alongside instrumental systematics and have the same freedom to
fit the light curve outside the eclipse.

  2. 2.

Once an eclipse mapping signal is identified, re-fit the map with as much
spatial freedom as possible and prevent overfitting by penalising the
likelihood by a factor proportional to the entropy of the map. The scale of
this factor is tuned to optimise the cross-validation score of the fitted map,
which identifies the correct penalty size that preserves an appropriate amount
of information content given the precision of the data. An eclipse map model
$\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$, oversampled in pixel space by a
factor of 3 (Luger et al., 2019), balances sampling time and spatial freedom.

We used k-fold cross-validation to calculate the cross-validation score, using
folds with a duration of one-fifth of the eclipse ingress or egress, to test
the ability of the fitted maps to predict an appropriate amount of spatial
information. This metric penalised both underfitting with too few degrees of
freedom and overfitting with too many degrees of freedom. To measure this
score in a reasonable time, we restricted the sampled folds to cover the
eclipse ingress and egress as well as an equally long section of the out-of-
eclipse phase curve on either side of the eclipse.

We conclude that this method is an improvement on fitting eclipse maps using
spherical harmonics alone without entropy-based regularisation. We showed that
the method correctly identified the presence of an eclipse mapping signal in
accurate and precise simulated datasets and the absence of an eclipse mapping
signal in inaccurate or imprecise simulated datasets. We then fitted a map of
appropriate precision to the datasets containing eclipse mapping signals.

We applied our method to three observational datasets – WASP-43b, WASP-18b,
and HD 189733b. The cross-validation test suggested that there was only a
robust eclipse mapping signal in the WASP-43b dataset. We used this to derive
a new eclipse map of WASP-43b, which makes use of the increased degrees of
freedom of a higher-order eclipse map model
$\mathcal{M}(\alpha,\ell_{\mathrm{max}}=4)$ to fit a flatter brightness
structure near the substellar point compared to previous low-order maps
(Hammond et al., 2024). This flat brightness structure could be consistent
with stationary wave structures seen in previous GCM simulations. However, the
lack of a double-peaked latitudinal structure is inconsistent with the
structure of these waves in the simulations. We conclude that there is
evidence for a small eastward dynamical hot-spot shift, a longitudinally broad
and uniform hot-spot, and a smoothly varying latitude structure with a
robustly detected equator-pole brightness difference.

The lower cross-validation scores for the eclipse mapping models for the
WASP-18b and HD 189733b datasets suggest that there is no eclipse mapping
signal in these datasets. Specifically, they show that the mapping models have
lower predictive ability for left-out data than the Fourier series model. We
suggest that the HD 189733b dataset has too low precision to resolve a mapping
signal. The WASP-18b dataset may have unresolved systematic effects that
prevent an accurate map being fitted to its sufficiently high-precision data,
or may have a complex day-side structure with largely uniform regions that are
not represented well by low-order spherical harmonics.

We conclude that eclipse mapping requires data of comparable quality to the
WASP-43b dataset and that care should be taken when fitting maps to future
datasets. Future work could improve the sampling method to allow fitting with
more pixels and more degrees of spatial freedom. An alternative model
comparison metric would be the Bayesian Evidence as estimated by nested
sampling algorithms, which would be prohibitively slow to calculate with our
method but could be used in future work. The k-fold cross-validation method
could also be varied to see if there is an alternative cross-validation metric
that is even better suited to comparing models of eclipse shape.

## Acknowledgements

We thank Louis-Philippe Coulombe, Adina Feinstein, and Nestor Espinoza for
providing access and insight to the WASP-18b NIRISS dataset. We thank Eric
Agol and Nicolas Cowan for providing access and insight to the HD 189733b
Spitzer dataset. M.H. acknowledges support from Christ Church, University of
Oxford. 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. Some of 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.

## Data Availability

The simulated and observational datasets, as well as the code used to
generate, fit, and plot them, are available at
https://zenodo.org/records/11394012.

## References

  * Akaike (1981) Akaike H., 1981, Journal of econometrics, 16, 3
  * Arlot & Lerasle (2016) Arlot S., Lerasle M., 2016, Journal of Machine Learning Research, 17, 1
  * Bell et al. (2024) Bell T. J., et al., 2024, arXiv preprint arXiv:2401.13027
  * Blecic et al. (2014) Blecic J., et al., 2014, ApJ, 781, 116
  * Boone et al. (2024) Boone S., Grant D., Hammond M., 2024, Monthly Notices of the Royal Astronomical Society, 528, 596
  * Buie et al. (1992) Buie M. W., Tholen D. J., Horne K., 1992, Icarus, 97, 211
  * Burnham & Anderson (2004) Burnham K. P., Anderson D. R., 2004, Sociological Methods & Research, 33, 261
  * Challener & Rauscher (2022) Challener R. C., Rauscher E., 2022, The Astronomical Journal, 163, 117
  * Challener & Rauscher (2023) Challener R. C., Rauscher E., 2023, The Astronomical Journal, 166, 176
  * Chen & Yang (2021) Chen Y., Yang Y., 2021, Stats, 4, 868
  * Chen et al. (2024) Chen X., et al., 2024, (submitted) MNRAS
  * Coulombe et al. (2023) Coulombe L.-P., et al., 2023, Nature, 620, 292
  * Cowan & Agol (2008) Cowan N. B., Agol E., 2008, The Astrophysical Journal Letters, 678, L129
  * Deitrick et al. (2020) Deitrick R., Mendonça J. M., Schroffenegger U., Grimm S. L., Tsai S.-M., Heng K., 2020, The Astrophysical Journal Supplement Series, 248, 30
  * Foreman-Mackey et al. (2020) Foreman-Mackey D., et al., 2020, exoplanet-dev/exoplanet v0.4.3, doi:10.5281/zenodo.1998447
  * Gillon et al. (2012) Gillon M., et al., 2012, A&A, 542, A4
  * Gull & Daniell (1978) Gull S. F., Daniell G. J., 1978, Nature, 272, 686
  * Hammond et al. (2024) Hammond M., et al., 2024, arXiv preprint arXiv:2404.16488
  * Hastie et al. (2009) Hastie T., Tibshirani R., Friedman J. H., Friedman J. H., 2009, The elements of statistical learning. Vol. 2, Springer
  * Hellier et al. (2011) Hellier C., et al., 2011, A&A, 535, L7
  * Horne (1985) Horne K., 1985, Monthly Notices of the Royal Astronomical Society, 213, 129
  * Kataria et al. (2015) Kataria T., Showman A. P., Fortney J. J., Stevenson K. B., Line M. R., Kreidberg L., Bean J. L., Désert J.-M., 2015, ApJ, 801, 86
  * Lewis & Hammond (2022) Lewis N. T., Hammond M., 2022, The Astrophysical Journal, 941, 171
  * Luger et al. (2019) Luger R., Agol E., Foreman-Mackey D., Fleming D. P., Lustig-Yaeger J., Deitrick R., 2019, The Astronomical Journal, 157, 64
  * Majeau et al. (2012) Majeau C., Agol E., Cowan N. B., 2012, The Astrophysical Journal Letters, 747, L20
  * Malik et al. (2017) Malik M., et al., 2017, AJ, 153, 56
  * Mansfield et al. (2020) Mansfield M., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 5151
  * Matsuno (1966) Matsuno T., 1966, Journal of the Meteorological Society of Japan Ser II, 44, 25
  * Mendonça et al. (2016) Mendonça J. M., Grimm S. L., Grosheintz L., Heng K., 2016, The Astrophysical Journal, 829, 115
  * Mendonça et al. (2018a) Mendonça J. M., Malik M., Demory B.-O., Heng K., 2018a, The Astronomical Journal, 155, 150
  * Mendonça et al. (2018b) Mendonça J. M., Tsai S.-m., Malik M., Grimm S. L., Heng K., 2018b, ApJ, 869, 107
  * Rauscher et al. (2007) Rauscher E., Menou K., Seager S., Deming D., Cho J. Y.-K., Hansen B. M., 2007, The Astrophysical Journal, 664, 1199
  * Rauscher et al. (2018) Rauscher E., Suri V., Cowan N. B., 2018, The Astronomical Journal, 156, 235
  * Salvatier et al. (2016) Salvatier J., Wiecki T. V., Fonnesbeck C., 2016, PeerJ Computer Science, 2, e55
  * Schwarz (1978) Schwarz G., 1978, Annals of Statistics, 6, 461
  * Snyder (1987) Snyder J. P., 1987, Map projections–A working manual. Vol. 1395, US Government Printing Office
  * Vehtari et al. (2017) Vehtari A., Gelman A., Gabry J., 2017, Statistics and computing, 27, 1413
  * Wahba (1990) Wahba G., 1990, Spline models for observational data. SIAM
  * Welbanks et al. (2023) Welbanks L., McGill P., Line M., Madhusudhan N., 2023, The Astronomical Journal, 165, 112
  * de Wit et al. (2012) de Wit J., Gillon M., Demory B.-O., Seager S., 2012, Astronomy & Astrophysics, 548, A128

Figure 15: First panel: the artificial emission map defined by Equation 16.
Second panel: the map represented by $\ell_{\mathrm{max}}=2$ spherical
harmonics, showing how this produces a slightly wider perturbation. Third
panel: the map represented by $\ell_{\mathrm{max}}=4$ spherical harmonics,
showing how they can capture most of its true spatial structure.

##  Appendix A Additional demonstration of optimising mapping information
content

Section 3.3 showed how to fit a map with optimised information content to a
simulated dataset from a numerical model of WASP-43b. In this appendix, we
demonstrate the model for a different input map, showing how it selects a
different spatial scale for the fitted map. We define an artificial flux map
$F_{p}/F_{S}$ that is entirely uniform except for a Gaussian perturbation:

$F_{p}/F_{S}=F_{0}(1+2e^{(-(r/r_{0})^{2}))},$ (16)

where $F_{0}=3000$ ppm is the magnitude of the map, $r$ is the angular
distance along the surface away from a point which we choose to be
$30^{\circ}$ east and $30^{\circ}$ north of the substellar point, and
$r_{0}=60^{\circ}$ is the scale of the perturbation. Figure 15 shows this map,
as well as $\ell_{\mathrm{max}}=2$ and $\ell_{\mathrm{max}}=4$ representations
of it. The $\ell_{\mathrm{max}}=2$ representation can capture the approximate
shape of the true map, but fails to capture the magnitude of the perturbation
correctly, while the $\ell_{\mathrm{max}}=4$ representation better matches the
magnitude.

We generate simulated light curves following the methodology in Section 2.1,
assuming the system parameters of WASP-43b again, and applying Gaussian noise
with standard deviations of 150 ppm and 250 ppm. Figure 16 shows our mapping
method applied to the two resulting datasets. The left-hand panel shows their
cross-validation scores as defined in Section 2.6 for eclipse map models
$\mathcal{M}(\ell_{\mathrm{max}}=2)$, relative to the scores for Fourier
series models $\mathcal{M}(n_{\mathrm{max}}=2)$. The eclipse mapping models
have significantly better cross-validation scores, confirming the strong
eclipse mapping signal in both datasets.

The right-hand panel of Figure 16 shows the cross-validation score for
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ eclipse map models fitted as described in
Section 2.7, as a function of the regularisation parameter $\alpha$ described
in Equation 15. The optimal value of $\alpha=4641.6$ is the same for both
datasets, and smaller than the optimal values in Figure 6 for the data
simulated from the GCM. The optimal value of $\alpha$ depends non-trivially on
the spatial scale and precision of the dataset. We found that simulating these
datasets again with different randomly generated Gaussian noise could result
in different optimal values of $\alpha$, as the mapping residuals have
relatively short durations and low signal to noise ratios. The optimised map
in the second row of Figure 17 is a much more precise fit to the true map than
the map in the first row which was fitted without optimising the information
content.

Figure 16: Left panel: The cross-validation scores (Equation 10) for eclipse
map models $\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the simulated
datasets defined by Equation 16, relative to Fourier series models
$\mathcal{M}(n_{\mathrm{max}}=2)$. Both datasets have strong eclipse mapping
signals, identified by the better cross-validation scores for the eclipse
mapping models $\mathcal{M}(\ell_{\mathrm{max}}=2)$. The relative scores are
normalised by their standard errors (Equation 12). Right panel: The relative
cross-validation scores for eclipse map models
$\mathcal{M}(\ell_{\alpha,\mathrm{max}}=4)$, as a function of $\alpha$
(Equation 15). The optimal value of $\alpha$ for both maps is smaller than in
Figure 6.

Figure 17: First row: an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=4)$ fitted to a light curve simulated from
Equation 16 with precision 150 ppm as described in Appendix A. Fitting the the
spherical harmonics alone results in a highly imprecise map. Second row: an
eclipse map model $\mathcal{M}(\alpha=4642,\ell_{\mathrm{max}}=4)$ with the
parameter $\alpha$ optimised according to the cross-validation score shown in
Figure 16. Figure 18: Supplementary Figure: an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the WASP-18b JWST NIRISS/SOSS
observations. Section 4.1 and Figure 12 suggest that this model is not
preferred over a simpler Fourier series model that does not fit an eclipse
map. This model appears to produce a well-constrained posterior distribution,
but this is not consistent with the eclipse shape in Figure 11. This shows why
statistical metrics like cross-validation are important, as it shows that this
apparently precise map model cannot predict out-of-sample sections any better
than a model assuming the eclipse shape of a uniform planet. Figure 19:
Supplementary Figure: an eclipse map model
$\mathcal{M}(\ell_{\mathrm{max}}=2)$ fitted to the HD 189733b Spitzer
observations. Table 2 shows how this model achieves a better $\chi^{2}$ value
than a Fourier series model that does not fit an eclipse map. However, Section
4.1 and Figure 12 use metrics that penalise overfitting to suggest that this
model is not in fact preferred over the Fourier series model.