An analytic theory for the resolution attainable using eclipse mapping of exoplanets

We present an analytic theory for the resolution attainable via eclipse mapping of exoplanets, based on the Fourier components of the brightness distribution on the planetary disk. We find that the impact parameter determines which features can and cannot be seen, via the angle of the stellar edge relative to the axis of the orbit during the eclipse. We estimate the signal-to-noise ratio as a function of mapping resolution, and use this to determine the attainable resolution for a given star-planet system. We test this theory against numerical simulations and find good agreement; in particular, our predictions for the resolution as a function of stellar edge angle are accurate to the simulated data to within 10% over a wide range of angles. Our prediction for the number of spatial modes that can be constrained given a light curve error is similarly accurate. Finally, we give a list of exoplanets with the best expected resolution for observations with the NIRISS SOSS, NIRSpec G395H, and MIRI LRS instruments on JWST.


INTRODUCTION
Understanding the temperature distributions of exoplanets is critical to testing models of heat redistribution, which in turn grant valuable insight into the mechanisms driving heat transport and the dynamic structure of the atmosphere (Lewis & Hammond 2022).While simulations and theoretical calculations (e.g.Perez-Becker & Showman (2013)) can inform us about what the temperature distributions of exoplanetary atmospheres might look like, these results cannot be tested without observational data.
Phase curves, which are observations of phase-resolved emission from planets, provide longitudinal information about their brightness temperatures.Their day-side and night-side temperatures, and the offsets of their hottest points from the substellar point, are especially valuable information.By decomposing the phase curve into its Fourier components, the longitudinal brightness distribution of the planet can be determined (Cowan & Agol 2008); this method was successfully applied to HD 189733 b (Knutson et al. 2007).Planets with axes of rotation inclined to the line of sight provide very limited latitudinal information, in that the spherical harmonics antisymmetric about the equator will also contribute to the light-curve (Cowan et al. 2013).
Eclipse mapping, a method first used on stars and accretion disks (Horne (1985), Baptista & Bortoletto (2004), Collier Cameron (1997)), is now available as a means to study exoplanets thanks to the precision of measurements with the James Webb Space Telescope (JWST) (Rauscher et al. 2007b).As a planet is eclipsed by its star, different parts of the planet are hidden at different times.Measuring the light curve therefore provides information about the brightness ★ E-mail: alexander.boone@hertford.ox.ac.uk distribution of the planet (Majeau et al. 2012).Eclipse mapping allows asymmetries about the equator to be constrained (as seen in simulations in e.g.Menou (2020), Cho et al. (2003)), and since secondary eclipses are short compared to the planet's orbital period, the method takes less telescope time and assumes less about the potential variability of brightness features on the planetary photosphere (Jackson et al. (2019), Rauscher et al. (2007a), Cho et al. (2003)).
While eclipse mapping was attempted by Majeau et al. (2012) using Spitzer data of HD 189733 b, they could only constrain the spherical harmonics up to  = 1, i.e. the dipole moment, and three of the four modes could be constrained with only the phase curve.Using JWST, it is now possible to apply eclipse mapping to some exoplanets and achieve reasonable resolution on the day side of the planet using a relatively short amount of telescope time; Coulombe et al. (2023) were able to constrain five spherical harmonics on WASP-18 b using a segment spanning only around a third of the full orbit, although there was a low degree of confidence in the fitted latitudinal structure.In general, eclipse mapping results cannot be trusted to very high resolution due to the presence of a "null space" of spatial patterns that do not contribute to a given eclipse signal (Challener & Rauscher 2023).
Hot Jupiters are gas giants in extremely short orbits around their stars ( <∼ 10 days).They have high equilibrium temperatures because of their proximity to their host stars, and are often inflated with radii moderately larger than Jupiter's.Their high temperatures and large sizes give them excellent SNRs and thus make them ideal targets for mapping.They are expected to be tidally locked due to their close-in orbits (Rasio et al. 1996;Guillot et al. 1996), and global circulation model (GCM) simulations predict them to have large day-night temperature contrasts due to the strong stellar heating (Perez-Becker & Showman 2013).Simulations also find a strong superrotating component at the equator, which in turn results in a shift of the hot-spot eastwards and a warming of the equator towards the terminator (Menou 2020).These features are validated by observations, with large phase curve amplitudes and phase offsets being ubiquitous among hot Jupiters (Parmentier & Crossfield 2017), though GCMs often have difficulty replicating the observed values (Showman et al. 2020).
The inclusion of drag (often in the form of Lorentz drag by the planetary magnetic field on a partially ionised atmosphere (Menou 2012)) can serve to hinder heat transport and thus results in a sharper drop in temperature at the terminator (Perez-Becker & Showman 2013).Clouds may also form on the night side and the cooler western hemisphere of hot Jupiters (Parmentier et al. 2013), and thereby create a larger thermal contrast and a reversed optical phase curve offset compared to the thermal phase curve (Roman et al. 2021).Clouds, photochemistry, and other processes may also create chemical gradients which will affect the temperature profile and spectrum emitted by the atmosphere (Parmentier et al. 2018).
While phase curves can probe some of these, eclipse mapping allows hot Jupiters to be observed in greater detail, and can access features which the phase curve alone cannot.It is clear that the resolution of an eclipse map can be improved by observing targets with better SNR and longer eclipses, and it might seem obvious that a planet with a perfectly edge-on orbit will not be able to provide much latitude information due to the system being symmetric along that axis.One might also guess that eclipse observations might be more sensitive to sharp gradients (e.g. as found in Parmentier et al. (2016)) than the phase curve, as the phase curve only gives hemisphereaveraged flux.That said, it is not immediately obvious how various parameters of the planet and its orbit affect the resulting eclipse map.We therefore derive a quantitative prediction of the expected resolution of an eclipse map derived from an observation of any planet given its properties, its orbital parameters, and the precision of the observation.
In Section 2, we provide an argument for the amount of degeneracy in eclipse mapping.We then derive the light-curve amplitude for sinusoidal surface modes, and finally derive the resolution achievable using eclipse mapping for surface modes up to an order .In Section 3, we discuss the methods used in our numerical tests.In Section 4, we present the results of those numerical tests.In Section 5, we apply our theory to real planets to predict which make the best targets for eclipse mapping.Finally, in Section 6, we give our conclusions.

