© The Authors 2023

1 Introduction

The study of hot Jupiters has opened up exciting opportunities and posed significant new challenges for the analysis of exo-planetary atmospheres. These gas giants provide an excellent opportunity for detailed observations of their atmospheres with current telescopes and instrumentation. Due to their large size and close proximity to the host star, they exhibit a relatively high planet-to-star flux ratio at near-infrared wavelengths, while the short orbital period allows for observations of full phase curves. Additionally, the high day-night temperature difference can cause fast winds on a global scale that redistribute heat and atmospheric constituents.

WASP-43b is a 2.0 M_Jup hot Jupiter that was discovered by Hellier et al. (2011). The planet orbits a young K7 dwarf with a gyrochronological age estimate of Myr (Hellier et al. 2011), which shows chromospheric and coronal emission consistent with a correspondingly high level of stellar activity (Czesla et al. 2013; Staab et al. 2017). Due to its short orbital period of only 0.81 d and the small orbital distance of 0.015 AU (see Table 1), WASP-43b is thought to be tidally locked, making it a suitable target to analyse variabilities across its atmosphere based on phase curve observations. Previous observational data include a low-resolution transmission spectrum and emission phase curves from the Hubble Wide Field Camera 3 (Bean 2013), as well as two datasets in the infrared region measured with the Spitzer InfraRed Array Camera (Blecic et al. 2014). Kreidberg et al. (2014) first found evidence for H₂O absorption and emission, while Chubb et al. (2020) recently reanalysed the data and reported evidence for the presence of AlO being favoured over all other investigated molecules. This points towards the existence of disequilibrium effects in the planet’s atmosphere causing the heavy AlO molecules to be detectable in the upper regions. While the presence of other molecules, such as CO and CH₄, has been suggested in previous studies (Kreidberg et al. 2014; Feng et al. 2016; Stevenson et al. 2017), the results are not conclusive.

Despite WASP-43b’s ultra-short period, its equilibrium temperature of T_eq ≈ 1400 K is considerably lower than that of many other comparable planets due to the low temperature of its host star (T_eff ≈ 4500 K). This opens up the possibility of cloud formation even on the dayside of WASP-43b, but the existence of such cloud layers is still under investigation. While Helling et al. (2020) and Venot et al. (2020) predict the presence of silicate and metal-oxide clouds based on 3D atmospheric modelling, Scandariato et al. (2022) analysed the geometric albedo and found no evidence for condensation of reflective clouds on the planetary dayside.

Exploring the C/O ratio of exoplanet atmospheres has been identified as a promising avenue towards unravelling the mysteries of planet formation and migration. Öberg et al. (2011) suggested that the C/O ratio could offer clues to the planet’s origin within the protoplanetary disk, and several subsequent studies have explored this hypothesis (e.g. Madhusudhan et al. 2016; Mordasini et al. 2016; Booth & Ilee 2019; Cridland et al. 2019). The composition of solid and gaseous material varies with distance from the host star, owing to differences in their condensation temperatures, and as the planet forms and undergoes migration its final composition and atmospheric makeup are shaped by the materials it accumulates along the way. Using simulations, Turrini et al. (2021) find a significant difference in C/O based on whether the migration is gas dominated or solid enriched, while the dependence on the formation distance is less pronounced.

Whether the atmosphere is above or below the tipping point of C/O = 1 has major implications for the entire composition, especially for hot and ultra-hot Jupiters. If there is more carbon than oxygen present, the entire oxygen supply is bound in CO and is not available to form other oxides. For more moderate atmospheric temperatures, however, the dominant carbon-bearing species is CH₄, and thus the oxygen is free to form H₂O (e.g. Giacobbe et al. 2021).

Additionally, the observable C/O is linked to the formation of clouds, as the condensation of certain elements removes them from the gaseous phase that is accessible with transmission and emission spectroscopy. According to Helling et al. (2021), cloud particles such as MgO and Al₂O₃ can remove significant amounts of oxygen from the local atmosphere, leading to an increase in C/O.

The previous observations of WASP-43b with the Hubble Space Telescope (HST) and Spitzer are suitable to measure the H₂O content in hot Jupiter atmospheres, but due to the limited wavelength coverage, the determination of other molecules with the majority of their strong emission lines in the infrared region is more difficult. This is the case for CO and CH₄, which emit most strongly at wavelengths longer than 2 µm. However, an estimation of the C/O ratio necessarily relies on robust estimates of the abundances of these main carbon-bearing molecules. Studies of one and the same planet that are based on different observations, reduction methods, or model assumptions can lead to vastly different C/O measurements, as was the case for example for WASP-127b (e.g. Spake et al. 2021; Boucher et al. 2023) or the ultra-hot Jupiter WASP-12b (e.g. Madhusudhan et al. 2011; Crossfield et al. 2012; Stevenson et al. 2014).

In the case of WASP-43b, the results are less controversial as they generally agree on a C/O < 1. Kreidberg et al. (2014) retrieved the abundances of the main carbon- and oxygen-bearing species using HST observations, but only the H₂O content was well constrained. Based on this result, they conclude that the C/O ratio is consistent with the solar value. Benneke (2015) found an upper limit of C/O < 0.87 in a reanalysis of the data, while Irwin et al. (2020) included the Spitzer/IRAC observations to retrieve C/O ~ 0.91. However, they caution that the abundances of CO and CH₄ effectively only depend on the two Spitzer datasets. The data reduction and assumed thermal profile thus have a strong influence on the retrieved C/O ratio, such that Changeat et al. (2021) recovered C/O = 0.68 with slightly different methods and assumptions.

