© The Authors 2023

1 Introduction

Shortly after the discovery of the first exoplanet orbiting a solar-type star, 51 Peg b (Mayor & Queloz 1995), we were able to observe the atmosphere around a similar exoplanet (HD 209458b) for the first time by detecting sodium through transmission spectroscopy using a space telescope, specifically the Hubble Space Telescope (HST; Charbonneau et al. 2002). This feature was later observed using data from a ground-based facility by the Subaru telescope/High Dispersion Spectrograph (HDS; Snellen et al. 2008). Since then, the field of exoplanet atmospheres has moved from the simple detection of chemical species to the characterization of atmospheric dynamics. This includes observations of wind patterns (e.g., Seidel et al. 2020b, 2021) or the asymmetry between the morning and evening terminators detected by Ehrenreich et al. (2020) and Kesseli & Snellen (2021). Furthermore, we can observe not only the strong resonant lines such as that of the sodium doublet (Seidel et al. 2019; Langeveld et al. 2022) but also target the cross-correlation function of the weak lines using spectra from high-resolution spectrographs such as High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003), Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO; Pepe et al. 2021), Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000), or HDS (Noguchi et al. 2002; Sato et al. 2002). Using the latter method allows for the creation of a “spectroscopic atlas,” as shown, for example, in Hoeijmakers et al. (2019, 2020).

While the field is rapidly gaining momentum through the use of larger telescopes, we have also started to realize the impact of stellar and telluric effects on the results of transmission spec-troscopy. The first sodium detection (Charbonneau et al. 2002) has been recently disputed by Casasayas-Barris et al. (2021); the reported signal seems to come from the Rossiter–McLaughlin (RM) effect (see also Sing et al. 2008; Casasayas-Barris et al. 2020). RM is an effect where we observe a radial velocity (RV) anomaly from the expected Keplerian value caused by the stellar rotation and transit of the planet. Due to the stellar rotation, the planetary disk occults regions of the stellar surface with different velocities and spectra properties. Since the RV is derived from the disk-integrated spectrum of the star, which is deformed by the occulted regions, we observe the velocity anomaly. In the case of ground-based observations, the telluric lines can also cause false positives and negatives (e.g., Seidel et al. 2020a,c; Langeveld et al. 2021; Orell-Miquel et al. 2022; Spake et al. 2022), and a careful correction of the Earth’s atmosphere is necessary.

The scientific exploration of exoplanetary atmospheres, namely the detection and characterization of chemical composition, is becoming more and more mainstream (e.g., Borsa et al. 2021; Chen et al. 2020; Deibert et al. 2021; Langeveld et al. 2022), such that we can start to perform population studies. The population sample size is mainly limited by the brightness of the host star and atmospheric scale height. There are also many potential exoplanets for planetary characterization not observed by any sufficient observatory yet.

Recently, several teams have attempted to summarize and characterize the detected species in exoplanetary atmospheres, including Guillot et al. (2022, Table 1) and Stangret et al. (2022, Fig. 7), for example. The former summarize chemical species identified in exoplanetary atmospheres, highlighting the size of the studied sample. The latter show species detected based on planets with T_eq > 1000 K, mainly focusing on (ultra)hot Jupiters. Atomic species, such as sodium or hydrogen, were observed only in several tens of exoplanetary atmospheres, sampling mainly hot Jupiters and warm Neptunes.

We could perform a population study to understand why some chemical species are (non)detectable in the exoplane-tary atmosphere. The goal of the Hot Exoplanet Atmospheres Resolved with Transit Spectroscopy (HEARTS) survey is to characterize a sample of hot Jupiters with the high-resolution spectrograph HARPS (Mayor et al. 2003) Fig. 1 (Wyttenbach et al. 2017; Seidel et al. 2019, 2020a,c; Bourrier et al. 2020; Hoeijmakers et al. 2020; Mounzer et al. 2022), and to serve as a pathfinder for VLT/ESPRESSO (Pepe et al. 2021). Data from HEARTS and other surveys are now used for such population studies, such as Langeveld et al. (2022). For the purpose of population study, a set of planets with very similar properties, such as planetary radius and stellar type, is interesting for atmospheric study. However, the sample of currently known exoplanets makes it difficult for such studies to be successful, as not many planets are similar enough and are orbiting stars bright enough to perform such a task.

The target of this paper is KELT-10 b (Kuhn et al. 2016). This hot Jupiter has many properties similar to the benchmark planet system HD 209458 b, providing a valuable opportunity to compare these two planets. It notably features a sodium detection obtained by McCloat et al. (2021) with VLT/UVES. We obtained two transits with HARPS as part of HEARTS, which is further discussed in Sect. 2, with simultaneous photometry observations discussed and analyzed in Sect. 3. The RM effect analysis is discussed in Sect. 4, followed by transmission spectroscopy in Sect. 5, searching for sodium (Na I) and hydrogen (H α). Finally, we discuss our results with concluding remarks in Sect. 6.

Fig. 1 Radius-insolation flux figure of known exoplanets. The color-coding used is as follows: yellow, KELT-10b; green, HEARTS and SPADES surveys’ sodium detections (SPADES is a twin survey using a similar setup as HARPS-N at the 3.58 m telescope at La Palma); and red, WASP-127, the only nondetection from HEARTS. WASP-127 b has not been detected in HARPS data (Seidel et al. 2020c). However, by using VLT/ESPRESSO, Allart et al. (2020) established a tenuous signal below the detection limit of Seidel et al. (2020c).

2 Observations

We observed two transits of KELT-10b on 26 June 2016 and 21 July 2016 from La Silla (Chile). We obtained 96 spectra¹ with HARPS (Mayor et al. 2003) at the 3.6 m telescope and two photometric series with EulerCAM (ECAM) at the 1.2 m Euler telescope² simultaneously. We refer the reader to Lendl et al. (2012) for more information about ECAM. We observed two new light curves that were analyzed jointly with already reduced light curves introduced by Kuhn et al. (2016) to refine the ephemeris system.

