Transiting exoplanets with the Mid-InfraRed Instrument on board the James Webb Space Telescope: From simulations to observations

The James Webb Space Telescope (JWST) has now started its exploration of exoplanetary worlds. In particular, the Mid-InfraRed Instrument (MIRI) with its Low-Resolution Spectrometer (LRS) carries out transit, eclipse, and phase-curve spectroscopy of exoplanetary atmospheres with unprecedented precision in a so far almost uncharted wavelength range. The precision and significance in the detection of molecules in exoplanetary atmospheres rely on a thorough understanding of the instrument itself and accurate data reduction methods. This paper aims to provide a clear description of the instrumental systematics that affect observations of transiting exoplanets through the use of simulations. We carried out realistic simulations of transiting-exoplanet observations with the MIRI LRS instrument that included the model of the exoplanet system, the optical path of the telescope, the MIRI detector performances, and instrumental systematics and drifts that could alter the atmospheric features we are meant to detect in the data. After introducing our pipeline, we show its performance on the transit of L168-9b, a super-Earth-sized exoplanet observed during the commissioning of the MIRI instrument. This paper provides a better understanding of the data themselves and of the best practices in terms of reduction and analysis through comparisons between simulations and real data. We show that simulations validate the current data-analysis methods. Simulations also highlight instrumental effects that impact the accuracy of our current spectral extraction techniques. These simulations are proven to be essential in the preparation of JWST observation programs and help us assess the detectability of various atmospheric and surface scenarios.


Introduction 1
The long-awaited James Webb Space Telescope (JWST) was 2 launched on 25 December 2021.Equipped with its four instru-3 ments NIRISS 1 , NIRCam 2 , NIRSpec 3 , and MIRI 4 it has now 4 started to provide its first observations in the infrared.Each in-5 strument has different modes for photometry or spectroscopy 6 and covers various regions of the spectra, from 0.6 to 28 µm.

7
In particular, the instrument that covers the longer wavelengths 8 is MIRI.

9
The demands for observations with MIRI are high, and a sig-

Simulations
To create time series of MIRI LRS spectra, we follow a fourstage process.First, we create the star-planet emission time series of 1D spectra with the exoNoodle package (Martin-Lagarde et al. 2020).Then, we use MIRISim (Klaassen et al. 2020) to convert the astrophysical signal into detector spectral images by considering both the telescope optical path and the MIRI instrument transmission coefficients.At this stage, we also simulate the detector behaviour and add effects at the pixel scale, but over only one non-destructive integration.However, transitingexoplanet simulations require another specific stage.The search for very faint flux variations requires considering faint detector persistence effects that may have an impact on the characterisation of the atmosphere.To include these features into the simulations, we created the MIRISim-TSO tool, which adds lowfrequency detector persistence effects to time-series observations (TSO).When our simulations are complete, we proceeded with the data reduction steps.

Star-planet system time series: exoNoodle
The first step is to simulate an astronomical scene which is the incoming light from a star-planet system as the exoplanet orbits its host star.Light-curve simulation codes are quite numerous within the community.Mostly based on the Mandel & Agol (2002) and the Giménez (2006) analytic light-curve and limbdarkening models, some tools compute exoplanetary transit light curves, such as PyTransit (Parviainen 2015), the TRIP module of ExoTETHyS (Morello et al. 2020), and even PyPplusS, which computes transiting exoplanets (Rein & Ofir 2019).Other packages are designed to perform fits of exoplanet transits and radial velocity variations, including TAP (Gazak et al. 2012), EXOFAST (Eastman et al. 2013), and JKTEBOP (Popper & Etzel 1981;Southworth et al. 2004), which was created to fit light curves of eclipsing binary stars.More recently, new Python frameworks have been released and offer comprehensive models and fitting toolkits, such as exoplanet (Foreman-Mackey et al. 2021) and starry (Luger et al. 2019).Both batman (Kreidberg 2015) and PyLightcurve (Tsiaras et al. 2016) provide exoplanet transit and occultation light curves, while SPIDERMAN (Louden & Kreidberg 2018) produces phase curves in addition to occultations.
Light-curve observations made with the aim to characterise atmospheres are highly sensitive to the limb-darkening effect, and all codes pay particular attention to this.We mention here specific limb-darkening computation codes that can be added to the long list given above: ExoTIC-LD (Laginja & Wakeford 2020), exoCTK (Bourque et al. 2021), ExoTETHyS.SAIL (Morello et al. 99 2020), and Limbdark.jl(Agol et al. 2020).
100 Each one of these simulation codes brings a different per-101 spective to the light-curve modelling: They either compute a 102 transit, an occultation, or a phase curve based on the star-103 planet emission and atmospheric transmission.However, most 104 of these codes simulate normalised light curves rather than abso-105 lute fluxes.In keeping with these models, Martin-Lagarde et al. 106 (2020) developed the exoNoodle 5 .This is a Python tool that 107 generates time series of spectra in absolute flux as the exoplanet 108 orbits the star.The implemented model is based on the Mandel 109 & Agol (2002) prescription.Simulating spectra over the orbital 110 phase allows us to prepare any MIRI LRS time-series observa-111 tions in absolute flux.

112
To create these simulations, exoNoodle requires several in-113 puts, which are listed below.3. The atmosphere transmission spectrum (i.e. in a transit ge-117 ometry) in units of (R p /R ) 2 .
118 4. The limb-darkening coefficients either as a quadratic or as 119 a four-coefficient law.The limb-darkening coefficients may 120 optionally depend on the wavelength.
121 5.The orbital parameters: the orbital period P, the semi-major-122 axis a, the distance to the star d, the inclination of the orbit 123 i, and the stellar and planetary masses and radii, M p , R p and 124 M , R , respectively.125 6. Phase values corresponding to the start and end of observa-126 tion (between -1 and 1).7. The sampling time that corresponds to the interval between 128 two spectrum computations in the time series.In the simu-129 lations presented in this paper, we decided to use the MIRI 130 detector integration time.We calculated the integration time 131 knowing the detector saturation level, the brightness of the 132 target, and the LRS slitless frame time, taken from Ressler 133 et al. (2015).We discuss the MIRI LRS readout pattern and 134 terminology in more detail in Sect.2.2.
135 8.A constant wavelength bin size ∆λ based on the MIRI LRS 136 wavelength dispersion model from 5 to 12 µm and the spec-137 tral resolution model R = λ ∆λ (Kendrew et al. 2015). 138 As an output, exoNoodle creates a time series of spectra in 139 µJy.The number of spectra we compute is based on the duration 140 of the time-series observation and the sampling we choose.After we modelled the time series of spectra with exoNoodle, 143 we created spectral images of the observed scene with an instru-144 ment simulator.There are several tools that simulate detector 145 images or spectra with a MIRI-like signal-to-noise ratio, such 146 as the PandExo tool from Batalha et al. (2017) and the Expo-147 sure Time Calculator for JWST (Pontoppidan et al. 2016) where P(λ) is given by s −1 m −2 µm −1 ), subtended by the telescope entrance pupil area where n is the number of pixels that form the diffraction pattern.

188
The pixel input flux F pixel,n (λ) is then converted into an elec-189 tronic signal S pixel,n (λ), and both quantities are related through where E n (λ) is given by S pixel,n (λ) is obtained by integrating photons during a given interval of time t frame (in s), within a range of wavelengths ∆λ (in µm).