THEORY
Eclipse mapping works by measuring the brightness of a system as a planet is eclipsed by its star.Since different parts of the planet get eclipsed at different points in time, this yields some information about the brightness distribution on the surface of the planet (see Figure 1A).
We are looking to determine the level of degeneracy in eclipse mapping, and for an estimate of the strength of the contribution of a given surface mode (e.g. a spherical harmonic or a Fourier mode) to the eclipse light curve.We begin by defining the differential light curve : where Φ is the light curve and  is a modified time coordinate such that  ranges from −1 to 1 as the planet goes from fully eclipsed to fully visible, with  = 0 halfway between first and second (or third and fourth) contacts.The differential light curve is easy to calculate assuming   <<   (where   and   are the planetary and stellar radii respectively) as the integral along a line across the planet (the line being the edge of the star at time ).We neglect the rotation of the planet between ingress and egress for simplicity.We also assume that 1 −  >>   /  so that the angle of the stellar edge does not change over the course of the ingress/egress.

Some geometry
We define an angle  0 such that the impact parameter  = sin 0 and thus We take the planet to be a circle of radius 1 centred on the origin, and call the coordinates parallel and perpendicular to the direction of motion  and  respectively.We find the stellar edge forms a line such that: (2) where ℎ is a coordinate parallel to the stellar edge, and can take values such that |ℎ| < √ 1 −  2 (see Figure 1B).For some brightness distribution (, ), the differential light-curve is: (5)

The degeneracy
Let us take a differential light-curve of the form where   is a polynomial of order .Consider a separable brightness pattern produced from the product of a polynomial of order  in , and another of order  in : which produces the differential light-curve .The result of this product will be some order  +  polynomial in ℎ and .Examining the contribution to the light curve of one term of this polynomial, and thus contributes a polynomial of order +  to   .As +  = +, we must set  +  =  if we want the light curve produced by  to be a polynomial of the same order as .For  to regenerate the differential light curve , we must match the coefficients of  and  such that they equal the coefficients on   .When both  and  are non-zero, the resulting equations are non-linear and involve products of the coefficients of  and .Nonetheless, we have  equations and  +  =  coefficients to solve for, and so can reduce ,  to at most a finite number of solutions.
While for a given , , we can solve for the coefficients of  and , we must recall that we do not know a priori which values of  and  to take, and so any ,  pair can satisfy the equations as long as  +  = : we instead have 1 2 ( + 2)( + 1) coefficients to solve for (the number of terms     where  +  ≤ ).Looking closely, we find another  equations given to us by the fact that the ingress and egress light curves will in general be different, and in The geometry we use in our derivation.We use Cartesian coordinates rather than spherical ones as they lend themselves more easily to analytic calculations.ℎ is a coordinate defined along the stellar edge, while  is a modified time coordinate perpendicular to the stellar edge.We also show the angle of the stellar edge  0 and the angle associated with the mode .The dashed red line shows the stellar edge for  0 = ; this clearly lines up with the features of ( , ).(C) The brightness contributions from points along the stellar edge for the mode shown in (B), when it is well-aligned (top) and poorly-aligned (bottom; see Equation 18).For simplicity, we rotate the coordinates from ( , ) to (ℎ,  ).The corresponding differential light curves  (  ) are also shown; the cancellation between positive and negative sections of (ℎ, 0) when the mode is poorly-aligned result in a very small differential light-curve compared to when the mode is well-aligned.
turn if   is odd, the ingress and egress light curves differ by a sign flip, while if   is even, they are identical.As there is one  = 0 mode, we must add that too.As such, we have 1 2 ( + 2) ( + 1) coefficients to solve for with 2 + 1 equations; for  ≥ 2, there is no unique solution.Equivalently, by subtracting solutions from one another, we get a set of brightness distributions which yield light curves equal to zero -the 'null space' (Challener & Rauscher 2023).While we cannot determine the true surface brightness distribution from the light curve due to these degeneracies, we can nonetheless make progress using some reasonable assumptions.