The observational conditions (Table 1) were not ideal either night, with high seeing for the first night and partial cloud cover for the second night. Due to the passage of clouds and exposure time (600 s) used on the second night, the observed spectra of this night suffer from a low signal-to-noise ratio (S/N, mean ~15). We ultimately rejected this night from our analysis due to the data quality.

The ECAM light curves are shown in Fig. 2. The entire transit duration is covered, with a short gap due to a telescope software issue during the beginning of transit on night 1. During the second night, a passage of clouds impacted the observations, leading to significant scatter during the transit.

Finally, the CORALIE (Queloz et al. 2000) RV follow-up from Kuhn et al. (2016) has also been used. CORALIE is a spectrograph located at the same telescope as ECAM, with RV measurement precision up to 3 m s⁻¹. CORALIE has been monitoring KELT-10 for potential companions and measuring the planetary mass of KELT-10b. We show all RV observations for both CORALIE and HARPS in Fig. 3.

Fig. 2 Light curves observed complementary to the HARPS dataset using the ECAM at the 1.2 m Euler Swiss telescope.

Fig. 3 RV observations from CORALIE (green) and HARPS (blue and orange for nights 1 and 2, respectively, analyzed dataset). The black line is a Keplerian model as reported in Kuhn et al. (2016). The inset plot in the upper right corner zooms in on the transit in detail. A strong and roughly symmetrical velocity anomaly due to the RM effect is observed during transit.

3 Photometry analysis

While the overall precision of ECAM light curves is better than those used in Kuhn et al. (2016), using only the two transits without the Kuhn et al. (2016) set of data provides compatible yet less precise system parameters. For that reason, we opted to analyze the entire collection of publicly available follow-up photometry (except for the discovery light curve from the KELT-South telescope), using the CONAN code (Lendl et al. 2017).

The first step in CONAN is to find a proper baseline model for the light curve. We searched for different trends, including airmass, sky, x and y coordinate shift, Full Width Half Maximum (FWHM), and sinusoidal pattern. Both ECAM light curves were detrended separately, and the appropriate model was selected using the difference in Bayesian Information Criterion (ΔBIC) value and checking by eye (essential to avoid over-fitting). The 26 June 2016 light curve has been detrended for time and FWHM trends, while the 21 July 2016 light curve has been detrended for time and sky trends.

For the set of limb-darkening coefficients, we used the publicly available code LDCU³ (Deline et al. 2022). The stellar parameters necessary for the limb-darkening computation come from the UVES analysis by Sousa et al. (2018).

Once we established a proper baseline, we ran CONAN on the entire set, using the differential evolution Markov chain Monte Carlo (DEMCMC) method described in Lendl et al. (2017). We ran the code multiple times, and it varied whether or not we fixed an initial parameter between the different usable jump parameters. In the end, we fixed semi-amplitude K to the value from Kuhn et al. (2016). We found better results with eccentricity set to 0, adopting these values. All resulting parameters are shown in Table 2. Finally, we searched for a wavelength dependence of the transit depth in different filters but measured consistent broadband R_p.

The resulting ephemeris is compatible with Kuhn et al. (2016) within 1σ, which can be attributed to different analysis methods. The precision of T₀ and P values has been improved significantly (by a factor of ~18, new precision of $T_0\text{\hspace{0.17em}is}{}_{-0.00029}^{+0.00028}$ and precision of $P\text{\hspace{0.17em}is}{}_{-0.0000032}^{+0.0000032}$) due to tue more extended baseline of observations.

4 Rossiter-McLaughlin effect analysis

4.1 Data reduction

We exploit one dataset obtained with HARPS during the transit of KELT-10 b on 26 June 2016. Forty exposures were acquired (with exposure times of 900 s), among which 14 are in transit, seven are before the transit, and 19 are after the transit. KELT-10 was observed on fiber A, while fiber B was on the sky. The spectra were extracted from the detector images, corrected for, and calibrated using version 3.5 of the Data Reduction Software (DRS) pipeline⁴. One of the DRS corrections concerns the color effect caused by the variability of extinction induced by Earth’s atmosphere (e.g., Bourrier & Hébrard 2014; Wehbe et al. 2020). The flux balance of the KELT-10 spectra was reset to a G2 stellar spectrum template, before they were passed through weighted cross-correlation (Baranne et al. 1996; Pepe et al. 2002) with a G2 numerical mask to compute cross-correlation functions (CCFs) with a step of 0.5 km s⁻¹. This mask is part of a new set that was built with weights more representative of the photonic error on the line positions, as described in Bourrier et al. (2021).

As stated in Sect. 2, another transit was recorded with HARPS on 21 July 2016. This visit was discarded from the analysis because of the poor S/N. In fact, including this transit introduces spurious line variations that bias the results and for which we could not correct. We refer to Appendix C for more details.

Throughout the RM analysis, the posterior probability distributions (PDFs) of free parameters describing models fitted to the data are sampled using emcee (Foreman-Mackey et al. 2013). The number of walkers and the burn-in phase are adjusted based on the degrees of freedom of the considered problem and the convergence of the chains. Best-fit values for these parameters are set to the median of their PDFs, and their 1σ uncertainty ranges are defined using highest density intervals (HDIs).

4.2 Extraction and analysis of individual exposures

We employed a novel technique, the Rossiter-McLaughlin Revolutions (RMR, Bourrier et al. 2021) technique, to analyze the data. This method fits the spectral line profiles from all planet-occulted regions at once with a joint model, which boosts the S/N by the number of in-transit exposures at first order, allowing access to fainter RM signals. The RMR is optimal to fully exploit a RM dataset, to assess biases due to line shape variations, and to measure obliquities from all types of systems in a more comprehensive way than traditional techniques. Indeed, a classical RM analysis only relies on the anomalous RVs of the disk-integrated lines, which condense the information contained in CCF profiles into a single measurement and limit our ability to detect the occultation of the stellar surface by the planet. This process can bias our interpretation of the RV anomaly if the occulted stellar line profile is not well modeled or varies along the transit chord. The RMR technique and its holistic analysis of the planet-occulted, rather than the disk-integrated, starlight, has been successfully used to detect the RM signal of the super-Earth HD 3167 b, the smallest exoplanet with a confirmed RM signal (Bourrier et al. 2021), and to refine the orbital architecture of the warm Neptune GJ 436 b (Bourrier et al. 2022).

