© The Authors 2025

1 Introduction

One of the most profound discoveries in the last decades in astronomy is that exoplanets are extremely common in the Universe. This discovery has not only introduced us to a range of strange and exciting worlds fundamentally different from the planets in our Solar System but has also allowed us to gain a deeper insight into their formation and evolution. Transit observations of exoplanets are crucial for comparative studies, and when combined with radial velocity data, these measurements allow the inference of the planet’s mass, radius, and composition. However, little is known about the processes that determine the atmospheric structure of these planets, which highlights the need for spectroscopic observations using space-based and ground-based instruments.

In particular, high-resolution cross-correlation spectroscopy (HRCCS) has emerged as a powerful technique for characterising exoplanetary atmospheres. With significant advancements in ground-based high-spectral resolution platforms (e.g. CRIRES+, ESPRESSO, IGRINS, IRD; Dorn et al. 2014, 2023; Pepe et al. 2021; Park et al. 2014; Mace et al. 2018; Tamura et al. 2012; Kotani et al. 2018), HRCCS has become instrumental in constraining basic atmospheric properties. First demonstrated by the detection of CO in the hot Jupiter HD 209458b (Snellen et al. 2010) using the CRyogenic high-resolution InfraRed Echelle Spectrograph (CRIRES; Kaeufl et al. 2004), HRCCS observations leverage the planetary Doppler shift of the dense forest of individually resolved spectral lines to detect both atomic and molecular species robustly. Since then, it has been used to detect several chemical species in the atmospheres of exoplanets, particularly those of close-in giant planets, such as hot and ultra-hot Jupiters (UHJs). Due to their proximity to their host stars, these tidally locked gas giants are a fascinating population without similar counterparts in our Solar System. Their high equilibrium temperatures (T_eq ≳ 2200 K) not only create conditions where major C- and O-bearing species (such as CO, OH and atomic O) are expected to exist in gaseous form but they also have a significant impact on their atmospheric structure and composition, such as the occurrence of inversion layers due to the presence of molecular compounds like TiO and VO (Prinoth et al. 2022; Pelletier et al. 2023; Maguire et al. 2024) high up in the atmospheres as well as due to absorption by atomic metals and metal hydrides (Lothringer et al. 2018; Gandhi & Madhusudhan 2019; Flagg et al. 2023). This makes UHJs ideal targets for atmospheric characterisation wherein they serve as a treasure trove for numerous chemical species (e.g. Fe I, Fe II, Ca I, Ca II, OH, CO, H₂O, etc.) that have been detected using high-resolution transmission and emission spectroscopy (e.g. Hoeijmakers et al. 2019; Gibson et al. 2020; Pino et al. 2020; Maguire et al. 2023; Nugroho et al. 2021; Brogi et al. 2023; Ramkumar et al. 2023).

The HRCCS technique is not limited to transmission spectroscopy; it is also powerful when applied to emission spectroscopy and phase curve studies of close-in gas giant planets. By analysing the thermal emission from the planet’s day-side and the variations in brightness throughout its orbit, we can obtain measurements of the temperature structure and chemical composition, which provides crucial information about day-night temperature differences, atmospheric dynamics, and heat redistribution. In addition, resolved phase-curve observations also enable us to obtain spectra as the planet rotates, allowing us to effectively spatially resolve the surface (e.g. Roman & Rauscher 2019; Parmentier et al. 2021). When combined with high-resolution Bayesian methods (e.g. Brogi & Line 2019; Gibson et al. 2020), these techniques allow us to recover vital constraints on the atmosphere, such as the abundance of chemical species, as well as the atmospheric C/O ratios and metallicities (e.g. Line et al. 2021; Brogi et al. 2023; Ramkumar et al. 2023; Smith et al. 2024a); quantities that can provide potential insights into the formation and subsequent evolution history of these extreme planets (e.g. Öberg et al. 2011; Madhusudhan et al. 2011).

Recently, CRIRES has been upgraded to a cross-dispersed echelle spectrograph (Dorn et al. 2014; Follert et al. 2014) and the improved spectrograph, now known as CRIRES+, offers a ten-fold improvement in spectral coverage, making it ideal for characterising the atmospheres of exoplanets. MASCARA-1b, an ultra-hot Jupiter orbiting a fast-rotating bright A8 star (V = 8.3, K = 7.7) with an orbital period of ≈2.15 days (Talens et al. 2017; Hooton et al. 2022), is a prime target for such studies due to its high equilibrium temperature (2594 K; Hooton et al. 2022). Analyses using high-resolution emission spectroscopy (e.g. CRIRES+, PEPSI, CARMENES) have detected both molecular (CO, H₂O; Holmberg & Madhusudhan 2022; Ramkumar et al. 2023) and atomic species (Fe I, Ti I, and Cr I; Scandariato et al. 2023; Ramkumar et al. 2023; Guo et al. 2024) in the day-side atmosphere of this gas giant, underscoring the planet’s suitability for atmospheric characterisation at high resolution.

We have recently analysed the CRIRES+ phase-curve observations of MASCARA-1b, taken during the SV run in September 2021 (program 107.22TQ.001; Ramkumar et al. 2023). The observations were taken in the K-band at the peak of the thermal emission from the planet just before the secondary eclipse (covering orbital phase ∼0.34–0.43) and have enabled high-significance detections of CO, H₂O, and Fe. The species are found as emission signatures, implying a temperature inversion in the atmosphere, and by using a Bayesian retrieval framework, we were able to place constraints on the C/O ratio and temperature-pressure (T –P) profile, confirming the consistency of C/O with the solar value. Our analysis also revealed intriguing evidence that we are sensitive to temporal and spatial variations in the planet’s atmosphere. For instance, our cross-correlation maps, as well as the retrievals, hint that the Fe feature originates in different regions of the atmosphere than the CO/H₂O, with the rotation (υ_eq~3.6 km s⁻¹ assuming tidal locking) and winds (>5 km s⁻¹, e.g. Kataria et al. 2016) separating the features in velocity space. To explore potential variations and understand the atmosphere’s chemical composition and temperature patterns, we observed MASCARA-1b using the same instrumental setup as the previous observations but after the secondary eclipse, thus extending the phase coverage. Analysing phase curves separately, for example, before and after the eclipse, can facilitate an immediate comparison of the T–P profiles, C/O ratios, and metallicities. Furthermore, incorporating velocity constraints will provide valuable information about how atmospheric properties vary as the planet rotates and will help explore the importance of 3D effects. When combined with an atmospheric retrieval, constraints on the parametric phase curve, as well as spatially varying T–P profile and abundances, can be placed, enabling us to pinpoint the chemical and temperature composition of the atmosphere.

Here, we present a complementary analysis of high-resolution secondary eclipse observations taken two years apart with the CRIRES+ spectrograph, covering the post-eclipse phases of MASCARA-1b’s orbit. The paper is organised as follows: in Section 2, we present our observations and data reduction; in Section 3, we describe our forward model atmosphere for emission and detail our methodology, including the cross-correlation technique and present our detection results. We then outline our atmospheric retrieval frameworks and present our retrieval results for the new data and the joint analysis in Section 4. Finally, we discuss our iron detection, examine signatures of atmospheric dynamics, and contextualise MASCARA-1b’s atmosphere in Section 5 before concluding the study in Section 6.

2 Observations and data reduction

Observations of MASCARA-1b were taken two years apart for two half nights in September 2021 and October 2023 as part of programs 107.22TQ.001 and 112.260X.001 (PI: Gibson, dPI: Ramkumar) using the upgraded CRIRES+ spectrograph (R ~ 100 000; λ-1921−2472 nm (K-band); Follert et al. 2014) installed on UT3 of the VLT. The data taken on 2021 September 16, covering the pre-eclipse phases, were previously presented and analysed in Ramkumar et al. (2023), and we refer the reader there for further details. The data taken on 2023 October 16, covering the post-eclipse phases, are presented here for the first time. This observing night was interrupted by thick clouds at intervals, resulting in the loss of the target for ∼30 minutes over the entire observing night (visible as horizontal bands in the cross-correlation map, e.g. Fig. 6). We obtained 256 exposures covering the orbital phase of MASCARA-1b from ϕ~0.54 to 0.63 (where ϕ = 0.5 corresponds to the secondary eclipse) and observed in the K-band using the K2166 wavelength setting with (interrupted) coverage from λ~1921 to 2472 nm. We calculated the orbital phase using the transit epoch from Hooton et al. (2022) and computed the barycentric velocity correction for each exposure using the online tool from Wright & Eastman (2014). An overview of the observations is shown in Table 1, and the observing conditions are presented in Fig. 1.

The basic calibration and extraction of the time-series spectra for each spectral order were performed with the ESO CR2RES data pipeline (version 1.1.5) and executed using the command-line interface ESOREX¹ (version 3.13.5). First, we combined each set of raw dark frames into master dark frames (which also produced a bad pixel map). Next, we computed master flat frames and performed trace detection. Following this, we performed wavelength calibration of the extracted spectra using the Fabry–Pérot Etalon (FPET) frame and checked for additional wavelength shifts during the observation by cross-correlating the data with Doppler-shifted telluric templates produced using ESO Sky Model Calculator (SkyCalc, Noll et al. 2012; Jones et al. 2013) over a velocity range of −300 to +300 km s⁻¹ in steps of 0.01 km s⁻¹ . However, we did not find any significant shifts (< 1 km s⁻¹); and therefore did not apply any additional refinement to the wavelength calibration. Finally, we calibrated and extracted the science spectra (1D spectra as a function of order). The extraction was done in pieces by dividing the image into several swaths with widths set to 400 pixels. In total, we obtained 7 echelle orders for the K2166 wavelength setting. Additionally, with CRIRES+, each spectral order is observed across 3 separate detectors (CHIP1, CHIP2, and CHIP3), so we treat each separate chip independently and hereafter refer to them as orders. We produced a 3D data cube (order × time and/or phase × wavelength) with 21 spectral orders ranging from λ ≈ 1920 to 2470 nm.