Using the large simultaneous wavelength coverage and high resolution (R ~ 100 000) of ground-based spectrographs such as CRIRES⁺ mounted on the ESO VLT/UT3 telescope, it becomes possible to circumvent dependencies on a limited number of measurements and the ensuing ambiguities. Observing the spectral bands of CO and CH₄ in the mid-infrared promises the tightest constrains on their abundances due to the large number and intensity of spectral lines, consequently allowing for a more precise determination of the C/O ratio.

In this work, we present high-resolution dayside observations of WASP-43b by employing the CRIRES⁺ (CRyogenic InfraRed Echelle Spectrograph ⁺, Dorn et al. 2014, 2023) instrument. This spectrograph is the upgraded version of the original CRIRES (Kaeufl et al. 2004), offering a significantly larger simultaneous wavelength coverage and new infrared detectors. It is mounted on an 8 m-class telescope, which allows for high S/N observations with sufficient time resolution. In contrast to previous observations using low-resolution instruments, the high resolution of CRIRES⁺ offers the opportunity to directly measure individual absorption and emission lines, while the capability to observe in the infrared region enables the study of species such as CO and CH₄, which do not produce significant spectral features in the visible wavelength range.

In Sect. 2 we describe the observations and our data reduction procedure, and in Sect. 3 we go into the details of the model generation, cross-correlation and retrieval analysis. We present our results and compare them to previous observations and simulations in Sect. 4, and conclude in Sect. 5.

Fig. 1 Orbital phase intervals of the two observations. The grey shaded area indicates the region in which the primary and secondary eclipses occur. The filled circles represent the planet WASP-43b, and the empty circle represents its host star. The semi-major axis and stellar radius are to scale.

2 Observation and data reduction

We observed WASP-43b on 11 March 2022 and 20 February 2023 with the CRIRES⁺ instrument as part of the guaranteed time observations of the CRIRES⁺ consortium. The two observing runs were conducted shortly before and shortly after the secondary eclipse, respectively, to observe the dayside of WASP-43b (see Fig. 1 for the phase coverage). We used the nodding mode technique to observe the target alternating between two different positions (A and B) along the slit to allow for the removal of detector artefacts and the sky background during the extraction. In our further analysis, we combined the A and B spectra to form a single time series for each night. We observed in the K2192 wavelength setting covering parts of the wavelength region from 1972 nm to 2486 nm with an exposure time of 300 s. The CRIRES⁺ science detector array consists of 3 chips, and thus each spectral order is separated into three segments. In the analysis, we treated each segment individually, resulting in 20 wavelength segments in the chosen wavelength setting.

We employed metrology to improve the wavelength calibration by fine-tuning the position of a defined set of emission lines from Kr- and Ne-lamps, and the 0.2″ slit was chosen to ensure a high spectral resolution. The first observing run was conducted under good weather with stable seeing and relative humidity (∼30%), and adaptive optics (AO) was employed. The second night was hampered by a varying coverage of thick and thin clouds, especially during the first half of the observation, which prohibited the use of AO and caused a substantial loss of flux. Although the seeing was similar in both nights, this led to a vastly increased point spread function (PSF) during the second night. The observational conditions that could influence the data quality are shown in Fig. A.1 for the two nights.

The raw spectra were reduced with the ESO CRIRES⁺ data reduction pipeline (version 1.2.3) using the EsoRex¹ tool. This pipeline includes the standard reduction steps for echelle time series observations, such as dark, bias and flat field correction, removal of bad pixels, and a calculation of the wavelength solution. In addition, we attempted to refine the wavelength solution by using Molecfit (Smette et al. 2015) to fit the telluric lines in the spectra, but we did not find any significant offset to the solution provided by the pipeline. In addition, by cross-correlating the telluric lines between individual spectra we measured the drift of the gratings, which manifested themselves in a slight shift of the spectra in time. The calculated drifts were consistently smaller than 0.1 pixels (corresponding to a velocity shift of ∼0.1 km s⁻¹), and therefore we did not apply any further correction.

In the case of good observing conditions that allow an excellent AO correction and coherence time, the PSF can be smaller than the minimum slit width of 0.2″, causing an increase in the spectral resolution (Dorn et al. 2023). The median FWHM of the PSF during night 1 was 2.8 pixel (corresponding to 0.16″), while in night 2 it was 5.8 pixel (0.33″). To determine the resolution, we assumed that the width of the PSF in dispersion direction is similar to the one in spatial direction. Then we converted the PSF from pixel space into wavelength space, which allows for the calculation of the resolution for each segment, and found an average resolution of R ∼ 120000 for night 1. The PSF in night 2 was significantly larger, and in this case the resolution was instead constrained by the slit width to R ∼ 97 000. This shows the effect of the cloud coverage and consequent lack of AO in the second observation, resulting in a significantly broader PSF. Details of the two observations are summarised in Table 2.

2.1 Normalisation

In the following, we describe the normalisation of the spectra, which consists of removing outliers, masking any sky emission lines, and removing the continuum variations. For the outlier removal, any data points deviating by more than 5σ in the time series of each pixel were removed and linearly interpolated. If more than three values of a pixel’s time series would have to be replaced, the entire pixel was masked for all spectra of the time series.

To mask the sky emission lines, a rough continuum approximation was computed and all regions with a flux above 103% of the continuum were masked. To compute such a continuum for each segment, a 90% percentile filter with a width of 100 pixel was applied to the spectra, before smoothing the result with a 50 pixel wide Gaussian filter.

To remove flux variations in the continuum, we determined some points following the shape of the continuum, to which we can fit a polynomial. To this end, we computed a master spectrum for each order as the time average of all exposures. We divided the master spectrum into 100 wavelength bins of equal size, and assigned to each bin the 90% percentile value of all flux values in the bin. Subsequently, the 100 points found in this way were binned again into only ten bins, and from each of these bins we selected the second largest value. The ten resulting points follow the overall shape of the blaze function while not being skewed by improperly removed emission lines. A third degree polynomial was fitted to these ten data points, and all spectra and errors were divided by the fit. Subsequently, each individual spectrum of each segment was corrected to the same continuum level by division of a linear fit to the entire spectrum. As the last normalisation step, deep telluric lines that drop below 40% of the continuum flux were masked.