The disk-integrated CCF_DI of each individual exposure as produced by DRS, corresponding to the light coming from the entire star, are aligned by correcting their velocity table for the Keplerian motion of the star. A master-out CCF, representative of the unocculted star, was built by coadding CCF_DI outside of the transit. The master-out CCF is fitted with a Gaussian model, and all CCF_DI are shifted to the star rest frame using the derived centroid. Aligned CCF_DI are then scaled to a common flux level outside of the transit and to the flux expected from the planetary disk absorption during transit using a light curve computed with the batman package (Kreidberg 2015). Cross-correlation functions from the planet-occulted regions are retrieved by subtracting the scaled CCF_DI from the master-out CCF and are then reset to a common flux level to yield intrinsic CCF_intr that allow for a more direct comparison of the local stellar lines (Fig. 4).

Following the RMR approach, we first analyzed individual exposures to evaluate the quality of each CCF_intr and assess the possibility for line variations along the transit chord. The various CCF_intr were fitted using a Gaussian model with uniform priors, which are broad enough to constrain the following parameters: 𝒰(–15, 15) km s⁻¹ on the RV centroid (i.e., about three times the maximum typical stellar surface RV for G-type stars); 𝒰(0,20) km s⁻¹ on the FWHM (i.e., about three times the width of the master-out CCF, which is assumed to be similar to that of the CCF_intr given the slow stellar rotation); and 𝒰(0,1) on the contrast (i.e., noninformative). A set of 100 walkers were run for 2000 steps, with a burn-in phase of 500 steps. The model was independently fitted to each CCF over [–150, 150] km s⁻¹ in the star rest frame. The stellar line is detected in all CCF_intr (Fig. 4) with well-defined PDFs for their model parameters (Fig. 5), except for the first and last exposures. This is a consequence of the darkened flux of the stellar limb and its partial occultation by the planet, which results in a lower S/N at ingress and egress. For this reason, we exclude these exposures from further analyses. We note that the derived contrasts and FWHMs are consistent between the different retained exposures (Fig. 5), as is expected from the stability of the G-type star. Surface RVs roughly go from −2 to +2 km s⁻¹ along the transit chord, meaning that the planet is successively blocking just as much blueshifted as red-shifted stellar surface regions, which already suggests an aligned orbit and a moderate stellar rotation.

Fig. 4 Map of the CCFintr during the transit of KELT-10 b. Transit contacts are shown as orange dashed lines. Values are colored as a function of their normalized flux, and plotted as a function of the RV in the star rest frame (in abscissa) and orbital phase (in ordinate). The core of the stellar line from the planet-occulted regions is clearly visible as a bright streak. The magenta solid line shows the stellar surface RV model from the RMR best fit. The small gap in the first part of the transit (at phase ~− 0.005) is due to a missing exposure.

4.3 Derivation of the RM parameters

The RMR technique exploits all of the information contained in the transit time series by simultaneously fitting a model of the stellar line to all CCFintr. The local stellar line is modeled as a Gaussian profile with a constant FWHM and contrast, consistently with the analysis of the individual exposures (Fig. 5). The time series of theoretical stellar lines was convolved with a Gaussian profile of width equivalent to the HARPS resolving power before being fitted to the CCF_intr. Stellar solid-body rotation is assumed. The MCMC jump parameters are hence the line contrast and FWHM, the projected obliquity λ, and the projected equatorial stellar rotational velocity υ_eq sin i_⋆, with the same priors for the contrast and FWHM as in Sect. 4.2, 𝒰(–180,180)° for λ and 𝒰(0,5) km s⁻¹ for υ_eq sin i_⋆.

The resulting best-fit line satisfyingly reproduces the local stellar lines along the transit chord, as can be seen in Fig. 4. The corner plot of the fit is shown in Fig. 6. The resulting PDFs are well defined and symmetrical, in a way that HDIs are here equivalent to quantile-based credible intervals. The correlation diagrams show no correlation between the various parameters, such that the system’s architecture does not depend on the line model. Consistently with the individual exposure analysis (Sect. 4.2), we find a low obliquity and stellar rotational velocity (λ = −5.2 ± 3.4°, υ_eqsin i_⋆ = 2.58 ± 0.12 km s⁻¹), and confirm that a constant stellar line model fits the data well.

Given the quality of the visit, we further tried to constrain finer stellar effects. We fit for stellar differential rotation (DR) employing a law derived from the Sun Ω(θ) = Ω_eq(1 - α sin² θ) (e.g., Cegla et al. 2016), where Ω_eq is the equatorial rotation rate, θ the latitude, and α ∈ [0,1] the unitless DR parameter. Compared to the solid-body fit, two more free parameters are considered here: α and sin i_⋆, for which we use noninformative priors. The resulting PDFs for α and sin i_⋆ are uniform within their definition intervals and the other parameters remain unchanged within 1σ, meaning that DR is not constrained. This is also supported by BIC comparison, which does not justify using the DR model over the solid body one.

For completeness, and to assess the stability of our results against the employed method, we also conducted a Reloaded Rossiter–McLaughlin (RRM, Cegla et al. 2016) fit of the surface RVs. The RRM technique is restricted to the interpretation of the RV centroids from the planet-occulted CCFs. It requires that the surface RVs have a strong enough S/N to be fitted with a stellar line model, which is the case here. We assume solid-body rotation, and the same priors as for the RMR fit were used. The resulting parameters are consistent within 1σ with those derived from the RMR analysis, indicating the robustness of our results. The various derived parameters can be found in Table 3.

