Two-Dimensional Eclipse Mapping of the Hot Jupiter WASP-43b with JWST MIRI/LRS
Mark Hammond Taylor J. Bell Ryan C. Challener Neil T. Lewis Megan Weiner Mansfield Isaac Malsky Emily Rauscher Jacob L. Bean Ludmila Carone João M. Mendonça Lucas Teinturier Xianyu Tan Nicolas Crouzet Laura Kreidberg Giuseppe Morello Vivien Parmentier Jasmina Blecic Jean-Michel Désert Christiane Helling Pierre-Olivier Lagage Karan Molaverdikhani Matthew C. Nixon Benjamin V. Rackham Jingxuan Yang






0000-0002-6893-522X]Mark Hammond
Atmospheric, Oceanic, and Planetary Physics, Department of Physics, University of Oxford, Oxford, UK
0000-0003-4177-2149]Taylor J. Bell
BAER Institute, NASA Ames Research Center, Moffet Field, CA, USA
Space Science and Astrobiology Division, NASA Ames Research Center, Moffett Field, CA, USA
0000-0002-8211-6538]Ryan C. Challener
Department of Astronomy, Cornell University, Ithaca, NY, USA
Department of Astronomy, University of Michigan, Ann Arbor, MI, USA

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

Steward Observatory, University of Arizona, Tucson, AZ, USA
NHFP Sagan Fellow

Department of Astronomy, University of Michigan, Ann Arbor, MI, USA
0000-0003-3963-9672]Emily Rauscher
Department of Astronomy, University of Michigan, Ann Arbor, MI, USA
0000-0003-4733-6532]Jacob L. Bean
Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL, USA

Space Research Institute, Austrian Academy of Sciences, Graz, Austria
0000-0002-6907-4476]João M. Mendonça
DTU Space, Technical University of Denmark, Kgs. Lyngby, Denmark
0000-0002-0797-5746]Lucas Teinturier
LESIA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CNRS, Meudon, France
Laboratoire de Météorologie Dynamique, IPSL, CNRS, Sorbonne Université, Ecole Normale Supérieure, Université PSL, Ecole Polytechnique, Institut Polytechnique de Paris, Paris, France


Tsung-Dao Lee Institute, Shanghai Jiao Tong University, 520 Shengrong Road, Shanghai, People’s Republic of China
School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai, People’s Republic of China

0000-0001-7866-8738]Nicolas Crouzet
Leiden Observatory, Leiden University, Leiden, The Netherlands

Max Planck Institute for Astronomy, Heidelberg, Germany
0000-0002-4262-5661]Giuseppe Morello
Department of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, Sweden
Instituto de Astrofísica de Canarias (IAC), La Laguna, Tenerife, Spain
Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain

Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, France


0000-0002-0769-9614]Jasmina Blecic
Department of Physics, New York University Abu Dhabi, Abu Dhabi, UAE
Center for Astro, Particle and Planetary Physics (CAP3), New York University Abu Dhabi, Abu Dhabi, UAE

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

Space Research Institute, Austrian Academy of Science, Schmiedlstrasse 6, 8042 Graz, Austria
Institute for Theoretical Physics and Computational Physics, Graz University of Technology, Petersgasse16/II, 8010 Graz, Austria

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

Universitäts-Sternwarte, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, D-81679 München, Germany
Exzellenzcluster Origins, Boltzmannstraße 2, 85748 Garching, Germany
0000-0001-8236-5553]Matthew C. Nixon
Department of Astronomy, University of Maryland, College Park, Maryland, USA

Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA, USA
0009-0006-2395-6197]Jingxuan Yang
Atmospheric, Oceanic, and Planetary Physics, Department of Physics, University of Oxford, Oxford, UK



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









§ INTRODUCTION
 

Eclipse mapping measures the two-dimensional spatial features of an object when it is eclipsed by another object. The eclipsed object is covered along one axis during the eclipse ingress, and is revealed along another axis during the eclipse egress. Combining the information from these two axes reveals two-dimensional information about the surface of the eclipsed object <cit.>.

Eclipse maps have been derived for objects like Pluto during its eclipse by Charon <cit.>, or the white dwarf BD +16^∘516B during its eclipse in a binary system <cit.>. <cit.> proposed applying this technique to exoplanets, and <cit.> showed that the James Webb Space Telescope (JWST) would have sufficient precision to accurately map exoplanets. Eclipse mapping is currently the only method by which 2D information can be measured for exoplanets, as out-of-eclipse “phase curves” provide low-resolution information as a function of longitude only. 2D spatial information is crucial for understanding atmospheric circulation <cit.>, chemical composition <cit.>, and cloud structure <cit.>.

<cit.> and <cit.> derived eclipse maps of the hot Jupiter HD 189733b using observations with the Spitzer Space Telescope. Their analyses were restricted to fitting large-scale shapes with a high degree of uncertainty due to the limited precision of the measurements, uncertainty in and degeneracy with orbital parameters, and the dependence of the result on the mapping model. <cit.> described the uncertainty in the mapping signal due to the impact parameter, stellar density, eccentricity, and argument of periastron in particular. <cit.>, <cit.>, and <cit.> developed more advanced methods to fit eclipse maps of 2D and 3D structure, based on fitting orthogonal phase curves rather than orthogonal surface maps, and applied these methods to the previous observations of HD 189733b. <cit.> presented the first eclipse map measured with JWST, deriving a map of thermal emission from observations of the hot Jupiter WASP-18b from 0.85 to 2.85 μm with the NIRISS/SOSS instrument. They found no longitudinal shift in the hot-spot, sharp gradients in brightness towards the terminator, and no clear measurement of latitudinal structure.

This study presents eclipse maps derived from broadband JWST MIRI/LRS <cit.> observations from 5 to 10.5 μm of a full phase curve of the hot Jupiter WASP-43b containing two eclipses and one transit <cit.>. WASP-43b is a “hot Jupiter” exoplanet with strong thermal emission from its permanent dayside, and a short, tidally-locked orbit around a K7 main sequence star that exhibits low variability <cit.>. These properties make it ideal for precise time-series observations of thermal emission. 

<cit.> analysed this MIRI/LRS phase curve, finding a large difference in dayside and nightside brightness temperatures and evidence for water absorption. We fit its inclination to be 82.11 ^+ 0.050  _-0.052, and the ratio of its semi-major axis to stellar radius to be 4.859^+0.013_-0.012, which corresponds to an impact parameter of 0.667^+0.006_-0.006. This means that the edge of the star crosses the planet at an angle of ∼42^∘ (the “stellar edge angle” defined in <cit.>). The longitudinal and latitudinal features of the day-side therefore affect the eclipse mapping signal almost equally, with a small bias towards longitudinal features. This geometry makes WASP-43b especially well suited for eclipse mapping.

We fit an eclipse map to this dataset, finding the clearest eclipse mapping signal to date and separating the latitudinal and longitudinal parts of this signal for the first time. In Section <ref> we describe the methods used to fit the eclipse maps. Section <ref> then shows the different maps we fit to the observations, and compares their structures and statistical evidence. These maps are then compared to numerical simulations in Section <ref>. In Section <ref> we conclude that the spherical harmonic components up to ℓ_ max=2 are constrained well, but that further observations are needed to precisely constrain higher-order mapping structures.