2.2 SYSREM

We used SYSREM (Tamuz et al. 2005; Birkby et al. 2013; see also Nugroho et al. 2017; Alonso-Floriano et al. 2019; Cabot et al. 2019; Gibson et al. 2020) for the removal of the telluric and stellar lines from the spectra. This principal component analysis algorithm takes into account the uncertainty of each data point to construct a model of the data consisting of linear components in wavelength and time. The stellar, telluric and planetary components of the observed spectra all have different Doppler shifts. While the stellar and telluric lines can be considered nearly static in velocity space over the course of the observation, the planetary radial velocity changes significantly of the order of several tens of kms⁻¹. The SYSREM algorithm iteratively subtracts linear trends in time from each spectral pixel, which primarily serves to remove the (quasi-)static components of the spectra. This leaves residuals consisting only of noise and the planetary signal, which shifts in wavelength over time due to the planetary motion and is therefore not modelled by SYSREM. During each iteration, the linear model created by SYSREM is continuously refined until a convergence criterion is met, after which the model is subtracted from the data and the computation of the next iteration starts. We assumed the model of each iteration to be converged when the average relative change is below 0.01.

The number of SYSREM iterations applied to the data can have a large influence on the recovered signal strength, and we explain our process for choosing this number in Sect. 3.3. The effect of the normalisation and application of SYSREM on the spectra is illustrated in Fig. 2.

3 Methods

Our analysis consisted of generating synthetic planet spectra with key molecular species of interest, cross-correlating these with our cleaned CRIRES⁺ spectra and subsequently applying a retrieval framework to constrain the atmospheric properties.

3.1 Generating planetary model spectra

The synthetic model spectra for the cross-correlation analysis were computed with the radiative transfer code petitRADTRANS (Mollière et al. 2019), using the CO, H₂O and CO₂ line lists of Rothman et al. (2010), the CH₄ line list of Hargreaves et al. (2020), the HCN line list of Harris et al. (2006) and Barber et al. (2014), and the NH₃ line list of Yurchenko et al. (2011). We assumed a non-inverted temperature-pressure (T-p) profile characterised by the parametric model from Guillot (2010). This model is determined by the irradiation temperature T_irr, the internal temperature T_int, the infrared opacity κ_IR and the ratio of visible to infrared opacity γ (see his Eq. (29)). Based on the assumed volume mixing ratios (VMRs), we calculated the mean molecular weight and subsequently the mass fractions of these molecules, assuming a solar-like VMR ratio between H₂ and He. We included collision-induced absorption (CIA) of H₂-H₂ and H₂-He (Borysow 2002; Richard et al. 2012, and the references therein), which we found to have a substantial effect on the line strength for this planet. For analysing the influence of clouds, we additionally added a grey cloud deck at varying pressures, which blocks the contribution of deeper layers. The resulting spectra were converted to a planet-to-star flux ratio. In addition we accounted for rotational broadening of the lines, which has a non-negligible effect for high-resolution spectroscopy. Following Eq. (3) in Díaz et al. (2011), a rotational profile was computed based on the equatorial rotation velocity υ_eq. Similar to Yan et al. (2023), we assumed a linear limb darkening law with a coefficient ϵ = 1 and an inclination angle of sin i ≈ 1 as the planet is expected to be tidally locked and thus the equatorial plane is aligned with the orbital plane. The model was then convolved with this profile in logarithmic wavelength space. Subsequently, the spectra were convolved to the previously determined resolution of CRIRES⁺ (Sect. 2) using a Gaussian instrumental profile, and continuum normalised.

We computed individual synthetic spectra for the species analysed with cross-correlation, and incorporated the model calculation in our retrieval framework. Figure 3 shows the CO, H₂O and CH₄ models used for cross-correlation, before rotational broadening is applied.

Fig. 2 Data reduction steps of a representative wavelength range. Top panel: unprocessed spectra as they are produced by the CRIRES+ pipeline. Middle panel: after continuum normalisation and masking of deep telluric lines. Bottom panel: after the removal of stellar and telluric lines with SYSREM.

3.2 Cross-correlation

We calculated the weighted cross-correlation function (CCF) from the residual spectra and the model following the work from Cont et al. (2022):

(1)

where t is the time index, i is the pixel index, R are the residual spectra, M is the model spectrum and σ is the uncertainty of R. The model spectrum was Doppler-shifted with velocities υ ranging from −500 km s⁻¹ to +500 km s⁻¹. This created a two-dimensional CCF map for each spectral segment. Combining the information of all segments into a single map was done by computing the mean CCF map. In this step we excluded some segments that were substantially contaminated by tellurics and in which no planetary signal could be recovered. We selected these segments by first visually inspecting the spectra, and subsequently calculating the signal strength of an injected model of CO and H₂O as outlined in Sect. 3.3 to determine whether the signal improves when a certain segment is excluded. Based on these criteria, we excluded segments 1 and 20 in night 1 and segments 1, 2 and 20 in night 2. An example of the CCF map for night 1 is shown in Fig. 4.

Following this, a K_p-υ_sys map was calculated from the CCFs. To this end, we considered a range of orbital semi-amplitude values (K_p), and for each value we shifted the CCFs into the corresponding planetary frame according to the Doppler velocity:

(2)