Sinusoidal brightness patterns
To examine the observable brightness of the various modes, we will take surface modes of the form: where  and  are the horizontal and vertical wavenumbers; we can take the overall order  such that  2 =  2 +  2 (for reasons that will become clear later on), though the exact combination is arbitrary.We include a phase factor along  and  for generality, though it should not have any effect on the ultimate result.This separable form results in analogous equations between the ingress and egress light curves, so we will only have to treat one and get the other for free (with a sign change in cos 0 ).

Calculating the differential phase curve
From Equation 4, the differential light curve can be calculated as: Using the identity and collecting like terms in ℎ and , we find We will discard the first term for two reasons: firstly, as ℎ is slowlyvarying around  = 0, the dominant frequency contributed by the term is ∼ ( cos 0 −  sin 0 ), which will end up being degenerate with the lower-order terms (which we expect to both be intrinsically brighter and to be damped by a much smaller denominator).Secondly, its denominator is always of order , whereas the second term has a denominator which vanishes for a specific (, ) pair (which we will call 'well-aligned').At this well-aligned (, ) pair, the the bright and dark areas of (, ) line up with the stellar edge, and so high-frequency (orthogonal) component of the resulting signal will be large (see Figure 1C, where  ≃  0 occurs when the mode becomes well-aligned).Note that neglecting this term does require us to assume that neither of (, ) is zero (i.e. that the mode varies along both axes); otherwise, the two terms are the same to within a phase difference.Thus we are left with: Expanding this using the angle summation formula yields where we can discard the term containing cos([...] ℎ) as it vanishes under the symmetric bounds.Evaluating the bounds, we reach: where ℎ 1 = √ 1 −  2 .Finally, we make the approximation that ℎ 1 is slowly-varying compared to  for most of the ingress/egress and thus its factor modulates  by a factor that is on average 1 √ 2 , giving: We can see that the first part forms a constant prefactor, and the second a (co)sinusoidal variation in time, with phase given by the phase factors  0 and  0 (as expected, they do not affect the amplitude of the signal).As the prefactor is a sinc, when its argument is near zero,  will be order-1 for all , however it will drop as ∼ 1/ for a poorly-aligned (, ) pair.Based on this, we can define  :  =  cos,  =  sin (see Figure 1B; this is why we picked our definition for  earlier) and our phase factor   = −  0 −   0 , and rewrite our formula as This form renders the notion of well-and poorly-aligned modes very clear (see Figure 1C).When sin( −  0 ) < 1/, the mode is well-aligned and the differential light curve is order 1.When sin( −  0 ) > 1/, the amplitude of the differential phase curve drops rapidly with worsening alignment, and takes a typical value of ∼ 1/.As there are ∼  modes of order , we can also see that there will always be a well-aligned mode, which will dominate the observed signal, especially for larger .Note that the zeroes of the sinc do not necessarily have any physical correspondence, they merely arise from our choice of surface modes.The scaling and width of the maximum, however, are physical; these tell us about how bright well-aligned features will be relative to their poorly-aligned counterparts, and what angular range a feature must be in to be well-aligned.Despite the notion of well-aligned/poorly-aligned modes being linked to the null space (in that well-aligned modes will dominate the non-null space if we choose to use the minimum (, ) that can create a given light curve), it is important to point out that neither the well-aligned modes nor the poorly-aligned modes are entirely in nor out of the null space; the null space will be some linear combination of well-aligned and poorly-aligned modes, if the contributions from well-aligned modes to  may be minor corrections.
While the argument from polynomial light curves only has two distinguishable modes per order (the inverting and non-inverting modes), the sinusoidal brightness distributions have four.This is because the even and odd polynomials are different orders, whereas their equivalents here -sine and cosine -are grouped into the same order.Hence the phase factor   links four separate modes; we could take them to be   = 0 or  2 , and either add 0 or  for the egress depending on whether the curve is inverting or not.

The full light curve
Up until now, we have been dealing with the differential light curve as it is easier to calculate.We now have a form which we can integrate: Since we are only interested in the variations in Φ (we can assume we know the overall brightness of the planet from the eclipse depth), we can drop the lower bound of the integral, and, by taking the rootmean-square average, we get rid of phase information.Thus the scale of the variation in Φ is For well-aligned modes where the sinc is order-1, ⟨Φ⟩ scales as 1/, while for poorly-aligned modes, it scales as 1/ 2 .One might notice the cos( −  0 ) in the denominator and suggest that perhaps some poorly-aligned modes will have significant contributions after all, but we should remember that 0 ≤  ≤  2 and 0 ≤  0 ≤  2 : only in rare cases when  ≃ 0 and  0 ≃  2 (or vice-versa) will this factor become significant; for most impact parameters, the term will remain relatively small.We also note that the divergence here only arises due to the frequency vanishing, in which case the mode becomes degenerate with (much brighter) low-order modes, and our approximations of ℎ 1 being slowly-varying break down.In short, we should not expect any poorly-aligned modes to become significant.
So while the brightness distributions are highly degenerate, with ∼  2 coefficients constrained by only ∼  equations, we find that only ∼  modes tend to dominate the light-curve.Because of this, we expect that features produced by these well-aligned modes can be resolved fairly easily, as they yield large signals, while features produced by the poorly-aligned modes yield very weak signals and cannot be detected.We have however made approximations to reach these results, and the accuracy of our predictions will determine whether these approximations are reasonable.