§ DATA AND METHODS
 

A full orbit of WASP-43b with the JWST MIRI/LRS slitless (SL) mode was observed as part of the JWST-ERS-1366 program, performing target acquisition with the F1500W filter and using the SLITLESSPRISM subarray for the science observation. The science observation lasted 26.5 hours, consisting of 9316 integrations lasting 10.34 s each. The two eclipses show a downwards systematic trend with time, similar to that identified in other MIRI/LRS time-series observations <cit.>. We discard the data beyond 10.5 μm due to the “shadowed region effect” described in <cit.> alongside a more detailed analysis of the data acquisition.

Figure <ref> shows the raw MIRI/LRS dataset (black points) in units of planetary flux divided by stellar flux. We use the  pipeline <cit.> to reduce the observed data, starting from the  Stage 4 output from the “Eureka v1” reduction in <cit.>. We trim the initial 780 integrations (not plotted) where instrumental systematic effects are the strongest. The only subsequent methodological difference to <cit.> is the method used to model the planetary emission in the phase curve — we used a  model of a planet with a 2D emission map over its surface <cit.>, instead of a Fourier series model of the phase curve.



 §.§ Astrophysical Model


We fit a 2D map of thermal emission to this phase curve, simultaneously with the orbital, planetary, stellar, and instrumental systematic parameters. Given the significant effect of the planetary emission map on the observed eclipse shape, this should derive more accurate parameters than those derived using a Fourier series model in <cit.>. This simultaneous fit also tests if there are any degeneracies between the system parameters and the eclipse map, such as between the eclipse timing and the longitudinal emission offset <cit.>. 

We fit the planet-to-star radius ratio, the linear ephemeris, the inclination, the ratio of the semi-major axis to the stellar radius, and the two parameterised stellar limb-darkening parameters in <cit.>. 


We set the orbital period constant at 0.813474 days as it is known with sufficient precision already <cit.>. We set the obliquity of the planet to be 0 as we assume a tidally locked orbit. We set the stellar radius to be constant at 0.665 R_⊙ <cit.>, although this is an arbitrary value that serves only to give our model a dimensional form, as the ratio F_p/F_S is only sensitive to the ratios R_p/R_* and a/R_* which we fit separately. The stellar mass does affect the signal <cit.> but is entirely determined in our model by the fitted value of a/R_* and the fixed value of the orbital period. 

<cit.> showed how orbital eccentricity can be degenerate with an eclipse mapping signal. <cit.> constrained the eccentricity of WASP-43b to be 0.0035_-0.0025^+0.0060, which has generally been used to assume zero eccentricity for this planet <cit.>. To check this assumption, we separately fitted the model including the eccentricity and periastron, finding an eccentricity of 0.001078_-0.00025^0.00185, and no meaningful constraint on periastron, as discussed in Appendix <ref>. We take this to verify the assumption of zero eccentricity, and proceed to fit the models with eccentricity set to zero. This is consistent with the circularization timescale of 3 Myr that we estimate using <cit.> (assuming Q_P=10^6), which is much shorter than the age of the system estimated by <cit.> to be 400^+200_-100 Myr.

For the instrumental systematic effects, we fit a uniform baseline and linear trend in time, the magnitude and timescale of an exponential ramp, and trends of the spatial position and width of the data on the detector. Table <ref> shows the new parameters fitted for the ℓ_ max=2 eclipse map model, compared to the parameters fitted with an n=2 Fourier series model (like in <cit.>).

The difference in the time it takes light to travel from the planet and the star has an important effect on this dataset, so we include it in our model. In eclipse, the light from the planet takes about 8.5 seconds longer to reach the observer than the light from the star. For the average gradient in eclipse ingress and egress of roughly 8 ppm per second (see Figure <ref>), this corresponds to a maximum effect of around 70 ppm on the eclipse shape. Figure <ref> shows that this corresponds to about 15% of the maximum deviation in eclipse shape due to the eclipse mapping signal itself.



 §.§ Eclipse Mapping


The eclipse mapping method works by constructing a map from a number of “basis maps”, which each have an associated basis phase curve. We match the observed light curve with a sum F(t) constructed from the the basis phase curves f_i(t) weighted by coefficients c_i:


    F(t) = ∑_i c_i f_i(t).


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


    Z(θ,ϕ) = ∑_i c_i z_i(θ,ϕ).


<cit.> shows the phase curves f_i(t) of each spherical harmonic z_i(θ,ϕ) (where θ and ϕ are longitude and latitude). We model the astrophysical signal using [<starry.readthedocs.io>] <cit.> and fit the coefficients c_i using PyMC3 <cit.>. We use “pixel sampling” to enforce a positive emission map globally (e.g. as used in <cit.>). This method samples the brightness of pixels distributed uniformly over the mapped surface and transforms them to spherical harmonic coefficients c_i to compute the actual phase curve. To fit the ℓ_ max=2 eclipse map alongside the orbital and systematic parameters, we sample 16 pixels evenly spaced on a Mollweide projection to represent the spherical harmonic space with 4ℓ^2 pixels <cit.>. We use a (natural) log-normal prior for the pixels with a mean magnitude of 6000 ppm and a standard deviation of 3000 ppm (transformed to log-space).

After fitting the pixels representing the ℓ_ max=2 eclipse map simultaneously with the orbital, planetary, stellar, and systematic parameters, we fix these additional parameters to their derived values, and re-fit the map with a variety of methods. We fit maps with:


    
  * ℓ_ max=3 spherical harmonics
    
  * An eigenmap as described in <cit.>
    
  * ℓ_ max=2 spherical harmonics to the first eclipse alone, the second eclipse alone, and both eclipses


The next section presents the results from each of these re-fitted maps, and describes how the details of each fitting procedure affect the derived map. It also presents the statistical evidence for an eclipse mapping signal, and compares the different mapping models.






§ RESULTS
 



 §.§ Eclipse map fitted with orbit and systematics


The red line in Figure <ref> shows the phase curve resulting from the ℓ_ max=2 eclipse map fitted alongside the orbital parameters, stellar parameters, and systematic model. The blue line in Figure <ref> shows the model of instrumental systematics fitted at the same time. We normalised the F_p/F_S values (the ratio of planetary to stellar flux) in Figure <ref> so that the systematic model has zero mean over the course of the fitted observations. 

Appendix <ref> shows the posterior distribution for the fitted orbital, planetary, stellar, and systematic parameters. Some of the parameters are consistent with those derived using a Fourier series fit like that used in <cit.>, but some are not consistent. There are small but statistically significant differences in the planetary radius, eclipse timing, and inclination, all of which we expect to be more accurately fitted by the more realistic eclipse shape in the eclipse mapping model. The posterior distribution of the stellar limb darkening parameters are different but we found that, when combined, these resulted in consistent limb darkening profiles. The stellar limb darkening is relatively weak at these long wavelengths, so is poorly constrained but is consistent between the two models. The systematic parameters are different to those produced by a Fourier series fit <cit.> but result in an almost identical model of systematics, as a short exponential ramp plus a linear trend is almost identical to a long exponential ramp.

