Toward Atmospheric Retrievals of Panchromatic Light Curves: ExPLOR-ing Generalized Inversion Techniques for Transiting Exoplanets with JWST and Ariel

Conventional atmospheric retrieval codes are designed to extract information, such as chemical abundances, thermal structures, and cloud properties, from fully “reduced” spectra obtained during transit or eclipse. Reduced spectra, however, are assembled by fitting a series of simplified light curves to time-series observations, wavelength by wavelength. Thus, spectra are postprocessed summary statistics of the original data, which by definition do not encode all the available information (i.e., astrophysical signal, model covariance, and instrumental noise). Here, we explore an alternative inversion strategy where the atmospheric retrieval is performed on the light curve directly, i.e., closer to the data. This method is implemented in EXoplanet Panchromatic Light curve Observation and Retrieval (ExPLOR), a novel atmospheric retrieval code inheriting from the TauREx project. By explicitly considering time in the model, ExPLOR naturally handles transits, eclipses, phase curves, and other complex geometries for transiting exoplanets. In this paper, we have validated this new technique by inverting simulated panchromatic light curves. The model was tested on realistic simulations of a WASP-43 b-like exoplanet as observed with the James Webb Space Telescope (JWST) and Ariel telescope. By comparing our panchromatic light-curve approach against conventional spectral retrievals on mock scenarios, we have identified key breaking points in information and noise propagation when employing past literature techniques. Throughout the paper, we discuss the importance of developing “closer-to-data” approaches such as the method presented in this work, and highlight the inevitable increase in model complexity and computing requirements associated with the recent JWST revolution.