Following the reduction and spectral extraction, we performed similar pre-processing steps as Ramkumar et al. (2023) to clean and blaze correct the extracted spectra. The noise estimates for each order were then extracted by assuming Poisson- dominated noise following Gibson et al. (2020, 2022), to which we refer the reader. However, while this approach is well-suited for optical data (e.g. Gibson et al. 2020), it may not fully account for the final noise amplitude in near-infrared data, particularly for poorly corrected telluric, stellar, or systematic effects. Therefore, as an alternative to the optimisation presented above, we also estimated the uncertainties of each order post-SYSREM by taking the outer product of the standard deviation of each wavelength channel and spectrum, normalised by the overall mean. Re-running SYSREM with these new uncertainties, we found that the choice of error estimation method did not have a discernible impact on our detections or retrieved values. Since the post-eclipse observations were affected by clouds, we note that applying the same treatment as the pre-eclipse data is not ideal, as the blaze correction does not perform optimally for these frames and affects the shape of the posteriors, though in small ways. Therefore, to ensure that the pre-processing steps did not significantly distort the underlying exoplanetary signal, we performed injection and recovery tests detailed in Sect. 3.2.

Before employing the cross-correlation technique to search for any exoplanetary signal, we0 removed the stellar and telluric lines using a de-trending algorithm, SysRem (Tamuz et al. 2005), which was initially adapted for high-resolution spectroscopy by Birkby et al. (2013) and has since been successfully utilised for transmission and emission spectra (e.g. Brogi et al. 2017; Nugroho et al. 2017, 2020b,a; Gibson et al. 2020, 2022; Ramkumar et al. 2023). Following the procedure outlined in Gibson et al. (2022), we first normalised the data by dividing through the median spectrum in each order. Then, we performed SYSREM on the normalised data and corresponding uncertainties for every iteration, generated a resultant model, subtracted it from the input data to get the processed (residual) matrix, and repeated the procedure for the subsequent iteration. Thus, after N iterations, the SYS REM model for a single order is: $$D={\displaystyle \sum _{i=1}^N}{u}_i{w}_i^{\text{T}}=U{W}^{\text{T}}$$(1)

where U and W are matrices containing column vectors u_i and w_i . U is then stored for processing the forward model (see Sect. 3). We applied 15 passes of SYSREM; while that is somewhat arbitrary, we reran the analysis in Sect. 4 with different values of SYSREM (N = 5, 10 and 20). A step-by-step overview of our pre-processing steps is shown in Fig. 2 for a single order.

Fig. 1 Left: phase coverage of our CRIRES+ pre- and post-eclipse observations. Middle and right: observing conditions showing the variation of airmass and mean signal-to-noise ratio (S/N) with orbital phase. Vertical grey bands denote cloud-affected regions.

3 A search for atmospheric species

3.1 Model spectra

To search for chemical species in the atmospheres of exoplanets and place quantitative constraints on atmospheric properties, both the cross-correlation and likelihood approaches require spectral templates. Therefore, for our analysis, we computed the 1D atmospheric models using IRRADIATOR (Gibson et al. 2022), where we defined 80 atmospheric layers equally spaced in log pressure from 10⁻⁸ bar to 100 bar and computed the T – P profile at each layer using the parametric model from Guillot (2010) with four input terms T_irr (irradiation temperature), κ_IR (mean infrared opacity), γ (ratio of visible-to-infrared opacity) and T_int (internal temperature). Next, we set the volume mixing ratios (VMRs), χ_species, for each species of interest at each layer. These are either given as free parameters, that is, the species are assumed to be well-mixed, or computed using FASTCHEM² (version 3.1.1, Stock et al. 2018, 2022; Kitzmann et al. 2024) using an additional free parameter to specify the metallicity (relative to solar) and C/O ratio, where we adjust both C and O abundances to have a desired ratio but preserve their sum. We compute the emission by integrating through the atmosphere for every point in the wavelength grid, assuming blackbody emission from the deepest layer. The spectra are calculated across a wavelength range of 1800–2850 nm at a constant resolution of R = 200 000. In our atmospheric models, we accounted for the opacity contributions from chemical species prominent in the atmospheres of UHJs at near-infrared wavelengths (e.g. Snellen et al. 2010; Birkby et al. 2013; Nugroho et al. 2020a; Brogi et al. 2023). For atomic and molecular absorption, we used the high-resolution pre-computed opacity grids provided by petitRADTRANS³ (version 2.7.7, Mollière et al. 2019, 2020) and also included collision-induced opacity due to H₂-H₂ and H₂- He and bound-free and free-free absorption from H⁻, using the cross-sections from Gray (2005) and Abel et al. (2011). Additionally, we incorporated a parameterised cloud deck pressure (P_cloud), assuming infinite opacity below the atmosphere. Lastly, since data pre-processing methods such as SYSREM can alter the underlying planet signal, applying these same processing steps to our forward models as we do to the data is necessary. Thus, we implemented the model filtering technique introduced by Gibson et al. (2022), which makes use of the output matrices U and W (Sect. 2), containing the column and row vectors u and w for each SYSREM iteration. Before applying this technique, we convolved our model emission spectra with a broadening kernel with a standard deviation of W_conv , followed by linear interpolation to match the 2D wavelength grid of our data. After this, we Doppler shifted this model to a planetary velocity, υ_p , for each order: $${\upsilon }_{\text{p}}=K_{\text{p}}\sin (2\pi \varphi )+{\upsilon }_{\text{sys}}+{\upsilon }_{\text{bary}}$$(2)

where Kp , ϕ, υ_sys and υ_bary are the radial velocity semi-amplitude of the planet’s orbit, orbital phase, the systemic velocity offset and barycentric velocity, respectively. This produced a 3D shifted forward model that was then filtered using the SYSREM basis vectors via the procedure outlined in Gibson et al. (2022).

In summary, we specify our 1D atmosphere using the parametric model from Guillot (2010) and two different chemical regimes, assuming either that the species of interest are well- mixed (constant VMRs with altitude) or by using FASTCHEM to compute the chemical profiles after specifying a C/O ratio and metallicity. We applied various combinations of these parametrisations to our cross-correlation and retrieval analyses, which we discuss in more detail in Sect. 4. The model spectrum (constant VMRs) and corresponding T–P profile are shown in Fig. 3.

Fig. 2 Data pre-processing steps, applied to a single order: (a) raw spectra, (b) cleaned and blaze corrected spectra, and (c) residual spectra after division by the median spectrum and subtraction of the SYSREM model (shown weighted by the uncertainties). We note that the blaze correction does not work well for cloud-affected frames, resulting in some horizontal structure (residuals) in these regions (the last two panels). However, they only have a minimal impact on our retrievals.

3.2 Injection and recovery tests

Before performing a retrieval analysis (outlined in Sect. 4) using our forward model, we conducted injection and recovery tests to ensure the pre-processing and filtering steps did not alter the underlying planet signal. We injected an atmospheric model containing CO, H₂O, and Fe with a negative value of the expected Kp and performed our retrieval analysis on the injected data. We found the retrieved model parameters and T–P profiles to be in agreement with the injected values (see Appendix A). Additionally, although the cloud-affected frames effectively contain weaker signals due to low S/N and could be filtered out, given the already limited phase coverage and data, we chose not to exclude these frames. To ensure their inclusion did not bias the results, we ran retrievals (Sect. 4) with and without the cloud-affected frames. We found the retrieved parameters and T – P profiles to be consistent within ≈1σ (see Appendix B). Although including these frames introduced small changes to the temperature gradient and irradiation temperature (T_irr), they did not impact the overall results.

3.3 Cross-correlation

After removing the stellar and telluric lines from our data using SYSREM, the remaining data is dominated by noise (panel c of Fig. 2). To extract the planetary signal, we employ the traditional cross-correlation technique (e.g. Snellen et al. 2010) and generate the cross-correlation function (CCF) by summing the product of the residual data (f_i) and the shifted (+filtered) model (m_i), weighted on the variance of the data (σ_i) over both spectral order and wavelength: $$\text{CCF}\left(\text{Δ}{\upsilon }_{sys}\right)={\displaystyle \sum _{i=1}^N}\frac{f_i{m}_i\left(\text{Δ}{\upsilon }_{sys}\right)}{{\sigma }_i^2}$$(3)

This produces a cross-correlation array of phase vs Δυ_sys, referred to as the cross-correlation map, where Δυ_sys represents the residual velocity shift relative to the stellar rest frame after correcting for systemic and barycentric velocities⁴. For strong planetary signals, the CCF will trace out the predicted radial velocity with phase in the CC-map, allowing us to confirm the presence of atmospheric species in exoplanets. However, the strength of the cross-correlation (CC) signal is often much weaker, and we need to integrate the cross-correlation map over a range of predicted planetary velocities (υ_p). Since the exact value for the radial velocity semi-amplitude (K_p) is unknown and to evaluate noise, we select a range of K_p values close to the predicted value from radial-velocity measurements (Table 2) and shift the CCF to a new planetary velocity (according to Eq. (2)) and then integrate over time to produce a K_p–∆υ_sys map. Integration helps isolate the source of the signal in velocity space, and we can be confident of the detection of species in the atmosphere.

The final step in securing the spectroscopic detection of the exoplanet is to determine the significance level of the signal, which measures the signal amplitude. There are several ways to do this (see, e.g. Snellen et al. 2010; Birkby et al. 2013). We subtract the K_p–Δυ_sys map by the mean in regions away from the peak before dividing through by the standard deviation (Brogi et al. 2012, 2018). Since it is difficult to interpret the significance of these values statistically due to the arbitrary selection of this region, we utilise the likelihood distribution (see Sect. 4.1) of our model scale factor, α, to compute the detection significance, similar to the approach of Gibson et al. (2020), which we discuss in the following sections.

Fig. 3 Left panel: best-fit model spectra (constant VMRs) of MASCARA-1b with contributions from the chemical species shown individually. Middle and right panels: the parametric T–P profile from Guillot (2010) and combined contribution function.

3.4 Cross-correlation results