Fig. 5 PDFs of the contrast (upper panels), FWHM (middle panels), and RVs (lower panels) of the Gaussian line model fitted to the individual CCFintr. Blue lines indicate the PDF median values with green dashed lines showing the 1σ highest density interval (HDIs). In-transit exposure indices are shown in each subplot.

Fig. 6 Correlation diagrams for the PDFs of the RMR model parameters for the KELT-10 b transit. Curved green and light blue lines show the 1 and 2σ simultaneous 2D confidence regions that contain 39.3% and 86.5% of the accepted steps, respectively. One-dimensional histograms correspond to the distributions projected on the space of each line parameter, with the green dashed lines limiting 68.3% of HDIs. The blue straight lines and squares show the median values.

5 Transmission spectroscopy

This section describes the transmission spectroscopy performed on the HARPS dataset, including all the corrections done before the calculation. We followed the methods by Mounzer et al. (2022), which further refine methods by Seager & Sasselov (2000), Redfield et al. (2008) and Wyttenbach et al. (2015).

We used the e2ds format of HARPS spectra (DRS version 3.5), which includes a 2D grid of 72 and 71 echelle orders (for fiber A and B, respectively), and each order has 4098 pixels. The first step in the analysis is deblazing, using the blaze calibration obtained for each night of observation. This study’s primary focus is orders 56 and 67 (using the order numbering by DRS, 55 and 66 for the fiber B data). The difference in order numbering is due to a gap in fiber B on the order of number 45 (at 5277 Å)⁵. These orders include Na I and H α.

We also excluded several additional spectra from the analysis observed at airmass >2. This filtering affected only the out-of-transit baseline and did not hamper our resulting master-out spectrum (the average out-of-transit spectrum used for stellar line correction). This filtering led to a final set of 37 spectra used in the analysis of transmission spectroscopy.

5.1 Correction for Earth’s atmosphere

After correcting for the instrumental effects, we proceeded to mitigate the impact of Earth’s atmosphere on the observed spectra. To do so, we used several corrections. First, we used molecfit (Smette et al. 2015; Kausch et al. 2015) to correct for the H₂O and O₂ lines from the atmosphere, following the procedures described in Allart et al. (2017) and Seidel et al. (2020c). In a few words, molecfit takes an input scientific spectrum and fits it with a synthetic high-resolution (R ~ 5 000 000) transmission spectrum in selected wavelength regions with prominent telluric absorption. It runs a line-by-line radiative transfer code using an extensive line list of atmospheric lines and calculates the transmission spectrum during the given observation. Finally, it corrects the scientific spectrum, significantly mitigating the effect of the telluric lines. A more detailed description of the program is in the references above.

After correcting for the H₂O and O₂ lines, we normalized the spectra using the median flux of each spectrum. This normalization still leaves a color effect in the spectrum due to the instrument that could be corrected by a better normalization approach or by subtracting a linear fit. Since this effect is stable for a night, it is unnecessary to rectify this trend since we will eventually divide by the master-out spectrum.

After normalizing the spectra, we corrected for the cosmic rays. We calculated a median spectrum and flag pixels as cosmic (and set them to the median spectrum value) if the flux difference from the median spectrum was beyond the 5σ level. We also manually checked the spectra for missed cosmics during the correction, finding that no cosmic rays hampered the region of Na I and H α.

Finally, molecfit does not correct for sodium emission in Earth’s atmosphere. We monitored the fiber B data (aimed at the sky background) and found significant sodium emission. Our attempts to subtract this emission were unsuccessful because of the fiber A/B efficiency ratio variability at night. Due to this variation, we observed the remainder of the emission feature at high airmass, creating a spurious signal in the transmission spectrum. As a result, we thus opted to mask a small part (around 8 pixels) of the spectrum.

Afterward, we shifted all the spectra into the rest frame of the star, using the BERV (approximately 2.7 km s⁻¹ for night 1) and systemic velocity value calculated from the analysis of the master-out CCF (Sect. 4; 31.96 ± 0.02 km s⁻¹). The value of systemic velocity calculated from HARPS is compatible with both the values from Kuhn et al. (2016, 31.9 ± 0.1 km s⁻¹) and the one used in McCloat et al. (2021, 31.61 ± 1.29 km s⁻¹).

5.2 Correction for stellar effects

Once we shifted to the stellar rest frame, we could correct for the stellar lines. At first order, we could calculate the average out-of-transit spectrum (M_out) and divide our dataset by it (Wyttenbach et al. 2015; Seidel et al. 2019). However, doing so disregards the RM and center-to-limb variation (CLV) effects that can affect our spectra. This approach is the one McCloat et al. (2021) took, as the lack of precise wavelength calibration in their UVES data disallowed the measurement of RVs from the spectra, preventing an RM correction. HARPS, however, is designated to provide accurate RVs, hence we could measure and correct for the RM effect and simulate its impact on the final transmission spectrum.

To correct our spectra for the RM effect, we adapted the formula by Mounzer et al. (2022): $$F_{\text{corr,i}}=1+\left(\frac{R_{\text{p}}^2}{R_{\text{s}}^2}\right)\left(1-\frac{M_{\text{out}}-F_i{\delta }_i}{M_{{\text{out,v}}_{\text{i}}}{\text{Δ}}_{\text{i}}}\right).$$

Here, $\frac{R_{\text{p}}^2}{R_{\text{s}}^2}$ corresponds to the planet/star radius ratio, F_i is the spectrum we wish to correct, M_out is the master-out spectrum, and $M_{{\text{out,v}}_{\text{i}}}$ is master-out shifted by the local stellar velocity. The δ_i and ∆_i correspond to the transit depth and white-light flux received during the individual exposure. The υ_i velocity was calculated using the λ = –5.2 ± 3.4° and υ_eq sin i_⋆ = 2.58 ± 0.12 km s⁻¹ values calculated in Sect. 4, using an analytical formula (for circular orbit) describing the RVs of the stellar surface along the transit chord of the planet for each exposure: $${\upsilon }_i\left(\varphi \right)=\left(\frac{a}{R_{\text{s}}}\right)\left(\sin \left(2\pi \varphi \right)\cos \left(\lambda \right)+\cos \left(2\pi \varphi \right)\cos \left(i\right)\sin \left(\lambda \right)\right){\upsilon }_{\text{eq}}\sin i_{★},$$