Introduction
In the field of exoplanetary atmospheres, the so-called atmospheric retrieval technique stands out as a paramount approach for extracting information from low-to mediumresolution (i.e., R < 5000) spectroscopic observations.Originally conceived to address the lack of prior knowledge and fully exploit the limited capabilities of past telescopes, this inversion technique is now considered as an unbiased method to explore the extensive and highly degenerate parameter space of exoplanet models (Madhusudhan & Seager 2009).Notably, improvements in theory (e.g., Line et al. 2013Line et al. , 2014;;Waldmann et al. 2015aWaldmann et al. , 2015b;;Madhusudhan 2018) as well as standardization of the core algorithms-facilitated by the availability of open-source software (Waldmann et al. 2015b;Lavie et al. 2017;Villanueva et al. 2018;Mollière et al. 2019;Zhang et al. 2019;Al-Refaie et al. 2021;Harrington et al. 2022)-have permitted the establishment of those sophisticated and robust tools across the exoplanet community.In parallel, benchmarking initiatives (Barstow et al. 2020) have played a crucial role in establishing some level of consensus for at least the simplest features.Nowadays, most quantitative observational studies of exo-atmospheres incorporate insights from atmospheric retrievals in one form or another (see, e.g., population studies and recent James Webb Space Telescope (JWST) works: Fisher & Heng 2018;Tsiaras et al. 2018;Pinhas et al. 2019;Welbanks et al. 2019;Min et al. 2020;Cubillos & Blecic 2021;Roudier et al. 2021;Changeat et al. 2022;Estrela et al. 2022;August et al. 2023;Coulombe et al. 2023;Edwards et al. 2023;Jiang et al. 2023;Kempton et al. 2023;Moran et al. 2023;Taylor et al. 2023;Dyrek et al. 2024;Edwards & Changeat 2024).
Many teams have engaged in developing exoplanet retrieval codes (Rengel & Adamczewski 2023).Despite those diverse efforts, there is a general trend toward standardizing the technique.Most of the codes are 1D, typically modeling a particular atmospheric region (e.g., the dayside or the terminator chord) and focusing on the inversion of eclipse or transit observations.However, most atmospheric studies today involve tidally locked exoplanets, which have thermally and chemically inhomogeneous atmospheres driven by the strong irradiation contrast experienced by their day-and nightsides-a phenomenon confirmed in theoretical studies (Showman & Guillot 2002;Choet al. 2003;Skinner et al. 2023).To characterize those processes properly, more constraining phase-resolved observations are typically necessary.A spectroscopic phase curve, however, requires significant telescope time, instrument stability, and involves more intricate analyses surpassing the capabilities of 1D models.Such difficulties could explain why only a handful of targets have been observed Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
by the Hubble Space Telescope (HST).So far, HST has only observed the phase curves of hot Jupiters, such as WASP-43 b (Stevenson et al. 2014), WASP-103 b (Kreidberg et al. 2018), WASP-18 b (Arcangeli et al. 2019), andWASP-121 b (Mikal-Evans et al. 2022).In the near future, the excellent pointing stability and moderate instrument systematics of JWST (Espinoza et al. 2023;Rigby et al. 2023) will allow us to exploit the richness of spectroscopic phase curves fully, which are now gaining popularity. 6n the HST era, initial approaches to interpret phase-resolved spectroscopic data consisted in performing 1D emission retrievals at discretized phases (Stevenson et al. 2014).However, this method falls short of maximizing the information content of those data sets and is susceptible to 3D biases (Changeat & Al-Refaie 2020).Recently, more sophisticated methodologies have emerged.For instance, Irwin et al. (2020), Yang et al. (2023) presented some of the first unified methods dedicated to spectroscopic phase-curve data.They fitted all the observed phases of the WASP-34 b HST data (Stevenson et al. 2014) using a single 2.5D model with optimal estimation (Irwin et al. 2020) and later nested sampling (Yang et al. 2023).Cubillos et al. (2021) developed a method to extract longitudinally resolved spectra, allowing more accurate 1D retrievals to be performed.Other works have successfully performed 3D atmospheric retrievals on phase-curve data (Chubb & Min 2022), highlighting the computational challenges arising from the augmented information content (e.g., transmission spectra versus panchromatic phase curves) and the need for more physically motivated prescriptions to manage the increased problem dimension.Finally, some works have adopted simplified geometry assumptions, for instance by dividing the planet into distinct atmospheric regions (Changeat & Al-Refaie 2020;Feng et al. 2020;Mikal-Evans et al. 2022), to accelerate the forward modeling and reduce the overall computational requirement.Despite their novelty, all those techniques act on transformed representations of the data (i.e., spectra).In this work, we discuss the advantages of performing atmospheric retrievals closer to the data, presenting a new code tailored to invert atmospheric information directly from panchromatic light curves of exoplanet observations.Below, we detail the reasons for conducting atmospheric retrievals at the light-curve level.Following correction of the detector images and processing the ramps, the current approach is to perform a preretrieval reduction stage.This stage involves fitting some simplified light-curve model convolved with an instrument systematic model.The goal is to remove the remaining instrument systematics and simplify the subsequent atmospheric retrieval step by eliminating the time component from the data.Despite its apparent advantage, this two-step process-i.e., performing a light-curve reduction before conducting atmospheric retrievals-can introduce many issues and biases, a few of which are noted here.
Wavelength binning.Wavelength binning at the light-curve stage is often required to boost the transit/eclipse signature and anchor the light-curve fit.Such binning dilutes the information contained in individual wavelength channels and leads to an overall loss of information content.Decorrelated fits.Typical reduction procedures are performed using wavelength-by-wavelength light-curve fits.Since each fit uses a completely independent Monte Carlo or nested sampling parameter exploration routine, this ignores the correlation between wavelength channels probing common processes.For instance, if a molecule has spectral features in two separated wavelength channels, a simultaneous fit of both channels would allow one to extract more information than performing two separate, independent, fits.As such, wavelength-by-wavelength light-curve fits leads to information loss.Noise properties.Spectra are often constructed from the mean and variance of fitted parameters (e.g., transit depth or eclipse depth) rather than their full distribution.Using those "summary statistics" as observables for retrievals is incorrect from a noise propagation perspective.Parameter conditioning.In current strategies, orbital and instrument systematics parameters are fitted separately (at the reduction stage) from the bulk planet and atmospheric parameters (at the retrieval stage).This leads to parameter conditioning and could bias our interpretations.Complex observations.Repeated observations of potentially variable atmospheres, phase-curve observations, eclipse mapping, and other complex events (e.g., multiplanet transit or planet + moon systems) cannot be naturally modeled using atmospheric retrievals codes that do not include time.Contamination.More generally, time-dependent astrophysical signals such as stellar activity or spot crossing events are difficult to disentangle when they are not simultaneously modeled with the planetary atmosphere.
Despite the obvious need for more fundamental and comprehensive analysis techniques of exo-atmospheric observations, no atmospheric retrieval strategy, except that of Yip et al. (2020), has attempted to include time (t) explicitly.Including time in atmospheric retrievals provides a natural solution to the above issues.Yip et al. (2020) pioneered atmospheric light-curve retrievals by demonstrating their relevance for transits.Specifically, they explored the impact of instrument systematics and orbital elements mismatch when fitting light-curve data from different telescopes, showing the precariousness of spectral retrievals combining data sets.
In this study, we generalize the work from Yip et al. (2020) and introduce a novel inversion approach tailored for interpreting spectroscopic observations of transiting exoplanets at the lightcurve level: panchromatic light-curve retrieval.Specifically, we leverage the TAUREX 3.1 atmospheric phase-curve model from Changeat et al. (2021), which was designed to handle phasedependent emission efficiently without undue computational scaling with the number of phases.This capability is integrated with the transit/eclipse light-curve models of PYLIGHTCURVE, resulting in a new atmospheric retrieval model EXoplanet Panchromatic Light curve Observation and Retrieval (EXPLOR).EXPLOR introduces time as a fundamental component of the atmospheric retrieval, allowing us for the first time to generalize exo-atmospheric fits to spectro-temporal flux data.In Section 2, we provide an overview of the EXPLOR code and outline the methodology employed in this study.In Section 3, we show examples of panchromatic light-curve retrievals with EXPLOR, with a focus on a WASP-43 b-like exoplanet.Section 4 discusses the advantages of our approach, emphasizing on a scenario where the information is diluted by the conventional reduction + retrieval strategy.Finally, in Section 5, we draw the conclusions of this work and highlight some potential next steps.