Following the procedure outlined in the previous section, we searched for individual species of interest by cross-correlating our data with a filtered model emission spectrum using the bestfitting model parameters from our atmospheric retrieval (see Sect. 4). We calculated the K_p–Δυ_sys maps by shifting the CCFs to the rest frame of the planet over a range of Kp and Δυ_sys from −300 to +300 km s⁻¹, in steps of 0.3 and 0.2 km s⁻¹, respectively, and summed over time. To estimate the detection significance (sometimes referred to as signal-to-noise, e.g. Brogi et al. 2018), we divided the K_p–Δυ_sys maps by their standard deviation calculated by avoiding the area ±50 km s⁻¹ from the expected planet signal.

In our search for chemical species in the day-side atmosphere of MASCARA-1b, we detected strong emission signatures of CO, H₂O and Fe in the pre-eclipse data and CO and H₂O in the post-eclipse data (second panel of Fig. 4), respectively. These emission features are also detected at high significance when combining the two nights (middle panel of Fig. 4). Despite the lower S/N of the post-eclipse data, CO and H₂O are detected at similar strengths as in the pre-eclipse data. We do not detect iron in the post-eclipse and the combined dataset. The species are detected roughly at the expected K_p and Δυ_sys, except for the pre-eclipse Fe signal that has a slight offset in both K_p and Δυ_sys. We discuss our iron detection in more detail in Sect. 5.2. The K_p–Δυ_sys maps of the individual species and the full atmosphere for the pre-, post- and pre+post datasets are shown in Figures 4 and 5, respectively. The colour bar shows the detection significance⁵ values computed from the maps and the strong negative values are deep side lobes near the signal. We note that the α method is less sensitive to these choices, offering a more objective and robust assessment of the detection significance. For both datasets, the emission signal for the entire atmosphere is strong and directly observable on the CCF map (see Fig. 6). In addition, we also applied the cross-correlation method to search for ¹³CO, Mg, CO₂, OH, HCN, Si and Ca but no significant detections are found in the K_p–Δυ_sys maps of the species (see Appendix C).

We emphasise that the detection significance values computed from the maps can change when cross-correlating with the same models due to the arbitrary selection of the noise region. Therefore, we followed the method outlined in Gibson et al. (2020) to evaluate the detection significance. We first computed the log-likelihood distribution (Sect. 4.1) based on the model scale factor, α, conditioned on the optimal values of all atmospheric parameters. Next, we calculated the mean of this distribution and divided it by the standard deviation to determine the “alpha” detection significance. This value indicates how many standard deviations the maximum signal deviates from an alpha of 0, corresponding to a null signal. Because the alpha detection significance is based on likelihood, the results are more objective, less arbitrary (i.e. require fewer assumptions) and provide a more conservative representation of the data. The penultimate panel of Fig. 4 shows a plot of the detection significance computed from the K_p–Δυ_sys maps and the alpha detection significance of the detected species for the datasets.

4 Atmospheric retrievals

While cross-correlation allows us to obtain detections of chemical species in the atmospheres of exoplanets, it does not provide quantitative information about the atmospheric properties. Additionally, this technique does not allow direct comparisons between observations and atmospheric models. However, recent advances in high-resolution Bayesian methods (e.g. Brogi & Line 2019; Gibson et al. 2020) have allowed us to “map” crosscorrelation values of a given atmospheric model to a likelihood value. This will enable us to perform atmospheric ‘retrievals’ that help constrain properties such as the temperature profile, chemical and/or elemental abundances, and subsequent products like the atmospheric C/O ratio and metallicity. We apply the same retrieval framework described in Ramkumar et al. (2023), which follows the approach of Gibson et al. (2020, 2022), to estimate the chemical composition and thermal structure of MASCARA-1b. We briefly outline our retrieval setup and present our results in the following sections.

4.1 Retrieval framework

To perform our retrieval analysis on the datasets, we computed our forward model, mi , by taking into account the atmospheric model parameters, θ = {κ_IR, γ, T_irr, T_int, P_cloud, χ_species×N_species}, the broadening kernel width {W_conv}, planetary velocities {K_p , Δυ_sys}, and the scaling factors for the model ({α}) and noise ({β}). The forward model was then filtered using the procedure outlined in Gibson et al. (2022). Finally, the log-likelihood for a given set of model parameters is computed as follows: $$\ln ℒ=-N\ln \beta -\frac{1}{2}{\chi }^2$$(4)

with $${\chi }^2={\displaystyle \sum _{i=1}^N}\frac{{\left(f_i-\alpha m_i\right)}^2}{{\left(\beta {\sigma }_i\right)}^2}$$(5)

where, f_i and σ_i are the data and the standard deviation, respectively. Assuming uniform priors (defined over linear or logarithmic parameter spacing, as listed in Tables 3 and 4), we compute the log-posterior by computing the log-prior and then adding the log-likelihood (i.e. Eq. (4)), which was folded into a Markov chain Monte Carlo (MCMC) framework to retrieve posterior distributions of the model parameters for each dataset. We use a Differential-Evolution Markov Chain (DEMC) (Ter Braak 2006; Eastman et al. 2013), running an MCMC chain with 192 walkers, with a burn-in length of 200 and a chain length of 300. This setup results in 96 000 samples, of which 38 400 are discarded. To assess convergence, we apply the Gelman-Rubin statistic (Gelman & Rubin 1992) after splitting the chains into four separate groups.

For our retrieval analysis, we focus on the chemical species detected in the day-side atmosphere of MASCARA-1b. Therefore, the number of atmospheric species under consideration, N_species , is three (i.e. CO, H₂O, and Fe). However, we note that this parameter may vary in different MCMC runs as we also include other C- and O-bearing species for the sake of completeness, in which case, the value of N_species is mentioned within parentheses (e.g. Sect. 5.3). We fixed T_int to be 100 K and α to be 1 in the model fits, and although α is fixed in the retrievals, it can be easily computed within the likelihood mapping, post-MCMC. That is, defining alpha as a free parameter in the likelihood rather than incorporating it into the model means we can quickly compute the conditional likelihood as a function of alpha from a fixed cross-correlation map. The last panel of Fig. 4 shows the conditional likelihood distribution of α of the detected species for the pre, post and pre+post datasets.

4.2 Retrieval results

We employ two approximations to describe the atmospheric chemical composition of MASCARA-1b in our retrieval setup. The first approach involves a free retrieval of the mixing ratios for the species, while the second approach uses FastChem (Stock et al. 2018, 2022; Kitzmann et al. 2024) to calculate the abundances self-consistently. The pre-eclipse data are analysed in Ramkumar et al. (2023) and we refer the reader there for a detailed discussion. The retrieved T–P profiles for the two approximations are over-plotted in Fig. 7 and the corresponding posterior distributions are plotted in Figures D.1 and D.2, respectively.

Fig. 4 Detection of individual species in the atmosphere of MASCARA-1b. Top, second and middle panels: Kp–∆υsys maps for the pre-eclipse, post-eclipse and pre+post datasets. Penultimate panel: detection significance computed using two ways: α-method is less sensitive to deep side lobes and is more objective (see text and Sect. 4). Bottom panel: conditional likelihood distribution of α.

Fig. 5 Kp–Δυsys maps of the combined CO, H2O and Fe detection in the atmosphere of MASCARA-1b for the pre-, post- and pre+post eclipse datasets. The atmospheric model used for the cross-correlation has combined contributions from the species and a thermal inversion as discussed in Sect. 3.1. The white dotted lines mark the detection peak, and the colour bar shows the detection significance computed from maps. The corresponding alpha detection significance and the conditional likelihood distribution of α are plotted in Fig. 4.

Fig. 6 Cross-correlation map for CO+H2O+Fe model, tracing out the planet signal throughout the observed phases. The dashed lines on either side of the CCF trail show the expected planet radial velocity (offset by ±25 km s−1), and the horizontal dashed lines indicate the secondary eclipse. The grey space between the end of our pre-eclipse observations and the start of the eclipse denotes phases at which no data were taken.

4.2.1 Free retrieval

We begin the retrieval process with a “free-retrieval” paradigm that assumes constant mixing ratios with altitude and uses the parametric model from Guillot (2010). The retrieval parameters and their prior ranges for the datasets are listed in Table 3. For the pre-eclipse, post-eclipse, and joint analysis (pre+post), we place bounded constraints on the abundance of H₂O while obtaining lower limits on CO and Fe as their retrieved volume mixing ratios approach the upper bounds of the prior distribution⁶. We did not detect Fe in the post-eclipse and combined datasets but obtained its 1σ lower limit as log₁₀(χ_Fe) > −4.97 for the latter and did not constrain it in the post-eclipse data. Despite the difference in signal-to-noise between the two datasets, the retrieved parameters for the post-eclipse (shown in blue in Fig. D.1) are consistent with the pre-eclipse values, though with slightly lower precision. When we combine the two nights of CRIRES+ data (shown in black in Fig. D.1), the abundance constraints are consistent with the pre- and post-eclipse sequences. However, they are not more precise than those provided by the pre-eclipse data alone, suggesting that the post-eclipse data contributed little new compositional information. We conclude that the precision is dominated by some degeneracies in the model that are not constrained at those wavelengths. Additionally, it is possible that the planet’s radial velocity trail is not perfectly sinusoidal for all species, which could introduce limitations when combining datasets if a suboptimal model is used. Nevertheless, combining the datasets enabled precise constraints on the planetary velocities (K_p, Δυ_sys). Overall, the abundances derived from our free-retrieval setup are consistent across the two datasets (pre-eclipse and post-eclipse) as well as the joint analysis (pre+post), with CO dominating relative to H₂O on the day-side of MASCARA-1b. However, this does not necessarily imply a carbon-rich atmosphere, as other oxygen-bearing species may influence the interpretation of the C/O ratio (see Sect. 5.1).

Lastly, as noted in Sect. 2, we use 15 SYSREM iterations for our cross-correlation and retrieval analysis presented in this work. While this number is arbitrary, we re-run the analysis with different values of SYSREM for the post-eclipse and combined dataset. The results of this analysis are shown in Figs. E.1 and E.2. We find that using N = 5,10 and 20 iterations give us consistent results for the post-eclipse dataset, whereas N = 20 filtered out some of our exoplanet signal for the pre-eclipse (Ramkumar et al. 2023, Fig. 10) and combined datasets.