194
The QE(λ) tr d (λ) product is called the photon-electron conver- According to Eq. 5, photons falling into a pixel are converted into electrons and are then read out by the detector proximity electronics to be converted into DN.A readout pattern is a complete scheme of integrating light while doing multiple readings of a detector subarray.In the case of the MIRI detector, this scheme is called the FASTR1 non-destructive readout: The detector is read out non-destructively at regular intervals while integrating light, until it reaches a signal level close to saturation.After integrating light non-destructively, two resets are performed.This description is displayed in Fig. 1.The regular time interval between two readouts is called the frame time, which is 0.159s for the LRS slitless subarray (Kendrew et al. 2015).This whole readout pattern composed of frames and two resets is called an integration, or a ramp.The integration time and therefore the total number of frames within an integration is chosen based on two parameters: the brightness of the observed target in e − s −1 and the detector pixel saturation level in e − .The whole observation is composed of several integrations and is called an exposure.
Although resets are performed to empty the pixel potential wells, there is a remaining number of electrons that causes an offset at the beginning of each ramp.Originally, only one reset was performed within the readout pattern (formerly called the FAST mode).Then, to minimise detector systematics linked to the reset step, a second reset was added to the readout pattern.As a consequence, the offset value was lowered from 10 000 DN to 3000 DN (Argyriou 2021).As MIRISim was coded based on the former FAST mode, we made the required changes in the MIRISim code to reproduce the FASTR1 mode.The offset value was changed to 3000 DN, and the timings were extended to include an additional reset.The whole reset pattern composed of two resets lasts 0.159 s, which is equivalent to a frame time.

The instrument settings
Simulations were made with the following MIRISim settings.As the observation mode is the LRS slitless mode, the focal plane subarray was set to SLITLESSPRISM.The filter parameter was set to P750L, which corresponds to the double-prism assembly.The telescope parameters were fixed to a beginning-of-life configuration, which is the telescope post-launch condition as determined during commissioning.The simulations included the following instrument settings: 1.A map of bad pixels (either hot or dead pixels) to be flagged as DO-NOT-USE in the simulations.
2. A dark current map in DN s −1 pixel −1 .The dark current is a random generation of electrons through heat in the depletion layer, when no photons enter the detector.Decreasing the detector temperature is a way to limit the dark current (Glasse et al. 2015).
3. The flat-field map, that is, the relative response of pixels illuminated with a uniform source (Glasse et al. 2015).
Some effects that are newly witnessed in the data were not added to the simulations because no model has been released so far from either ground-based or in-flight testing.This is the case for the reset-switch charge decay (RSCD) (Ressler et al. 2023;Morrison et al. 2023), caused by resetting the detector.This was not added in our simulations.The detector reset is accompanied by a transient effect that impacts the first frames of the subsequent integration, creating non-linearities at the start of the ramp (Rauscher et al. 2007;Ressler et al. 2008Ressler et al. , 2015)), which are only measured when the detector is in the dark (i.e.not illuminated).
The reset anomaly is caused by trapped charges in the detector (possibly on the surface, at the indium bumps; see Rauscher et al. (2007)), even after reset and without illumination.When the detector is illuminated, the problems associated with the reset are known as the RSCD, which causes a greater increase in the ramp at the start of an integration (from the second integration onwards) and hence non-linearities.Finally, the reset causes a drop in offset at the start of the ramp, so that the starting point of the ramps would be different for the first and subsequent integrations.The addition of a second reset in the MIRIm ramps is partly linked to this problem and corrects offset errors at the start of the ramp after two resets (Argyriou 2021).
No background or noise are included at this stage of simulations.The readout noise in e − , caused by fluctuations in the readout amplifiers (McMurtry et al. 2005), is added to the simulations at this stage.

The calibration data products
MIRISim is based on the use of calibration data products (CDPs) that were initially created from ground-based testing campaigns.
The in-flight calibration of the instrument made during commissioning provided a new set of files.These files are meant to be used by the Space Telescope Science Institute (STScI) jwst reduction pipeline6 .To create up-to-date simulations, CDPs were adapted to be compatible with MIRISim.They were taken from the jwst STScI pipeline calibration reference data system7 (CRDS).CDPs of MIRI LRS taken from commissioning include non-linearity coefficients, dark map, readnoise map, bad-pixel mask, flat-field map, absolute flux calibration coefficients, spectral dispersion coefficients, gain value, and pixel area value.Only the monochromatic PSF file comes from former optical modelling and ground-based calibration because no in-flight calibration file is available.Commissioning revealed that the electronic gain value that converts electrons into DN is lower than expected.The simulations therefore include a gain of 3.1 e − DN −1 instead of 5.5 e − DN −1 (private communication, S. Kendrew, 2023)8 .

Detector persistence effects: MIRISim-TSO
In this section, we discuss the impact of detector persistence effects on the flux level over the whole exposure and how we include them in our simulations.Detector effects that were previously added in our simulations such as the dark current or the readout noise are additive and only depend on the pixel location on the detector.Non-linearities also affect the signal level, but only over one integration.In the end, MIRISim applies these effects independently on each integration.However, detectors feature low-frequency drifts that evolve over several integrations and might affect the whole exposure (Teams 2021).These drifts Fig. 1.Non-destructive readout pattern for the Mid-InfraRed Instrument Low Resolution Spectrometer (called FASTR1).Two consecutive integrations are displayed here.Each integration is composed of two resets and a series of frames.The number of frames is determined with the magnitude of the source and the detector pixel saturation level.An exposure is a set of several integrations.Adapted from Ressler et al. (2015).
account for 1 to 2 % of the absolute flux, and have the same order 314 of magnitude or might even exceed the atmospheric feature am-315 plitudes that we are meant to detect in our observations.These 316 temporal drifts may have different origins.An origin might ei-317 ther be a deviation in the temperature stability of the proximity 318 electronics or a telescope-pointing deviation, known as jitter.Jit-319 ter is known to have a low impact on the MIRI detector outputs 320 as commissioning confirmed a pointing stability better than ex-321 pected of 1 mas (Rigby et al. 2023).As the pointing stability 322 is excellent and the PSF is well sampled from 5 to 12 µm and 323 above, no intra-pixel noise is measured.This intra-pixel noise 324 was known to be one of the most problematic issues for Spitzer 325 data reduction (Ingalls et al. 2012;Morello et al. 2016).In the 326 particular case of MIRI, detector persistence effects may also af-327 fect the flux stability.Some tests performed on the MIRI detector 328 test model in 2018 by the Jet Propulsion Laboratory (JPL) that 329 were presented in Martin-Lagarde (2020) show different types of 330 persistence effects over an exposure.These data consist of pho-331 tometric time series obtained by illuminating the detector with a 332 set of blackbodies at different temperatures.The use of different 333 blackbodies provided a set of light curves at different flux levels, 334 expressed in DN s −1 .The independent fit of these light curves re-335 veals that the persistence effects are flux-dependent.Using these 336 test data, we identified and modelled these persistence effects 337 and included them in our simulations.Analysis of the commis-338 sioning data indeed revealed these effects.Sect. 5 broadly dis-339 cusses their nature and the quantification of their impact on time-340 series observations.Persistence effects are induced by previous uses of the detector.343 They therefore depend on its history.These previous operations 344 tend to modify the detector behaviour at the beginning of an ex-345 posure.The direct consequence of these modifications is that the 346 output flux level is altered and no longer corresponds to its ex-347 pected value.The JPL tests of 2018 revealed three persistence 348 effects in the data: the response drift, the anneal recovery, and the 349 idle recovery (Martin-Lagarde 2020).Each one of these effects is 350 described more extensively in Appendix A. The MIRISim-TSO 351 tool was created to add these effects to the simulations9 .
All three persistence effects are described using an exponential model and come from Martin-Lagarde (2020).They are implemented in this form in MIRISim-TSO.For all effects, the time variable t varies between t 0 set to 0, the beginning of an observation, and t 0 + n int ∆t int , n int ∈ N, where n int is the number of integrations, and ∆t int is the integration time.The response drift effect was tested for fluxes within the [0DN/s; 5000DN/s] range.
It can be expressed as follows: where S 0 is the expected flux level in DN s −1 , a 1 (S 0 ) and a 2 (S 0 ) are the amplitudes of the two exponentials in DN s −1 , and α 1 (S 0 ) and α 2 (S 0 ) are the time constants in seconds.
Anneal recovery has a similar aspect as the idle recovery.
As soon as the annealing process is stopped and the detector is cooled down again, anneal recovery starts.This may be prior to the beginning of an observation.The time variable is therefore expressed as follows: t + t A , where t A is a negative value that refers to the starting time of the anneal recovery.Anneal recovery is expressed as follows: where S 0 is the expected flux level in DN s −1 , b 1 and b 2 are the amplitudes of the two exponentials in DN s −1 , and β 1 and β 2 are the time constants in seconds.
The idle recovery effect is also expressed with one exponential.Its amplitude depends on the time spent resetting before the observation.Idle recovery is expressed as follows: where S 0 is the expected flux level in DN s −1 , c(S 0 , ∆t I ) is the amplitude of the exponential in DN s −1 , ∆t I is the time spent idling before the observation, and γ(S 0 ) is the time constant in seconds.