The final form
We are now ready to estimate the signal-to-noise ratio expected for a given mode.The SNR is: where   is the error in , and  is the signal from the particular mode we are interested in, of order .Assuming the measurements are independent and have the same error,  2  =  2 1   where   is the number of measurements integrated together and  1 is the error per point. can also be estimated as  =   ⟨Φ⟩  , where   =    2  ;   is the typical coefficient of modes at order  in the expansion of (, ), our brightness distribution.Hence: For regular measurements separated by a period   , where Δ is the duration of the ingress or egress, and the factor of 2 comes from a full eclipse including both an ingress and egress. is the orbital period.Defining  points = /  , the number of measurements that would be made in a full orbit, and approximating Δ as We can use this approximation for Δ as long as 1 −  >>   /  , which we assumed earlier so we could take  0 to be constant.We can define  =  1 / √︁  points , the error in the flux from the planet as integrated over an entire orbit.Plugging in our value of ⟨Φ⟩ from 20 and considering only the contribution from the well-aligned modes (where  ≈  0 ), We will assume the features to have a spatial power spectrum of   ∼  − where  is the gradient of the power spectrum.Our analytic formulation is agnostic to the value of , but to calculate estimates of the mapping resolution of real planets we assume  ≈ 2 based on the feature scale in simulations of hot Jupiters (Cho et al. 2003).By taking SNR ≃ 1 at the limit of detectability,   ≃ 2 0  − , and then inverting this equation, we can calculate  max , the highest order detectable: where  0 is the phase-averaged flux from the planet and  is the error on the phase-averaged flux as integrated over an entire orbit.The factor of 2 in   arises from the day-side flux being typically substantially larger than the night side, and thus the phase-averaged flux is around half the total day-side flux.Note that this falloff is better than inverting the phase curve, where the exponent is 1 +2 (see Cowan & Agol (2008)).This is because the brightness contributed by a given patch of planetary surface drops discontinuously to zero as it passes behind the star, whereas it falls off linearly as the planet rotates.If we follow Cowan et al. (2013) and look at the effect of convolving our brightness distributions with these (multiplying the power spectra), the step (for eclipse maps) scales as the inverse of wavenumber, while the linear increase (for phase curves) scales as the inverse square of the wavenumber.While the exponent is better for eclipse mapping, phase curve inversion benefits from a longer integration time (the whole orbit compared to just the eclipse ingress/egress), so may outperform eclipse mapping when  0 / is lower.The threshold is around where Δ is the duration of the ingress/egress and  the orbital period.This is only accurate to within an order-unity factor, since the modes used here and in phase curve inversion are not directly comparable.Despite each individual mode's strong  0 dependence,  max only depends on  0 insofar as it affects the duration of the ingress/egress; as long as there is a well-aligned mode, it does not matter which mode that is.

Resolution at a given order
Consider a brightness distribution dominated by a single sharp peak which integrates to 1.The corresponding differential light curve must also share this peak.Without loss of generality, we will place it at (0, 0), the centre of the planetary disk, which also places it at  = 0 in ().Fitting this differential light curve using the well-aligned modes ( =  0 ) up to order , we have: where we are approximating the sharp peak as a Dirac delta function.
The coefficients   are the Fourier coefficients of the delta function,  0 = 1 2 ,  ≠0 = 1.Let us estimate the size of the resulting brightness peak using the second derivatives of  with respect to  and  around the origin: and similar for  with sin 0 instead of cos 0 .This sum evaluates to: The height of the peak is  + 1 2 , and approximating the peak using a Gaussian and matching the second derivative at the origin and the height of the peak, we find:

𝑟
where Similarly, we can define   as follows: For a Gaussian of standard deviation   , the full width at half maximum is   = 2 √ 2 ln 2   .Then the resolved spot covers an angular range of   = 2 arcsin(  /2).Note that this estimate for   is only accurate for a spot in the centre of the planetary disk (at the substellar point), and the angular resolution worsens for features towards the edge of the disk as the planetary surface curves away from the line of sight.In the limit of large  ( >∼ 4),   scales as ( cos 0 ) −1 with a prefactor of √ 3 .For the ideal resolution,  =  max , with  max from Equation 26.
Thus, based on the eclipse light curve alone, we cannot resolve features smaller than ∼   in longitude and ∼   in latitude.Since   diverges as  0 →  2 and   diverges as  0 → 0, we have a tradeoff between longitude and latitude resolution; unless we have an excellent SNR, we cannot expect to have good resolution in both, and picking a target with good latitude resolution will mean compromising on longitude resolution.
Unfortunately, while Equations 26 and 33 tell us that a high stellar edge angle  0 grants us both a better SNR via longer ingress/egress durations and a better resolution along the latitude direction (which cannot be constrained via phase curves), planets are roughly evenly spaced in impact parameter  = sin 0 and thus planets with good latitude resolution will be few and far between compared to those with good longitude resolution.This makes the few bright planets with high impact parameters ( >∼ 0.7) very valuable for fitting maps with good two-dimensional resolution.