where ϕ is the orbital phase at which we observed the planet. Other parameters follow the same notation as in previous sections (cf. Table 2 and Sect. 4). The shift is between –2.6 to +2.6 km s⁻¹. This roughly corresponds to 3–4 pixels (1 pixel ~820 ms⁻¹) in HARPS spectra from the unshifted M_out.

To estimate the impact of the RM effect on the transmission spectrum of KELT-10 b (for comparison with the McCloat et al. 2021 detection), we used the EVaporating Exoplanets (EVE) code (e.g., Bourrier & Lecavelier des Etangs 2013; Allart et al. 2018). This code allows for simulations of the 3D atmosphere of a transiting planet and calculations of synthetic transmission spectra accounting for the geometry of the transit and the properties of the stellar surface. First, we accounted for CLV of the local stellar spectrum by tiling the EVE stellar grid using synthetic line profiles calculated at different limb darkening angles with the Turbospectrum⁶ code (Plez 2012). Then the RM effect was naturally accounted for by shifting the local stellar spectra according to the rotational velocity of KELT-10. Using these settings, we simulated a transit of KELT-10b without including any atmosphere to assess the distortion induced purely by the RM effect and CLV. We set the code to return a time series of KELT-10 spectra as they would be observed with HARPS and we used them to compute an average transmission spectrum comparable to our observations.

Figure 7 shows the comparison between our EVE simulation and the best-fit model to the UVES transmission spectrum from McCloat et al. (2021) in the region of the sodium doublet. We are unable to explain the observed signal with the pure RM and CLV effects. The amplitude of the effect is on the order of 0.2%, and opposite to an atmospheric absorption signal, impacting the transmission spectrum significantly. It is even possible that the RM effect is obscuring part of the putative atmospheric signal; however, given the issues with the wavelength calibration in the UVES data (Sect. 5.3 in McCloat et al. 2021), it is difficult to draw a definite conclusion.

If the observed UVES signal is truly redshifted, we could explain the “emission-like” shoulders apparent in the UVES transmission spectrum (Fig. 4 in McCloat et al. 2021). However, only the D1 signal there seems to be as broad as our simulation, while the D2 signal appears to be due to a few pixels. The amplitude of the signature would be slightly reduced, but not enough to compromise its detection.

On the other hand, if the redshift is purely of nonastro-physical origin, the RM effect would significantly decrease the observed signal, meaning the true amplitude of the sodium signal is higher. In that case, KELT-10 b would have one of the strongest sodium signals detected yet, which, considering its location in the hot Jupiters’ population (Fig. 1), seems curious, as we would expect a stronger signal in hotter and puffier planets (as long as the sodium itself is not ionized). It is also possible that due to the lower resolution of UVES spectra, the RM effect is smeared out and hence weaker.

After stellar line correction, including the RM effect (CLV is negligible compared to RM), a visible Na l remnant is still present. This remnant results from very low S/N in the core of the stellar Na I and should be removed, as it can create spurious signals. We masked the core region with NaN values, following the procedures presented in Seidel et al. (2020a, c). This way, we could mask the sodium remnant, but we unfortunately also masked part of the planetary signal.

To estimate how much of the planetary signal was removed, we used the value of the mean FWHM of the Na I signal from Seidel et al. (2019) (FWHM = 0.65 Å). The area within the FWHM marked the region for the planetary signal, and we assumed no shift of the sodium line (in the planetary rest frame). The fraction of information lost is $\frac{FWH{M}_{\text{mask}}}{FWH{M}_{\text{signal}}}=28\%$. After we shifted the spectra to the rest frame of the planet, we calculated an average transmission spectrum over full in-transit exposures, weighting them by S/N.

5.3 Transmission spectroscopy

The final transmission spectrum is shown in Figs. 8 and 9 for Na I and H α, respectively. Although we corrected our spectrum for the RM effect, this effect is negligible within our detection limits. We also checked the RM uncorrected spectrum, finding no significant changes in the notable features of the spectrum.

Neither of the sodium lines were detected. As part of the comparison with McCloat et al. (2021), a Gaussian model with their reported parameters has been over-plotted. To calculate the upper limit, we took the mean uncertainty in the ±1.5 Å region around the lines, marked as a yellow area, and used it as a 1σ estimate. We used a broader region to estimate the uncertainty since we do not know of any potential shift of the possible signal. These kinds of shifts can occur due to wind patterns and this is quite common (e.g., Snellen et al. 2010; Seidel et al. 2020b, 2021; Mounzer et al. 2022). The 3σ upper limits are 3.34% for Na I and 3.97% for H α. These upper limits are likely underestimated due to the distortion of the putative atmospheric signal by the planet-occulted stellar lines.

In the region of Na I, we see a few bins that show absorption. However, the location of these bins in each of the doublet lines is different, with –30 km s⁻¹ for the D2 feature and –10 km s⁻¹ for the D1 feature. This suggests that these features are of nonastro-physical origin. In fact, the D2 feature can be attributed to telluric sodium. For the D1 feature, we attribute this to a low S/N remnant (the count in the line core is below 100) and systematics in the measurements – several features of a similar strength have been observed (in the red part of the spectrum) as well. These features also hint at potential correlated noise. Finally, this kind of signal, if real, would be visible in the UVES data. Based on these arguments, we, therefore, claim nondetection.