where ф is the orbital phase, υ_sys and υ_bary are the systemic and barycentric velocities and υ_max is a potential deviation from the planetary rest frame. Subsequently, the shifted CCFs were collapsed into a one-dimensional vector by averaging along the time axis. Stacking these vectors for every trial K_p value resulted in a two-dimensional K_p-υ_sys map, which was converted into units of S/N by dividing by the standard deviation of the map in regions far away from the expected signal position (|υ| > 50km s⁻¹).

Fig. 3 Synthetic model spectra of CO (top panel), H2O (middle panel) and CH4 (bottom panel) in the wavelength range of the observations, calculated using petitRADTRANS. The shown models are convolved to a resolution of R = 100 000 and are not rotationally broadened. The grey regions indicate the wavelength ranges covered by the K2192 setting.

Fig. 4 Cross-correlation function of CO for night 1 after five sysrem iterations. The dashed lines indicate the expected planetary trail after the eclipse (which is located between the two lines). The noise structure at υ = 0 kms−1 is caused by tellurics that could not be completely removed.

3.3 Selecting the number of SYSREM iterations

The removal of stellar and telluric lines with SYSREM is an iterative process, where the choice of the number of iterations influences the final result. Too few iterations lead to strong residuals in the K_p−υ_sys map which may hinder the detection of the planetary spectral signature. With too many iterations, SYSREM begins to affect and remove the planetary trail as well. In order to identify the optimal number of iterations in an objective way, we followed the method proposed by Cheverall et al. (2023, see their Sect. 3.8). We first applied SYSREM to the normalised data and calculated the cross-correlation function CCF_obs. Then we repeated these analysis steps on the same data, but with a planetary model spectrum injected at the expected Doppler velocity. This yielded the signal-injected cross-correlation function CCF_inj. Subsequently, we derived the differential cross-correlation function ΔCCF = CCF_inj − CCF_obs, and calculated the K_p−υ_sys map and the S/N for each iteration. Finally, we chose the SYSREM iteration that maximises the S/N of ΔCCF, and applied the same number of iterations to the non-injected data. The results of this selection process are summarised in Table 3 and Fig. 5.

3.4 Retrieval framework

The properties of the atmosphere of WASP-43b were analysed by performing an atmospheric retrieval, in which the observed data is compared to a large number of different models. Using a sampling algorithm to determine the posterior distributions of the model parameters allows for the statistical measurement of atmospheric characteristics. In this way, any properties that have an effect on the observed emission spectra can be analysed. Both the T-p profile and the abundances influence the line strengths, while the position of the lines over the course of the observing run depends on the K_p and systemic velocity. Additionally, the equatorial velocity υ_eq can be determined as a rotation increases the line width by means of rotational broadening. The existence of clouds at a certain pressure P_cioud truncates the spectral lines across the entire wavelength range as contributions from regions below the cloud deck are blocked.

For the retrieval, we assumed that the abundances are iso-baric and that the temperature is described along the parameterisation of Guillot (2010). While this profile does not have enough degrees of freedom to capture the temperature variations across the entire height of the atmosphere, we expect it to be sufficient to model the region with the steepest temperature gradient, from which the majority of the signal originates. We did not include Rayleigh scattering as this process did not have any appreciable effect on the models.

We followed the method of Yan et al. (2023), which is based on the work of Yan et al. (2020) and inspired by Brogi & Line (2019), Shulyak et al. (2019) and Gibson et al. (2020). The forward models were computed as described in Sect. 3.1. To be able to compare these models to the data, it is important to consider any factors that impact the planetary signal. Therefore we modelled the effect of smeared spectral lines and applied corrections to account for distortions introduced by the application of SYSREM.

Fig. 5 S/N detection strength as a function of SYSREM iterations for the data of night 1. The measured signal strengths are a result of using a pure CO model (blue colour), a pure H2O model (orange colour), and a pure CH4 model (black colour). The top panel shows the signal strength for the real data, while the bottom panel shows the strength of an injected signal using the differential ΔCCF according to the method described in Sect. 3.3. The diamond shape indicates the SYSREM iteration with the strongest detection of the injected signal.

Fig. 6 Profile of a single CO line for models with different combinations of smearing and rotational broadening. All models have the instrumental broadening applied, and were additionally smeared with a velocity change of 5kms−1 (orange), rotationally broadened with υeq = 11 km s−1 (purple), or had both effects combined (brown).

3.4.1 Accounting for smearing of the signal

Due to its short orbital period, the radial velocity of WASP-43b changes significantly during a single exposure of 300 s. For night 1, the velocity difference between start and end of an exposure was between 2 and 5 km s⁻¹, causing a substantial smearing of the spectral lines. To account for this effect, the model was interpolated to the same wavelength grid as the 2D data (wavelength × time), resulting in a model matrix M. The model could then be shifted according to the planetary Doppler velocity (see Eq. (2)). For each exposure, we calculated these radial velocities at the beginning and end time, shifted the model to ten evenly spaced planetary rest-frames between these two velocities, and computed the smeared model as the mean over the ten sub-exposures. The effect of smearing on the shape of the spectral lines is shown in Fig. 6. For a model without rotational broadening, the smearing changes the line shape significantly, while in the case of a strongly rotationally broadened model, the smearing is barely noticeable.

3.4.2 Filtering the models

Applying SYSREM not only removes static signals in the data, but also distorts the remaining planetary signal. Therefore, it is important to account for this distortion effect when comparing the planetary signal with synthetic model spectra via Bayesian inference algorithms to retrieve precise atmospheric properties. One could apply the same detrending procedure to both data and model, but this is computationally expensive. Instead, Gibson et al. (2022) proposed a fast alternative filtering technique using just the output of SYSREM, which is briefly summarised in the following:

During each iteration, SYSREM calculates column and row vectors u and w as a best approximation of the data. Combining all column vectors into a matrix U (with the dimensions N_exposures × N_iterations), the correction for a model M_unflltered can be computed:

(3)

where A^† = (A^TA)⁻¹ A^T is the Moore-Penrose inverse of a matrix A. Each entry in the diagonal matrix Λ represents the reciprocal mean uncertainty in wavelength direction of the corresponding exposure. We precomputed this filter and applied it to every model during the retrieval, yielding the filtered model M_unfiltered − M′ which was then compared to the data.

3.4.3 Likelihood function

To retrieve the atmospheric properties, we applied a Markov chain Monte Carlo (MCMC) algorithm using the EMCEE package (Foreman-Mackey et al. 2013) to the logarithm of the likelihood function (Hogg et al. 2010; Yan et al. 2020):

(4)

where R_i,j is the residual spectrum at pixel index i and time j, M_i,i is the 2D matrix of a filtered model spectrum, σ_i,j are the uncertainties of the residuals and β is a scaling factor to σ_i,j. For the MCMC, we used 32 walkers with 15 000 steps each. We chose uniform priors for all input parameters, which are summarised in Table 4.

4 Results

4.1 Cross-correlation results

We calculated K_p-υ_sys maps for CO, H₂O and CH₄ as detailed in Sect. 3.2. The models used for the cross-correlation were computed assuming abundances of log₁₀ CO = −4, log₁₀ H₂O = −4 and log₁₀ CH₄ = −4, respectively. We calculated the temperature-pressure profile using the parameterisation of Guillot (2010) with T_irr = 1400 K, T_int = 800 K, log₁₀(κ_ir)= −2 and log ₁₀(γ) = − 0. 4, and calculated individual model spectra for these species. We only included instrumental broadening and did not account for any other broadening effects for these models. The cross-correlation analysis was conducted separately for each night. Table 3 summarises the resulting S/N of the signal in the K_p−υ_sys maps, as well as the number of SYSREM iterations, determined as outlined in Sect. 3.3. Figures 7 and B.1 show the K_p−υ_sys maps for CO and H₂O of night 1 and night 2, respectively, while the results for CH₄ are shown in Fig. B.2.

We detected the presence of CO and evidence for H₂O in the first night. However, we were not able to find either of the species in night 2. This does not come as a surprise, because the data quality of night 2 is considerably worse (see Sect. 2). We assessed the quality of night 2 by injecting a signal with the expected strength, and were not able to recover it. Therefore the non-detection in night 2 is most likely due to the poor data quality. Pino et al. (2022) showed that principal component analysis algorithms remove much of the planetary signal at phases close to quadrature, as there is only a small shift in the line positions. This could affect especially the signal in night 2, as the observation started at a phase of 0.28. Using an injection-recovery test, we determined that SYSREM removes a substantial amount of the planet signal at phases up to 0.4. We ran the analysis again and excluded any exposures prior to this point in orbit, but we were still not able to detect CO or H₂O in night 2.

We found no evidence for CH₄ in either of the nights. In addition we searched for other molecules that are expected to be present based on equilibrium chemistry calculations (CO₂, HCN and NH₃), but we found no evidence for any of these species in either night. We note that there are spurious peaks which have a S/N of 3–4 visible in some of the K_p−υ_sys maps. However, these are far away from the expected position of a real signal, and are generally not stable against small changes in the model. We tested different models with variations in the T−p profile or the inclusion of rotational broadening, and found that the strength and position of these spurious peaks is not consistent. On the other hand, the CO and H₂O signal in night 1 remain at the same position and at a stable S/N for all of the tested models. This shows that one has to be careful in the interpretation of signals with S/N ∼ 4, as these can in principle be created by spurious noise structures, and the stability of such signals against changes in the model should be investigated.

The results of night 1 agree with earlier studies of WASP-43b using HST and Spitzer observations. Since the first detection of H₂O by Kreidberg et al. (2014), water features have been continuously reconfirmed in subsequent analyses (Stevenson et al. 2017; Chubb et al. 2020). While these studies found no conclusive evidence for CO absorption or emission lines, this is not surprising as the K band, which hosts the most prominent spectral bands of CO, was not covered by these previous observations.

Fig. 7 Kp-υsys map for the CO signal after five SYSREM iterations (top panel) and H2O signal after six iterations (bottom panel) in night 1. The white dotted lines indicate the expected position of a signal at the literature Kp value in the planetary rest-frame. The CO signal has a S/N of 5.02, and the H2O signal a S/N of 4.03.

4.2 Retrieval results

While the cross-correlation method can confirm the presence of a species, it does not provide robust estimates of the abundances. We therefore applied the retrieval framework detailed in Sect. 3.4 to constrain the atmospheric properties. As we did not find evidence for a signal of CO or H₂O (or any of the other species) in night 2, we applied the retrieval only to the data of night 1.

Following a 15 000 step MCMC sampling, we discarded the initial 5000 steps as burn-in, and subsequently calculated the posterior distributions of the atmospheric and orbital parameters, which are illustrated in Fig. 8. The retrieved parameters and their uncertainties are given in Table 4. We retrieved a non-inverted T−p profile (see Fig. 8), which agrees in the upper atmospheric layers with the equilibrium temperature of T_eq = 1426.7 ± 8.5 K. We found a planet velocity semi-amplitude and an offset in velocity from the planetary rest frame of υ_max = −2 ± 3kms⁻¹. The retrieved K_p value is slightly below the literature value of K_p,lit = 202 ± 4kms⁻¹, which can be explained by the elongated diagonal structure of the H₂O detection peak in the K_p−υ_sys map and the dependency of K_p,lit on the semi-major axis, which is not precisely known.