4.2.2 Chemical equilibrium

For our retrievals assuming chemical equilibrium, we find that the H₂O, H₂ and H abundances change quite drastically with altitude in this temperature regime (e.g. right panel of Fig. 7). Therefore, instead of retrieving the individual gas volume mixing ratios, we used FASTCHEM (version 3.1.1, Stock et al. 2018, 2022; Kitzmann et al. 2024) to calculate the abundances of chemical species. This chemical regime is preferred over a well- mixed model, as the latter is known to break down for ultra-hot Jupiters, where strong altitude-dependent variations in composition occur. The results of this analysis (i.e. assuming chemical equilibrium) are shown in Fig. D.2, and the marginalised distributions for each parameter and dataset are summarised in Table 4. We placed bounded constraints on the planetary velocities and the T–P profile parameters for the two datasets and the joint analysis. The retrieved temperature profiles under this setup are over-plotted in Fig. 7 (middle panel). Similar to the free-retrieval setup, constraints on the retrieved parameters for the post-eclipse data remain consistent with the pre-eclipse values, although less precise. In addition, combining the pre-and post-eclipse data provides consistent constraints, albeit weighted more in favour of the pre-eclipse data that appears to be driving these inferences. The retrieved values for the T–P profile parameters between the two retrieval frameworks remain consistent for the individual datasets and the joint analysis (pre+post) within ≈1.1σ.

Lastly, in our retrieval framework, H⁻ bound-free and free-free opacities are included in the chemical equilibrium retrievals but not the free-chemistry retrievals. While the former is less relevant in the K-band, omitting free-free opacity in the free-chemistry retrieval does not significantly affect our results. To verify this, we re-computed our best-fit model emission spectra (assuming chemical equilibrium) with and without the H⁻ free-free opacity (and CIA) and found that the relative flux difference is ∼10⁻⁵. This indicates that although H⁻ free-free opacity is present, its impact on the overall emission spectra in the observed K-band wavelength range (1921–2472 nm) is minimal, suggesting that the omission of H⁻ opacities in the free retrievals does not significantly affect the model outcomes for this wavelength range.

Fig. 7 Left and middle panel: retrieved T–P profiles (Guillot 2010) for the pre-eclipse (red), post-eclipse (blue) and pre+post (black) datasets under the free-retrieval and chemical equilibrium setups. The red, blue and grey shading marks the 1σ and 2σ recovered distribution computed from 10 000 samples from the MCMC for both regimes. Right panel: the atmospheric structure for the post-eclipse data from the best-fitting chemical model. The volume mixing ratio profiles for the continuum and detected species are shown as dashed lines (calculated using FASTCHEM), with the parametric T–P profile shown as a solid blue line.

5 Discussion

5.1 Chemistry of MASCARA-1b

The composition (elemental abundances) of exoplanets, especially of highly irradiated gas giants (e.g. UHJs), can provide clues about the planetary formation and migration histories within the protoplanetary disk (Öberg et al. 2011) via diagnostics like the bulk metallicity and the C/O ratio. Here, we infer the chemistry of MASCARA-1b in terms of the atmospheric C/O and metallicity obtained from our retrieval frameworks. These quantities can be derived from the inferred gas abundances or fitted directly from chemical equilibrium model fits. The former becomes invalid for UHJs due to significant altitude-dependent variations in chemical composition (e.g. Fig. 7), highlighting the need to incorporate vertical chemical equilibrium when interpreting atmospheric properties.

Nonetheless, to assess potential biases introduced by the well-mixed assumption, we compute the free-retrieval based C/O and metallicity of MASCARA-lb derived from our retrieved gas abundance as follows: $$\text{C}/\text{O}=\frac{n_{\text{C}}}{n_{\text{O}}}=\frac{n_{\text{CO}}}{n_{\text{CO}}+n_{{\text{H}}_2\text{O}}}$$(6)

and $$[\text{M}/\text{H}]={\log }_{10}\left[\frac{2n_{\text{CO}}+n_{{\text{H}}_2\text{O}}+n_{\text{Fe}}}{2n_{{\text{H}}_2}{\left[\left(n_{\text{O}}+n_{\text{C}}+n_{\text{Fe}}\right)/n_{\text{H}}\right]}_{\odot }}\right]$$(7)

We assume that CO and H₂O are the dominant carbon- and oxygen-bearing molecules in the atmosphere of MASCARA-1b and compute the planetary metallicity by normalising the elemental abundances relative to hydrogen (the numerator in Eq. (7)), relative to that in the Sun (the denominator in Eq. (7)). The solar elemental abundances are taken from Asplund et al. (2009) and the median values of these derived quantities for the datasets are listed in Table 3. The derived metallicities for the pre-eclipse and combined datasets are consistent with the solar value within 1σ and are broadly consistent with the stellar value within their uncertainties. For the post-eclipse data, despite the strong CO detection, the derived [M/H] appears lower than that of the pre-eclipse. This difference is likely driven by the absence of detectable iron in the post-eclipse data, combined with the slightly weaker H₂O detection and/or constraint. However, the pre- and post-eclipse metallicities remain statistically consistent.

Looking at the derived C/O ratios, we obtain very precise constraints with $\text{C}/\text{O}={0.98}_{-0.02}^{+0.01}$, ${0.97}_{-0.02}^{+0.01}$ and ${0.98}_{-0.01}^{+0.01}$ for the pre-, post- and combined datasets. While there are no reported/measured values of the C/O ratio for the host star, MASCARA-1, the derived values are all markedly super-solar and are only constrained by our detection of CO and H₂O, assuming that every bit of C and O is within these molecules. However, this assumption could be incomplete. Qualitatively, if a significant portion of the oxygen is present in other species, such as OH or atomic oxygen, due to thermal dissociation of H₂O, it could influence the C/O ratio we are measuring. OH is expected to be present on the day side of MASCARA-1b, as it has been detected in colder planets such as WASP-121b and WASP-76b (Smith et al. 2024b; Landman et al. 2021; Gandhi et al. 2024; Weiner Mansfield et al. 2024). However, since K-band CRIRES+observations are not sensitive to OH lines and atomic oxygen is not observable in the near-infrared, our current retrievals may not fully account for these species. Therefore, we included potential C- and O-bearing species that exhibit sharp spectral lines in this wavelength setting but were not detected (e.g. OH, ¹³CO, HCN and CO₂; here N_species = 7), in our free-chemistry retrieval and find that the inclusion of these species does not change the derived C/O, which remains super-solar. Additionally, similar to Brogi et al. (2023), if we account for missing oxygen by assuming atomic O with a VMR of 10⁻³ we indeed obtain a lower C/O≈0.6, although we are unable to provide a precise estimate of uncertainty for this value. Overall, our free-retrieval-based results strongly favoured a H₂O-depleted model over one with a similar abundance of CO and H₂O, thereby driving the C/O towards 1 resulting in precise constraints.

While the free-chemistry model offers flexibility by allowing the gas volume mixing ratios to be retrieved individually, the assumption of a well-mixed atmosphere does not accurately reflect the vertical variation in chemical composition, particularly for highly irradiated gas giants. This can lead to biases in the elemental abundance estimates (e.g. Brogi et al. 2023; Ramkumar et al. 2023). Therefore, we fit directly for the C/O ratio and metallicity derived from chemical equilibrium fits (see Table 4 and Fig. D.2) using FastChem. This results in a $\text{C}/\text{O}={0.68}_{-0.22}^{+0.12}$, ${0.75}_{-0.17}^{+0.11}$ and ${0.74}_{-0.14}^{+0.10}$ for the pre-, post-, and pre+post eclipse datasets, respectively. The retrieved values are all consistent with the solar value within ≈1.1–2σ and the combined dataset provides slightly tighter constraints on the retrieved values than the individual datasets. In addition, the post-eclipse seems to constrain the C/O ratio better than the pre-eclipse, though the difference is negligible. This could be because the uncertainties remain consistent across all datasets, as suggested by our injection-recovery tests (see Appendix A) where the injected and recovered values for both a constant-with-altitude abundance and equilibrium chemistry modelling setups are in good agreement. Nevertheless, we speculate that the uncertainties in these parameters are likely dominated by degeneracies in the model and the wavelength range of the observations, as well as the limitations of our 1D models that do not account for 3D effects such as spatial inhomogeneities in the planet’s atmosphere. These factors may contribute to the consistency of the C/O constraints despite the lower data quality of the post-eclipse sequence. Moreover, although the metallicity prior in our chemical equilibrium retrieval is relatively narrow, once we broaden the metallicity constraints (from 0.001 to 1000x solar), the retrieved C/O ratio remains consistent for the post-eclipse data and the joint analysis. However, we find a strong correlation between C/O and metallicity for the pre-eclipse data (see Fig. D.3), leading to inconsistencies in the retrieved value. While this does not alter our overall conclusions, it highlights an important degeneracy we aim to explore further. Lastly, similar to our free-chemistry analysis, for completeness, we also included opacities for OH, CO₂ and HCN in our equilibrium chemistry retrievals and find that the inclusion of these species did not change the retrieved C/O value, which remains consistent with the solar value within ≈1σ.

Overall, our results indicate that the atmosphere of MASCARA-1b has a C/O ratio consistent with the solar value, with the chemical equilibrium retrieval providing more realistic and conservative constraints. While it would be interesting to frame the planet’s chemistry in terms of the chemistry of its parent star, the C/O ratio of MASCARA-lb alone is inconclusive in determining its formation as there are no reported or measured values of the C/O ratio for the host star, MASCARA-1. Figure 8 highlights the diversity of C/O ratios among hot and ultra-hot Jupiters using ground-based high-resolution observations. Most planets in the sample exhibit C/O ratios near or above the solar value, with the observed spread in the reported values reflecting the diversity of formation and migration pathways. MASCARA-lb adds to this population with ratios retrieved from pre-, post-, and combined datasets, all consistent with the solar value, making it a valuable addition to the sample of UHJs with well-constrained C/O ratios. While future observations could further refine these constraints, offering new insights into the planet’s chemical inventory and its formation and migration history, it also underscores the challenges in interpreting these ratios in isolation since carbon and oxygen measurements alone are insufficient to determine the origin of a gas-giant planet due to the intrinsic degeneracies in the planet-building gas and solid compositions (Mordasini et al. 20l6).