Fig. 7 Modeling of the RM + CLV effect on the transmission spectrum in the rest frame of the planet, specifically focused on the sodium detection presented in McCloat et al. (2021, AD2 = 0.66%, σD2 = 0.06 Å, ΔλD2 = +0.23 Å; AD1 = 0.43%, and σD1 = 0.06 Å, ΔλD1 = +0.23 Å). Top panel: Gaussian fit to the lines attributed by McCloat et al. (2021) to planetary sodium (their Fig. 4, blue) and simulated RM and CLV effects using EVE (red). We note that compared to McCloat et al. (2021), we use a reversed y-axis (due to different units used), that is absorption is above continuum and emission is below it. The 3σ and 5σ level reported by McCloat et al. (2021) is shown as dark blue dashed lines. Bottom panel: O–C values (UVES–EVE simulation, blue), and 3σ and 5σ UVES detection levels (blue).

Fig. 8 Transmission spectroscopy of KELT-10 b in the region of the sodium doublet. In gray, the entire transmission spectrum is shown (corrected for RM effect, as described in Sect. 5.2), with 20x bins in dark blue. The detection of McCloat et al. (2021) is overplotted (dark red) as a visual cue for comparison. The black dashed line corresponds to the 3σ and 5σ limit above the continuum $\left(R_{\text{p}}^2/R_{\text{s}}^2\sim 1.44\%\right)$. The yellow area is the ±1.5 Å region around each of the Na I lines.

Fig. 9 Same as Fig. 8, but in the region of H α. The color-coding is the same as in Fig. 8. The additional red line is a continuum of one equivalent height. No detection is visible in this region either, with corresponding 3σ and 5σ limits on the plot.

6 Conclusion

6.1 Bulk planetary and orbital properties

We have carried out an independent photometry analysis, which shows consistent results with parameters from Kuhn et al. (2016) within 1σ. We improved the precision of the P and T₁₄ parameters by a factor of ~18.

We were able to study the RM effect with great precision using the RMR method (Bourrier et al. 2021), finding the orbit of KELT-10 b to be aligned with the stellar rotation plane, in projection onto the plane of the sky λ = –5.2 ± 3.4°. The system is likely truly aligned, as previous results show that hot Jupiters around late-type stars (T_eff < 6250 K) tend to live on aligned orbits (Winn et al. 2010; Albrecht et al. 2022). Indeed, late-type stars are expected to have substantial convective envelopes, which trigger stronger tidal interactions, thus realigning their systems more efficiently. This argument is supported by the proximity of KELT-10 b to its host star (a/R_⋆ < 10, Table 2), and by the relatively old age of the system (4.5 ± 0.7 Gyr, Kuhn et al. 2016).

6.2 Transmission spectroscopy

The nondetection of 3σ_lim = l.8R_p(Na I) and 3σ_lim = 1.9R_p (H α) corresponds to no extreme extension of upper atmospheric layers, as expected from such a “well-behaved” system. The RM + CLV effects create a spurious signal, very similar to the one detected by Casasayas-Barris et al. (2021) on the example of HD 209458 b. In future studies of this system, these effects should be corrected for in the transmission spectra or included in EVE-like simulations to derive (a) confident atmospheric detection(s).

6.3 Comparison with previous sodium detection

Recently, McCloat et al. (2021) have claimed sodium detection from a VLT/UVES single transit dataset. However, these data could not allow these authors to obtain RV measurements precise enough to correct for the RM effect (which was attributed to molecf it being unable to fine-tune the wavelength grid – see Sect. 4.6 in McCloat et al. 2021).

Since the authors could not measure the RVs during the transit, a potential RM effect could have created a significant signal, as in the case of Casasayas-Barris et al. (2021). Furthermore, the tentative sodium detection had a significant redshift, likely due to the wavelength calibration issues mentioned above.

We modeled the effect of RM and CLV on the transmission spectrum, finding that the distortion by planet-occulted stellar lines cannot explain the measured sodium detection. This discrepancy with the detection of McCloat et al. (2021) is surprising, especially if the observed redshift is instrumental, as the RM effect would significantly strengthen the sodium detection, making it one of the strongest detections of the Na I doublet to date.

Our transmission spectrum cannot confirm the UVES sodium detection (McCloat et al. 2021), but we also cannot rule it out, as the lower quality of the final transmission spectrum does not allow us to probe to such detection limits. There are several potential reasons why our results do not fully match.

First, our RM correction method might introduce several biases. In this work, we did not thoroughly investigate different stellar models to use in the RM+CLV modeling. Incorrect RM correction can then introduce spurious effects that blur the signal (under-correction) or create a fake one (over-correction).

Furthermore, molecfit struggles with the correction of low S/N spectra, which is the case locally in the core region of the sodium doublet. In addition to that, due to low counts on the detector in the stellar core (<100), we are no longer photon-noise dominated but red-noise dominated (Seidel et al. 2022; Pepe et al. 2021, see Sect. 5 and Appendix D). This is the case for all observations when the observed planetary signal is on the order of (or less than) a single photon on the detector and needs to be treated carefully.

Finally, it is possible that either this work or McCloat et al. (2021) is further affected by some unaccounted effect (e.g., stellar variability, instrument issue), which could cause a discrepancy between the observations. In fact, the peaks in our transmission spectrum further in the red hint at the possibility of correlated noise and systematic errors, which are also on the same order as the expected signal based on UVES (McCloat et al. 2021). On the other hand, atmospheric variability can significantly affect the dataset, especially for slit-based spectrographs such as UVES.

We conclude that to properly address the question on the presence of chemical species in the atmosphere of KELT-10b (and confirmation of the previous Na I detection), we need additional data from larger telescopes than the 3.6m telescope at La Silla to obtain a good S/N and temporal sampling and to correct for the stellar effects such as RM that cause a significant impact (~0.2%) on the transmission spectrum.

Acknowledgements

We would like to thank the anonymous reviewer for the insightful comments. We would like to thank Sean McCloat, who was kind enough to provide reduced UVES data for comparison purposes with the analysis presented in this paper. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Four Aces; grant agreement No 724427). It has also been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). DE and MS acknowledge financial support from the SNSF for project 200021_200726. M.A. and V.B.’s work has been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). M.A. and V.B. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Spice Dune; grant agreement No 947634). This project has also been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the SNSF. M.L. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 742095; SPIDI: Star-Planets-Inner Disk-Interactions; http://www.spidi-eu.org). S.G.S. acknowledges the support from FCT through the contract nr.CEECIND/00826/2018 and POPH/FSE (EC). Python packages used: numpy, scipy, matplotlib, astropy (Astropy Collaboration 2022), specutils, batman.