Figure <ref> shows that this ℓ_ max=2 eclipse map model fits the observed phase curve well. The information about the 2D eclipse map structure is contained within the shapes of the ingress and egress of each eclipse. Figure <ref> shows the residual difference in ingress and egress between the observations (shown as black, green, and purple points), and a “control” model representing the phase curve as a Fourier series up to order n=2. This Fourier series was fitted with orbital and systematic parameters fixed to the values derived when fitting the ℓ_ max=2 eclipse map, listed in Table <ref>, so that the residual is only due to the different emission models. 

The green and purple points in Figure <ref> show the residual of the observed data in the first and second eclipses, binned every 8 points, and the black points show the average of the two eclipses. The ℓ_ max=2 eclipse mapping model and the n=2 Fourier series model fit the out-of-eclipse phase curve as well as each other, so the residual outside the range plotted in Figure <ref> is determined by the uncertainty in the data. This is not surprising, as both models have access to almost identical fitting functions for the out-of-eclipse phase curve – the ℓ_ max=2 harmonics produce the same phase curves as the n=2 Fourier series, apart from a small modification from the orbital obliquity.

The size of the residual between the observed data and the n=2 control model shows the size of the eclipse mapping signal, which is the effect of partial stellar coverage of non-uniform emission from the planet <cit.>. For example, when the star covers an area that is emitting more than the average of the planetary disk, the observer measures less flux than would be expected for a uniform disk, so the residual signal is negative. 

The blue shaded area shows the range of phase curves fitted using the ℓ_ max=2 eclipse map model, with the Fourier series control model subtracted. The two shaded regions show the first and second quantiles, containing 68.27% and 95.45% of the posterior distribution respectively. The eclipse map model matches the residual signal in ingress and egress, with a uncertainty comparable to the error bars of the observed data. This implies a good fit to a robust eclipse mapping signal, which we quantify in the next section.

Figure <ref> shows the ℓ_ max=2 eclipse map itself. The plotted 2D map is constructed from the median values of the posterior distribution of each of the spherical harmonic coefficients. This median map has a hot-spot near the substellar point, with a small meridionally averaged shift of (7.75^+0.36_-0.36)^∘ eastwards and a shift of ( 10.72 ^+ 4.14  _- 4.68 )^∘ southwards[The observations cannot distinguish between an orbit about a vector pointed “up” or “down” in the sky, so our choice of north and south poles is made arbitrarily and kept consistent for plotting. The eastward direction corresponds to the direction of rotation of the planet.].

We present the maps as a ratio of planetary flux to stellar flux, as a conversion to temperature is more complex and requires modelling or assumptions about the temperature structure. We make a simple estimate of the brightness temperature assuming planetary black-body emission B_T,λ in the bandpass λ=5 to 10.5 μm, and a PHOENIX <cit.> stellar model spectrum with effective temperature T_ eff = 4300 K and surface gravity log g = 4.50 as used in <cit.>. The stellar radius that we assumed to be 0.665 R_⊙ does affect the value of the brightness temperature. We do not include effect of the uncertainty of this and other stellar parameters on the derived brightness temperature, in order to highlight the uncertainty due to mapping alone. This determines that this ℓ_ max=2 map corresponds to a brightness temperature of ( 1790.0 ^+ 23.0  _- 29.0 ) K at the substellar point, ( 1293.0 ^+ 27.0  _- 34.0 ) K at the equatorial east terminator, ( 1114.0 ^+ 30.0  _- 36.0 ) K at the equatorial west terminator, ( 1100.0 ^+ 128.0  _- 109.0 ) K at the south pole, and ( 1011.0 ^+ 111.0  _- 108.0 ) K at the north pole (where the distinction between the poles is our arbitrary choice).

The posterior distribution of the fitted maps is shown east-west along the equator and north-south through the substellar point. This range of fitted maps includes any degeneracies with the orbital or systematic parameters as these are fitted simultaneously. We confirmed that the “median map” plotted in 2D corresponds closely to the center of the east-west and north-south posterior distributions (not shown). The maximum likelihood map is very similar to the median map, except near the poles where its latitudinal structure deviates slightly from the median but remains within the first quantile.

Figure <ref> shows an ℓ_ max=2 eclipse map fitted to the same dataset, but with all of its fitting basis maps averaged in latitude, removing the ability of the map to fit latitudinal structure. This makes some of the basis maps uniform everywhere, so we remove them from the fitting process (reducing the number of parameters). The resulting map is very tightly constrained by its need to fit the out-of-eclipse phase curve. Its residual signal in ingress and egress is shown by the red region in Figure <ref>, which fails to match the observed residual points. This shows the presence of a latitudinal mapping signal of around 200 ppm (the difference between this longitude-only fit, and the observed residual), compared to the longitudinal mapping signal of around 250 ppm (the magnitude of the residual of this longitude-only fit).






 §.§ Eclipse Mapping Signal


We can quantify the significance of the signal in Figure <ref> by comparing the Bayesian evidence of the ℓ_ max=2 eclipse map model (the blue shaded region) and the n=2 Fourier series control model (the zero line) in Figure <ref>. As applied in <cit.> to a similar analysis of exoplanet phase curves, the log-likelihood function is:


    log L=-1/2 σ^2∑_i=1^N(F(t_i)-D_i)^2-N/2log 2 πσ^2,


for data D_i at N times t_i fitted by a model F(t) with variance σ^2 . The log-odds ratio ln O = ln L_1 - ln L_2 then represents the relative performance of models 1 and 2, with a value of ln O above 1 corresponding to a better fit to model 1. The log-odds ratio for the ℓ_ max=2 map model compared to the Fourier series model (with the eclipse shape of a uniformly emitting disk) is 52.6, showing that the eclipse model is overwhelmingly preferred according to the criterion in <cit.>. Both models fit the out-of-eclipse data very similarly as they have access to almost identical out-of-eclipse fitting functions, so this difference in likelihood is entirely due to the improved fit in the ingress and egress of the eclipses shown in Figure <ref>. Later, we use the Bayesian Information Criterion (BIC) to compare the fit quality of the eclipse map models while taking the number of fitted parameters into account. We cannot use the BIC to compare a map model to the Fourier series model as they are not nested.

Notably, the log-odds ratio of an n=2 Fourier series fit with free orbital, planetary, stellar, and systematic parameters (used to derive the parameters in Table <ref>) is 20.7 compared to the ℓ_ max=2 map model. This shows how the non-mapping parameters are modified by the Fourier series model to slightly better fit the ingress and egress shape. However, this makes a minor difference overall and the eclipse map model is still very strongly preferred.

Comparing the eclipse map and Fourier series models confirms that there is a strong eclipse mapping signal in this dataset. This does not immediately imply the presence of 2D information in the dataset, as <cit.> only found evidence for longitudinal information in their eclipse map of WASP-18b. We can show that latitudinal information is present by comparing the 2D ℓ_ max=2 map to the “1D” ℓ_ max=2 map where all of the basis maps are averaged in latitude, shown in Figure <ref>. The red region in Figure <ref> shows the residual eclipse mapping signal for this 1D map, which does not fit the real dataset as well as the 2D map.  

The log-odds ratio of the 2D map to the 1D map fits in Figure <ref> is 9.8, showing that the 2D map with variable latitudinal structure is greatly preferred.  Notably, the uncertainty on the residual of the 1D map in Figure <ref> is very small, as the longitudinal structure is well constrained by the additional information from the out-of-eclipse phase curve. This implies that most of the uncertainty on the 2D model in Figure <ref> is related to uncertainty about the latitudinal structure of the map.