Eclipse Mapping Metric
Based on our theory, we propose a simple metric which approximates (and is monotonic in) the resolution attainable using eclipse mapping.We combine Equation 26 with  = 2, Equation 32, and subsequent calculations to give: EMM is a combined resolution which estimates the diagonal extent of a resolution element.Other parametrisations for EMM(EMM  , EMM  ) are possible and may be useful for certain observations.We remind the reader that in these expressions,  0 is the flux from the planet averaged over one orbital phase, and  is the corresponding error, i.e.  =  1 / √︁  points where  1 is the error for a single data point and  points is the number of points in a complete orbit.As above,  is the impact parameter,   is the planetary radius, and  is the semi-major axis.

METHODS
We implement a system to fit a brightness distribution from a secondary eclipse light curve.Our aim is to produce a numerical fit to RMS light curve brightness (arbitrary units) Raw light curve Orthogonal component Figure 2. The RMS brightness for the  = 5 light-curves and their orthogonal components as a function of stellar edge angle.The orthogonal components are calculated by orthogonalising each curve with respect to the lightcurves with  < 5, so they are not orthogonal to each other.While the non-orthogonalised light curves have a high baseline intensity, this is contributed by low-frequency components which are degenerate with the lower-order modes.The peaks are from the high-frequency component; since this component cannot be found in the lower-order modes, it also remains in the orthogonalised curves.We also note that the peaks in the orthogonal components are several times higher than the typical background, which matches the theory which suggests that the well-aligned modes should be ∼  times brighter than the background.
a light curve that matches the analytic process described above, in order to test our analytic predictions of mapping resolution.We first generate the light curves corresponding to the well-aligned modes up to an order  (defined here as  +  = , for easier comparisons between modes), as well as the poorly-aligned modes up to the same order.There are only two modes whenever one of ,  is zero (in this situation, the phase factor controls the amplitude and the phase-shifted curves are zero), and so we borrow the phase-shifted modes from the next ,  up by adding 1 to whichever of ,  is zero, and these modes are then removed from the set of poorlyaligned modes.These should still be on the central peak of the sinc, if lower down, so they should still be reasonably bright.We orthogonalise all the modes using the Gram-Schmidt procedure, running all the well-aligned modes first, in order of increasing .The wellaligned modes do not overlap significantly so their orthogonalised versions are mostly unchanged; the orthogonalised poorly-aligned modes vanish due to the degeneracy, so we are left with a basis for the light curve, and a 'null space' of brightness distributions which contribute zero to the light curve.
We use a Markov Chain Monte-Carlo (MCMC) to fit the coefficients of the basis light curves to the input curve, and generate a posterior distribution of coefficients.To fit the poorly-aligned modes, we borrow from the 'slice mapping' method of Majeau et al. (2012) and assume that the true distribution is the smoothest possible distribution for that particular combination of well-aligned modes; we use a least squares solver to maximise the smoothness of the map.While smoothness maximisation is not generally realistic, it does have the result of removing features that are not constrained by the light curve (thus probably spurious) while retaining those which are well-constrained.We estimate the Gaussian widths using the second derivative of the brightness at the origin.The error bars are calculated using the posterior distribution from the MCMC fit.The input brightness distribution is a radially symmetric Gaussian with width 0.1   .Because the  = 6 simulations had resolutions comparable to the size of the Gaussian input, the predicted curves are calculated as √︁   (  0 ) 2 + 0.1 2 where  = ,  (see Equation 32for the definitions of   and   ).We note that the simulated data oscillate around the predicted values with frequency proportional to  ; we do not as yet have an explanation for this.
When fitting the light-curves for Figure 6, as we are not interested in the distributions of parameters that correspond to any given lightcurve (or the resulting maps), we calculated the coefficients for the modes by taking the scalar product of the input light-curve with the orthogonalised well-aligned mode light curves.
It is worth noting two things.First, that while we use the wellaligned modes for fitting the light curve, the final map's light curve will draw from some mixture of well-and poorly-aligned modes, and the same map should be reached via any method which breaks the degeneracy by optimising smoothness.Second, the map will be incorrect, as the real map will not in general maximise smoothness.Smoothness maximisation does however provide us with a way to break the degeneracy while yielding reasonable values; it mostly preserves large features, while removing usually spurious smaller oscillations.
We use RK4 integration assuming   <<   (straight line stellar edge) to generate the differential light-curve, which we then integrate with RK4 again.We do not include the planet's rotation from ingress to egress.To impose smoothness, we calculate the fitted brightness distribution on a 50 × 50 grid, and define smoothness per grid point as the sum of the squared differences between the grid point and its four neighbours.The total smoothness is then the sum over the smoothness of all the grid points in the planetary disk.Note that the smoothness score defined here is minimised when the solution is most smooth.