Appendix A Observational log

Detailed observational log for both nights of HARPS. The first night (Figure A.l) suffered from high seeing, which varied during the night. The second night (Figure A.2) shows a low S/N, leading to removal from the dataset.

Fig. A.1 Observational log for the first night (26 June 2016). Top panel: Airmass (blue) and seeing (red). Bottom panel: S/N for order number 56 (blue) and 67 (red).

Fig. A.2 Same as Figure A.1, but for the second night (21 July 2016).

Appendix B Fiber B sodium emission

As discussed in Section 5.1, we detected a significant telluric emission feature from sodium. This is shown in Figure B.1, on a total master spectrum of the first night of fiber B.

Fig. B.1 Sodium emission observed on fiber B during the first night of observation. Not correcting for it causes a significant absorption feature. No features were found on fiber B for the wavelength region of Hα. The red dashed line signifies the masked part of the emission.

Appendix C RM analysis of the discarded visit

An additional transit was recorded with HARPS on 21 July 2016. The S/N of this visit is substantially lower than the 26 June 2016 one (16 on average for the former versus 31 for the latter, see Fig. Figure A.1 and Figure A.2). Furthermore, the 21 July 2016 spectra were Moon-contaminated shortly after the end of the transit, resulting in spurious features in the affected CCFs. This is visible in the disk-integrated contrast fluctuations starting after the transit (Fig. C.1, upper panel). We show the contrasts of the 26 June 2016 CCFs for a visual comparison of the stability of the two nights.

Fig. C.1 Contrast of the disk-integrated CCFs as a function of orbital phase (top) and S/N (bottom) for the 26 June 2016 visit (blue points) and the 21 July 2016 visit (red points). The horizontal dashed lines indicate the mean value for each visit. The vertical dashed lines represent the transit contacts. Out-of-transit (in-transit, resp.) exposures are shown as squares (circles, resp.). The full-line curve in the lower plot represents the best polynomial fit.

The 21 July 2016 fluctuations correlate with the S/N, as seen in Fig. C.1 (lower panel). We have captured this correlation by fitting a third order polynomial to the affected visit, based on BIC comparison. We then corrected the disk-integrated (DI) CCFs of this visit according to the fit result, which reduced the dispersion around the mean contrast value by a factor of ~ 3. Ultimately, this correction turned out to be detrimental to the intrinsic CCFs. Indeed, Fig. C.2 shows a large dispersion of the derived intrinsic RVs after correction, which results in very large error bars for the derived RM parameters. We have thus discarded the correction.

We further tried to circumvent the Moon contamination by excluding the post-transit baseline from the construction of the master out. Expectedly, the systemic velocity derived from the master out now coincides 26 times better between the two nights. Yet, the RMR analysis after this correction yields a very large difference in the derived parameters between the two nights. The same settings for the MCMC exploration as in Sect. 4 were used. The best-fit results for the affected visit are υ_eq sin i_⋆ = 1.88 ± 0.21 km s⁻¹ and $\lambda ={23.0}_{-9.0}^{+8.7}{}^{\circ }$. Due to the 2.9-σ discrepancy between the two nights for the obliquity value, we have also discarded the correction.

Fig. C.2 Intrinsic RVs as a function of orbital phase before (top) and after (bottom) the correction for the contrast–S/N correlation, for the 26 June 2016 visit (blue points) and the 21 July 2016 visit (red points). The horizontal dashed line indicates the zero RV level. The vertical dashed lines represent the transit contacts. The solid black line is the analytical RV curve for the adopted values of Table 3.

As a last attempt, we excluded the Moon-affected regions from the DI CCFs (roughly [-20, 0] km s⁻¹), which should represent a smoother correction than excluding the full post-transit baseline. Still, a significant discrepancy in the RMR results between the two nights persists, even if it is lower than with the last correction: ${\upsilon }_{\text{eq}}\sin i_{\star }={2.43}_{-0.22}^{+0.21}\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$ and $\lambda ={9.1}_{-9.1}^{+8.8}{}^{\circ }$. Even if the stellar equatorial velocities are now compatible, the obliquities have a 1.5-σ discrepancy between the two nights. In the absence of agreement between the two visits and despite our tentative corrections, we have discarded the Moon-contaminated night.

Appendix D Masking stellar sodium remnant

After correcting for master out, we observed a low S/N remnant in the sodium doublet due to line saturation. To remove this effect, we masked the region, as shown in Figure D.1

Fig. D.1 Used mask in the steps as described in Section 5. Top panel: Colormap of all spectra during the first night in the rest frame of the star. The red dashed lines correspond to transit contact points. The white mask is the mask used for telluric sodium, located at roughly –30 km s−1. Middle panel: Master-out spectrum, with the following colored lines: white dash-dot line, masked region; and green, sodium line $\pm \frac{FWHM}{2}$ region as reported in Seidel et al. (2019). Bottom panel: Masked stellar signal based on the mask shown in the middle panel (white lines), which is still in the rest frame of the star. The black line corresponds to the sodium line in the planetary rest frame, considering no additional shifts (e.g., due to winds). Green lines are offset by $\pm \frac{FWHM}{2}$ as used in the middle panel, and they correspond to the area from which we estimated our lost signal. The noisy spectra at ingress and egress (around red dashed line) are due to RM correction and light-curve weighting. These two spectra were not used in the final transmission spectrum.

References

All Tables

All Figures

Fig. 1 Radius-insolation flux figure of known exoplanets. The color-coding used is as follows: yellow, KELT-10b; green, HEARTS and SPADES surveys’ sodium detections (SPADES is a twin survey using a similar setup as HARPS-N at the 3.58 m telescope at La Palma); and red, WASP-127, the only nondetection from HEARTS. WASP-127 b has not been detected in HARPS data (Seidel et al. 2020c). However, by using VLT/ESPRESSO, Allart et al. (2020) established a tenuous signal below the detection limit of Seidel et al. (2020c).