Model
For this project, we have developed a new model, EXPLOR, to invert exoplanet atmospheric properties from panchromatic light-curve data.The code employs a similar class structure to TAUREX 3.1 7 (Al-Refaie et al. 2021, 2022a).It combines the recently developed 1.5D phase-curve plugin (Changeat & Al-Refaie 2020;Changeat et al. 2021;Changeat 2022) with the open-source light-curve package PYLIGHTCURVE8 (Tsiaras et al. 2016a).By combining the capabilities of those two libraries, EXPLOR considers time (t) as a new component in the atmospheric model, enabling full Bayesian retrieval of panchromatic light-curve observations (including events such as limb-darkened transits and eclipses).The role of each components of the code is described below.
Modeling the atmospheric flux along the orbit.The planetary emission is modeled along t using the TAUREX 3 1.5D phasecurve model (Changeat et al. 2021).This model simulates an atmosphere as three separated regions (hot spot, dayside, and nightside) of homogeneous properties and computes the emission contribution of each region at each phase using a quadrature integration scheme.Interested readers can refer to Changeat & Al-Refaie (2020), Changeat et al. (2021), and Changeat (2022), but we here summarize the main features of the model.As described in the aforementioned works, each region has distinctive properties-thermal structure, chemical profiles, and cloud properties-which can be set separately or coupled in a flexible manner.While relatively simple compared to the complex 3D structures of exo-atmospheres (Choet al. 2003;Kane 2022;Skinner & Cho 2022;Skinner et al. 2023), this approach has many advantages.In terms of computational requirements, it scales broadly linearly with the number of regions (here fixed to three), and has no scaling with planetary phase (i.e., time), making it ideal for light-curve modeling.By design, the model already includes the main phase-curve observables (i.e., hot-spot offset and hot-spot shape).Additionally, since the model is built using the TAUREX 3.1 framework, EXPLOR is automatically compatible with all the TAUREX 3 plugins.EXPLOR can therefore use GPUaccelerated emission and transmission (Al-Refaie et al. 2022a), chemical codes (ACE, GGCHEM, FASTCHEM, and FRECKLL; Agúndez et al. 2012;Stock et al. 2018;Woitke et al. 2018;Al-Refaie et al. 2022b), cloud models (YUNMA; Ma et al. 2023), stellar activity models (ASTERA; Thompson et al. 2024), and optimizers (MULTINEST AND ULTRANEST; Feroz et al. 2009;Buchner 2021).EXPLOR uses an implementation of the 1.5D phase-curve model to calculate the spectroscopic thermal emission of the planet and its atmosphere along t, denoted F PC (λ, t), and the spectroscopic contribution to the planetary radius in transit, denoted Δ T (λ).Here, λ is the wavelength.The radiative transfer calculations from this module do not assume any particular spectral resolution, and are typically ran at the native resolution of the cross sections; in this work R = 15,000 with ExoMol9 cross sections (Tennyson et al. 2016).EXPLOR also accounts for most of the 3D biases arising from the 3D nature of the exoplanet (Feng et al. 2016;Caldas et al. 2019;Changeat & Al-Refaie 2020;Taylor et al. 2020;MacDonald & Lewis 2022;Nixon & Madhusudhan 2022;Zingales et al. 2022).
Modeling transits and eclipses.To construct the light-curve part of the EXPLOR model (i.e., time/phase conversion and transit and eclipse events), we employ the open-source package PYLIGHTCURVE (Tsiaras et al. 2016a).PYLIGHTCURVE allows us to imprint the transit and eclipse events to the planetary emission calculated by the 1.5D TAUREX phase-curve model.Combining the time-dependent transit light-curve model (F T ) and the normalized eclipse light-curve model ( ¢ F E ) from PYLIGHTCURVE with the time-dependent planetary emission model (F PC ), the final panchromatic light-curve model F(λ, t) is given by: , . 1 In this equation, ¢ F E is normalized from zero to one, as done in Dang et al. (2018), since the 1.5D phase-curve model already accounts for the wavelength-dependent emission during this event.To compute the transit light-curve drop at each wavelength, we provide the planet-to-star radius ratio (including atmospheric contribution) calculated with the 1.5D phasecurve model (see point above) to PYLIGHTCURVE.This means that EXPLOR automatically accounts for the changing apparent size of the planet with wavelength in both transit and eclipse events (Taylor 2022).In PYLIGHTCURVE, limb darkening is modeled by the Claret (2000) law, for which the coefficients are precomputed at each wavelength, using for instance an ATLAS grid of stellar models.Note that future instances of EXPLOR could consider retrieval strategies for the limbdarkening coefficients.Here we obtain those coefficients using an automated call to the EXOTETHYS10 package (Morello et al. 2020(Morello et al. , 2021)).
The forward model in EXPLOR, when used for a light-curve simulation containing 200 integrations and at R = 15,000 between λ ä [0.3, 50] μm, runs on the order of 10 s on a single core, meaning that it can be used for sampling with full Bayesian optimizers.

Applications to a WASP-43 b-like Planet
To validate the EXPLOR panchromatic light-curve strategy, we simulate examples of a WASP-43 b-like exoplanet, as observed by NASA/ESA/CSA's JWST (Greene et al. 2016) and ESA's Ariel telescope (Tinetti et al. 2021).Our overall methodology is similar to Changeat et al. (2019): (1) EXPLOR is used in forward model mode to obtain a high-resolution spectro-temporal map of a phase-curve observation, knowing the true input parameters; (2) the simulated map is convolved to an instrument noise instance for a WASP-43 b-like exoplanet (we consider JWST and Ariel); and (3) EXPLOR is used in retrieval mode to recover the free parameters of the model.