Recent studies by Lothringer et al. (2021) and Chachan et al. (2023) elucidated the diagnostic power of refractory elements such as Fe, Mg, Si, and Ca, and introduced the concept of the refractory-to-volatile ratio – the abundance of refractory elements like iron relative to volatile species such as CO and H₂O – as a valuable tool for tracing the formation and migration pathways of hot and ultra-hot Jupiters. This ratio provides a framework for predicting how a planet’s solid-to-gas ratio and accretion location within the protoplanetary disk influence its atmospheric composition and has very recently been applied to high-resolution observations of the ultra-hot Jupiter WASP-121b, enabling the exploration of potential formation and migration pathways for such planets (Pelletier et al. 2025; Smith et al. 2024b). In the case of MASCARA-lb, our pre-eclipse detection of iron allows us to estimate the refractory-to-volatile ratio within this framework. Assuming that the Fe signal is not spurious (see Sect. 5.2), we derive the refractory-to-volatile ratio to be $[\text{R}/\text{V}]={0.32}_{-1.49}^{+1.17}$ from our free-chemistry retrievals. However, it is important to note that attributing the entire refractory inventory to Fe alone or the volatile inventory solely to CO and H₂O likely introduces significant bias, as it excludes contributions from other important rock- and ice-forming elements. The non-detections of other refractory (e.g. Mg, Si) and potential volatile (e.g. OH and atomic O) species, partly due to the scarcity of strong spectral lines in the observed wavelength range, could lead to an incomplete representation of refractory and volatile content. Furthermore, our retrieved ratios are derived assuming a constant-with-altitude modelling regime. For these reasons, we choose not to over-interpret the derived refractory-to-volatile ratio and emphasise the need for future observations to better constrain the inventory of both refractory and volatile species in MASCARA-1b.

Nevertheless, UHJs present a tremendous opportunity to diagnose planet formation histories by revealing a wealth of volatile and refractory species in their atmospheres. Combining optical and near-infrared high-resolution emission spectroscopy datasets should provide interesting constraints on refractory-to-volatile elemental ratios for irradiated gas giants.

Fig. 8 Planetary mass versus C/O ratio for hot and/or ultra-hot Jupiters with constrained abundances of C- and O-bearing species from ground-based high-resolution observations. Our retrieved C/O values for MASCARA-1b are shown as diamonds (labelled M-1b). Values for other planets are from previous studies using either free-chemistry retrievals or chemical equilibrium models (see Table F.1). The dashed line indicates the solar C/O value, and Jupiter is included for reference (from Weiner Mansfield et al. 2024). To reduce clutter, the labels for planets excluding MASCARA-1b and WASP-18b are shown in the inset, which provides a zoomed-in view. For clarity, small mass offsets are applied to WASP-121b, WASP-76b, WASP-189b and MASCARA-1b.

Fig. 9 Kp−Δυsys map of the pre-eclipse Fe signal, with the detection significance computed using a region where KP>210 km s−1 and Δυsys< − 40 km s−1, avoiding the planetary signal near the expected Kp and Δυsys. While choosing a different region reduced the detection significance to ≈5σ, Fe remains significantly detected, with the alpha detection significance constrained to be non-zero at more than 4.2σ, supporting the robustness of the Fe signal.

5.2 Revisiting the iron detection

The detection of iron in the pre-eclipse dataset (Ramkumar et al. 2023) presents an interesting, yet somewhat ambiguous, result. While a prominent peak is found around the consistent K_p and Δυ_sys in the detection map (see top panel of Fig. 4), albeit a little shifted, we cannot definitively conclude that it is a strong detection of Fe based on our CCF analysis alone. Specifically, the presence of significant negative peaks (e.g. the −6σ feature at a Δυ_sys≈23 km s⁻¹ and ≈90 km s⁻¹) in the CCF map raises reasons to question the robustness of this detection. We also note that the reported detection significance does not have to be 8σ because the choice of the region used to calculate this is arbitrary and can result in varied values (e.g. Fig. 9). Furthermore, when we perform a free retrieval analysis with a model containing Fe alone, we cannot constrain it, likely because the approach does not account for other species (e.g. CO, H₂O) that can mask certain spectral lines or alter the effective line strengths that could potentially weaken the Fe signal. Nonetheless, the presence of a positive CCF peak in our K_p−Δυ_sys map combined with constraints in our retrieval analysis with all the species seems to indicate that Fe signal in the pre-eclipse data is genuine. In addition, the detection significance derived from the likelihood (alpha value) is constrained to be non-zero at more than 4.2σ (see penultimate panel of Fig. 4), which further boosts our confidence in its detection.

Assuming the pre-eclipse detection is real, the absence of iron in the post-eclipse data becomes particularly intriguing. To explore this, we performed injection-recovery tests with the best-fitting Fe model, injected at a negative K_p. For the pre-eclipse data, the injected and recovered values are in agreement (see Appendix G), with the detection significance computed from the K_p−Δυ_sys map and the alpha detection significance being ≈10.92σ and ≈8.31σ, respectively. However, we note that the ‘actual’ signal – corresponding to the Fe detection originally observed in the CCF analysis – is not visible in the recovered map. This absence may result from the injected signal adding noise, which can weaken the effective signal, particularly pronounced for weaker detections, such as iron, which has an alpha detection significance of only ≈4.4σ. In contrast, we do not detect iron in the post-eclipse dataset and performing an injection and recovery test yields detection significance and alpha detection significance values of ≈7.73σ and ≈7.53σ, respectively. These values are not significantly different from those of the pre-eclipse dataset; however, they are slightly lower, indicating that the data quality is somewhat reduced, and the sensitivity to iron is consequently diminished. This suggests that the non-detection in the post-eclipse phase could result from the weaker Fe signal falling below the detection threshold due to the slightly poorer signal-to-noise ratio. The non-detection could, therefore, be attributed to the data quality rather than an inherent absence of iron in this phase sequence.

Lastly, to assess whether the inclusion of Fe impacts our inference of other species, we performed retrievals identical to those described in Sect. 4.2, but excluding Fe from the model (i.e. assuming CO and H₂O as the only opacity sources). The results are shown in Appendix H. The retrieved temperature profiles with and without Fe are consistent (Fig. 10), suggesting that including iron does not significantly alter the atmospheric structure of MASCARA-1b. This consistency indicates that the model is not artificially forcing iron into the atmosphere by substantially affecting the T–P profile, which might occur if the iron detection were spurious. Furthermore, the absolute abundances of CO and H₂O remain unchanged across both models. Our chemical model also predicts iron to exist at these temperatures (right panel of Fig. 10), and the corresponding retrieved T–P profiles remain consistent. Additionally, the abundance retrieved from our free-chemistry retrieval analysis agrees with the predictions from our chemical model. Thus, even if we were to ignore iron, we would still expect to see it in the atmosphere at these temperatures.

5.3 Signatures of atmospheric dynamics

We looked for shifts in the retrieved K_p and Δυ_sys as proxies for possible atmospheric dynamics, searching for any velocity offsets between the detected species (e.g. Brogi et al. 2023; Cont et al. 2021) or between the pre- and post-eclipse datasets (e.g. Pino et al. 2022). Figure 11 illustrates the 2D posteriors of the detected species for the pre-, post- and pre+post-eclipse datasets. We first search for shifts between the detected species in the individual datasets. For the post-eclipse and pre+post data, CO and the combined species (CO+H₂O+Fe) are approximately centred around the expected rest-frame velocity (Δυ_sys). However, H₂O appears slightly shifted in the 2D parameter space but when we examine the 1D marginalised distributions for H₂O across the datasets, the results remain consistent with no significant shift. We note that correlations between K_p and Δυ_sys in the 2D posteriors can be obscured when marginalised, making species-specific shifts challenging to interpret. While H₂O seems to be shifted in the 2D space, we have not performed statistical tests to quantify its significance. Therefore, we remain cautious in over-interpreting these shifts based purely on posterior distributions. However, in our analysis of the pre-eclipse data (Ramkumar et al. 2023), tentative evidence for shifts is seen for the Fe feature in both K_p and Δυ_sys. Specifically, while the cross-correlation analysis shows that the peak (mean) of the map is red-shifted in Δυ_sys for Fe, the 2D posteriors from the retrieved samples suggest that Fe in general is blue-shifted relative to CO and H₂O. If the iron signal is real, this blue shift may indicate that iron is tracing a different dynamical region of the atmosphere, which may also be linked to its non-detection post-eclipse. In addition, to investigate potential differences in dynamics between species, we examine the planetary velocity (υ_p) as a function of the orbital phase for the detected species in the pre-, post-, and combined datasets (Fig. 12). This plot shows that while the CO and H₂O signals exhibit consistent trails across the observed phases, the iron signal appears slightly offset during the earlier phases of the pre-eclipse dataset (see also Fig. I.2 for a zoomed-in version), with the trace overlapping towards the end (i.e. close to the secondary eclipse), where it appears to be more strongly detected.

To further explore these trends and potential longitudinal variations in the atmospheric parameters, we split the pre- and post-eclipse datasets into two subsets based on the orbital phase: phases near quadrature and phases close to the secondary eclipse (see Appendix I). A joint retrieval performed on these subsets indicates that iron is better constrained during the phases close to the secondary eclipse, as corroborated by the K_p−Δυ_sys maps (Appendix I), and aligns with the trace plot. This raises two possibilities: (i) the stronger detection of iron near the eclipse may be due to localised atmospheric phenomena on the planet, or (ii) when combined with the already weak detection (based on the alpha-detection significance), pre-processing steps like SYSREM is likely getting rid of the signal more efficiently at the lower and/or shorter phases when the planet is more static that could be weakening the iron signal at this phase subset. Both scenarios merit further investigation to disentangle their contributions. Nonetheless, it is intriguing that Fig. 12 combined with the phase-resolved retrieval shows a slightly different trail for iron, which might imply that there are some 3D effects at play that may lead to the non-detection or could in principle lead to a lower signal post-eclipse. That said, we exercise caution in interpreting these features without detailed statistical tests and, therefore, only speculate – underscoring the need for follow-up studies using high-resolution spectroscopy and 3D atmospheric models to explore these phenomena fully.