The input format of MIRISim-TSO is compatible with
MIRISim outputs, that is, a series of integrations of 3D arrays in DN.The first two dimensions are the dimensions of the LRS slitless subarray (72 x 416 pixels), and the third dimension is the number of frames within one integration.To add persistence effects to the input data, Eq. 6, 7 and 8 are integrated between t and t + ∆t frame , where ∆t frame is the frame time.Effects are therefore added frame per frame in DN.This operation is applied to each pixel independently.In this way, the flux dependence of each frame is taken into account at the pixel scale.

Adding the background to the simulations
To include the background in the simulations, we used the background observation acquired on 26 May 2022 of the calibration target BD+60-1753.It consists of four LRS slitless integrations of 125 frames each.The flux level in DN s −1 was then computed by fitting the slope over the integration time.To add the background in our simulations, we created a unique background image from all four integrations by taking the median value of the four images.Bad pixels were masked using the data-quality flags (DQ)10 .These flags are meant to report any pixel issue that 400 could be related to an unreliable behaviour of the detector.For 401 example, flags report bad, hot, or saturated pixels.The slope val-402 ues in DN s −1 were integrated over the frame time to match the 403 ramp structure of the data at a pixel scale.For each simulated tar-404 get, the background image was scaled to match the observational 405 background.Photon noise was applied to the simulations in two steps.First, 408 the signal S at the frame level in DN was converted into elec-409 trons using the electronic gain.Then, samples were drawn from 410 a Poisson distribution and applied to the difference between pairs 411 of frames.The number N of photons (or electrons in our case) 412 received by a detector over a given time interval is described by 413 the standard Poisson distribution, where the S g product is the expected number of electrons, g be-415 ing the electronic gain.The variance Var[N] of this distribution 416 is 417 Photon noise in DN therefore is the standard deviation of the 418 Poisson distribution and varies as the square root of the signal, where S diff is the signal value in DN of the difference between 420 pairs of frames at the pixel scale.

421
The output files of MIRISim-TSO were then combined into 422 files of 2 Go, which correspond to the official segmented raw 423 data products of JWST provided by STScI11 .These segmented 424 files are 4D datasets with the same structure as the raw data prod-425 ucts.The first two dimensions are the dimensions of the LRS slit-426 less subarray (72 x 416 pixels), the third dimension is the number 427 of integrations, and the last dimension is the number of frames 428 within one integration.In this way, the simulation output files 429 are fully compatible with the STScI jwst reduction pipeline.We chose to reproduce the observation of L168-9b with the MIRI LRS slitless mode to fulfil two main purposes.First, our objective is to improve our simulations and make them as realistic as possible in order to prepare appropriately for future observation cycles.Then, our goal is to understand the origin of the systematics residuals that are witnessed in the data after reduction with the jwst pipeline (Bouwman et al. 2022).The reduction pipelines are still under optimisation for time-series observations and play a part in the spectrophotometric precision that we obtain.Reliable synthetic data can therefore be used to improve the reduction methods.In this section, we present our simulations of L168-9b and their comparison to real data.

Building the simulations
To be as realistic as possible, we used exactly the same setup for our simulations as in the real observations.We simulated the observation of a transit of L168-9b following the approach detailed in Sect.For the sake of consistency with real data, the same integration time was used to sample the simulated time series with exoNoodle.Fig. 3 shows the associated white light-curve in absolute flux computed with exoNoodle.This white light-curve as well as the whole spectral time series was then provided to MIRISim.MIRISim was run using the LRS slitless mode, with nine frames and two resets per integration.The total number of integrations matching the data sampling on purpose, we obtained 9371 integrations.The MIRISim setup is detailed in Table 2. Using MIRISim-TSO, we added the background, the photon noise, and the persistence effects at the pixel and frame scale.
As mentioned in Sect.2.3, we found evidence for persistence effects in the data that show an amplitude up to 1 to 2 % of the absolute flux.Fig 4 shows the overall aspect of the white light-curve of the time-series observation of L168-9b.We note a  strong presence of an exponential decay in flux at the beginning 488 of the observation.Based on the telemetry (see Appendix B), 489 several resets were performed prior to the observation, which 490 means that the persistence effect visible in the white-light data 491 is likely to be the idle recovery.To add the idle recovery to the 492 simulations that comply with the real data, we fit the real white 493 light-curve to derive the idle time parameter depicted in Eq. 8. 494 To do this, we selected a subarray centred on the spectral trace, 495 which is two times the half-width source aperture (4 pixels).This 496 corresponds to the pixels that are the most strongly impacted by 497 the idle recovery.Then, we extracted more than 35 spectral light 498 curves (following the process described in Sect.3.2.2).For each 499 extracted light curve, we fitted a flux-independent Eq. 8, and we 500 derived the amplitude c and time constant γ(S 0 ), knowing the 501 ∆t I parameter, which is given by the telemetry in Fig. B.1 of 502 Appendix B. Based on the amplitude and time-constant values 503 obtained for the spectroscopic light curves, we linearly interpo-504 lated the values to create flux-dependent parameters that can be 505 injected into Eq.8. Then we applied this idle formula to all pixels 506 at different signal levels.We discuss the persistence effects wit-507 nessed in the real data in more detail in Sect. 5.As mentioned 508 in Sect.2.3, no jitter was added to the simulations because the 509 stability of JWST is confirmed to be better than expected (Rigby 510 et al. 2023).This hypothesis is indeed verified in the real data, 511 for which no changes in position caused by the pointing stabil-512 ity are reported.The results of our simulations are displayed in 513 Fig. 5.The left panel shows a spectral image produced by the 514 simulation, and the right panel shows its comparison to a real 515 uncalibrated image.Even though the two are really similar, the 516 real image displays some features that do not appear in the sim-517