RESULTS
Our analytic theory makes three main predictions.Firstly, the space of brightness distributions is highly degenerate with respect to light curves.Secondly, that the SNR at a given order  goes as  −1 multiplied by the intrinsic brightness of the modes of that order.Finally, the achievable resolutions along the horizontal and vertical axes are related to the angle of the stellar edge with respect to the planet's motion as calculated in Equation 32.
When orthogonalising the light-curves, we find that for  = 4, the average RMS for the orthogonalised well-aligned modes is usually a factor of a hundred or more larger than the RMS of the brightest orthogonalised poorly-aligned mode.The reason the null space curves are not reduced all the way to zero is that each light-curve is composed of many higher frequencies (especially at the beginning and end of the ingress and egress curves where there is a discontinuity in the second derivative), and these components, while small, are different for each light-curve.Figure 2 shows the  = 5 light curves and their components orthogonal to all the  < 5 modes.While the low-frequency term dropped in Equation 13 contributes a roughly constant baseline to each mode which dominates the signal (solid lines), the high-frequency component (not found in the  < 5 modes) dominates the orthogonal components (dashed lines).We see the peaks predicted by Equation 18, and that there is always a well-aligned mode.

Original map
Well-aligned modes only After smoothing Figure 4.The method described in 3, as applied to the Gaussian map used for Figure 3.We have intentionally picked a very small hot-spot in the input map so as to approximate a delta-function, so the results can be more easily compared with our theory in Section 2.7.When the well-aligned maps are already fairly smooth (top), there is little change during the smoothing procedure.However, if the well-aligned map is less smooth, or if one of the modes selected as well-aligned is on the edge of being poorly-aligned (bottom), the smoothing procedure has a larger effect and generally does well in retrieving a reasonable result.We used  = 4 for these tests.
To demonstrate that the system described in 3 works reliably, we ran several tests giving the system known maps and seeing how well it could recreate them from the light curve.As expected, the method is able to fairly reliably retrieve features which are reasonably wellaligned with the stellar edge during the eclipse.
In Figure 4, we test our method's ability to retrieve what is effectively a delta-function map (we use  = 4, which is unable to resolve the small hot-spot), and hence is directly comparable to our theory in Section 2.7.As expected, the resulting maps have large brightness peaks at the correct position, but whose sizes along  and  depend strongly on the stellar edge angle  0 , being stretched perpendicular to the stellar edge.
Figure 5 shows the input map and maps retrieved from lightcurves generated with three different stellar edge angles (marked by the green lines).When the stellar edge is almost vertical, the light curve does not distinguish between a narrow equatorial band and a uniformly warm planet, and so the smoothing routine prefers the lessdetailed map with no equatorial band.When the stellar edge is more favourable for resolving horizontal features, the equatorial band is correctly retrieved, as is the presence of a hot-spot, though the exact position is slightly offset.For an almost horizontal stellar edge, the horizontal band is easily retrieved but the position of the hot-spot in longitude is lost.Note however that including phase curve data will give some longitude information (though as found in Equation 26compared with Cowan & Agol (2008), the phase curve is much less sensitive to small features); combining these is beyond the scope of this paper.
Figure 3 shows how we applied the method to fit Gaussian maps using different SNRs and correspondingly different values of fitting order .We find that the widths of the fitted Gaussians along the horizontal and vertical axes trace the values predicted by Equation 32 very closely, with RMS errors of 5, 3, 2 • over the | −  0 | < 60 • range for the  = 2, 4, 6 tests respectively, each corresponding to an average error of 10% compared to the predicted values.
We also ran several tests fitting randomly-generated maps where the coefficient of each mode is calculated from a standard normally- At  0 = 0.1, the horizontal band is very poorly-aligned and so is not detected; the smoothing process then prefers a smooth map with few features.At  0 = 0.9, the horizontal band can be detected, alongside variations along .Finally, at  0 = 1.4, the horizontal band is well-aligned but variations in the  direction cannot be detected and are smeared out.Around  =  4 , the smoothing cannot distinguish between horizontal and vertical features as long as they are symmetric across both axes, so turns the horizontal band into a plus shape; in real situations, these features could be distinguished from one another using the phase curve, so we picked a stellar edge angle reasonably close by ( 0 = 0.9) which does not suffer from this problem.Best fitting order  as a function of .We find  (  ) numerically by fitting a known light-curve with noise added; we calculate the value of  2 for the fit as compared to the noiseless curve. (  ) is then the value of  for which  2  is minimised.We chose this score as the mean squared error  2  2 is sensitive to overfitting but not underfitting, while  2 is sensitive to underfitting but not overfitting.Their product then is sensitive to both.From Equation 26, we expect  (  ) ∼  +1 if the brightness of the modes of order  follows  − .In this case we used  = 2 and find our best fit from  (  ) at  +1 = 3.37±0.05,reasonably close to the predicted value.The best fit gradient also drops closer to the prediction if we exclude the lower-order points.For phase curve inversion, the power law would have an exponent of  + 2 and thus the same number of data points would constrain fewer surface features.The dark and light grey bands around the best fit line represent the 1 and 3 bounds respectively.We used  0 =  4 for these runs.distributed variable multiplied by  − , where we used  = 2. From the theory (Equation 25), we expect that the value of  required to be able to justify adding another order should go as  −−1 ; in this case  −3 .Using scipy.optimize'scurve-fitting routine, we find a best fit value of  + 1 = 3.37 ± 0.05 (see Figure 6), reasonably close to the predicted value.Most of the deviation comes from low orders which are more significantly discrepant from  + 1 = 3, and fitting only the higher orders ( ≥ 10) yields  + 1 = 3.02 ± 0.01.While we use the well-aligned modes to perform our fitting and one might suggest that this would alter the results (since Equation 26 only considers the well-aligned modes), all bases with identical span would yield the same light-curves and thus the same errors and power law.