The inclusion of a cloud deck leads to a muting of the signal, as contributions of layers below the cloud deck pressure P_cloud become inaccessible to observations. This introduces a degeneracy between the cloud height and abundances, as a muting of the signal due to clouds can be counterbalanced by increased volume mixing ratios. Figure 9 shows the correlation between P_cloud and the VMRs of CO and H₂O, and the strong degeneracy for cloud deck pressures below one bar is clearly visible. Above this threshold, the VMRs remain constant with decreasing cloud height, because collision-induced absorption (CIA) of H₂−H₂ and H₂−He acts as an additional continuum opacity source, preventing us from probing deeper into the atmosphere. This effect is illustrated in Fig. 10, showing the relation between line strength and cloud deck pressure for models with and without CIA. At pressures of one bar, the line strength of the model with CIA levels off, and a further increase of P_cloud does not affect the model.

Previous retrievals of low-resolution observations found no evidence for the existence of clouds on the dayside of WASP-43b. Fraine et al. (2021) used HST WFC3/UVIS to study the planet’s reflected light and were not able to detect the eclipse, leading them to conclude that a significant cloud coverage is unlikely to exist above the probed pressures of around one bar. Phase curve observations with CHEOPS, TESS and HST constrained the maximum cloud deck pressure to a similar value (Scandariato et al. 2022), and the retrieval of Chubb & Min (2022) using Spitzer data also found a good fit using cloud-free models. Murphy et al. (2023) found that, when comparing their models to the Spitzer phase curve, either thin or very thick clouds were favoured. Although they are unable to reject either of the two scenarios due to a limited model resolution and degeneracies, they suggest a relatively cloud-free dayside to be most likely.

Contrary to these results, pure general circulation and cloud formation models of WASP-43b’s atmosphere by Venot et al. (2020) and Helling et al. (2020) predict several species of clouds at pressures above the detection limit of the aforementioned observational studies, although these results could strongly depend on the initial assumptions and fidelity of the simulations.

Our retrieval results favour an atmosphere free of clouds down to pressures of at least one bar. In addition, due to the unphysically large CO volume mixing ratios when assuming larger cloud deck pressures, and the absence of observational evidence for clouds in the low-resolution data, we deem it most likely that clouds are either not present or located below this threshold pressure on the dayside. Constraining the posterior distributions to include only MCMC steps with log₁₀ P_cloud < 0, we retrieve a super-solar carbon monoxide abundance of , in agreement with the findings of Irwin et al. (2020) and Chubb & Min (2022). If the H2O signal is real and not created by spurious noise, the detected feature corresponds to a water abundance of log₁₀(H₂O) = −4.4± 0.3. This H₂O abundance is consistent within ∼1σ with the results of Stevenson et al. (2017) and Welbanks & Madhusudhan (2019, 2022), but lower than the detections of Kreidberg et al. (2014) and Irwin et al. (2020).

We found a broadening of the spectral lines corresponding to an equatorial velocity , which is larger than the rotational velocity υ_rot = 2πR_p/P = 6.29 ± 0.11 kms⁻¹ assuming the planet is tidally locked. The difference in velocity can be explained by a high velocity jet causing equatorial super-rotation with a wind speed of . For the rotational profile, we assumed a limb darkening coefficient of ϵ = 1, which is justified by the existence of a hot spot creating a strong contrast between the centre and limb region. Setting this parameter to significantly smaller values does not influence the results by more than 1–2 km s⁻¹.

The existence of such a jet is supported both by previous observations and simulations. Multiple studies of Spitzer observations found an eastward shift of the hot spot by ∼10° (e.g. Stevenson et al. 2017; Chubb & Min 2022; Murphy et al. 2023), while Scandariato et al. (2022) analysed TESS phase curves and found marginal evidence for a ()° eastward phase offset. The formation of prograde winds at the equator is generally predicted by 3D general circulation models of tidally locked hot Jupiters (e.g. Showman et al. 2009; Rauscher & Menou 2010; Heng et al. 2011; Mayne et al. 2014). In simulations specifically of WASP-43b, Kataria et al. (2015) find that the jet can in principle reach velocities upwards of 4kms⁻¹, while Zhang et al. (2017) report a wind speed of 6.3 km s⁻¹ at pressure levels associated with CO lines. There has been the notion that both prograde and retrograde streams exist in different regions and heights, and that the retrograde equatorial flow could dominate the overall dynamics (Carone et al. 2020). However, the follow-up study by Schneider et al. (2022) suggests that the retrograde flow may be an artefact of imposed temperature forcing. In their refined model, the prograde jet emerges as the dominant structure with a velocity of ~5 km s⁻¹. Our retrieved equatorial wind speed thus agrees well with previous observations and various simulations of the atmospheric dynamics.

We investigated the influence of smeared spectral lines (see Sect. 3.4.1) on the retrieved υ_eq, but we found no significant difference between retrievals with and without smearing. As shown in Fig. 6, the addition of smearing is relevant in the absence of rotational broadening, but only has a marginal effect on the shape of the lines when a strong rotational broadening is present.

Fig. 8 Posterior distributions of the atmospheric and orbital parameters for WASP-43b using the data of night 1. The top right panel shows the median T−p profile, which was computed from 10 000 random samples of the MCMC walkers. The blue shaded region indicates the 1σ uncertainty interval, while the vertical dotted line shows the equilibrium temperature of WASP-43b (Teq = 1426.7 K). The mean contribution function for the wavelength range of our observation is shown in grey.

Fig. 9 Posterior distributions of the CO abundance (top panels) and H2O abundance (bottom panels). The left panels show the correlation between the VMRs and the cloud deck pressure, and the right panels show the total distribution of the VMRs divided into two subsets: with log10 Pcloud > 0 (blue and red) and with log10 Pcloud < 0 (green and purple). The threshold (indicated by the dashed line) is caused by collision-induced absorption, which blocks any signals originating in the deeper layers.