Following this, we also checked for offsets between the pre-and post-eclipse sequences for each species and found no signatures of atmospheric dynamics (e.g. right panel of Fig. 11). While it is known that different line lists can influence our ability to detect planetary signals, discrepancies in line strengths and positions might affect the S/N (e.g. Webb et al. 2020; Nugroho et al. 2021; Rafi et al. 2024) or retrieved atmospheric properties. Nevertheless, previous analyses, such as that of Gandhi et al. (2020), confirm that the line lists for CO and H₂O are appropriate for high-resolution studies up to spectral resolutions of R= 100 000. To assess the robustness of our detections and to determine whether the choice of line list impacts our retrieval results, we performed cross-correlation analysis using the POKAZATEL line list from Exomol for H₂O (Polyansky et al. 2018). We found no significant differences between the retrieved parameters or the detection significance computed from the conditional likelihood distribution (see Appendix J) when comparing the POKAZATEL line list with HITEMP 2010 (Rothman et al. 2010), which we use for the results presented in Section 3.3. Additionally, we looked for shifts in the retrieved K_p and Δυ_sys for H₂O alone using two line lists and did not observe any significant shifts in velocity, with all the datasets remaining centred around the expected systemic velocity (see Fig. 13).

Fig. 10 Left and middle panel: retrieved T–P profiles (Guillot 2010) for the pre-eclipse data with the inclusion of iron (i.e. CO+H2O+Fe; red) and without iron (i.e. CO+H2O; black) under the two retrieval setups outlined in Sect. 4.1. The red and grey shading marks the 1σ and 2σ recovered distribution computed from 10 000 samples from the MCMC. Right panel: the atmospheric structure for the pre-eclipse dataset from the best-fitting CO+H2O only model. The volume mixing ratio profiles for the continuum and detected species are shown as dashed lines (calculated using FASTCHEM), with the parametric T–P profile shown as a solid blue line.

Fig. 11 Searching for signatures of atmospheric dynamics. Left and middle panel: posterior distributions for the planetary velocities of CO, H2O and the combined spectrum (i.e. CO+H2O+Fe) for the post-eclipse and pre+post datasets. Right panel: marginalised posterior distributions of the combined atmosphere for the pre- and post-eclipse datasets.

Fig. 12 Retrieved planetary orbital velocity (υp) for the detected species from our observations, computed from 10 000 random samples of the MCMC. Dashed lines mark the median, with shaded regions showing 1σ and 2σ contours. For clarity, only the median values (dotted lines) are plotted for the pre+post data.

5.4 Phase and/or time dependence

Exoplanet atmospheres are inherently 3D structures characterised by dynamic processes, including circulation, chemical mixing, and temperature variations. This is especially true for tidally locked, highly irradiated planets, such as ultra-hot Jupiters, where significant day-to-night temperature gradients can occur. However, our retrieval frameworks currently utilise 1D forward models that assume a globally uniform atmosphere where the atmospheric signal is constant over time. Therefore, it is necessary to exercise caution when interpreting the atmospheric properties retrieved from ID retrievals, as they may not accurately represent the global chemistry of the planet, particularly in cases where phase-dependent variations are expected. For example, recent emission spectroscopy analyses of the UHJ, WASP-33b, detected a phase dependence found via the model scaling parameter, α (e.g. Herman et al. 2022; van Sluijs et al. 2023), and reported that a larger scaling is required to best model the observations after the secondary eclipse. Similarly, Pino et al. (2022) explored atmospheric dynamics in KELT-9b by tracking the planet’s velocity as a function of time using phase-resolved CCFs. Such analyses have demonstrated the power of high-resolution spectroscopy phase curves in probing atmospheric dynamics, temperature gradients, and even potential “3D-ness”.

Our analysis of MASCARA-1b suggests possible variations in the strength of the cross-correlation signal over time, particularly in the pre-eclipse data (see Fig. 6). Notably, the strength of the CCF trail changes at the beginning and end of the observing sequences, with the lower half of the plot (i.e. the start of the observations) showing a diminished trail. This is most likely attributed to the SYSREM pre-processing since de-trending algorithms can remove significant parts of the planet signal. While our filtering technique (Gibson et al. 2022) accounts for these alterations by applying the same pre-processing to the forward model as we do to the data, we will always reduce the signal a little in a phase-dependent way with SYSREM and/or PCA. Nonetheless, the filtering approach works well within our retrieval frameworks, allowing us to recover consistent posteriors with reduced sensitivity to how aggressively the data is filtered until we reach a higher number of SYSREM passes (e.g. Ramkumar et al. 2023, Fig. 10). Additionally, there appears to be a slight dependence on the signal-to-noise, which may also contribute to the variation seen in the CCF. For instance, the pre-eclipse signal disappears towards the end of the sequence, coinciding with a significant increase in airmass and a decrease in S/N. This drop in signal at phases furthest from the eclipse is also expected due to the reduced visibility of the hot day side from the field of view. However, it does not drop off rapidly (see right panel of Fig. 1).

To determine whether the variation in the strength of the CCF trail is solely due to fluctuations in signal-to-noise or if it indicates a physical change in the planet’s signal, we chose not to fix the model scaling α, at 1. Instead, we treated it as a free parameter in our combined (pre+post) retrieval fits. This approach allowed us to examine potential time and/or phase dependence before considering a phase-dependent alpha implementation within the likelihood equation. Upon doing so, we find the retrieved values for α before and after the eclipse to be consistent, with ${\alpha }_{\text{pre }}={0.90}_{-0.08}^{+0.08}$ and ${\alpha }_{\text{post }}={0.95}_{-0.09}^{+0.09}$. Building on this, we further investigated potential phase-dependent trends by computing a likelihood map as a function of α using a combination of Equations (4) and (5). We began with our cross-correlation map (ϕ vs Δυ_sys) and calculated it as a function of α (expansion of Eq. (5)). Next, we Doppler-shifted it to the planet rest frame using Eq. (2) (conditioned on the optimal values of all atmospheric parameters). Instead of summing over time, we calculated it for each phase point. This was then used to compute a log-likelihood map as a function of alpha for every frame, which was converted into a likelihood to obtain a conditional likelihood distribution of α, similar to the method outlined in Sects. 3.3 and 4.1 (see also Fig. 4), but for every frame. Lastly, we determined α and its corresponding uncertainties for each phase by calculating the mean and standard deviation of the distribution. A plot of α versus orbital phase for the pre- and post-eclipse datasets is shown in Figure 14. The α values are centred around 1, with no significant trend seen between the two datasets, suggesting that the line contrast is constant over time. However, given the degeneracy between the temperature gradient and chemical abundances (+ continuum), it does not necessarily mean the atmosphere is constant over time. Nonetheless, more sophisticated methods and high-resolution observations may be required to better constrain and probe the phase dependence of the model scaling factor.

Fig. 13 2D posteriors of H2O for the two datasets using the HITEMP and POKAZATEL line lists (see text).

Fig. 14 Model scaling factor (α) versus the orbital phase for MASCARA-1b. For visualisation purposes, frames from the individual dataset have been binned together. The α values (blue dots) and their uncertainties for each frame are derived from the conditional likelihood distribution (see text). The black markers show the retrieved median α for the individual datasets. The dark blue points denote outliers corresponding to cloud-affected frames, located within the regions marked by vertical grey bands.

5.5 Comparison with phase curve measurements

Recent analyses of Spitzer’s 4.5 µm phase curve of MASCARA-1b (Bell et al. 2021; Dang et al. 2025) have measured day-side and night-side temperatures of 3000±100 K and 1300±300 K, respectively, showing a fairly inefficient recirculation, which is consistent with expectations for ultra-hot Jupiters. Building on these findings, we calculated equilibrium temperatures for different heat redistribution scenarios given by the heat redistribution factor, f, as follows: $$T_{\text{p},\text{eq}}=T_{\ast }{\left(\frac{R_{\ast }}{a}\right)}^{\frac{1}{2}}{\left[\frac{f}{4}\left(1-A_{\text{B}}\right)\right]}^{\frac{1}{4}}$$(8)

where T_*, R_*, a and A_B denote the stellar effective temperature, stellar radius, semi-major axis and the bond albedo, respectively. The heat redistribution factor, f, accounts for the fact that the absorbed energy may not all be re-radiated uniformly over the planet’s surface, with f=1 corresponding to uniform distribution (e.g. full heat redistribution across the planet), f=2 for re-radiation over half the surface of the planet (e.g. day-side only re-radiation) and $f=\frac{8}{3}$ corresponding to instant re-emission (Lambertian sphere). For an $A_{\text{B}}={0.32}_{-0.11}^{+0.09}$ (Dang et al. 2025), we calculate the equilibrium temperature to be ≈2376±35 K, 2826±41 K and 3036±44 K for f= 1, f=2 and $f=\frac{8}{3}$, respectively. On comparing these with the irradiation temperature from our retrieval analysis, we find that while the results are broadly consistent with f=1, particularly for the free-chemistry setup, given the uncertainties, our results also overlap with f=2, allowing for consistency with any of the heat redistribution scenarios within ≈1σ. However, our uncertainties on the T–P profile remain consistent with the Spitzer phase curve analysis. Additionally, in the context of high-resolution spectroscopy, our analysis of the phase-dependent scaling factor α – which represents variations in line contrast as a function of the orbital phase – reveals no significant trends. This further suggests that we do not observe any substantial changes in the atmosphere of MASCARA-1b as it rotates, likely explaining the consistency in our retrieved parameters between the pre- and post-eclipse datasets. It is also worth noting that the near-symmetric phase coverage of the pre- and post-eclipse datasets may further contribute to the observed consistency. For instance, if the planet’s hotspot exhibits minimal longitudinal offset or its spatial variations are symmetric relative to the substellar point, retrievals from both phases could appear similar even with strong spatial variation. This consistency between pre- and post-eclipse also aligns with the Spitzer phase curve observations, which do not show any significant phase offset (Dang et al. 2025).