<cit.> conducted the same comparison between 1D and 2D eclipse maps of WASP-18b and found no evidence that the 2D map was preferred over the 1D map. They therefore concluded that the eclipse mapping signal revealed longitudinal information only. In contrast, the eclipse mapping signal in this study of WASP-43b is sensitive to both longitudinal and latitudinal information, as would be expected from the higher impact parameter of the orbit. This is therefore the first time that latitudinal information has been shown to be detected on an exoplanet.




 §.§ Eclipse map fitted with ℓ_ max=3 spherical harmonics


The spherical harmonic order of the ℓ_ max=2 map fit to the dataset was limited by the number of parameters that could be fitted simultaneously with the orbital, stellar, and systematic parameters. Figure <ref> in Appendix <ref> shows the posterior distribution for these parameters, demonstrating that there are no significant degeneracies. We therefore fix these parameters to the values (listed in Table <ref>) derived using the ℓ_ max=2 eclipse map model and re-fit the map with a variety of methods.

<cit.> found that an n=1 Fourier series fitted the out-of-eclipse phase curve very poorly, so we do not consider any ℓ_ max=1 spherical harmonic fits. The top row of Figure <ref> shows the eclipse map re-fitted using spherical harmonics up to order ℓ_ max=3.  We fix the orbital, planetary, stellar, and systematic parameters to the median values derived from the fit in Figure <ref>. We found that the ℓ_ max=3 maps required too many pixels to sample effectively with PyMC3, so in this case we sampled the spherical harmonic coefficients directly with mc3 <cit.>. This method allows us to impose positivity by excluding negative maps from the fitting process, instead of sampling pixels with positive priors. We set the spherical harmonic coefficients to have Gaussian priors with a mean of 0 and a standard deviation equal to the mean of the observed planetary flux.

The resulting map in the top row of Figure <ref> shows a similar median structure to the map in Figure <ref>, with more uncertainty due to the increased degrees of freedom. There is a longitudinal offset of (0.50^+14.79_-8.04)^∘, which  is consistent with the ℓ_ max=2 fit value of (7.75^+0.36_-0.36)^∘. Its median value is almost zero, which may seem unusual as the overall phase curve has an offset of (7.3^+0.4_-0.4)^∘. This apparent discrepancy is due to the more complex shapes allowed by the ℓ_ max=3 spherical harmonics, which produce a dayside hot-spot which peaks at the substellar point but extends further east than west (as shown in Figure <ref>) resulting in an overall shift in the phase curve. We discuss this distinction between hot-spot position and phase curve offset in Section <ref>, showing how the phase curve offsets and hot-spot positions can differ in numerical simulations. 

The latitudinal structure of the ℓ_ max=3 map is somewhat different to the ℓ_ max=2 map, being flatter near the equator with no significant latitudinal hot-spot shift. We will also discuss in Section <ref> how latitudinal structure could be degraded by using low-order basis maps to fit the data. For all these reasons, the spatial resolution of the map used to fit the data is a key question. In Section <ref>, we will therefore discuss how many parameters  are justified to fit the dataset.




 §.§ Eclipse map fitted with eigenmapping


The bottom row of Figure <ref> shows a map fitted using “eigenmapping” <cit.> with ThERESA <cit.>, with the orbital and systematic parameters again fixed to the median values derived from the ℓ_ max=2 fit in Figures <ref> and <ref>.

Briefly, eigenmapping starts with spherical harmonic phase curves up to some degree ℓ_ max, orthogonalizes them with principal component analysis to create “eigencurves”, and fits the observed phase curve as a sum of N_E highest-variance eigencurves. ℓ_ max and N_E are selected to minimize the BIC. Each eigencurve has a corresponding eigenmap, and the fitted map is the corresponding sum of these eigenmaps. We again enforce a positive-flux constraint on the fitted map globally. Eigenmapping has been used to map  observations of HD 189733b with the Spitzer Space Telescope <cit.>, and observations of WASP-18b with JWST NIRISS/SOSS <cit.>. 

We find the best BIC for ℓ_ max = 2 and N_E = 6, a larger number of eigencurves than used for WASP-18b in <cit.>, showing the improved precision of this WASP-43b dataset. We find a slight eastward hot-spot offset of (7.5^+0.5_-0.5)^∘. This eigenmap achieves the best BIC score in Table <ref>; its fitted map is almost the same as the ℓ_ max=2 map in Figure <ref>, suggesting that it achieves this improved BIC score by discarding the map components that contribute the least to the observed phase curve. This means that it produces the same fit quality using fewer parameters (see Section <ref> for a related discussion of the mapping “null space”).






 §.§ Eclipse maps fitted with single eclipses


The maps in Figures <ref> and <ref> are fitted to the entire dataset, including 2D information from both eclipses as well as 1D information from the rest of the phase curve. To isolate the 2D mapping information in each eclipse, Figure <ref> shows the result of re-fitting the ℓ_ max=2 eclipse map using the first and second eclipses only. As in Figure <ref>, we fix the orbital and systematic parameters to the median values derived with the model in Figures <ref>, listed in Table <ref>. The exceptions are the magnitude of the exponential ramp and the magnitude of the linear trend, which we re-fit as these require different values given the different starting points in time

The first row in Figure <ref> shows an ℓ_ max=2 eclipse map fitted using a ∼2.2 hour section of data containing the first eclipse only with a small section of phase curve either side. We only show the map on the dayside of the planet, as there is no information about the nightside contained in this limited dataset. The map is similar to Figure <ref> and <ref>, but has a larger latitudinal hot-spot shift than the map in Figure <ref>. The latitudinal peak of the fit to the first eclipse only is weakly constrained, showing how an observation of only a single eclipse is likely to be insufficient to derive an accurate eclipse map.  There is a greater range of fitted maps at -90^∘ latitude than +90^∘ latitude because the inclination of the orbit angles this pole away from the observer.

The second row in Figure <ref> shows an eclipse map fitted to the second eclipse only, with a smaller latitudinal hot-spot shift than the fit to the first eclipse only. The third row in Figure <ref> shows an eclipse map fitted to both eclipses together. We suggest that the original fit to the entire dataset in Figure <ref> which found a latitudinal shift of ( -10.72 ^+ 4.14  _- 4.68 )^∘ is better constrained in latitude than this two-eclipse fit (despite containing no additional latitudinal information) because the additional measurement of the rest of the phase curve provides independent longitudinal information. This breaks degeneracies between the longitudinal and latitudinal information in the eclipse mapping signal <cit.>. 

The maps in Figure <ref> show that the eclipse mapping signals in both eclipses are mostly consistent with each other, as would be expected from the similar residual signals in Figure <ref>. However, they have different degrees of latitudinal asymmetry, which also manifests in the overall fit in Figure <ref>. Figure <ref> shows that the first eclipse implies more latitudinal asymmetry. This may be a real effect, but it may also be due to the increased effect of instrumental systematics on this eclipse. The blue line in Figure <ref> shows that the systematic signal is still relatively large during the first eclipse, but that it is almost negligible during the second eclipse. This shows the utility of observing a full phase curve where the periodicity of the astrophysical signal isolates it from the systematic signal. 





 §.§ Model Selection