REAL PLANETS
Our analytic result in Equation 32 predicts the angular size of features resolvable by eclipse mapping of a planet given the parameters of its system and the signal to noise ratio of the observations.Given this, we can derive the expected resolution of an eclipse map of any planet with any instrument.We apply this to all known exoplanets for three instruments: MIRI LRS, NIRSpec G395H, and NIRISS SOSS.We do not include grazing orbits as they are not described by our theory, but some of these may be good targets for partial mapping.
We select the top 100 targets for each instrument using the more approximate and non-absolute eclipse mapping metric ∼     √︃ Δ Δ 0 10 0.4( − 0 ) (Mullally et al. 2019), where Δ is the duration of eclipse ingress and egress and  is the K-band magnitude of the observed star (Δ 0 and  0 are the reference values for the planet HD 209458 b).We calculate the average flux ratio in the relevant wavelength band, and calculate the cadence and the error per point in a phase curve using Pandexo (Batalha et al. 2017).
We then calculate the eclipse mapping metrics in Equation 34 for the top 100 targets selected in this way for each instrument.We estimate  0 as the average flux ratio in the relevant wavelength band, where the planetary flux is calculated from its equilibrium temperature.Figures 7, 8, and 9 show the resulting longitudinal, The dashed line in each figure shows the effect of varying the stellar edge angle on the expected mapping resolution for a planet with the SNR of HD 189733 b, where we include the duration of its eclipse in the SNR, defining it as SNR  0 .Figures 7, 8, and 9 show a skew towards better longitudinal resolution EMM  and worse latitudinal resolution EMM  , with planets clustering in the top left of the plot.This is because if planets are uniformly distributed in impact parameter , their stellar edge angle  0 = arcsin() will be distributed non-uniformly in the range 0 to 90 • , with more small values than large values, resulting in EMM  being typically larger than EMM  .
Tables 1, 2, and 3 list the top 15 targets for each instrument, ordered by the best overall resolution EMM.As would be expected, the shorter-wavelength NIRISS SOSS is better suited to mapping higher-temperature planets, while the longer-wavelength MIRI LRS is better suited to mapping lower-temperature planets.In general, an impact parameter around 0.5 to 0.8 (giving an intermediate stellar edge angle) is key for a good eclipse mapping signal.We have not included the effect of information from the out-of-eclipse phase curve when deriving an eclipse map.This provides more longitudinal information (albeit with less spatial resolution, as discussed above) than the eclipse only.Including phase curve information will therefore make planets with high  where EMM  < EMM  better targets for eclipse mapping.
Table 1 also lists the expected resolution of WASP-18 b as Coulombe et al. ( 2023) derived an eclipse map of WASP-18 b from an observed secondary eclipse with NIRISS SOSS, finding good constraints on its horizontal (longitudinal) structure but poor constraints on its vertical (latitudinal) structure.This is consistent with the analytic prediction of Table 1 of EMM  = 19.8 and EMM  = 52.6 for this planet and instrument.