6 Conclusions

In this work, we have presented high-resolution emission spectroscopy observations of the UHJ MASCARA-1b with CRIRES+ in the K-band, covering the post-eclipse phases of the planet’s orbit. We applied the standard cross-correlation technique and a Bayesian retrieval framework that includes both pre-eclipse and post-eclipse datasets and learned the following about the thermal and chemical properties of the planet:

• In addition to the CRIRES+ pre-eclipse phase-curve observations presented in Ramkumar et al. (2023), we present here an additional night taken two years apart, capturing the planet’s day-side emission in the post-eclipse phases. Using the traditional cross-correlation method, we detected the presence of CO and H₂O in the post-eclipse data, with the features also identified at high significance in the joint (pre+post) analysis (Sect. 3.3);

• The molecular species are seen as emission signatures through which we also confirmed the presence of a thermal inversion layer in the day-side atmosphere of MASCARA-1b. No significant signature of the chemical species CO₂, ¹³CO, OH, Mg, HCN and CH₄ could be detected in our pre-, post-, and pre+post-eclipse data (Sect. 3.3);

• The likelihood framework introduced in Gibson et al. (2020, 2022) was applied to obtain quantitative constraints on the atmospheric properties, such as the parametric T–P profile, planetary velocities (K_p, Δυ_sys), abundances of the species as well as the atmospheric C/O ratio while simultaneously marginalising over the noise properties of the dataset (Sect. 4);

• We implemented two different approximations in our retrieval framework to describe the chemical composition of MASCARA-1b’s atmosphere (Sect. 4.1). While our free-chemistry retrievals that assume a well-mixed atmosphere, indicate a super-solar C/O ratio, incorporating a chemical equilibrium model results in a C/O of ${0.75}_{-0.17}^{+0.11}$ for the post-eclipse data and ${0.74}_{-0.14}^{+0.10}$ for the joint analysis, both consistent with solar values within ≈1.1–2σ (Sect. 4.2). The chemical equilibrium model, which accounts for altitude-dependent variations in composition, is the preferred model for ultra-hot Jupiters, given their complex atmospheric dynamics;

• Despite the difference in signal-to-noise between the pre-and post-eclipse datasets, the retrieved model parameters, under both free-chemistry and chemical equilibrium setups, are consistent with the pre-eclipse values, though less precise. Combining the two nights of the CRIRES+ data enabled slightly tighter constraints on the retrieved parameters than the individual datasets. Overall, the retrieved values for the T–P profile parameters between the two retrieval frameworks remain consistent across the datasets, falling within ≈1.1σ (Sect. 4.2);

• We revisit the detection of the Fe signal in our pre-eclipse data (Sect. 5.2) and found no signatures of atmospheric dynamics in the post-eclipse and combined datasets (Sect. 5.3);

• A phase-dependent analysis of the scaling factor α showed no significant variation between the pre- and post-eclipse phases, with the retrieved median a values centred around ≈1 (Sect. 5.4);

• We find no strong evidence for 3D effects in the atmosphere of MASCARA-1b, and given the large uncertainties on the T–P profile, our findings can be explained by any of the f values and are not in conflict with the Spitzer phase curve (Sect. 5.4). Nonetheless, we note that our observations cover only a subset of the planet’s rotation, which may limit our sensitivity to detecting longitudinal temperature variations or other 3D effects.

While the consistency between pre- and post-eclipse retrievals may not be particularly intriguing from a 3D perspective, it remains interesting that this agreement persists even though the observations were taken two years apart, with the post-eclipse data probing the atmosphere of MASCARA-lb after the planet had rotated by approximately ≈106 degrees compared to the pre-eclipse observations. Future high-resolution observations with broader phase coverage and/or complementary wavelengths could further enhance our understanding of longitudinal and temporal variability in ultra-hot Jupiter atmospheres. Such observations, in synergy with space-based facilities such as the JWST, could help provide critical insights into the chemical and thermal dynamics shaping these extreme planetary environments.

Acknowledgements

We thank the anonymous referee for careful reading of the manuscript and providing helpful comments. This work relied on the observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programs 107.22TQ.001 (PI: Gibson) and 112.260X.001 (PI: Gibson, dPI: Ramkumar). We are extremely grateful to the CRIRES+ instrument teams and observatory staff who made these observations possible. S.R. gratefully acknowledges support from the Provost’s PhD Award from Trinity College Dublin. N.P.G is supported by Science Foundation Ireland and the Royal Society in the form of a University Research Fellowship and S.K.N. is supported by JSPS KAKENHI grant No. 22Kl4092. We are grateful to the developers of the NumPy, SciPy, Matplotlib, iPython, corner, petitRADTRANS, FASTCHEM, and Astropy packages, which were used extensively in this work (Harris et al. 2020; Virtanen et al. 2020; Hunter 2007; Perez & Granger 2007; Foreman-Mackey 20l6; Mollière et al. 20l9; Stock et al. 2018, 2022; Astropy Collaboration 2022).

Appendix A Injection-recovery tests

Fig. A.1 Summary of our injection tests for pre-, post- and pre+post-eclipse datasets outlined in Sect. 3.2, with the 1D and 2D marginalised posterior distributions of each of our model parameters displayed. The red, blue and black posterior distributions represent the datasets, with the injected values highlighted by horizontal and vertical solid purple lines. The dotted lines show the retrieved median values. Upper right: the retrieved T–P profiles are computed from 10,000 random samples of the MCMC; the solid curves show the median profiles and the shaded regions show the 1σ and 2σ contours. The dashed black curve shows the injected T–P profile.

Fig. A.2 Same as Fig. A.1 but for a chemical model.

Appendix B Additional retrieval analyses including the cloud-affected frames in the post-eclipse data

Fig. B.1 Free-retrieval results for the post-eclipse data with (red) and without (blue) the inclusion of cloud-affected frames (Sect. 3.2). While excluding these frames slightly improved the Fe constraint, it is still a long-tailed distribution and did not impact the overall conclusions (Fe remains a non-detection).

Fig. B.2 Same as Fig. B.1 but assuming chemical equilibrium.

Appendix C Additional K_p−Δυ_sys maps

Fig. C.1 Kp−Δυsγs maps for the non-detected species. Left panels: emission model of each species with the wavelengths covered by CRIRES+ K2166 setting used in this work shown as vertical grey bands. Middle and right panels: the detection maps for the pre-ecl., post-ecl., and combined nights of the CRIRES+ data (Sect. 3.4). The high-resolution opacities were obtained from the following sources: 13CO, CO2 and OH from the HITEMP line list (Rothman et al. 2010), metal opacities from (Kurucz 2018) and HCN opacities from ExoMol database (Harris et al. 2006; Barber et al. 2014).

Appendix D Atmospheric retrievals of MASCARA-1b

Fig. D.1 Summary of our free-retrieval results for the combined atmosphere with the 1D and 2D marginalised posterior distributions of each model parameter displayed within a corner plot. The red, blue and black distributions represent the pre-, post- and pre+post eclipse datasets (Sect. 4.2).

Fig. D.2 Same as Figure D.1, but for the chemical retrieval described in Section 4.2.2.

Fig. D.3 Same as Figure D.2, but with a broader metallicity prior (from 0.00l to l000x solar).

Appendix E Consistency between SYSREM iterations on the post-eclipse and pre+post datasets

Fig. E.1 Summary of our free-retrieval results with the 1D and 2D marginalised posterior distributions of the parameters from the MCMC fit. The different colours show the samples obtained using four values for the SYSREM iterations (Sect. 4.2.1).

Fig. E.2 Same as Fig. E.1 but for the pre+post dataset.

Appendix F Ground-based high-resolution measurements of well-constrained C/O ratios

Appendix G Injection and recovery test for Fe

Fig. G.1 Results of the injection and recovery test for the pre-eclipse Fe (Sect. 5.2). The red and black posterior distributions represent independent sub-chains of the same MCMC chain, both converging to similar distributions. Upper right: the retrieved T–P profile computed from 10,000 random samples of the MCMC; the solid red curve shows the median profile and the shaded regions show the 1σ and 2σ contours. The dashed black curve shows the injected T–P profile.

Fig. G.2 Recovered Kp−Δυsys maps of iron for the pre- and post-eclipse datasets. The white dotted lines mark the peak of detection.

Appendix H Revisiting the pre-eclipse iron detection

Fig. H.1 Results of the free-retrieval for the pre-eclipse dataset with (red) and without (black) the inclusion of Fe as described in Sect. 5.2.

Fig. H.2 Same as Fig. H.1 but for a chemical model.

Appendix I Phase-dependent retrieval

Fig. I.1 Results of the joint (pre+post) phase-dependent free-retrieval as described in Sect. 5.2. The subsets of phases close to the quadrature are shown in red and phases close to the eclipse are shown in black. Top right: retrieved T–P profile.

Fig. I.2 Left: Kp−Δυsys maps for the iron signal using the phase-dependent retrieval setup. Right: The retrieved planetary orbital velocity (υp) of the detected species for the pre-eclipse data alone, computed from 10,000 random samples of the MCMC. The dashed lines show the median value, and the shaded regions show the 1σ and 2σ contours (see Sect. 5.3).

Appendix J Plot of alpha detection significance using different line lists for the pre-eclipse, post-eclipse and pre+post datasets

Fig. J.1 Plot of alpha detection significance computed from the conditional likelihood distribution (detailed in Sects. 3.3 and 4) using different line lists (see Sect. 5.3) shown for the pre-eclipse, post-eclipse and pre+post datasets. The labels ‘H10’ and ‘H19’ correspond to the HITEMP 2010 and HITEMP 2019 line lists, respectively.