Creation of the Simulated Observations
In our example, the planet and star have the same properties as WASP-43 b (see Hellier et al. 2011), which is a benchmark target for both observatories (Bean et al. 2018;Edwards et al. 2019;Edwards & Tinetti 2022).Note, however, that we do not attempt to simulate the real atmosphere of WASP-43 b.For the forward model, we employ the same setup as in Changeat et al. (2021), with the thermal structure of each region being the same as in their Figure 1.More specifically, the thermal structure is parameterized by five temperature-pressure (T-p) nodes for the hot spot and dayside regions, and three T-p nodes for the nightside.The hot-spot free parameters (i.e., offset, Δ HS , and size, α HS ) are also set to the same values, respectively: Δ HS = −12°.3and α HS = 40°.0.For the chemistry, we investigate two scenarios showcasing the flexibility of EXPLOR.
Scenario 1.We couple all the regions of the model and assume a hydrogen/helium-dominated atmosphere with a solar He/H 2 ratio.Trace water and carbon dioxide are then added with respective volume mixing ratios (VMRs) of log(H 2 O) = −3 and log(CO 2 ) = −4.
Scenario 2. We simulate the atmosphere at chemical equilibrium using the TAUREX plugin ACE (Agúndez et al. 2012;Al-Refaie et al. 2022a).The metallicity (Z) and the C/O ratio in ACE are coupled between the model's regions, but the chemical profiles are not (i.e., each region has a different chemical abundance profile for each molecule).To make this case more realistic, we also include gray cloud cover on the nightside using a cloud top pressure of p = 1000 Pa.As fully opaque clouds create continuum emission, a thermal gradient is required to probe their altitude.To facilitate this, we increase the temperature at the bottom of the atmosphere compared to the previous scenario to T = 800 K.This scenario is a more challenging and realistic case.
The radiative transfer calculations include opacities from the molecules H 2 O (Polyansky et al. 2018), CO 2 (Yurchenko et al. 2020;Chubb et al. 2021), CO (Li et al. 2015), and CH 4 (Yurchenko et al. 2017).Rayleigh scattering (Cox 2015) and collision-induced absorption by H 2 -H 2 and H 2 -He pairs (Abel et al. 2011(Abel et al. , 2012;;Fletcher et al. 2018) are also considered.The forward model is ran with 200 phases of equal integration time (here 7 minutes) and at the native resolution of the cross section (i.e., R = 15,000).The constructed spectro-temporal flux maps cover two eclipses and one transit, including pre-and posteclipse baselines.This is the preferred observing phase-curve strategy for most HST and JWST approved programs, and the baseline option for Ariel phase curves (Tinetti et al. 2021).

JWST and Ariel Instrument Noise Models
As described previously, the forward model is convolved with an instrument noise model (i.e., that of JWST or Ariel) for WASP-43 b to simulate an observation.For JWST-MIRI, the noise in an observation is obtained using the EXOCTK PANDEXO tool (Batalha et al. 2017).For Ariel, the official ARIELRAD radiometric model (Mugnai et al. 2020) is employed.In the case of JWST, we focus on MIRI with the Low Resolution Spectrometer as this observation was successfully executed as part of the Early Release Science transiting exoplanet program (Bean et al. 2018).As the telescope noise models were designed with transit and eclipse simulations in mind, the default outputs are the uncertainties in the transit and eclipse depths (δT), assuming half a baseline before and after the event.For independent and normally distributed noise (i.e., ideal instrument behavior), δT is related to the noise of individual exposures (δe) by: where t 14 is the transit duration of the exoplanet and t e is the integration time of the exposure.This equation is used to obtain the noise in the final spectro-temporal flux maps.Finally, the simulated observations are fed back as input to conduct the panchromatic light-curve retrievals, using EXPLOR in retrieval mode.

Retrieval Model Setup
Scenario 1: self-retrieval.We perform a self-retrieval where the chosen retrieval model has the same assumptions as the forward model.In total, this run has 18 free parameters.We include parameters describing the planet's orbit (mid time, t mid ), the bulk planetary properties (planetary radius, R p ), the atmospheric dynamics (hot-spot offset, Δ HS ), the atmospheric thermal structure (13 temperature points at fixed pressure labeled with T), and atmospheric chemistry.The pressure values from the T-p nodes are fixed (for this case to their true values).Changeat et al. (2021) and Rowland et al. (2023) have demonstrated the relevance and computational advantage of this approach.For this scenario, the chemical parameters are the global VMR of water, log(H 2 O), and carbon dioxide, log(CO 2 ), fitted from log(VMR) ä [−12, 0].Scenario 2: uninformed retrieval.Contrary to Scenario 1, we voluntarily introduce differences between the forward and retrieval models, using an uninformed thermal profiles, to explore a more challenging case.The thermal structure is still parameterized using nodes, but in this case, we arbitrarily select the pressure nodes.We use seven fixed-pressure nodes for the hot-spot and day regions with p ä {10 6 , 10 5 , 10 4 , 1000, 100, 10, 0.1} Pa, and five nodes for the nightside with p ä {10 6 , 10 5 , 1000, 10, 0.1} Pa.Fewer nodes can be used for the nightside since less information can be extracted from the weaker planetary emission.The chemical model is also ACE, with the corresponding free parameters being the metallicity (Z) and the C/O ratio, fitted from log(Z) ä [−2, 2] and C/O ä [0.01, 2].In addition we retrieve the cloud top pressure for the nightside region.In total, this is a more extensive retrieval containing 25 free parameters.
In both scenarios, we explore the parameter space using noninformative uniform priors and the MULTINEST (Feroz et al. 2009;Buchner et al. 2014) optimizer.We use 250 live points and an evidence tolerance of 0.5.

Results
With the JWST and Ariel mock simulations, we demonstrate the capabilities of EXPLOR.This section focuses on full phasecurve retrievals.More specific examples of eclipse only are discussed in Section 4.

Scenario 1: Self-retrieval
For Scenario 1, we show the simulated white light curves for JWST/MIRI and Ariel and four example spectra at key orbital phases in Figure 1.As WASP-43 b is a benchmark target, these simulations are representative of typical phase curves for both observatories.The white light curves imprint the main dynamical properties of the atmosphere, as shown by the day-night flux contrast and the hot-spot offset.In EXPLOR, flux is a function of wavelength and time, so each data point on the white light curve has an underlying spectrum.The spectra (labeled preeclipse, quadrature, pretransit, and in transit) illustrate the spectro-temporal variations of the planet's emission.Entire panchromatic phase curves are best depicted in spectro-temporal flux maps, as shown in the top panel of Figure 2 for JWST/MIRI and in the top panel of Figure 3 for Ariel.In those maps, horizontal lines of excess/reduced flux correspond to the contribution of radiatively active sources.For instance, H 2 O induces oscillating flux dimming around λ ä {2, 3, 7} μm via its absorbing properties, while CO 2 is responsible for the much narrower absorption band at λ = 4.5 μm.Such maps are the input for EXPLOR retrievals.In our tests, the simulated observations are not binned in wavelength and are passed to EXPLOR at the native instrument resolution.While some of the JWST/MIRI detector systematics appear to be mitigated by spectral binning (Bouwman et al. 2023), noise properties are generally best handled at native resolution due to correlated noise (Espinoza et al. 2023;Holmberg & Madhusudhan 2023).In practice, the simulations have 200 time steps with the JWST/MIRI observation having 247 wavelength channels11 and the Ariel observation having 52 wavelength channels.
Performing the self-retrievals as described in Section 2, we obtain the best-fit spectro-temporal flux maps shown in the middle panels of Figure 2 (JWST/MIRI) and Figure 3 (Ariel).The residuals (bottom panels) are consistent with white noise, as expected in a self-retrieval exercise.We show the extracted T-p structures for the three simulated regions (hot spot, dayside, and nightside) in Figure A1, which are consistent with the input's true profiles.The 1.5D model in EXPLOR, combined with the highly informative panchromatic light-curve approach, allows us to infer detailed thermal structures in those WASP-43 b simulations.The chemistry is also accurately constrained as shown by the full posterior distributions in Figure A2.Note that such a high level of information extraction is made possible by bypassing the construction of reduced spectra, which is required for conventional atmospheric retrievals (see later discussion).From the posterior distributions, we also demonstrate that orbital, bulk, and atmospheric (dynamical, thermal, and chemical) information can be efficiently extracted from the panchromatic light curves, thus avoiding the parameter conditioning paradigm emphasized in Yip et al. (2020).In this WASP-43 b simulation, the observations do not probe well the fully opaque deep layers of the atmospheres, leading to larger uncertainties on the values of T surf , representing the temperatures at the surface (p = 10 bar), as expected.Overall, this example demonstrates the relevance of panchromatic lightcurve retrievals.

Scenario 2: Uninformed Retrieval
A similar exercise is reproduced considering a more complex cloudy atmosphere at chemical equilibrium.For simplicity, only the JWST/MIRI case is explored but similar results should be achievable for Ariel.As explained, the retrieval model is slightly different to the forward model (i.e., this is not a self-retrieval), making this case more realistic.
As with Scenario 1, we demonstrate that the retrieval recovers a relevant interpretation of the simulated atmosphere.The spectro-temporal flux maps for this scenario are shown in Figure A3 and the corresponding retrieved thermal structure is shown in Figure A1 (right panel).Due to the added complexity of the model (i.e., higher dimensionality), this case has a much higher computing requirement than the self-retrieval described above.This is also because the plugin ACE (Agúndez et al. 2014) is not GPU enabled, thus reducing the performance gain from the GPU-enabled phase-curve emission model.The recovered posterior distributions for this scenario are shown in Figure A4, which show the correlation between the 25 free parameters of our problem.Importantly, note that there are no true values for the T-p nodes of this retrieval since the node pressures are voluntarily made different between the forward and retrieval models.Reasonable reference values are instead provided by interpolating the temperature from the forward model at the relevant pressures.Importantly, we find that the introduction of slight modeling differences between the forward and retrieval models would not affect the interpretation of this atmosphere, which show the robustness and stability of our approach.
As opposed to more standard atmospheric retrievals on reduced spectra, fitting directly the spectral light curves marginalizes over the full set of relevant parameters, allowing us to ensure a consistent propagation of the covariance between each parameter.By interpreting the "unprocessed" data directly, our panchromatic light-curve retrieval also allows us to maximize the information content and get more precise estimates of the planetary atmosphere.This is shown by the particularly tight constraints that can be obtained on the free parameters by using this method.

Why Should We Care about Time in Spectral Retrievals?
The atmosphere of a transiting exoplanet is primarily characterized using spectroscopic time series obtained around transit and eclipse events.Those events are important because they isolate the planetary signal from the stellar signal.Time is a key element in those observations, but it is rarely employed at the retrieval stage.As previously said, the time series are usually first cleaned of instrument systematics and compressed into a planetary spectrum, using data reduction pipelines.For instance the open-source codes IRACLIS (Tsiaras et al. 2016b), CASCADE (Lahuis et al. 2020), and EUREKA!(Bell et al. 2022) transform raw HST/JWST images into a single processed spectrum by fitting a simplified light-curve model (usually convolved with some instrument systematics model).The processed spectrum, however, remains in the form of summary statistics or real-value statistics (i.e., mean, median, and variance).In other words, the raw observation is compressed to a series of Gaussian distributions with a vector of means and standard deviations (i.e., a multivariate Gaussian distribution with diagonal elements only).Here, we discuss how the astrophysical signal (i.e., the exo-atmospheric information) is propagated during those reduction steps.We illustrate this using a controlled scenario restricted to the eclipse of a WASP-43 b-like planet.We briefly describe our approach here, but more details are available in Appendix B. Ultimately, our goal is to compare retrievals performed at the light-curve level (i.e., referred as EXPLOR retrievals) against retrievals performed at the spectrum level (referred as TAUREX retrievals).Following the standard practice of the literature, we obtain an eclipse spectrum from the EXPLOR spectro-temporal flux map by fitting a simplified light-curve eclipse model with POP (Changeat et al. 2024).The resulting spectrum is then interpreted using the standard 1D retrievals of TAUREX 3 and compared against the EXPLOR retrieval.Our analysis is done for an unscattered spectro-temporal flux map, and reproduced for five independently scattered instances (i.e., the spectro-temporal flux map is normally scattered according to JWST/MIRI instrument noise).To ensure fair and realistic comparison, we produce the spectro-temporal map of a 1D atmosphere (i.e., spatial variation are not included) using our Scenario 1 slightly modified.We modify the T-p structure of Scenario 1 to the hot-spot profile everywhere on the planet.A more subtle case with spatial variations is presented in Appendix B, which shows similar conclusions.For our 1D case, the input and EXPLOR best-fit maps are shown in Figure B1 while the reduced spectra and best-fit TAUREX 3 retrievals are shown in Figure 4.The comparison demonstrates the relevance of the literature approach (i.e., conventional spectral retrievals can interpret the input data), but also illustrates the compression of information from the 2D spectro-temporal maps to the 1D spectra.The posterior distributions in Figure 5 further reinforce this point, quantifying the imperfect propagation of the astrophysical signal and the noise properties as shown by the much broader distributions.CO 2 for instance, which is consistently (i.e., in all the tested scattered instances) retrieved with EXPLOR is often missed by the TAUREX retrievals.
From a theoretical point of view, working on summary statistics has major drawbacks, which are summarized here.
1. Information dilution.As the extracted spectrum is constructed from summary statistics, it does not encode the full information content of the observations.This is shown explicitly by the posterior distributions (Figure 5), which have larger uncertainties and lower levels of detection significance (see the CO 2 posteriors) when performing spectral retrievals (TAUREX retrievals).An obvious example of such information dilution is also the time-dependent information available during eclipse, which allows us to perform eclipse mapping (see additional case in Figure B2).2. Noise propagation.The underlying assumption used to construct the observed spectra is to assume Gaussian error propagation.Observational noise and light-curve models have some non-Gaussian components that cannot be well handled by summary statistics (Yip et al. 2020).For instance, when reduction pipelines retrieve the planetary radius from each light curve independently, the posterior distribution is not always Gaussian, but often spectral retrievals only handle the mean and variance for each wavelength.Additionally, when performing data reduction, each light curve is considered as an independent measurement, meaning that their covariance is assumed to be zero.This is not expected to be true, especially when instrument systematics are present.Noise propagation along the reduction chain is imperfect in our example of scattered data, as shown by Figure 5. 3. Parameter conditioning.Since the standard retrieval techniques do not marginalize over all the parameters of the problem (i.e., they do not explicitly include instrument systematics and orbital configuration), the extracted parameters could be biased or their uncertainties underestimated.This is because the solution is conditioned over a particular instance (rather than the full probability distribution) as the fitting of the instrument and the orbital dynamics are decoupled at the reduction stage.We do not explore this point more, but the panchromatic light-curve approach should naturally avoid parameter conditioning by exploring the relevant parameter space.
In addition to solving those issues, the introduction of time to retrievals offers several advantages.1.It is more flexible.This method allows one to model timevarying phenomena naturally.By resolving the ingress and egress of the transit and eclipse, this allows eclipse mapping and could also be used to constrain patchy clouds (i.e., when the east and west terminator limbs are different).Inspection of residuals could help uncover and characterize stellar activity, spot crossing events, and chemical inhomogeneities (i.e., disequilibrium chemistry).2. It is a more natural approach to model complex observations.It can handle complex situations with minimum modeling efforts.For instance, modeling and the retrieval of transits, eclipses, and phase curves involving multiple planets of the same system, or from observations of the same planet using different techniques, are trivial.

Note on Computational Requirements
Recently, the analysis of revolutionizing JWST data has challenged the scalability and efficiency of atmospheric retrieval methods, which by including more accurate and complex processes require significantly more computational resources.An important barrier encountered by modern retrievals of the JWST era-on top of the increase in the forward model compute time-regards the wider prior space.When explored using mainstream sampling-based techniques, the prior space can easily span 30+ free parameters, reaching the "curse of dimensionality" (Buchner 2021).This is particularly the case for "more-than-1D" retrievals (i.e., like EXPLOR), which require additional parameterizations for the spatial and temporal dimensions.Additional limitations can arise when moving away from clean simulated data sets.Real observations display strong spectro-temporal instrument systematics (see, e.g., Bell et al. 2023;Kempton et al. 2023) or require improved stellar properties (i.e., accurate limb-darkening properties; Rustamkulov et al. 2023) that need to be modeled simultaneously, potentially adding several degrees of freedom to the problem.On one side, potential avenues to overcome those issues include the development of more physically motivated retrievals (i.e., with fewer degrees of freedom) based on our continuously evolving understanding of exo-atmospheres and instruments.On the other side, novel inference techniques, either based on more scalable samplingbased methods-such as dynamic nested sampling (Higson et al. 2019) or Phantom-powered nested sampling (Albert 2023) -or alternative techniques-such as variational inference (Yip et al. 2024), simulation-based inference (Gebhard et al. 2023;Vasist et al. 2023) or machine-learning-accelerated surrogate modeling (Himes et al. 2022)-need to be developed.
For reference, our EXPLOR retrieval of the unscattered eclipse case shown in Figure 5 utilizes 230.4 CPUh (about 111 mn for 61,000 evaluations) using two AMD EPYC 7543 nodes (a total of 128 cores) when a TAUREX retrieval of the derived spectrum uses 9.2 CPUh (about 5 mn for 17,000 evaluations) on the same machines.While direct comparisons are difficult since the EXPLOR run includes an additional free parameter for the midtransit time and models the planet in 1.5D, this difference clearly demonstrates the extent of modern computational challenges for retrievals.

Conclusion
We introduce an innovative retrieval method designed to infer atmospheric properties from panchromatic light curves directly.This method, embedded in the novel EXPLOR model, operates closer to the data as compared to conventional retrieval strategies acting on spectra, providing several advantages.In particular, the noise properties of the observations can be more directly integrated into the retrieval process, avoiding issues such as information dilution, loss of noise characteristics, and parameter conditioning.By treating time in the atmospheric model, this generalized approach also offers a more flexible and natural strategy for modeling observations of transiting exoplanets.In particular, EXPLOR can model transit and eclipse situations, but also phase curves, eclipse mapping, and other more complex situations.In this study, we conduct mock simulations to validate our model implementation and showcase its capabilities.Overall, this work lays the groundwork for the development of approaches that are more intimately connected to observational data, which will ultimately lead to a more efficient use of telescope data.In future investigations, we will more directly compare the performances of panchromatic light-curve retrievals with the standard approach and apply our model to real observations obtained by next-generation space-based facilities.perfect case where the input spectro-temporal flux map is not scattered according to the uncertainties and a realistic scenario where we repeat our retrievals for five normally scattered instances of the input spectro-temporal flux maps.As done previously, the inputs fed to EXPLOR are at the native instrument resolution.We illustrate our panchromatic lightcurve approach for the eclipse case in Figure B1, showing a scattered instance in the top row and an unscattered case in the bottom row.
Standard reduction + retrieval case (TAUREX retrieval).To illustrate the standard practice of the literature, we start from the same spectro-temporal maps (left column of Figure B1), but now perform a light-curve reduction step to produce an eclipse spectrum.The spectrum is then fed to a 1D spectral retrieval.We first bin the simulated observation (i.e., the spectrotemporal maps) into 60 linearly separated wavelength channels -this is typically done to increase the signal-to-noise ratio in each light curve.We then employ the POP (Changeat et al. 2024) tool to fit each light curve independently.In this fit, the orbital parameters are kept fixed (i.e., we assume those parameters can be obtained with a prior fit of the white light curve, or via external measurements), and we model the astrophysical signal with a single free parameter: the eclipse depth.We add an instrument systematics model, composed of a single vertical offset, as a free parameter.This offset parameter illustrates the addition of an independent wavelength-bywavelength component, as would be the case in a real data reduction scenario (Tsiaras et al. 2016b).However, we here note that the addition of this parameter does not change our results, probably because the EXPLOR model did not include any instrument systematics commonly found in real data (i.e., short exponential ramps and long-term quadratic ramps).The inclusion of such effects can be explored in future works.Once the light curves are fitted, we construct the eclipse spectrum by taking the median and variance of the dayside flux, as done by many reduction pipelines (Tsiaras et al. 2016b;Bell et al. 2022).The information from this reduced spectrum is then inverted with the TAUREX 3 retrieval framework using the base 1D eclipse model (Al-Refaie et al. 2021, 2022a).As with the EXPLOR retrievals, we retrieve a 1D thermal profile using the same pressure nodes (i.e., with five free parameters), the abundances of H 2 O and CO 2 , and the planetary radius.
In addition to those tests, we performed eclipse retrievals on the data from Scenario 1 where spacial inhomogeneities are present (see main text).Since the T-p profile is different for each region, traditional 1D spectral retrievals should struggle more (phase-dependent behavior is visible even when only considering the eclipse only).To account partially for this-as would be done with real data-the light-curve stage is performed using a sine phase-curve model composed of three free parameters: the nightside flux, the dayside flux, and the hot-spot phase offset.The EXPLOR retrieval is performed using independent T-p profiles for each regions.The rest of the method proceeds as described with the previous eclipse tests.Figure B2 shows the posterior distributions and retrieved profiles for this case, highlighting that (1) the chemical parameters of the 1D retrieval are not biased in this case, (2) information dilution also occurs (similar conclusion to the 1D case), and (3) the T-p profile in the eclipse retrieval converges toward a mean value between the dayside and hot-spot profiles.5 are found but the EXPLOR retrieval is allowed to perform eclipse mapping while the TAUREX retrieval is 1D.The retrieval TAUREX T-p profile can only obtain an average of the dayside and hot-spot profiles, but it is reassuring to find that no significant bias is introduced in this case.By design, the EXPLOR retrieval accounts for 3D biases in eclipse and transit (as described in, e.g., Feng et al. 2016;Caldas et al. 2019;Taylor et al. 2020;MacDonald & Lewis 2022), allowing one to extract more realistic information from the data than transitional 1D retrievals.

Figure 1 .
Figure1.Simulated light curves for a WASP-43 b-like exoplanet using our panchromatic light-curve model EXPLOR.The top panel shows the white light curve convolved with the observational noise from JWST's /Mid Infrared Instrument (MIRI), and the middle panel shows the same but convolved with the noise profile from Ariel.The bottom two panels show the underlying simulated spectra (in full resolution) at phases color marked in the white light curves (top and middle panels).Note that the top and middle panels are zoomed to illustrate the phase variations and eclipses, while the insets show the full light curve with the transit.Each data point corresponds to Δt = 7 mn integrations.By definition of the instruments, the JWST/MIRI phase curve covers λ ä [5, 12] μm, while the Ariel phase curve covers λ ä [0.5, 7.8] μm.

Figure 2 .
Figure 2. Simulated panchromatic light-curve observation of a WASP-43 b exoplanet with JWST/MIRI (top) for our Scenario 1. Middle: the corresponding best-fit solution from the EXPLOR atmospheric retrieval.Bottom: the residuals between the simulated observations and the best-fit solution, normalized by the observational uncertainties σ.The simulation is created with the same model as shown in Figure 1.

Figure 3 .
Figure 3. Simulated panchromatic light-curve observation of a WASP-43 b exoplanet with Ariel (top) for our Scenario 1. Middle: the corresponding best-fit solution from the EXPLOR atmospheric retrieval.Bottom: residuals between the simulated observations and the best-fit solution, normalized by the observational uncertainties σ.The simulation is created with the same model as shown in Figure 1.

Figure 4 .
Figure4.Observed spectra obtained when the EXPLOR simulations are reduced using POP (data points), and best-fit spectra (solid lines) from corresponding 1D TAUREX 3 retrievals.The reduction was performed on simulated eclipses of WASP-43 b that are similar to Scenario 1 (the raw data shown in FigureB1).The black case corresponds to an unscattered spectrotemporal flux map, while the colored cases are five different instances of scattered spectro-temporal flux maps.Note that the spectra are offset for better clarity.Overall the 1D TAUREX 3 retrievals are able to explain the data.

Figure 5 .
Figure 5. Corner plots (left: unscattered input; middle: combination of five independently scattered inputs) from the EXPLOR (black) and the TAUREX (blue) retrievals obtained for the 1D (i.e., spatially homogeneous atmosphere) JWST eclipse examples, and retrieved T-p structure (right) for the unscattered case.In the scattered case the final posteriors are obtained by adding the resampled weighted traces from the five different scattered instances.In the T-p plot, the solid lines and shaded regions are the median and 1σ and 3σ regions, respectively, while the dashed lines are the true values from the forward model.Widening of the posteriors and the retrieved T-p profiles for the TAUREX retrievals, even when the maps are unscattered, shows the effect of information dilution arising from the reduction process.In the scattered case, the TAUREX solutions are less stable to the noise properties, most likely due to incorrect propagation of non-Gaussian components arising during the reduction process.

Figure A2 .
Figure A2.Full corner plot of the EXPLOR panchromatic light-curve retrieval for the WASP-43 b-like case for Scenario 1.This shows simulations for JWST/MIRI (black) and Ariel (blue).The true values of this simulation are indicated by the red crosses.

Figure A3 .
Figure A3.Simulated panchromatic light-curve observation of a WASP-43 b exoplanet with JWST/MIRI (top) for our Scenario 2. Middle: corresponding best-fit solution from the EXPLOR atmospheric retrieval.Bottom: residuals between the simulated observations and the best-fit solution, normalized by the observational uncertainties σ.

Figure B1 .
Figure B1.Panchromatic light-curve retrievals (ExPLOR) for an eclipse of a WASP-43 b-like planet when the data are scattered (top row) vs. unscattered (bottom row).A simulated spectro-temporal flux map convolved with a JWST/MIRI instrument model (left) serves as input to the retrievals.The retrieved best-fit maps obtained with EXPLOR (middle) and the residual maps (right) show the ability of the model to interpret such observations.

Figure B2 .
Figure B2.Corner plots (left: unscattered input; middle: combination of five independently scattered inputs) from the EXPLOR (black) and the TAUREX (blue) retrievals obtained for the 1.5D (i.e., spatially varying atmosphere) JWST eclipse examples, and retrieved T-p structure (right) for the unscattered case.Similar resultsto the 1D case shown in Figure5are found but the EXPLOR retrieval is allowed to perform eclipse mapping while the TAUREX retrieval is 1D.The retrieval TAUREX T-p profile can only obtain an average of the dayside and hot-spot profiles, but it is reassuring to find that no significant bias is introduced in this case.By design, the EXPLOR retrieval accounts for 3D biases in eclipse and transit (as described in, e.g.,Feng et al. 2016;Caldas et al. 2019;Taylor et al. 2020;MacDonald & Lewis 2022), allowing one to extract more realistic information from the data than transitional 1D retrievals.