In the text

	Fig. 2 Light curves observed complementary to the HARPS dataset using the ECAM at the 1.2 m Euler Swiss telescope.
In the text

	Fig. 3 RV observations from CORALIE (green) and HARPS (blue and orange for nights 1 and 2, respectively, analyzed dataset). The black line is a Keplerian model as reported in Kuhn et al. (2016). The inset plot in the upper right corner zooms in on the transit in detail. A strong and roughly symmetrical velocity anomaly due to the RM effect is observed during transit.
In the text

Fig. 4 Map of the CCFintr during the transit of KELT-10 b. Transit contacts are shown as orange dashed lines. Values are colored as a function of their normalized flux, and plotted as a function of the RV in the star rest frame (in abscissa) and orbital phase (in ordinate). The core of the stellar line from the planet-occulted regions is clearly visible as a bright streak. The magenta solid line shows the stellar surface RV model from the RMR best fit. The small gap in the first part of the transit (at phase ~− 0.005) is due to a missing exposure.

In the text

	Fig. 5 PDFs of the contrast (upper panels), FWHM (middle panels), and RVs (lower panels) of the Gaussian line model fitted to the individual CCFintr. Blue lines indicate the PDF median values with green dashed lines showing the 1σ highest density interval (HDIs). In-transit exposure indices are shown in each subplot.
In the text

Fig. 6 Correlation diagrams for the PDFs of the RMR model parameters for the KELT-10 b transit. Curved green and light blue lines show the 1 and 2σ simultaneous 2D confidence regions that contain 39.3% and 86.5% of the accepted steps, respectively. One-dimensional histograms correspond to the distributions projected on the space of each line parameter, with the green dashed lines limiting 68.3% of HDIs. The blue straight lines and squares show the median values.

In the text

Fig. 7 Modeling of the RM + CLV effect on the transmission spectrum in the rest frame of the planet, specifically focused on the sodium detection presented in McCloat et al. (2021, AD2 = 0.66%, σD2 = 0.06 Å, ΔλD2 = +0.23 Å; AD1 = 0.43%, and σD1 = 0.06 Å, ΔλD1 = +0.23 Å). Top panel: Gaussian fit to the lines attributed by McCloat et al. (2021) to planetary sodium (their Fig. 4, blue) and simulated RM and CLV effects using EVE (red). We note that compared to McCloat et al. (2021), we use a reversed y-axis (due to different units used), that is absorption is above continuum and emission is below it. The 3σ and 5σ level reported by McCloat et al. (2021) is shown as dark blue dashed lines. Bottom panel: O–C values (UVES–EVE simulation, blue), and 3σ and 5σ UVES detection levels (blue).

In the text

Fig. 8 Transmission spectroscopy of KELT-10 b in the region of the sodium doublet. In gray, the entire transmission spectrum is shown (corrected for RM effect, as described in Sect. 5.2), with 20x bins in dark blue. The detection of McCloat et al. (2021) is overplotted (dark red) as a visual cue for comparison. The black dashed line corresponds to the 3σ and 5σ limit above the continuum $\left(R_{\text{p}}^2/R_{\text{s}}^2\sim 1.44\%\right)$. The yellow area is the ±1.5 Å region around each of the Na I lines.

In the text

	Fig. 9 Same as Fig. 8, but in the region of H α. The color-coding is the same as in Fig. 8. The additional red line is a continuum of one equivalent height. No detection is visible in this region either, with corresponding 3σ and 5σ limits on the plot.
In the text

	Fig. A.1 Observational log for the first night (26 June 2016). Top panel: Airmass (blue) and seeing (red). Bottom panel: S/N for order number 56 (blue) and 67 (red).
In the text

	Fig. A.2 Same as Figure A.1, but for the second night (21 July 2016).
In the text

	Fig. B.1 Sodium emission observed on fiber B during the first night of observation. Not correcting for it causes a significant absorption feature. No features were found on fiber B for the wavelength region of Hα. The red dashed line signifies the masked part of the emission.
In the text

Fig. C.1 Contrast of the disk-integrated CCFs as a function of orbital phase (top) and S/N (bottom) for the 26 June 2016 visit (blue points) and the 21 July 2016 visit (red points). The horizontal dashed lines indicate the mean value for each visit. The vertical dashed lines represent the transit contacts. Out-of-transit (in-transit, resp.) exposures are shown as squares (circles, resp.). The full-line curve in the lower plot represents the best polynomial fit.

In the text

Fig. C.2 Intrinsic RVs as a function of orbital phase before (top) and after (bottom) the correction for the contrast–S/N correlation, for the 26 June 2016 visit (blue points) and the 21 July 2016 visit (red points). The horizontal dashed line indicates the zero RV level. The vertical dashed lines represent the transit contacts. The solid black line is the analytical RV curve for the adopted values of Table 3.

In the text

Fig. D.1 Used mask in the steps as described in Section 5. Top panel: Colormap of all spectra during the first night in the rest frame of the star. The red dashed lines correspond to transit contact points. The white mask is the mask used for telluric sodium, located at roughly –30 km s−1. Middle panel: Master-out spectrum, with the following colored lines: white dash-dot line, masked region; and green, sodium line $\pm \frac{FWHM}{2}$ region as reported in Seidel et al. (2019). Bottom panel: Masked stellar signal based on the mask shown in the middle panel (white lines), which is still in the rest frame of the star. The black line corresponds to the sodium line in the planetary rest frame, considering no additional shifts (e.g., due to winds). Green lines are offset by $\pm \frac{FWHM}{2}$ as used in the middle panel, and they correspond to the area from which we estimated our lost signal. The noisy spectra at ingress and egress (around red dashed line) are due to RM correction and light-curve weighting. These two spectra were not used in the final transmission spectrum.

In the text