Data reduction and analysis
After the spectra were simulated, we proceeded to reduce them with the exact same methods as for the real data.In order to compare the results of our simulations, we also reprocessed the real data to verify our approach.For the sake of consistency, we used the same tools and compared the outcomes at different stages of the data reduction and analysis.The two main tools are the jwst pipeline (Bushouse et al. 2023), version 1.8.5 under CRDS context 1075 for the data reduction steps and the Eureka!Python package (Bell et al. 2022) for the analysis.This process was divided into five stages.The first two stages focus on detector-level corrections and calibration.The second stage includes background subtraction.Spectral extraction is performed in stage 3, and spectroscopic light curves are extracted in stage 4. The last stage is the light curve fitting with an astrophysical model and systematics detrending.Stages 1 and 2 were run with a slightly modified version of the default setup of the 1.8.5 version of the jwst pipeline, and stages 3 to 5 followed the steps defined in the Eureka!pipeline.In this section, we describe our approach for both the data reduction (stages 1 and 2) and the data analysis (stages 3 to 5).

Stages 1 and 2
Stages 1 and 2 from the jwst pipeline are meant to remove and correct for detector systematics.Stage 1 operates at the pixel and frame level to extract the mean count rate out of the nondestructive readouts of the detector.To do this, stage 1 is initialised using the DQ flags.Then, flagged pixels are used to generate a mask that is applied to the whole subarray and ensures that no unreliable pixels are included in further calculations.The next step is to apply the dark current correction by subtracting a 4D map of the dark current, interpolated towards the dimensions of the dataset.Some effects do not yet benefit from corrections: the first frame effect, the last frame effect, the RSCD for at least the first four frames (Argyriou 2021) and saturated pixels at the end of the ramp.The first and last frame effects are a direct consequence of the detector reset and present parity effects.As the name suggests, these effects impact the first and last frames of an integration, which are systematically lower than the empirical values evaluated by linear ramp extrapolations.The value of the last frame is even lower for odd lines, whose reset voltage is influenced by that of the adjacent line (Ressler et al. 2015;Argyriou 2021).Whenever a pixel is impacted by these effects, it is just flagged as a DO-NOT-USE pixel.The first and last frames are systematically rejected instead of being corrected.After the correction steps and flagging steps are applied, cosmic rays are detected and replaced using the two-point difference method (Anderson & Gordon 2011).In the first version of the pipeline described here (version 1.5.3 and CRDS context 0916), the RSCD effect did not benefit from a correction or a dedicated step.In the next versions, the RSCD step was not activated because its impact is not known when going through the first reductions.
The main function of stage 1 is to fit the ramp in order to derive a mean count rate, which is the slope of the ramp.A least-squares minimisation method is applied to fit the ramp.The default ramp-fitting algorithm uses an optimal weighting of the ramp that gives an additional weight to the first and last frames of the ramp, based on the number of frames (Robberto 2013).
As non-linearities may distort the ramp, a non-linearity correction is applied before fitting the ramp.The output of stage 1 is a 3D time series of the spectral images.The resulting quantity is the mean count rate, which is a flux of DN that is a number of DN per second crossing the section of a pixel (DN s −1 pixel −1 ).
This physical quantity can be directly related to a flux of photons after applying absolute flux calibration, which is the main goal of stage 2. The jwst pipeline relies on the use of reference and cali-622 bration files managed by the CRDS.To reduce our simulations 623 and data, we used the latest in-flight version of these files.Only 624 the electronic gain file was modified from 5.5 e − DN −1 to 3.1 625 e − DN −1 to comply with the value inferred during commission-626 ing.The gain reference file has not yet been updated in the CRDS 627 system, and the error arrays returned by the calibration pipeline 628 therefore currently underestimate the true noise.The transit is well centred for simulations, but slightly offset in 647 time for real data as some small uncertainties in the ephemeris 648 of the planet existed when the commissioning observation were 649 planned (in May 2022).In addition, the flux scatter is higher at 650 larger wavelengths in general and at the beginning of the expo-651 sure at short wavelengths in particular.This is discussed in more 652 detail in Sect. 4.

653
The main purpose of stage 4 is to produce a set of light 654 curves that stem from spectral binning, which increases the spec-655 tral signal-to-noise ratio.To do this, we divided our dataset into 656 25 spectral bins, from 5 to 12.2 µm with a step of 0.145 µm.A 657 subsample of 11 normalised light curves obtained from stage 4 658 of the 25 extracted bins is shown in the centre panel of Fig. 6 659 for the simulated and real data.Each light curve corresponds to 660 the time series of a given spectral bin.The second step of stage 661 4 is to apply a temporal sigma clipping of each light curve to re-662 move the outliers.The jwst pipeline currently did not correctly 663 handle dark current subtraction for segmented files.We therefore 664 rejected the first integrations of each segment from the dataset.665 In the particular case of the L168-9b observations, some data 666 points are also affected by the High Gain Antenna (HGA) move.667 Fig. 7 shows the outliers witnessed in the real data.Each vertical 668 black line marks the first integration of a segmented file, and the 669 red line shows the HGA move during the exposure.We clipped 670 the light curves using a box-car filter with a width of 100 integra-671 tions, with a maximum of ten iterations and a rejection threshold 672 of 5σ to reject these outliers.