To compare our fitted eclipse maps, we calculate the Bayesian Information Criterion (BIC) <cit.> and the Akaike Information Criterion (AIC) <cit.>. The BIC is:


    BIC=k ln (N)- 2 ln (L),


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


    ln (L) = -1/2∑_i^N( M_i - D_i/σ_i)^2,


for a sum over all N data points D_i with uncertainty σ_i, fitted by a model M_i. This uncertainty is derived from the residuals of the originally fitted phase curve, by scaling the expected error by a factor referred to as the “Uncertainty Scaling Factor” in Figure <ref> (see <cit.>).

A smaller BIC implies a better model, with a better fit or fewer fitting parameters. The difference in BIC is the relevant quantity for model comparison, so Table <ref> shows the ΔBIC for each model compared to Model M5, the eigenmap in Figure <ref>. This criterion allows comparison of nested models where it is assumed that the true model is inside the set of tested models <cit.>. 

We also consider the AIC of each model, which is:


    AIC= 2 k - 2 ln (L),


Due to the large number of points N in this dataset, the BIC applies a stronger penalty for each parameter than the AIC. We also present the “weights” of each model, wBIC and wAIC, which are the relative probabilities of each compared model in a set, based on their ΔBIC or ΔAIC. The wBIC for a particular model is:


    wBIC=exp(-ΔBIC / 2)/∑_r=1^R exp(-ΔBIC_r / 2)


where ΔBIC is compared to the best-performing model, and there are R total models being compared with a score of ΔBIC_r each. wAIC depends on ΔAIC in the same way. Table <ref> summarises our compared models, showing their χ^2 values, their ΔBIC and ΔAIC relative to the model with the best BIC and AIC, and the implied probability of each model.

Model M1 is the Fourier series model used in <cit.>. It achieves a significantly worse χ^2 value than all the eclipse mapping models, due to a worse fit to the eclipse shape. Model M2 is the eclipse map model shown in Figure <ref>, where the latitudinal structure is fixed to be flat. This achieves a better χ^2 value than the Fourier series model, as it includes the effect of longitudinal map structure on the eclipse shape. However, its χ^2 value is worse than the fully two-dimensional eclipse map fits, as shown by its poor fit in Figure <ref>. Interestingly, due to its smaller number of parameters, it achieves the second-best BIC score of all the eclipse mapping models due to the strong penalty applied to the number of parameters by the BIC. However, it achieves the worst AIC score of the eclipse map models, where the number of parameters is penalised less strongly. We can confidently reject models M1 and M2 compared to the 2D eclipse mapping models.

Model M3 is the ℓ_ max=2 eclipse mapping model in Figure <ref>. It achieves a better χ^2 value than models M1 and M2, as shown by its good fit in Figure <ref>. It has a worse BIC score than the eigenmap (model M5), as they have similar χ^2 values, but model M3 uses two more parameters. It has a more similar AIC score to the eigenmap, where these additional parameters are penalized less heavily.

Model M4 is the ℓ_ max=3 eclipse mapping model in Figure <ref>. It achieves the best χ^2 value due to its increased degrees of mapping freedom, but has poor BIC and AIC scores due to its increased number of parameters. These metrics imply that the data quality does not justify this number of degrees of freedom, which is also implied by the large uncertainty on the ℓ_ max=3 map in Figure <ref>.

Model M5 is the ℓ_ max=2 eigenmap using 6 basis maps. It achieves a similar χ^2 to the other eclipse mapping models, and has the best BIC and AIC scores due to its reduced number of parameters. It achieves these by discarding mapping structures that contribute weakly to the observed phase curve, so is able to match the observations with fewer degrees of freedom.

It is not clear whether the BIC or AIC is the better metric for model comparison, or if another metric would be more appropriate. The BIC assumes that the set of fitted models includes the “true” physical system, while the AIC assumes that all of the tested models are inexact representations of this system <cit.>. It could be argued that eclipse maps fitted to real planets do include the “true” map, as spherical harmonics form a complete basis set on the sphere. On the other hand, information about the real map is inevitably lost in the mapping process due to the need to truncate fits to low-order harmonics, and due to the presence of a “null space” (see Section <ref> for discussion of both issues). 

To summarise our model comparison,


    
  * The Fourier series model (M1), and the eclipse map with flat latitudinal structure (M2) perform badly on all statistical metrics, so we confidently reject them.
    
  * The latitude-longitude ℓ_ max=2 eclipse map (M3) is favoured over the longitude-only eclipse map (M2), showing the presence of latitudinal information in the data.
    
  * The ℓ_ max=3 eclipse map achieves a better χ^2 value than model M3, but at the cost of many more parameters, so has worse BIC and AIC scores.
    
  * The ℓ_ max=2 eigenmap achieves the best BIC and AIC scores, producing a fit of comparable quality to model M3 with two fewer parameters


Therefore, the eigenmap (model M5) is the best performing model on these metrics. However, it requires model M3 to fit the orbital, stellar, and systematic parameters simultaneously first, and then both models produce a similar result in the end. The eigenmap achieves a better AIC and BIC by discarding structures that do not contribute strongly to the map, producing a very similar map to model M3 overall. We suggest that more work can be done on the process of fitting eclipse maps and comparing levels of model complexity. Cross-validation may provide a more practical metric for eclipse mapping model comparison, such as the leave-one-out cross-validation applied by <cit.>, which is asymptotically equivalent to the AIC but provides advantages such as an absolute measurement of model predictive power <cit.>.






§ COMPARISON TO GENERAL CIRCULATION MODELS


In this section, we interpret the eclipse maps by comparing them to three-dimensional General Circulation Models (GCMs). We use some of the GCMs presented in <cit.>, selecting those that matched the out-of-eclipse phase curve well. We include a simulation from four models: THOR (; Simulation 31 in ), expeRT/MITgcm (; Simulation 20 in ), the Generic Planetary Climate Model (PCM) (; Simulation 7 in ), and the RM-GCM (; Simulation 24 in ). 

The THOR and RM-GCM simulations use semi-grey radiative transfer, while the PCM and expeRT/MITgcm simulations use multi-band correlated-k schemes. The THOR, RM-GCM, and PCM simulations feature clouds (on the nightside only for THOR), whereas the expeRT/MITgcm simulation is free of clouds. Output from each model was post-processed using several multi-band radiative transfer codes to calculate spectrally resolved emission from each column, which was then integrated from 5 to 10.5 μm, weighted by the MIRI/LRS throughput. Appendix <ref> and <cit.> give more detail on each model.



 §.§ Observable Features of GCMs