References

All Tables

All Figures

	Fig. 1 Left: phase coverage of our CRIRES+ pre- and post-eclipse observations. Middle and right: observing conditions showing the variation of airmass and mean signal-to-noise ratio (S/N) with orbital phase. Vertical grey bands denote cloud-affected regions.
In the text

Fig. 2 Data pre-processing steps, applied to a single order: (a) raw spectra, (b) cleaned and blaze corrected spectra, and (c) residual spectra after division by the median spectrum and subtraction of the SYSREM model (shown weighted by the uncertainties). We note that the blaze correction does not work well for cloud-affected frames, resulting in some horizontal structure (residuals) in these regions (the last two panels). However, they only have a minimal impact on our retrievals.

In the text

	Fig. 3 Left panel: best-fit model spectra (constant VMRs) of MASCARA-1b with contributions from the chemical species shown individually. Middle and right panels: the parametric T–P profile from Guillot (2010) and combined contribution function.
In the text

	Fig. 4 Detection of individual species in the atmosphere of MASCARA-1b. Top, second and middle panels: Kp–∆υsys maps for the pre-eclipse, post-eclipse and pre+post datasets. Penultimate panel: detection significance computed using two ways: α-method is less sensitive to deep side lobes and is more objective (see text and Sect. 4). Bottom panel: conditional likelihood distribution of α.
In the text

Fig. 5 Kp–Δυsys maps of the combined CO, H2O and Fe detection in the atmosphere of MASCARA-1b for the pre-, post- and pre+post eclipse datasets. The atmospheric model used for the cross-correlation has combined contributions from the species and a thermal inversion as discussed in Sect. 3.1. The white dotted lines mark the detection peak, and the colour bar shows the detection significance computed from maps. The corresponding alpha detection significance and the conditional likelihood distribution of α are plotted in Fig. 4.

In the text

Fig. 6 Cross-correlation map for CO+H2O+Fe model, tracing out the planet signal throughout the observed phases. The dashed lines on either side of the CCF trail show the expected planet radial velocity (offset by ±25 km s−1), and the horizontal dashed lines indicate the secondary eclipse. The grey space between the end of our pre-eclipse observations and the start of the eclipse denotes phases at which no data were taken.

In the text

Fig. 7 Left and middle panel: retrieved T–P profiles (Guillot 2010) for the pre-eclipse (red), post-eclipse (blue) and pre+post (black) datasets under the free-retrieval and chemical equilibrium setups. The red, blue and grey shading marks the 1σ and 2σ recovered distribution computed from 10 000 samples from the MCMC for both regimes. Right panel: the atmospheric structure for the post-eclipse data from the best-fitting chemical model. The volume mixing ratio profiles for the continuum and detected species are shown as dashed lines (calculated using FASTCHEM), with the parametric T–P profile shown as a solid blue line.

In the text

Fig. 8 Planetary mass versus C/O ratio for hot and/or ultra-hot Jupiters with constrained abundances of C- and O-bearing species from ground-based high-resolution observations. Our retrieved C/O values for MASCARA-1b are shown as diamonds (labelled M-1b). Values for other planets are from previous studies using either free-chemistry retrievals or chemical equilibrium models (see Table F.1). The dashed line indicates the solar C/O value, and Jupiter is included for reference (from Weiner Mansfield et al. 2024). To reduce clutter, the labels for planets excluding MASCARA-1b and WASP-18b are shown in the inset, which provides a zoomed-in view. For clarity, small mass offsets are applied to WASP-121b, WASP-76b, WASP-189b and MASCARA-1b.

In the text

Fig. 9 Kp−Δυsys map of the pre-eclipse Fe signal, with the detection significance computed using a region where KP>210 km s−1 and Δυsys< − 40 km s−1, avoiding the planetary signal near the expected Kp and Δυsys. While choosing a different region reduced the detection significance to ≈5σ, Fe remains significantly detected, with the alpha detection significance constrained to be non-zero at more than 4.2σ, supporting the robustness of the Fe signal.

In the text

Fig. 10 Left and middle panel: retrieved T–P profiles (Guillot 2010) for the pre-eclipse data with the inclusion of iron (i.e. CO+H2O+Fe; red) and without iron (i.e. CO+H2O; black) under the two retrieval setups outlined in Sect. 4.1. The red and grey shading marks the 1σ and 2σ recovered distribution computed from 10 000 samples from the MCMC. Right panel: the atmospheric structure for the pre-eclipse dataset from the best-fitting CO+H2O only model. The volume mixing ratio profiles for the continuum and detected species are shown as dashed lines (calculated using FASTCHEM), with the parametric T–P profile shown as a solid blue line.

In the text

	Fig. 11 Searching for signatures of atmospheric dynamics. Left and middle panel: posterior distributions for the planetary velocities of CO, H2O and the combined spectrum (i.e. CO+H2O+Fe) for the post-eclipse and pre+post datasets. Right panel: marginalised posterior distributions of the combined atmosphere for the pre- and post-eclipse datasets.
In the text

	Fig. 12 Retrieved planetary orbital velocity (υp) for the detected species from our observations, computed from 10 000 random samples of the MCMC. Dashed lines mark the median, with shaded regions showing 1σ and 2σ contours. For clarity, only the median values (dotted lines) are plotted for the pre+post data.
In the text

	Fig. 13 2D posteriors of H2O for the two datasets using the HITEMP and POKAZATEL line lists (see text).
In the text

Fig. 14 Model scaling factor (α) versus the orbital phase for MASCARA-1b. For visualisation purposes, frames from the individual dataset have been binned together. The α values (blue dots) and their uncertainties for each frame are derived from the conditional likelihood distribution (see text). The black markers show the retrieved median α for the individual datasets. The dark blue points denote outliers corresponding to cloud-affected frames, located within the regions marked by vertical grey bands.

In the text

Fig. A.1 Summary of our injection tests for pre-, post- and pre+post-eclipse datasets outlined in Sect. 3.2, with the 1D and 2D marginalised posterior distributions of each of our model parameters displayed. The red, blue and black posterior distributions represent the datasets, with the injected values highlighted by horizontal and vertical solid purple lines. The dotted lines show the retrieved median values. Upper right: the retrieved T–P profiles are computed from 10,000 random samples of the MCMC; the solid curves show the median profiles and the shaded regions show the 1σ and 2σ contours. The dashed black curve shows the injected T–P profile.

In the text

	Fig. A.2 Same as Fig. A.1 but for a chemical model.
In the text

	Fig. B.1 Free-retrieval results for the post-eclipse data with (red) and without (blue) the inclusion of cloud-affected frames (Sect. 3.2). While excluding these frames slightly improved the Fe constraint, it is still a long-tailed distribution and did not impact the overall conclusions (Fe remains a non-detection).
In the text

	Fig. B.2 Same as Fig. B.1 but assuming chemical equilibrium.
In the text

Fig. C.1 Kp−Δυsγs maps for the non-detected species. Left panels: emission model of each species with the wavelengths covered by CRIRES+ K2166 setting used in this work shown as vertical grey bands. Middle and right panels: the detection maps for the pre-ecl., post-ecl., and combined nights of the CRIRES+ data (Sect. 3.4). The high-resolution opacities were obtained from the following sources: 13CO, CO2 and OH from the HITEMP line list (Rothman et al. 2010), metal opacities from (Kurucz 2018) and HCN opacities from ExoMol database (Harris et al. 2006; Barber et al. 2014).

In the text

	Fig. D.1 Summary of our free-retrieval results for the combined atmosphere with the 1D and 2D marginalised posterior distributions of each model parameter displayed within a corner plot. The red, blue and black distributions represent the pre-, post- and pre+post eclipse datasets (Sect. 4.2).
In the text

	Fig. D.2 Same as Figure D.1, but for the chemical retrieval described in Section 4.2.2.
In the text

	Fig. D.3 Same as Figure D.2, but with a broader metallicity prior (from 0.00l to l000x solar).
In the text

	Fig. E.1 Summary of our free-retrieval results with the 1D and 2D marginalised posterior distributions of the parameters from the MCMC fit. The different colours show the samples obtained using four values for the SYSREM iterations (Sect. 4.2.1).
In the text

	Fig. E.2 Same as Fig. E.1 but for the pre+post dataset.
In the text

Fig. G.1 Results of the injection and recovery test for the pre-eclipse Fe (Sect. 5.2). The red and black posterior distributions represent independent sub-chains of the same MCMC chain, both converging to similar distributions. Upper right: the retrieved T–P profile computed from 10,000 random samples of the MCMC; the solid red curve shows the median profile and the shaded regions show the 1σ and 2σ contours. The dashed black curve shows the injected T–P profile.

In the text

	Fig. G.2 Recovered Kp−Δυsys maps of iron for the pre- and post-eclipse datasets. The white dotted lines mark the peak of detection.
In the text

	Fig. H.1 Results of the free-retrieval for the pre-eclipse dataset with (red) and without (black) the inclusion of Fe as described in Sect. 5.2.
In the text

	Fig. H.2 Same as Fig. H.1 but for a chemical model.
In the text

	Fig. I.1 Results of the joint (pre+post) phase-dependent free-retrieval as described in Sect. 5.2. The subsets of phases close to the quadrature are shown in red and phases close to the eclipse are shown in black. Top right: retrieved T–P profile.
In the text

	Fig. I.2 Left: Kp−Δυsys maps for the iron signal using the phase-dependent retrieval setup. Right: The retrieved planetary orbital velocity (υp) of the detected species for the pre-eclipse data alone, computed from 10,000 random samples of the MCMC. The dashed lines show the median value, and the shaded regions show the 1σ and 2σ contours (see Sect. 5.3).
In the text

	Fig. J.1 Plot of alpha detection significance computed from the conditional likelihood distribution (detailed in Sects. 3.3 and 4) using different line lists (see Sect. 5.3) shown for the pre-eclipse, post-eclipse and pre+post datasets. The labels ‘H10’ and ‘H19’ correspond to the HITEMP 2010 and HITEMP 2019 line lists, respectively.
In the text