673
Stage 5 fits both astrophysical and systematics models to in-674 fer the distributions of the transit parameters.Our fit was made 675 using the no-U-turn sampler (NUTS) described in Hoffman & 676 Gelman (2011), which is an extension of a Hamiltonian Monte 677 Carlo algorithm that automatically tunes the step size and the 678 number of steps per sample to avoid a random-walk behaviour.679 The transit was modelled using the starry package (Luger et al. 680   2. The reparametrised quadratic limb-darkening coefficients (q 1 , q 2 ) defined by Kipping & Sandford (2016), with uniform priors.This allowed them to vary between 0 and 1, U(0, 1). 3. The parameters of the one decaying exponential we used to correct for the persistence effect r 0 and r 1 .The model of systematics is expressed as where t is the time.We used loose normal priors derived from a preliminary least-squares fit on the white light-curve for c 0 , r 0 , r 1 .These priors are reported in Table 3.

Simulated data
Table 3. Priors used for the systematics detrending at stage 5 of Eureka!, derived from a preliminary least-squares optimisation.

Parameter
Prior distribution c 0 (DN s −1 ) N(1, 0.1 2 ) r 0 (DN s −1 ) N(−0.01, 0.005 2 ) r 1 (s −1 ) N(75, 30 2 ) The following parameters were first fitted on the white lightcurve and then fixed for the analysis per wavelength to avoid variations that would be non-physical.4 shows the comparison between the input parameters and the retrieved parameters for the simulations.Finally, stage 6 produces a transmission spectrum for the simulated and real data.Fig. 8 shows the transmission spectra obtained from the analysis of the simulated and real data.We demonstrate 722 that the simulations replicate the data with similar scatter and er-723 ror bars, thus leading to a mean transit depth of 569 ± 92 ppm 724 for the real data and 531 ± 90 ppm for the simulated data.To confirm our result, we performed a parallel analysis with 726 two other sampling methods: emcee (Foreman-Mackey et al. 727 2013), and dynesty (Speagle 2020).We obtained very consis-728 tent results, meaning that the transmission spectrum is not de-729 pendent on the sampling method.The N/S was calculated based on the standard deviation of an 734 out-of-transit part of the normalised light curves (between in-735 tegrations 3500 and 4000).This N/S computation method is 736 equivalent to the method used in Bouwman et al. (2022), which 737 plots the error bars obtained after the spectral fits because there 738 is no correlated red noise in the residuals.The Allan plots of 739 the residuals for spectral bins between 5 and 7 µm are shown in 740 Appendix C. In Fig. 9, the simulation mimics the real data very 741 well, except at short wavelengths, where the real data show an 742 excess of N/S up to 30% between 5 and 7 µm.We do not ex-743 pect a difference like this between the two datasets.This same 744 behaviour of the real data noise estimate is seen in Bouwman 745 et al. ( 2022).This additional N/S at short wavelengths does not 746 have any known physical interpretation, but it results in larger 747 error bars for the transmission spectrum at shorter wavelengths 748 and may consequently weaken future atmospheric retrievals.In 749 tion for it.

The noise-to-signal ratio
In this section, we investigate the additional N/S based on the tions and simulations.No excess in the N/S between 5 and 7 µm is visible any longer.
The reason for this increase in N/S at short wavelengths is complex.To provide a clear explanation of it, we start by recalling the cosmic-ray detection process.This process is based on the two-point difference method (Anderson & Gordon 2011), which rejects outlier frames that deviate from the mean value of the successive frame differences.When no cosmic ray or other non-linearity occurs in the ramp, the shape of the successive differences is expected to be a horizontal straight line.In reality, the ramps display residual non-linearities that remain at the end of stage 1, just before fitting the ramps.Fig. 11 shows the first four integrations of the data after correcting from ramp nonlinearities in Stage 1.The figure shows the ramps of the brightest pixel [389,36] (top panel) and those of a pixel with a medium signal [300,36] (third panel from the top).Panels 2 and 4 show the successive frame differences of these first four integrations.We observe non-linearity residuals of up to 4% for high signals.This value meets the 5% requirement made during ground-based testing, but is still not sufficient.
Because of these residual non-linearities, the cosmic-ray rejection stage flags frames that are not impacted by cosmic rays, but by these non-linearity residuals.These particular frames are mostly between frames 1 and 5.However, the first frame being affected by the first frame effect is never excluded from the dataset because it is not taken into account in the two-point difference method.Fig. 12 shows the flagging status of all nine frames for the first five integrations for all pixels of column 36 (rows ranging from 300 to 390).Fig. 12 displays a fluctuation in the number of frames flagged by the cosmic-ray detection step for different integrations (depicted in yellow), which creates a variability in the ramp adjustments as the number of frames considered for the adjustment varies significantly only at high fluxes, that is, at short wavelengths.The resultant slope values are therefore highly scattered at short wavelengths, thus creating the excess of N/S observed in the data.At medium and low flux level, we do not observe any additional frame rejection by the cosmic-ray detection stage.
The residual non-linearities may indeed be due to the RSCD effect or to any other detector effect (Argyriou 2021;Argyriou et al. 2023).As mentioned in Sect.2.2.2, RSCD is caused by resetting the detector (Ressler et al. 2023) and generates nonlinearities at the beginning of the ramp, which appear starting from the second integration.After commissioning, Morrison et al. (2023) investigated this effect on dark exposures and concluded that the RSCD displays two features.The fast decay appears in the second integration, where an additional signal at the beginning of the ramp decays exponentially.The first integration does not exhibit this decay as the preceding detector idling prevents any signal from accumulating in the detector traps.The second decay is a slope difference between integrations and depends on the number of frames.It appears that shorter integrations display larger differences that become even larger for the LRS subarray.The RSCD effect does not benefit from any correction in the pipeline, but the corresponding RSCD step flags the impacted four first frames.In other words, by activating the RSCD step, we remove the first four frames from the dataset.When the RSCD step is activated, which is applied upstream of the cosmic-ray detection step, the first four frames of the dataset are systematically excluded from the ramp-fitting step, thus limiting the variability of the number of frames used in the adjustments and limiting the increase of N/S at short wavelengths.
The RSCD step was not originally included in the pipeline.The first data reduction applied to the L168-9b data on 7 July  or a surplus of signal at the beginning of a time-series observation, and they are due to previous operations of the detector.In order to quantify the impact of these effects on exoplanet time series observed with the MIRI LRS slitless mode, we computed the light curve of each pixel of the slitless subarray.Depending on the location of the pixel in the subarray and therefore on the amount of flux it receives, we witness different types of persistence effects in the data.Fig. 13 shows the light curves of a set of pixels taken from the slitless subarray.The top left panel shows a rectangular selection of nine pixels located in the spectral trace, at the highest flux levels.The top right panel shows the corresponding time series.Pixel 5, and therefore, the column amid the selection (depicted in yellow), is the brightest pixel in the subarray, with a flux level up to 32 000 DN s −1 .The columns on the right (in light purple) and on the left (in dark purple) receive 38 % and 57 % less flux on average, respectively, than the columns in the middle.The difference in flux between the column on the left and on the right points out an asymmetry of the PSF, which is centred on columns 37 and 38.The loss of signal at 1.2 hours of observation, which is visible in the light curves, is due to the HGA move during the exposure.The bottom left panel shows a selection of pixels located in the background.All pixels display a very low, but still similar amount of flux between 30 and 50 DN s −1 .The corresponding time series are displayed in the bottom right panel.Remarkably, no significant persistence effect is visible for these pixels.
Persistence effects are only visible in the top right panel at the beginning of each time series, for less than approximately 30 minutes.Each column shows a different behaviour of these persistence effects.The light curves of pixels located in the brightest column display a decay of flux, whereas those from the left column display an increase in flux over time.Those from the right column display no change in flux at the beginning of the observation.Those from the background do not seem to show any persistence effect at all.The discrepancies between the persistence effects show that they are highly flux dependent.We therefore point out several flux regimes.First, high fluxes above 30 000 DN s −1 are impacted by the idle effect, introduced in Sect.2.3.Then, as the LRS slitless subarray receives background flux even when no target is observed, the response drift has less impact than during a first complete illumination.Therefore, fluxes be-Fig.11.First four integrations of L168-9b target data after application of the non-linearity correction step.Figures 1 and 3 (starting from the top) show the ramps of pixels [389,36] and [300,36], where 36 is the brightest column.The two-pixel rows (389 and 300) are representative of a strong and medium signal, respectively.Figures 2 and 4 show the successive frame differences, where non-linearities still remain even after the correction.
tween 15 000 and 30 000 DN s −1 must be in an equilibrium state  The overall aspect of the light curves after spectral binning 939 between 5 and 12 µm, displayed in Fig. 6, mainly shows a decay 940 rather than an increase in flux, except at the highest wavelengths.

941
The idle effect therefore dominates the response drift in the trace.

942
As discussed in Sect.3.1, this justifies the use of the idle effect 943 alone in our simulations.944 6. Perspectives

Test of several atmospheric models
In Sect.3.2 we have demonstrated that we are able to produce realistic MIRI LRS slitless simulations.In this section, we simulate MIRI LRS data for two atmospheric scenarios for L168-9b: -Scenario 1: a hydrogen-and helium-rich atmosphere.
-Scenario 2: a thick Venus-like atmosphere with a metallicity X 1000, 50% CO 2 , and 15% CO The atmospheric absorption model was generated using ATMO, a 1D-2D radiative-convective equilibrium model for planetary atmospheres (Tremblin et al. 2015(Tremblin et al. , 2017)).It solves the radiative transfer equation for a given set of opacities and computes a P-T profile that satisfies hydrostatic equilibrium and conservation of energy.It can compute equilibrium and nonequilibrium chemical abundances with the kinetic network of Venot et al. (2012).Here, for the two atmospheric scenarios, we assumed an equilibrium chemistry.Fig. 12. Flagging status at the frame level for the first five integrations (y-axis) of the L168-9b data after applying the cosmic-ray detection step (also known as the jump detection step) using version 1.8.5 of the pipeline.The flagging status is shown for all pixels located in the brightest column, 36 (x-axis).The colour code shows the status of each frame.Frames in blue are valid and therefore kept in the dataset for the ramp adjustment step.Frames in grey are rejected from the dataset by the last frame-effect detection step.Frames in yellow are rejected by the jump detection step and are mostly between number 1 and 5.They are located at short wavelengths (rows of pixels ranging from ∼ 340 to 390, mainly focused on pixels with high fluxes).come more critical as the fit is performed on only a few frames.Moreover, if the signal level in a pixel exceeds 10 000 DN, the last frames are likely to be impacted by the brighter-fatter effect that introduces non-linearities as well, which alters not only the signal level, but also the size of the PSF.The brighter-fatter effect (Argyriou et al. 2023) is caused by charge migration from the central brightest pixel of the PSF into the surrounding pixels as charge accumulates towards saturation.The effect causes a broadening of the PSF, and it is particularly marked at short wavelengths, where the contrast between the central and nearby pixels is stronger; a secondary manifestation is a steepening of the ramp in the pixels into which additional charge is migrating, which adds to the residual non-linearities in their ramps.As presented by Morrison et al. (2023), the RSCD effect also includes non-linearities at the beginning of the ramp.A correction for the RSCD effect may therefore be needed to remove the resulting non-linearity behaviour and to keep the impacted frames in the dataset to provide more frames for the fit.Thus, correcting for all non-linearities in general may become crucial to provide a linear ramp to fit and to keep all the frames within the dataset.Finally, the next step is to examine all the LRS slitless observations that are available to date, including calibration data acquired during commissioning on HD167060, HD180609, and HD37962, to derive correlations between the non-linearity behaviour and the number of frames, as well as correlations between non-linearities and the signal and flux levels in a pixel.We expect the outcome of this work to provide key insights into the extent of non-linearities and correction methods that are currently at stake.

Future work on persistence effects
As presented in Sec. 5, the persistence effects witnessed in the L168-9b data are flux-dependent and affect the spectral time series.The only way to remove them from the dataset is to fit them with a model made of one or two exponentials or even with a polynomial function.As persistence effects occur at the curve modulations are strongly correlated to persistence effects, and fitting them is therefore highly degenerate.Bell et al. (2023) showed the strong impact of persistence effects in the MIRI LRS slitless phase-curve observation of the ERS target WASP-43b.
Not only are they strongly correlated to the phase-curve parameters, they also follow the same regimes of the idle effect and response drift as found in the L168-9b data in Sect. 5. Furthermore, Bell et al. (2023) showed that persistence effects are not only flux dependent, but also result from the pixel location in the subarray.As the LRS slitless subarray overlaps the MIRI coronographs (Boccaletti et al. 2015) and constantly receives observational background, the persistence effect may depend on the fil-ter position prior to the observation.Long-wavelength filters are 1048 likely to increase the observational background.To characterise 1049 the persistence effects, one way is to pull down the telemetry of 1050 the MIRI LRS slitless subarray to find correlations between these 1051 effects and the filter wheel position and the time spent idling 1052 prior to the observation.We introduced realistic simulations of transiting exoplanets with 1055 the MIRI LRS slitless mode.In particular, we simulated the ob-1056 servation of L168-9b, a super-Earth-sized exoplanet chosen to 1057 meet the requirements of the instrumental stability calibration 1058 program conducted during commissioning in the LRS mode of 1059 the MIRI instrument.To ensure that our simulations complied 1060 with real data, we refined and adapted the detector set-up, the 1061

10
nificant part of the observations is dedicated to exoplanet obser-11 vations, either directly or indirectly.When considering all Cycle 12 1 programs proposed for Early Release Science (ERS), Guaran-13 teed Time Observation (GTO) and General Observations (GO), 14 115 distinct transiting exoplanets are being observed with the 15 JWST.Twenty-one of these 115 planets have been observed with 16 the MIRI instrument, 9 of them with the Low Resolution Spec-

114 1 .
The emission spectrum of the star.

115 2 .
The day-and nightside emission spectrum of the planet. 116 127

181A
(in m 2 ), and T T is the dimensionless telescope transmission 182 function.183 As the diffraction pattern is sampled into a given number of 184 pixels, each pixel receives a sub-amount of the overall absolute 185 PSFs flux.The absolute PSFs flux can be described as the sum 186 of the flux that arrives at each pixel, 187

430 3 .
The case of L168-9b 431 In this section, we perform simulations of the transiting exo-432 planet L168-9b.JWST observed L168-9 as part of the MIRI LRS 433 slitless commissioning under program ID 1033.The aim of this 434 program was to test the time-series observation mode including 435 the spectrophotometric stability.L168-9 is a bright M1V star lo-436 cated 25 pc away.It is orbited by a warm super-Earth that was 437 initially announced by Astudillo-Defru et al. (2020).L168-9b 438 has a radius of R p = 1.39 R ⊕ and a mass of M p = 4.6 M ⊕ .439 This target was selected to meet the objectives of this calibra-440 tion program because it probably shows no strong atmospheric 441 features.The planet is expected to have an equilibrium temper-442 ature between 668 K and 965 K, and it has been found to be 443 free from any primordial hydrogen-helium envelope (Astudillo-Defru et al. 2020).The outcomes of this program are presented in detail in Bouwman et al. (2022).L168-9b was observed on 29 May 2022 for 4.2 hours.The full observation is composed of 9371 integrations and represents 6 Gb of data.
2. As inputs for exoNoodle, we used a synthetic PHOENIX stellar spectrum (Husser et al. 2013) interpolated to the appropriate temperature and metallicity.Fig. 2 shows the stellar emission spectrum of L168-9.The planet parameters are derived from the literature (Astudillo-Defru et al. 2020; Patel & Espinoza 2022), and
(2020) a (AU) 0.02091 ± 0.00024 Astudillo-Defru et al. (2020) ulation.In the real data, we note a hot pixel column on the left 518 that is fully saturated.The top part shows spectral contamination 519 (Bouwman et al. 2022; Bouchet et al. 2022) that remains under 520 investigation.The corresponding pixels are below the MIRI LRS 521 dispersion range, between 4 and 5 µm.We chose not to replicate 522 these features as they are systematically removed during the re-523 duction steps.

Fig. 5 .
Fig.5.MIRI LRS slitless subarray 416 x 72 pixels, cut at pixel 155 at the bottom of the spectra.L168-9b raw simulations and comparison to uncalibrated real data.The left panel shows the simulated spectral image, and the right panel is taken from the uncalibrated data.Both images are the last (ninth) frame of the first integration.The vertical and horizontal directions are pixels.The vertical axis is the spectral dispersion direction and is related to the wavelengths.Thus, pixels with high fluxes are located in the top part of the vertical axis, at short wavelengths.Then, the flux level decreases as the wavelength increases along the vertical axis.The real data display a hot pixel column on the left that is fully saturated (a small part of it is encircled in orange).The upper part of the trace shows two distinct zones.The top part, indicated with the blue arrow, shows a spectral contamination(Bouwman et al. 2022;Bouchet et al. 2022).

Stage 2
uses as input a set of mean count rates in DN s −1 at the pixel level and converts them into MJy sr −1 , applying a calibration factor determined during commissioning.In the peculiar case of transiting exoplanets, absolute flux calibration is skipped.As planetary flux variations are always relative to the stellar flux, the light curves are normalised, and therefore, no absolute calibration is required.Before calibration, each spectral image is divided by the flat-field reference image.The last step of stage 2 is to subtract the background from each spectral image.Because there is no slit to isolate the point source, slitless spectroscopy also integrates background flux over time, which can be removed.Background is subtracted following the method explained inBouwman et al. (2022).Ten columns on the left and ten on the right side of the trace (column 36) are selected.Then, the median value is used as a reference to create a background image.This image is then subtracted from the spectral image.As each spectral image has its own background, this process has the benefit of removing any time-variable instrumental features from the data (such as the 390 Hz noise; private communication, M.Regan, 2023).Persistence effects are flux-dependent, and therefore, the background pixels corresponding to a low level of flux are less impacted by these effects.Removing the background does not allow us to correct for these effects in the pixels located in the spectral trace.We refer to Sect. 5 for further explanations of persistence effects.Then, a spatial filter of outlier detection is applied to remove any hot pixels that would have been left in the subarray.This technique consists of applying a running median on each column vertically to detect and reject any outlier (i.e. a hot pixel or a saturated pixel) left in the subarray based on a 5σ rejection threshold.Finally, the spectrum is cut at pixel 395, around 4.5 µm to avoid scattered light at short wavelengths.
5, we used the Eureka!pipeline(Bell et al. 631  2022).The input of stage 3 are .calintsfiles, which are spec-632 tral images produced by the jwst pipeline stage 2. In stage 3 the 633 1D spectrum is extracted from each image.Similarly to Bouw-634 man et al. (2022), we chose a rectangular selection from pixel 635 155 to pixel 385 (height) and from pixel 13 to pixel 64 (width).636 We used a half-width source aperture of 4 pixels.The result-637 ing waterfall plots from stage 3 for both real and simulated data 638 are shown in the top panel of Fig 6.The top panel shows the 639 temporal evolution over the whole exposure represented on the 640 y-axis in days, of each 1D spectrum displayed along the x-axis 641 in µm.The colour bar refers to the level of normalised flux as a 642 function of wavelength and time.The horizontal strip between 643 0.06 and 0.11 days in the left panel and that between 0.44 and 644 0.48 days in the right panel is darker than the rest of the time 645 series.This deficit of flux corresponds to the planetary transit.646

Fig. 6 .
Fig. 6.Comparison between real and simulated data at different stages of the reduction and analysis using Eureka!(Bell et al. 2022).Top panel: Waterfall plots from stage 3. Centre panel: Transit light curve as a function of wavelength, offset for clarity and produced with the chromatic 13 visualisation tool.Bottom panel: Raw and detrended white light-curve with the best-fit model resulting from the optimisation performed in stage 5 (red curves) and the corresponding Allan plot, showing the evolution of the root-mean-square (RMS) of the light curve as a function of the binning size.The photonnoise limit is almost reached in both cases, attesting to the Gaussianity of the residuals. 685

Fig. 7 .
Fig. 7. White light-curve over the whole MIRI LRS slitless subarray of the real observation of L168-9b.Outliers are detected in the whole time-series data points.The black lines mark the position of the first integration of each segmented file, and the red line shows the timing of the High Gain Antenna (HGA) move during the exposure.

1.
The time of mid-transit t 0 , with a normal prior of N(59728.4603,0.005 2 ) BMJD TDB computed for the nearest epoch based on the best-fit ephemeris derived from the analysis of the most recent TESS observations (courtesy of Billy Edwards).2. The orbital period P, with a normal prior based on the value reported by Astudillo-Defru et al. (2020) N(1.40150, 0.00018 2 ). 3. The inclination i, with a normal prior based on the value reported by Astudillo-Defru et al. (2020) N(85.5, 0.7 2 ) 4. The eccentricity e, with a normal prior N(0., 0.1 2 ) taken arbitrarily from the unique known constraint, that is, e < 0.21 (Astudillo-Defru et al. 2020).Using the probabilistic programming framework pymc3, which features a NUTS sampler, we ran three chains with a number of iterations to tune the set to 2000 and a number of draws set to 2000.In the bottom panel of Fig.6, we show the resulting best fit of the white light-curve from real and simulated data using NUTS.The root-mean-square (RMS) of the detrended white light-curve is 501 ppm and 519 ppm for simulated and real data, respectively.Table Fig. 8. Transmission spectrum of L168-9b obtained from simulated (red dots) and real data (grey dots).The blue dots are the points from Bouwman et al. (2022), and the purple curve is the reference value from Patel & Espinoza (2022) along with the 95% confidence interval.

Fig. 9
Fig.9shows the evolution of the noise-over-signal (N/S ) es-732 timate of the spectral light curves for simulations and real data.733 The N/S was calculated based on the standard deviation of an 734 out-of-transit part of the normalised light curves (between in-735 tegrations 3500 and 4000).This N/S computation method is 736 equivalent to the method used inBouwman et al. (2022), which 737 plots the error bars obtained after the spectral fits because there 738 is no correlated red noise in the residuals.The Allan plots of 739 the residuals for spectral bins between 5 and 7 µm are shown in 740 Appendix C. In Fig.9, the simulation mimics the real data very 741 well, except at short wavelengths, where the real data show an 742 excess of N/S up to 30% between 5 and 7 µm.We do not ex-743 pect a difference like this between the two datasets.This same 744 behaviour of the real data noise estimate is seen in Bouwman 745 et al. (2022).This additional N/S at short wavelengths does not 746 have any known physical interpretation, but it results in larger 747 error bars for the transmission spectrum at shorter wavelengths 748 and may consequently weaken future atmospheric retrievals.In 749

Fig. 9 .
Fig. 9. Noise-over-signal (N/S ) estimates of the spectral light curves of L168-9b for the simulations and real data.The blue curve shows the estimate for the N/S of the data, and the orange curve shows the estimate for the N/S of the simulation.

Fig. 10 .
Fig. 10.Noise-over-signal (N/S ) estimates of the L168-9b observation and simulation after removing the first four frames of each integration of the exposure before fitting the ramp using the jwst pipeline Stage 1.No additional N/S at short wavelengths is witnessed any longer.

884 5 .
Investigating persistence effects 885 In this section, we present a comprehensive investigation of the 886 persistence effects found in the L168-9b data based on a quanti-887 tative analysis of their amplitude, time constant, and flux depen-888 dence.As presented in Sect.2.3, persistence effects are a deficit 889 932 between idle and response drift.Fluxes lower than 15 000 DNs −1 933 seem to be affected by the response drift effect, which occurs 934 when the subarray is illuminated after a pause in the observa-935 tion.Finally, fluxes with an order of magnitude of a few tens, 936 generally coming from the observational background, are only 937 affected by a very low amplitude response drift. 938

Fig. 13 .
Fig. 13.Depiction of persistence effects in the L168-9b data.Upper part: Full light curves.Bottom part: Zoom into the beginning of the light curves up to half an hour.Top left panel: Rectangular selection of nine pixels located in the spectral trace, at the highest flux levels.Top right panel: Corresponding time series.Persistence effects are visible at the beginning of each time series for less than approximately 30 minutes.Bottom left panel: Selection of pixels located in the background.All pixels display a really low but yet similar amount of flux between 30 and 50 DN s −1 .Bottom right panel: Corresponding time series in which no persistence effect is visible for these pixels.

Fig. 14 .
Fig. 14.Transmission spectrum of L168-9b from the analysis of simulated data with MIRISim-TSO.Top: Transmission spectra for two atmospheric scenarios, Venus-like (red dots) and H/He-dominated (blue dots), as well as from real commissioning data (grey dots).Middle: Differences between the observed transit depth and the simulated transit depth for a hydrogen-helium-dominated atmosphere, where σ = (depth obs [λ]−depth HHe [λ]) err 2 depth,obs +err 2 depth,HHe

Fig. B. 1 .
Fig. B.1.MIRI telemetry during the 48 hours preceding the observation of L168-9b that provides the detector status, either in exposure (observing) mode, or in clocking mode (in repeated idle procedure).The observation is preceded by a series of idles lasting 9h13.

Fig. C. 1 .
Fig. C.1.Allan plots of the first nine spectral bins of the L168-9b data between 5 and 7 µm.All plots show the minimum correlated noise in the residuals of the fit.

Fig. D. 1 .
Fig. D.1.L168-9b spectrum in DN extracted using only the penultimate frame of the ramp.No ramp-fitting is applied.

Fig. D. 2 .
Fig. D.2.Noise estimate in DN of the L168-9b observation extracted using only the penultimate frame of the ramp.

Fig. D. 3 .
Fig. D.3.Top: Noise-over-signal (N/S ) estimate for the L168-9b data between 5 and 12 µm, using only the penultimate frame of the ramp to compute the signal, instead of fitting the ramp.The blue curve shows the data N/S estimate, and the orange curve shows the simulation N/S estimate.Bottom: Zoom-in between 2500 and 20000 ppm to focus on the short wavelength N/S values.
. To 148 create realistic simulations of time-series observations, we use 149 the stable version 2.4.2 of the MIRI official simulator (Klaassen 150 et al. 2020, MIRISim).MIRISim simulates almost all MIRI ob-151 servation modes (except for coronagraphy), including imaging, 152 the LRS and Medium-Resolution Spectrometer (MRS) modes.153 The code itself is based on the use of calibration data products 154

Table 1 .
Input parameters for exoNoodle to simulate a transit observation of L168-9b with the MIRI LRS slitless mode.

Table 2 .
L168-9b observation parameters used in MIRISim.All parameters come from the observational setup of L168-9b validated by the Astronomer's proposal tool 12 .

Table 4 .
Retrieved parameters from the white light-curve fit and comparison to the exoNoodle input parameters were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.These observations are associated with calibration program 1033.This work is based [in part] on observations made with the NASA/ESA/CSA James Webb Space Telescope.The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST.These observations are associated with program 1280.MIRI draws on the scientific and technical expertise of the following organisations: Ames Research Center, USA; Airbus Defence and Space, UK; CEA-Irfu, Saclay, France; Centre Spatial de Liège, Belgium; Consejo Superior de Investigaciones Científicas, Spain; Carl Zeiss Optronics, Germany; Chalmers University of Technology, Sweden; Danish Space Research Institute, Denmark; Dublin Institute for Advanced Studies, Ireland; European Space Agency, Netherlands; ETCA, Belgium; ETH Zurich, Switzerland; Goddard Space Flight Center, USA; Institut d'Astrophysique Spatiale, France; Instituto Nacional de Técnica Aeroespacial, Spain; Institute for Astronomy, Edinburgh, UK; Jet Propulsion Laboratory, USA; Laboratoire d'Astrophysique de Marseille (LAM), France; Leiden University, Netherlands; Lockheed Advanced Technology Center (USA); NOVA Opt-IR group at Dwingeloo, Netherlands; Northrop Grumman, USA; Max-Planck Institut für Astronomie (MPIA), Heidelberg, Germany; Laboratoire d'Études Spatiales et d'Instrumentation en Astrophysique (LESIA), France; Paul Scherrer Institut, Switzerland; Raytheon Vision Systems, USA; RUAG Aerospace, Switzerland; Rutherford Appleton Laboratory (RAL Space), UK; Space Telescope Science Institute, USA; Toegepast-Natuurwetenschappelijk Onderzoek (TNO-TPD), Netherlands; UK Astronomy Technology Centre, UK; University College London, UK; University of Amsterdam, Netherlands; University of Arizona, USA; University of Bern, Switzerland; University of Cardiff, UK; University of Cologne, Germany; University of Ghent; University of Groningen, Netherlands; University of Leicester, UK; University of Leuven, Belgium; University of Stockholm, Sweden; Utah.We would like to thank the following National and International Funding Agencies for their support of the MIRI development: NASA; ESA; Belgian Science Policy Office; Centre Nationale d'Études Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luft-und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economiá y Competividad; Netherlands Research School for Astronomy (NOVA); Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Board; and UK Space Agency.Software: numpy (Harris et al.2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), astropy (Collaboration et al. 2022), nuts (Hoffman & Gelman 2011), emcee (Foreman-Mackey et al. 2013), dynesty (Speagle 2020), starry (Luger et al. 2019), exoNoodle (Martin-Lagarde et al. 2020), MIRISim (Klaassen et al. 2020), eureka!(Bell et al. 2022), pipeline_parallel (https://gitlab.com/jwst_fr/pipeline_parallel/-/tree/master/) and jwst (https: //jwst-pipeline.readthedocs.io/en/latest/).
). 1099 P.O.L. acknowledges funding support from CNES.This work uses observa-1100 tions made with the NASA/ESA/CSA James Webb Space Telescope.The data 1101