Fig. 10 Line depth of the deepest line in a model with CO and H2O as a function of cloud deck pressure, for a model with collision-induced absorption (orange) and without this effect (blue).

4.3 Determining the C/O ratio

We conducted a second retrieval, in which the C/O ratio and metallicity [M/H] were set as free parameters instead of the abundances of CO and H₂O. We used the poor_mans_nonequ_chem subpackage from petitRADTRANS to calculate chemical equilibrium abundances from C/O ratios and metallicities. Because we concluded previously that there is no evidence for clouds in the probed region of the atmosphere, we did not include a cloud deck in this second retrieval.

We found C/O = 0.78 ± 0.09 and , while the other parameters (T−p profile, K_p, υ_max, υ_eq, β) were consistent with ourinitial retrieval with free abundances (see Fig. C.1). Calculating the abundances of CO and H₂O based on this C/O ratio and metallicity also gave VMRs close to the results of our initial retrieval. According to these calculations, there could be a significant amount of CH₄ present (log₁₀(CH₄) ≈ −3.2) given the retrieved atmospheric conditions, which we were not able to detect with our data. We conducted an injection-recovery test on the data of night 1 with a pure CH₄ model of that abundance, and were not able to recover any signal. The line strength of the model had to be increased by a factor of 3 to be able to tentatively detect the injected signal in the data. Therefore CH₄ could be present in the atmosphere with an abundance consistent with our retrieved C/O ratio, without being detectable in the data.

Our result is in agreement with the general consensus that the gaseous C/O ratio is significantly smaller than unity (Benneke 2015; Changeat et al. 2021; Chubb & Min 2022). The retrieved super-solar metallicity matches with findings of Kreidberg et al. (2014, [M/H] = 0.4−3.5 × solar) and Stevenson et al. (2017, 0.4 − 1.7 × solar), while it is lower than the estimates of Kataria et al. (2015, 5 × solar) and Changeat et al. (2021, [M/H] = 1.81).

Many previous studies used the Spitzer datasets to constrain the C/O ratio, and these results show a strong dependency on the reduction method and assumed thermal profiles (Changeat et al. 2021). The observations presented in this work cover the CO band head in the K band with the high resolution of CRIRES . Therefore the determination of the CO abundance does not rely on individual data points, as it is informed by a large number of CO absorption lines. Nevertheless, similar to previous works the retrieved C/O ratio is only based on two molecular species, as we were not able to detect other carbon- and oxygen-bearing molecules. This could be either due to small overall atmospheric abundances, or because the observations were not sensitive to the relevant wavelength or pressure regions for these species.

WASP-43b is a target of the JWST Cycle 1 ERS Program (DD-ERS 136, PIs: N. Batalha, J. Bean, K. Stevenson; Stevenson et al. 2016; Bean et al. 2018), which observed a full phase curve of the planet with JWST/MIRI in the infrared 512 µm range (Venot et al. 2020). This observation will be used to study the planetary night side, and to constrain abundances and properties of clouds. A joint retrieval using these not yet published observations with JWST will offer the best bet for an accurate determination of the abundances and the C/O ratio.

5 Conclusions

We analysed the thermal emission spectra of WASP-43b observed with the high-resolution spectrograph CRIRES by implementing the cross-correlation method. We found a robust CO absorption signal, the first direct detection of this species in WASP-43b’s atmosphere, as well as evidence for H₂O absorption lines. For this, we removed stellar and telluric contributions using SYSREM, and optimised the number of iterations based on an injection-recovery test. Subsequently, forward-model retrievals were conducted to constrain the atmospheric and orbital parameters of WASP-43b. Our retrieval accounts for the smearing of spectral lines due to the planetary motion during exposures, applies a filter to the model to account for distortion effects introduced by SYSREM, and includes rotational broadening due to the tidally locked rotation and equatorial winds.

We recovered a non-inverted T−p profile that agrees well with the equilibrium temperature in the upper regions of the atmosphere. The retrieved CO and H₂O abundances are generally consistent with results of previous studies, although there are some studies that report a notably greater H₂O VMR. Assuming equilibrium chemistry in the atmosphere, we found a super-solar metallicity and a C/O ratio of 0.78 ± 0.09. This result is based on the detected signal of CO and H₂O only, and does not account for any other carbon- and oxygen-bearing species that were not detectable with our data. In agreement with simulations of tidally locked hot Jupiters, we found a rotational broadening consistent with a fast super-rotating equatorial jet with a speed of . Furthermore, there was no evidence for a significant cloud deck coverage on the dayside above pressures of one bar based on our data. A combination of these results with the JWST observations of WASP-43b will lead to a larger wavelength coverage and more robust abundance measurements not only of CO and H₂O, but potentially also of other carbon- and oxygen-bearing species.

Acknowledgements

CRIRES⁺ is an ESO upgrade project carried out by Thüringer Landessternwarte Tautenburg, Georg-August Universität Göttingen, and Uppsala University. The project is funded by the Federal Ministry of Education and Research (Germany) through Grants 05A11MG3, 05A14MG4, 05A17MG2 and the Knut and Alice Wallenberg Foundation. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 108.22PH.004 and 110.249Y.004. F.L. acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project number 314665159. D.C. is supported by the LMU-Munich Fraunhofer-Schwarzschild Fellowship and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2094 – 390783311. E.N. acknowledges the support by the DFG Research Unit FOR2544 “Blue Planets around Red Stars”. M.R. and S.C. acknowledge the support by the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (DFG PR 36 24602/41 and CZ 222/5-1, respectively). D.S. acknowledges funding from project PID2021-126365NB-C21(MCI/AEI/FEDER, UE) and financial support from the grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033.