Figure <ref> shows thermal emission maps post-processed from each GCM simulation. The top row shows the thermal emission integrated from 5 to 10.5 μm. These maps contain small-scale structure that is not present in the observed eclipse maps in Section <ref>. A chevron-like structure is present in all the GCM maps, which is commonly attributed to the temperature structure associated with equatorial Kelvin and Rossby waves <cit.>. THOR, expeRT/MITgcm, and the RM-GCM each show cold `lobes' on the nightside, which are likely associated with stationary equatorial Rossby waves <cit.>. A Rossby wave-like structure also appears to dominate the dayside emission in the RM-GCM simulation. 








  §.§.§ Observable and Null Maps


Not all spatial brightness patterns produce a signal in phase curve space <cit.>. For instance, rotational phase curves of spatially unresolved objects (i.e., exoplanet phase curves) are insensitive to latitudinal structures. If the object is eclipsed, as is the case for WASP-43b, then the shape of eclipse ingress and egress breaks many of these degeneracies, but still leaves a “null space” of unobservable patterns <cit.>. This null space means that retrieved eclipse maps will not always match thermal emission output from GCM simulations. This effect may account for some of the discrepancy between the GCM maps shown in the top row of Figure <ref> and the observed eclipse maps in Section <ref>. 

To illustrate this effect, we separate our GCM maps into their observable and null components following <cit.>. The second row of Figure <ref> shows the observable component for each GCM. The most significant difference between the observable GCM maps and the original GCM maps is the lack of spatial variation on the nightside of the observable maps. In this region, we only have information from the phase curve, so our mapping capabilities are limited to large-scale longitudinal variation. This means that the fine structure associated with, for example, nightside equatorial Rossby waves, cannot be constrained by observations. By contrast, the dayside is scanned by the eclipse, which means that some features remain observable, such as the broad shape and location of the hot-spot. While the observable maps appear to have some new structures compared to the original GCM maps, this is simply the result of removing the high-order structures of the null space. 

There are a number of differences between the observable parts of the GCMs, and the observed eclipse maps. Each observable GCM map has a large hot-spot that is offset eastwards from the substellar point, while the eclipse maps <ref> have smaller or negligible hot-spot sifts. In addition, the dayside emission in the GCMs varies less with latitude near the equator than the emission in any of the ℓ_ max=2 eclipse maps (Figures <ref> and <ref>), or the eclipse map derived using the eigenmapping method (bottom row of Figure <ref>).

None of the GCM maps in the top row of Figure <ref> show a significant latitudinal hot-spot offset, which is present in the observed ℓ_ max=2 full phase curve map, the ℓ_ max=2 map derived from the first eclipse only, and the eigenmap. It is unsurprising that the GCMs do not show a latitudinal offset, given that the forcing and boundary conditions for each model are hemispherically symmetric. However, the observable maps for each do have small latitudinal offsets due to asymmetries introduced by the inclined viewing angle. This implies that asymmetries could be introduced by mapping at high orders.





For the purposes of comparison to the ℓ_ max=2 and ℓ_ max=3 maps we fit to the observations, it is important to note that ℓ_ max=2 and ℓ_ max=3 maps have no null space for this observation. The non-zero null spaces identified in Figure <ref> are due to the high-order spherical harmonic bases used to represent the small-scale structure of the GCM results. Therefore, while the presence of a null space places a theoretical upper limit on the accuracy of eclipse mapping it should not have a direct effect on our fitted ℓ_ max=2 and ℓ_ max=3 maps, and should not introduce a latitudinal asymmetry at these low orders. The eigenmap basis may have some latitudinal asymmetry introduced by its calculation of the structures which contribute most strongly to the light curve.



  §.§.§ Effect of Mapping Order


The eclipse maps in Section <ref> are fitted with low-order spherical harmonics up to ℓ_ max=2 and ℓ_ max=3, as including higher-order harmonics could lead to overfitting. This means that some realistic structures cannot be fitted accurately — for example, a sharp brightness temperature gradient caused by the onset of cloud formation needs high-order harmonics to fit it accurately <cit.>. To illustrate this effect, the third and fourth rows of Figure <ref> show emission maps from each GCM, truncated to use information from harmonics up to ℓ_ max=2 and ℓ_ max=3 only (computed using ). 

Latitudinal and longitudinal cross-sections of the truncated maps (taken at the sub-stellar longitude, and along the equator, respectively) are shown in Figure <ref>, where they are compared to the relevant observed eclipse maps. The eigenmaps fitted in Figure <ref> are also compared to each GCM expressed using a basis of the eigenmaps from this particular fit. 

Truncation of the GCM emission to ℓ_ max=2 forces the dayside emission to be more strongly peaked on the equator than in the higher-order representations. This difference is consistent with the different latitudinal structure suggested by the ℓ_ max=2 and ℓ_ max=3 eclipse maps from Section <ref> (comparing, e.g., Figures <ref> and the top row of Figure <ref>). The ℓ_ max=2 fit is therefore consistent with the real WASP-43b either having i) a peaked latitudinal structure in reality, or ii) a flatter latitudinal structure in reality (like the GCMs), which is being masked by the spherical harmonic truncation. While the flatter fits in the ℓ_ max=3 map in Figure <ref> are in better agreement with the GCM simulations, the ℓ=3 modes are not constrained well enough to conclude that the flatter structure is a better representation of the “real” emission from WASP-43b. 

Turning to longitudinal structure, Figure <ref> shows that restricting the emission to ℓ_ max=2 broadens the longitudinal structure of the hot-spot, most notably for the THOR and expeRT/MITgcm simulations. It also removes sharp gradients in emission at the terminators of the THOR and RM-GCM simulations, which are associated with the formation of clouds on the nightside. As with the latitudinal structure, these differences between ℓ_ max=2 and ℓ_ max=3 are consistent with those present in the eclipse maps presented in Section <ref> and Figure <ref>.







 §.§ Longitudinal Offsets


The eastwards hot-spot shift of the temperature structure of tidally locked planets is one of their key observable features, traditionally derived from the phase offset of the maximum of their phase curve. This phase curve offset is not, however, identical to the actual shift of the 2D emission structure. Instead, it represents the integration of the emission over the observed hemisphere, weighted by viewing angle. Figure <ref> shows how the simulated phase curves for each GCM have different phase curve offsets than the offsets of their emission map, which we define as the maximum longitude θ_max of ∫^+π/2_-π/2 F_p(θ,ϕ) cosϕ dϕ. We refer to this averaging as “meridional averaging” from now on, and always consider this to be the “longitudinal offset” of a 2D map, for comparison to the “peak offset” of a phase curve.

By this measure, the THOR GCM has a phase curve offset of 11.1^∘ but an emission map offset of 17.5^∘ (meridionally averaged). Figure <ref> shows how this difference is especially pronounced for this simulation, due to its sharp gradients shown in the top row of Figure <ref>. These gradients are poorly captured by the ℓ_ max=2 representation, so the phase curve has a very different offset, as it is dominated by low-order modes <cit.>. This also results in different emission map offsets in its ℓ_ max=2 and ℓ_ max=3 representations, which are 13.2^∘ and 26.5^∘ respectively. The ℓ_ max=2 emission map offset is similar to the phase curve offset, as expected. This shows how low-order representations of emission maps can bias measurements of the large-scale offsets, as well as discard small-scale information. 

Figure <ref> shows that the meridionally averaged emission map offset for the ℓ_ max=2 map in Figure <ref> is (7.75^+0.36_-0.36)^∘. This is consistent with the emission map offset for the first, second, and combined eclipses in Figure <ref>, which are ( 7.03 ^+ 7.57  _- 5.41 )^∘, ( 6.67 ^+ 5.77  _- 4.32 )^∘, and ( 4.86 ^+ 3.6  _- 2.88 )^∘ respectively. The fit to the full phase curve is more precise due to the increased longitudinal information present in the out-of-eclipse phase curve, which constrains the low-order longitudinal structure very well and break degeneracies with the latitudinal structure <cit.>. The ℓ_ max=3 map fit to the whole dataset shown in Figure <ref> has a map offset of ( 0.50 ^+ 14.79  _- 8.04 )^∘ that is consistent with the ℓ_ max=3 eclipse-only fits, as well as with the ℓ_ max=2 map fits. As discussed in Section <ref>, this map has a hot-spot offset with a median position on the substellar point, but a non-zero phase offset of around 7 degrees, consistent with the overall phase curve, showing the difference between these two metrics. The ℓ_ max=3 fit also has much more uncertainty due to its increased degrees of mapping freedom, resulting in the large range of fitted maps in Figure <ref>.

The longitudinal offsets of the n=2 Fourier series fit, and the ℓ_ max=2 map fit to the full dataset, appear very precise compared to the ℓ_ max=3 map fit to the full dataset. This reflects the choice of fitting functions rather than a truly increased precision. The out-of-eclipse phase curve in Figure <ref> is very well constrained, so the n=2 Fourier series model has a very precise peak as shown in Figure <ref>. Similarly, most of the ℓ_ max=2 modes in the map in Figure <ref> are tightly constrained by the out-of-eclipse phase curve, resulting in an apparently more precise value in Figure <ref>. This precision is due to the limited fitting functions rather than actual statistical certainty. This is shown by the greater uncertainty on the ℓ_ max=3 map in Figure <ref>, which has more degrees of mapping freedom. It can be consistent with the out-of-eclipse phase curve with many different maps, resulting in the larger uncertainty on its meridionally averaged longitudinal shift in Figure <ref>.

In summary, the derived emission map offset varies between around 0 and 20 degrees based on the mapping method used. The map offsets in the GCMs used for comparison are larger, varying between around 10 and 30 degrees depending on the mapping basis. We conclude that while the ℓ_ max=2 map and eigenmap are the best constrained maps, they may be discarding important structures from the real map. The ℓ_ max=3 map has access to more structures, but has too large an uncertainty to provide precise conclusions. Observations of more eclipses would provide the information needed to better constrain the ℓ=3 modes.




§ CONCLUSIONS


In this study, we have presented eclipse maps of the thermal emission of WASP-43b derived from a JWST MIRI/LRS phase curve. We fitted a map using ℓ_ max=2 spherical harmonics simultaneously with the parameters of the system and the systematic parameters of the instrument. There is a clear residual mapping signal in the ingress and egress of the eclipses, which a 2D map model fits much better than a Fourier series model. This residual signal is at most 450 ppm, composed of a signal ∼250 ppm due to longitudinal structure and ∼200 ppm due to latitudinal structure. This is the first time that the magnitude of the signal of the latitudinal structure has been statistically identified in an eclipse map of an exoplanet.

Fitting this dataset with an eclipse map model derived statistically significantly different parameters to those derived with a Fourier series model as in <cit.>. We suggest these updated parameters are more accurate due to a more accurate fit to the eclipse shapes, but note that they have little effect on the eclipse map itself, as shown in Figure <ref>. There were no significant degeneracies between the fitted parameters.

Figure <ref> shows that the ℓ_ max=2 map found a small eastward hot-spot shift of (7.75^+0.36_-0.36)^∘ (defined using a meridional average), compared to an eastward phase shift of (7.3^+0.4_-0.4)^∘ derived from fitting a Fourier series to the phase curve. This ℓ_ max=2 fit finds a sharply peaked latitudinal structure, with a small latitudinal offset of ( -10.72 ^+ 4.14  _- 4.68 )^∘.

Figure <ref> shows an ℓ_ max=3 eclipse map fitted with the non-mapping parameters fixed to those derived from the ℓ_ max=2 map fit. This found a smaller but more uncertain hot-spot shift of ( 0.5 ^+ 14.79  _- 8.04 )^∘ degrees east, and a flatter latitudinal structure near the equator than the ℓ_ max=2 map. It achieved a better χ^2 than the ℓ_ max=2 model, but a worse BIC and AIC due to its higher number of parameters.

In addition, we fitted an “eigenmap” (with the orbital and systematic parameters fixed), which fitted the data with six eigenmaps with orthogonal phase curves from an ℓ_ max=2 basis. Figure <ref> shows that this derived a map very similar to the original ℓ_ max=2  eclipse map in Figure <ref>, but with fewer parameters as it discarded the mapping structures that contribute least to the phase curve. Table <ref> shows how this model achieved the best BIC score, with a similar fit quality but fewer parameters.

As described in Section <ref>, a simple estimate of brightness temperature from the observed broadband flux corresponds to a temperature at the substellar point of ( 1790.0 ^+ 23.0  _- 29.0 ) K for the ℓ_ max=2 map. This is consistent with the brightness temperature of the eigenmap at the substellar point of ( 1783.0 ^+ 51.0  _- 70.0 ) K, as well as with the brightness temperature of the ℓ_ max=3 map which is ( 1822.0 ^+ 75.0  _- 92.0 ) K at the substellar point.


We tested the mapping information in each eclipse separately by re-fitting a map to each eclipse individually. Figure <ref> shows that the second eclipse produced a dayside map similar to that derived from the full phase curve, but the first eclipse produced one with a larger latitudinal offset. This may be due to increased systematic effects on the first eclipse; it could also be a real difference in thermal emission structure but we suggest that such large-scale variability over one orbit is unlikely. We recommend that those looking to use MIRI/LRS for eclipse mapping observations allow for significant settling time at the start of their observations to reduce the impact of the large instrumental systematics.

Figures <ref> and <ref> compare the fitted eclipse maps to four GCM simulations from <cit.>. In general, the GCMs match the eclipse maps well given the effect of truncating the spherical harmonic order of the maps. The ℓ_ max=2 map is largely consistent with the ℓ_ max=2 GCM representations. The ℓ_ max=3 map has almost no hot-spot shift, different to the GCMs. It has a flatter latitudinal structure near the equator, which appears more consistent with the structure of the GCMs, but the higher uncertainty on the ℓ_ max=3 map makes it difficult to compare exactly.

We conclude that there is a strong eclipse mapping signal in this MIRI/LRS observation of WASP-43b, which can strongly constrain the ℓ_ max=2 spherical harmonic components of the planetary emission. Our fiducial map is the ℓ_ max=2 eclipse map in Figure <ref>, fitted simultaneously with the orbital, stellar, and systematic parameters. This fiducial map finds a (meridionally averaged) hot-spot shift of (7.75^+0.36_-0.36)^∘  eastward (shown in Figure <ref>), compared to a shift in the peak of the phase curve of (7.3^+0.4_-0.4)^∘ eastward.

We also highlight our statistically preferred map (according to the BIC and AIC), which is the eigenmap fitted using ℓ_ max=3 and N_E=6 basis maps; this produces a very similar map to the fiducial ℓ_ max=2 map using fewer fitting parameters. However, the limited spatial order of the fitting functions means we may be missing important structure as shown in Figure <ref>. The limitation of our maps in general to ℓ_ max=2 structures is a very model-dependent result; there may be information about smaller-scale structures in the data that the spherical harmonics are not well suited to capture. Future studies could investigate mapping methods with more degrees of freedom to relax the dependence of our conclusions on the form of the low-order spherical harmonics <cit.>. Measuring further eclipses would allow better constraints on the smaller-scale structures. In general, we suggest that observing multiple eclipses, and ideally observing a full phase curve, is necessary to obtain reliable eclipse maps for even the best targets like WASP-43b. 

In summary, the dayside eclipse maps fitted to this dataset are generally symmetric about the equator, have a small hot-spot shift eastward from the substellar point, and vary smoothly away from the substellar point. This structure implies a weak eastward heat transport by atmospheric dynamics <cit.>, no strong latitudinal asymmetries driven by magnetic fields or atmospheric dynamics <cit.>, and no dayside homogeneity driven by magnetic fields <cit.>. They are generally consistent with our four numerical simulations of WASP-43b, which have small eastwards hemispheric hot-spot shifts, and are symmetric about the equator. However, the eastward hot-spot shift measured by the eclipse maps is generally smaller than that predicted by the simulations. This dataset contains spectroscopic information that we averaged out in this broadband analysis; future studies could derive maps in particular wavelength bands or regions of molecular absorption, in order to measure three-dimensional temperature structure or the distribution of particular chemical species. However, the accompanying decrease in precision may weaken the mapping signal and the accuracy of the resulting maps. More precise mapping may also reveal subtler spatial features that are not present in this analysis, and mapping of planets with different properties may reveal different spatial structures on their daysides.


 
This work is based on observations made with the NASA/ESA/CSA JWST. The data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. The specific observations analyzed are associated with program JWST-ERS-01366 and can be accessed via [DOI: 10.17909/kj9a-8d81]https://doi.org/10.17909/kj9a-8d81 Support for this program was provided by NASA through a grant from the Space Telescope Science Institute. The results reported herein benefited during the design phase from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate.

M.H. acknowledges support from Christ Church, University of Oxford.
L.T. acknowledges access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
T.J.B. acknowledges funding support from the NASA Next Generation Space Telescope Flight Investigations program (now JWST) via WBS 411672.07.05.05.03.02. 
 








§ GCMS AND POST-PROCESSING


We compare four GCM simulations from <cit.> to our derived eclipse maps. The first is THOR <cit.> (Simulation 31 in ), with the same model configuration as was used to simulate WASP-43b previously <cit.>. The simulation presented here uses semi-grey radiative transfer, and a simple parameterized cloud scheme on the nightside of the planet <cit.>. It was run for roughly 9400 orbits, and the output data was averaged over the last 500 days. The emission from each column was calculated by post-processing the results with a multiwavelength radiative transfer model <cit.>, assuming 1× solar metallicity and equilibrium chemical species concentration calculated with the FastChem model <cit.>.

The second is expeRT/MITgcm (Simulation 20 in ). This model uses the dynamical core of the MITgcm <cit.> on a C32 cubed-sphere grid, coupled to a non-grey radiative transfer scheme based on petitRADTRANS <cit.>. It follows a setup used in <cit.> and <cit.> to investigate the deep dynamics of hot Jupiters, with radiative transfer as described in <cit.>. The only difference between the configuration used here and that in <cit.> is the omission of TiO and VO. Eleven frequency bins are used, with 16 k-coefficients for each. The simulation shown here has 47 vertical levels and  10× solar metallicity, and was run for 1500 days, with the results shown here averaged over the last 100 days. The spectrally resolved emission was post-processed using petitRADTRANS <cit.> and prt_phasecurve <cit.> using a spectral resolution of R=100.

The third is the Generic Planetary Climate Model () (Simulation 7 in ), which has been used to model exoplanets <cit.> and the planets of the Solar System <cit.>. The simulation presented here uses a horizontal resolution of 64×48 with 40 vertical levels. The  treats clouds as radiatively active tracers of fixed radii. The correlated-k radiative transfer uses 27 shortwave and 26 longwave bins, each with 16 k-coefficients. The model is run for 2000 orbits, and the data presented is an average of the model output over the last 100 days. The simulation presented here has 1× solar metallicity and models Mg_2SiO_4 clouds with radius 1μm.  The simulation was post-processed with the Pytmosph3R code <cit.> to calculate the spectrally resolved emission from each column, which was then integrated from 5 to 10.5 μm.

The fourth is the RM-GCM (Simulation 26 in ), which was adapted from the GCM of <cit.> by <cit.>, <cit.>, and <cit.>, and has been applied to investigations of exoplanets using semi-grey radiative transfer and aerosol scattering <cit.>. The simulation presented here has 1x solar metallicity and condensate clouds represented by aerosols as described in <cit.>, with 50 vertical levels, and was run for over 3500 orbits. It was post-processed following <cit.> and <cit.> to produce the spectrally resolved emission.










§ ORBITAL AND SYSTEMATIC PARAMETERS


Figure <ref> shows the orbital, planetary, stellar, and systematic parameters fitted with the n=2 Fourier series model described in Section <ref>, compared to the parameters fitted with the ℓ_ max=2 eclipse map model shown in Figure <ref>. We find no significant degeneracies between these parameters, or between these parameters and the Fourier series coefficient or map pixels used to fit the light curve shape as described in Section <ref> (not shown). There are some degeneracies between the parameters of the systematic model, due to inherent degeneracies between a linear slope and a long-timescale exponential ramp.

As described in Section <ref>, the eclipse map model produces a better fit to the data (the smaller χ^2 value listed in Table <ref>). This is as expected, as it should provide a more realistic model of the eclipse shape. We therefore assume that the parameters derived by the eclipse map model (coloured blue in Figure <ref>) are more accurate. The eclipse timing is very well constrained despite the possibility for degeneracy between eclipse timing and eclipse mapping <cit.>. We suggest this is due to the high precision and cadence of the observations, the independent measurement of longitudinal structure from the phase curve, and the presence of a transit in the dataset.

The Fourier series model finds a statistically significantly smaller planetary radius, higher orbital inclination, and larger semi-major axis, than the eclipse map model. The Fourier series model also finds different stellar limb darkening parameters to the eclipse mapping model, but the resulting limb darkening profile posteriors (not shown) are almost the same. This reflects degeneracies in the limb darkening model, and the limited limb darkening information available at these longer wavelengths <cit.>. The two fits in Figure <ref> also find different systematic parameters, but these result in almost exactly the same systematic model. The Fourier series model fits a weak exponential ramp with a long timescale, plus a strong decreasing linear trend. The eclipse map model fits a strong exponential ramp with a short timescale, plus a weak decreasing linear trend. These two systematic models are almost identical within observational uncertainty when the linear and exponential trends are combined.

These statistically significant differences in the parameters fitted by the two models do not produce significant differences in the resulting eclipse map. Figure <ref> shows an ℓ_ max=2 eclipse map fitted using orbital, planetary, stellar, and systematic parameters fixed to those derived using the n=2 Fourier series model (listed in Table <ref>). The resulting map is very similar to that fitted simultaneously with the additional parameters in Figure <ref>, showing that the differences in parameters in Figure <ref> do not have a significant effect on the eclipse map. We suggest that it is still best to fit the orbital and systematic parameters simultaneously with an eclipse map model in general, to derive more accurate parameters and to search for degeneracies with the fitted map.



aasjournal





