CONCLUSIONS
We derived an analytic theory for the resolution achievable using eclipse mapping on an exoplanet.We find that: (i) The space of surface brightness patterns is highly degenerate in the space of eclipse light-curves; that is, there are many brightness patterns which all result in the same light-curve.Equivalently, there are very few brightness patterns which yield non-zero light-curves compared to the number of brightness patterns whose light-curves are zero (the "null space").
(ii) Brightness patterns whose features align with the edge of the star during ingress and egress yield large signals, while others experience cancellation between bright and dark regions and thus yield small signals.As such, we can constrain the well-aligned brightness patterns well, and have to use other methods (e.g.imposing non-negativity or maximising smoothness) to constrain the others.
(iii) Because of this, our resolutions in the latitude and longitude directions are affected strongly by the impact parameter: at high impact parameters, the longitudinal resolution is worsened and small features are not detectable, while the latitudinal resolution is good.Thus the choice of target planet depends strongly on what features one intends to observe: a planet with high impact parameter is ideal for measuring the width of the equatorial jet, while for measuring steep of east-west brightness temperature gradients (e.g.where cloud formation starts), one should opt for a planet with an impact parameter near to zero.
(iv) We derive an Eclipse Mapping Metric (Equation 34) which estimates the achievable resolution in latitude and longitude; we calculate this for the MIRI LRS, NIRSpec G395H, and NIRISS SOSS instruments on JWST.We find our values to be consistent with JWST observations of WASP-18 b (Coulombe et al. 2023) which constrain well features along the longitude axis, but where features in latitude were almost unconstrained.We find our theoretical predictions validated by our numerical tests, which yield the correct resolutions with an error of 10% from the predicted values, as well as the correct scaling (to a similar accuracy) for the number of resolvable modes as a function of light curve error.
Given our theory, we present a list of exoplanets with particularly good resolution for eclipse mapping.34.HD 189733 b is marked with an asterisk as it is expected to partially saturate the detector but we leave it here for reference.
Figure 1.(A) The premise of eclipse mapping.As a planet passes behind its star, different slices of the planet are hidden at different times.The change in sign of the stellar edge angle from ingress to egress (and to a lesser extent, the curvature of the star) gives 2D information about the brightness distribution.During some time d, an area of the planet d is hidden.The flux contributed by d is dΦ, equal to the integral of the brightness distribution ( , ) over the slice.(B)The geometry we use in our derivation.We use Cartesian coordinates rather than spherical ones as they lend themselves more easily to analytic calculations.ℎ is a coordinate defined along the stellar edge, while  is a modified time coordinate perpendicular to the stellar edge.We also show the angle of the stellar edge  0 and the angle associated with the mode .The dashed red line shows the stellar edge for  0 = ; this clearly lines up with the features of ( , ).(C) The brightness contributions from points along the stellar edge for the mode shown in (B), when it is well-aligned (top) and poorly-aligned (bottom; see Equation18).For simplicity, we rotate the coordinates from ( , ) to (ℎ,  ).The corresponding differential light curves  (  ) are also shown; the cancellation between positive and negative sections of (ℎ, 0) when the mode is poorly-aligned result in a very small differential light-curve compared to when the mode is well-aligned.

Figure 3 .
Figure3.Horizontal and vertical resolutions as a function of stellar edge angle  0 .We estimate the Gaussian widths using the second derivative of the brightness at the origin.The error bars are calculated using the posterior distribution from the MCMC fit.The input brightness distribution is a radially symmetric Gaussian with width 0.1   .Because the  = 6 simulations had resolutions comparable to the size of the Gaussian input, the predicted curves are calculated as √︁   (  0 ) 2 + 0.1 2 where  = ,  (see Equation 32 for the definitions of   and   ).We note that the simulated data oscillate around the predicted values with frequency proportional to  ; we do not as yet have an explanation for this.

Figure 5 .
Figure5.A randomly generated map (left) and the retrieved maps for  = 4 at three different values of  0 (right).At  0 = 0.1, the horizontal band is very poorly-aligned and so is not detected; the smoothing process then prefers a smooth map with few features.At  0 = 0.9, the horizontal band can be detected, alongside variations along .Finally, at  0 = 1.4, the horizontal band is well-aligned but variations in the  direction cannot be detected and are smeared out.Around  =  4 , the smoothing cannot distinguish between horizontal and vertical features as long as they are symmetric across both axes, so turns the horizontal band into a plus shape; in real situations, these features could be distinguished from one another using the phase curve, so we picked a stellar edge angle reasonably close by ( 0 = 0.9) which does not suffer from this problem.
Figure6.Best fitting order  as a function of .We find  (  ) numerically by fitting a known light-curve with noise added; we calculate the value of  2 for the fit as compared to the noiseless curve. (  ) is then the value of  for which  2  is minimised.We chose this score as the mean squared error  2  2 is sensitive to overfitting but not underfitting, while  2 is sensitive to underfitting but not overfitting.Their product then is sensitive to both.From Equation 26, we expect  (  ) ∼  +1 if the brightness of the modes of order  follows  − .In this case we used  = 2 and find our best fit from  (  ) at  +1 = 3.37±0.05,reasonably close to the predicted value.The best fit gradient also drops closer to the prediction if we exclude the lower-order points.For phase curve inversion, the power law would have an exponent of  + 2 and thus the same number of data points would constrain fewer surface features.The dark and light grey bands around the best fit line represent the 1 and 3 bounds respectively.We used  0 =  4 for these runs.

Figure 7 .
Figure 7.The best 100 targets for eclipse mapping with NIRISS SOSS, selected as described in Section 5.The five best targets by single-eclipse overall resolution EMM according to Equation 34 are labelled.

Figure 8 .Figure 9 .
Figure 8.The best 100 targets for eclipse mapping with NIRSpec G395H, selected as described in Section 5.The five best targets by single-eclipse overall resolution EMM according to Equation 34 are labelled.

Table 1 .
Coulombe et al. (2023)r eclipse mapping with NIRISS SOSS, ordered by their single-eclipse overall resolution EMM according to Equation34.Additionally, the predicted resolution of WASP-18 b is listed for comparison with the real eclipse map derived byCoulombe et al. (2023)with NIRISS SOSS.The targets marked with asterisks are expected to partially saturate the detector but we leave them here for reference.

Table 2 .
The best 15 targets for eclipse mapping with NIRSpec G395H, ordered by their single-eclipse overall resolution EMM according to Equation