Appendix A Observational conditions

Fig. A.1 figureObservational conditions for night 1 (blue) and night 2 (orange) as a function of phase. Shown are the seeing (first panel), relative humidity (second panel), airmass (third panel), total flux per exposure (fourth panel) and median FWHM of the PSF in spatial direction (fifth panel). The worse data quality of night 2 can be contributed to clouds (as seen by the drop in total flux) and the lack of AO causing a significantly larger PSF.

Appendix B Additional K_p-υ_sys maps

Fig. B.1 Kp-υsys map for the CO signal (top panel) and H2O signal (bottom panel) in night 2 after eight SYSREM iterations. The white dotted lines indicate the expected position of a signal at the literature Kp value in the planetary rest-frame. The absence of a significant signal in either of these species can be attributed to the poor data quality of night 2.

Fig. B.2 Kp-υsys, map for the CH4 signal in night 1 after ten SYSREM iterations (top panel) and night 2 after eight iterations(bottora panel). The white dotted lines indicate the expected position of a signal at the literature Kp value in the planetary rest-frame.

Appendix C Posterior distribution of C/O retrieval

Fig. C.1 Posterior distribution of the retrieval using C/O and [M/H] as free parameters instead of the CO and H2O abundances.

References

All Tables

All Figures

	Fig. 1 Orbital phase intervals of the two observations. The grey shaded area indicates the region in which the primary and secondary eclipses occur. The filled circles represent the planet WASP-43b, and the empty circle represents its host star. The semi-major axis and stellar radius are to scale.
In the text

	Fig. 2 Data reduction steps of a representative wavelength range. Top panel: unprocessed spectra as they are produced by the CRIRES+ pipeline. Middle panel: after continuum normalisation and masking of deep telluric lines. Bottom panel: after the removal of stellar and telluric lines with SYSREM.
In the text

	Fig. 3 Synthetic model spectra of CO (top panel), H2O (middle panel) and CH4 (bottom panel) in the wavelength range of the observations, calculated using petitRADTRANS. The shown models are convolved to a resolution of R = 100 000 and are not rotationally broadened. The grey regions indicate the wavelength ranges covered by the K2192 setting.
In the text

	Fig. 4 Cross-correlation function of CO for night 1 after five sysrem iterations. The dashed lines indicate the expected planetary trail after the eclipse (which is located between the two lines). The noise structure at υ = 0 kms−1 is caused by tellurics that could not be completely removed.
In the text

Fig. 5 S/N detection strength as a function of SYSREM iterations for the data of night 1. The measured signal strengths are a result of using a pure CO model (blue colour), a pure H2O model (orange colour), and a pure CH4 model (black colour). The top panel shows the signal strength for the real data, while the bottom panel shows the strength of an injected signal using the differential ΔCCF according to the method described in Sect. 3.3. The diamond shape indicates the SYSREM iteration with the strongest detection of the injected signal.

In the text

	Fig. 6 Profile of a single CO line for models with different combinations of smearing and rotational broadening. All models have the instrumental broadening applied, and were additionally smeared with a velocity change of 5kms−1 (orange), rotationally broadened with υeq = 11 km s−1 (purple), or had both effects combined (brown).
In the text

	Fig. 7 Kp-υsys map for the CO signal after five SYSREM iterations (top panel) and H2O signal after six iterations (bottom panel) in night 1. The white dotted lines indicate the expected position of a signal at the literature Kp value in the planetary rest-frame. The CO signal has a S/N of 5.02, and the H2O signal a S/N of 4.03.
In the text

Fig. 8 Posterior distributions of the atmospheric and orbital parameters for WASP-43b using the data of night 1. The top right panel shows the median T−p profile, which was computed from 10 000 random samples of the MCMC walkers. The blue shaded region indicates the 1σ uncertainty interval, while the vertical dotted line shows the equilibrium temperature of WASP-43b (Teq = 1426.7 K). The mean contribution function for the wavelength range of our observation is shown in grey.

In the text

Fig. 9 Posterior distributions of the CO abundance (top panels) and H2O abundance (bottom panels). The left panels show the correlation between the VMRs and the cloud deck pressure, and the right panels show the total distribution of the VMRs divided into two subsets: with log10 Pcloud > 0 (blue and red) and with log10 Pcloud < 0 (green and purple). The threshold (indicated by the dashed line) is caused by collision-induced absorption, which blocks any signals originating in the deeper layers.

In the text

	Fig. 10 Line depth of the deepest line in a model with CO and H2O as a function of cloud deck pressure, for a model with collision-induced absorption (orange) and without this effect (blue).
In the text

Fig. A.1 figureObservational conditions for night 1 (blue) and night 2 (orange) as a function of phase. Shown are the seeing (first panel), relative humidity (second panel), airmass (third panel), total flux per exposure (fourth panel) and median FWHM of the PSF in spatial direction (fifth panel). The worse data quality of night 2 can be contributed to clouds (as seen by the drop in total flux) and the lack of AO causing a significantly larger PSF.

In the text

	Fig. B.1 Kp-υsys map for the CO signal (top panel) and H2O signal (bottom panel) in night 2 after eight SYSREM iterations. The white dotted lines indicate the expected position of a signal at the literature Kp value in the planetary rest-frame. The absence of a significant signal in either of these species can be attributed to the poor data quality of night 2.
In the text

	Fig. B.2 Kp-υsys, map for the CH4 signal in night 1 after ten SYSREM iterations (top panel) and night 2 after eight iterations(bottora panel). The white dotted lines indicate the expected position of a signal at the literature Kp value in the planetary rest-frame.
In the text

	Fig. C.1 Posterior distribution of the retrieval using C/O and [M/H] as free parameters instead of the CO and H2O abundances.
In the text