© ESO 2021

1 Introduction

The transmission spectroscopy method is applicable when an exoplanet crosses the disk of its host star as seen from Earth. Within this geometry, a part of the stellar light passes through the exoplanet atmosphere, and the atmosphere imprints its signature onthe spectrum. Through low- and high-resolution spectroscopy, one can access embedded spectral information from different pressure levels in the atmosphere. Low-resolution transmission spectra probe deeper levels of the atmosphere; they usually contain information regarding the atmospheric temperature, abundances, the presence or absence of haze and clouds, and the general shape of the atomic and molecular continuum features (e.g., Sing et al. 2008; Vidal-Madjar et al. 2011; Fischer et al. 2016). On the other hand, high opacities at lower atmospheric pressures can imprint their spectral signatures on high-resolution spectra, resulting in unambiguous detection of atmospheric atomic and molecular features (e.g., Hoeijmakers et al. 2019; Yan et al. 2019; Sánchez-López et al. 2019; Nortmann et al. 2018; Khalafinejad et al. 2017). Therefore, acombined or even parallel analysis of both low- and high-resolution spectra provides the complementary information (e.g., Brogi et al. 2017; Pino et al. 2018; Carone et al. 2021) required to fully unveil exoplanet atmospheres.

WASP-69 b is a hot, inflated Saturn-mass exoplanet orbiting the active K-type star BD–05 5432 with an orbital period of ~3.87 d (Anderson et al. 2014). The most relevant parameters of this star–planet system are listed in Table 1. Following the discovery of the planet, there have been several attempts to characterize the atmosphere of WASP-69 b using both high- and low-resolution instruments. Casasayas-Barris et al. (2017) investigated the sodium absorption using high-resolution spectra collected with the High Accuracy Radial velocity Planet Searcher-North (HARPS-N) at the 3.56 m Telescopio Nazionale Galileo. Through comparison between the in-transit and out-of-transit spectra, these latter authors detected an excess absorption only at the Na D₂ line in the passband of 1.5 Å, at the level of 5σ. Later, Deibert et al. (2019) searched for sodium, potassium, and water in the atmosphere of WASP-69 b by cross-correlating model spectra with the observed transit spectra taken with the Gemini Remote Access to CFHT ESPaDOnS Spectrograph (GRACES). Their analysis did not show any significant signal (> 3σ) of these species. Additionally, using the high-resolution transit spectra taken with the CARMENES instrument, Nortmann et al. (2018) detected excess absorption of the helium triplet in the near-infrared regime (He I λ10 830 Å) with a signal-to-noise ratio (S/N) of 18. These authors also observed a post-transit absorption together with a line blueshift of several kilometers per second, which they interpreted as escape from the upper atmospheric layers. Recently, Murgas et al. (2020) obtained optical broad-band spectroscopy of WASP-69 b by conducting transit observations with the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) at the 10.4 m Gran Telescopio Canarias (GTC). Their transmission spectrum did not show any signs of excess absorption at the sodium doublet wavelengths. However, the slope of the spectrum obtained by these latter authors was consistent with a hazy atmosphere that causes atmospheric scattering. In addition to these ground-based studies, Tsiaras et al. (2018) used the Wide Field Camera 3 (WFC3) at the Hubble Space Telescope (HST) in the near-infrared and homogeneously studied a sample of giant exoplanets, including WASP-69 b. Their low-resolution transmission spectra revealed water vapor features.

In the present paper, we use high-resolution transit spectra of CARMENES (the same observations as Nortmann et al. 2018 plus a new observation from August 2020) and mainly search for exoplanet atmospheric features in the visible region of the spectrum (e.g., sodium, potassium, and excited atomic hydrogen). To increase the S/N of the transmission spectrum, we combine our final results from CARMENES with those of HARPS-N (Casasayas-Barris et al. 2017). We also compare the three individual nights of CARMENES observations and evaluate the level of stellar activity. To obtain a more complete picture of the atmosphere, we perform atmospheric retrieval modeling on both low- and high-resolution transmission spectra. For the low-resolution part, we make use of the transmission spectra of WFC3/HST (Tsiaras et al. 2018) and OSIRIS/GTC (Murgas et al. 2020) already available in the literature, and look for water features in infrared (IR).

The structure of the paper is as follows. In Sect. 2, we present the high-resolution transit observations with CARMENES and describe our initial data-reduction procedure. We also describe our long-term photometric observations and analysis with the STELLA telescope. In Sect. 3, we analyze the high-resolution transmission spectra around prominent atomic features in the optical and near-infrared spectral regions and we discuss the stellar chromospheric activity. At the end of this section, we discuss our observational results and a measure of the statistical robustness of our signal. In Sect. 4, we perform atmospheric retrievals on both low- and high-resolution data and discuss our findings. Finally, in Sect. 5 we summarize our findings and provide our conclusions.

2 Observations and data reduction

2.1 Transit spectroscopy with CARMENES

Three transits of WASP-69 b starting on 2017 August 22, 2017 September 22, and 2020 August 13 were observed with CARMENES (Calar Alto high-Resolution search for M dwarfs with Exoearths with near-infrared and optical Echelle Spectrographs; Quirrenbach et al. 2014, 2018) in the two visible (VIS) and near-infrared (NIR) channels. The details of the first two observations and the data reduction are explained by Nortmann et al. (2018). The third night also followed the same observing strategy.

Below, we first summarize the information on the observations and the initial data reduction, and then explain our own data analysis of the CARMENES VIS data in detail. The target star was placed in fiber A, while the sky was monitored with fiber B at 88 arcsec to the west. The observation from the first night consisted of 35 spectra: 18 taken before and after the transit (out-of-transit) and 17 during the transit (in-transit). On the second night, we obtained 31 spectra, consisting of 14 out-of-transit and 17 in-transit exposures. Finally, on the third night we captured 34 exposures, of which 17 were out-of-transit and 17 were in-transit. The seeing averages on the first, second, and third nights were 0.63, 0.69, and 1.11 arcsec, respectively. The exposure time for each frame was 398 s in all cases, which was short enough to keep the change in planetary velocity to less than 1 km s⁻¹. We used a mean CARMENES pixel equivalent-velocity size of 1.3 km s⁻¹ so as to not smear the signal. The airmass change on the first observing night is plotted in the top panel of Fig. 1 (airmass changes on the other two nights are shown in Fig. B.1). The initial raw data frames were processed through the CARMENES pipeline, caracal (CARMENES Reduction And Calibration; cf. Caballero et al. 2016). The spectra were alignedwith the radial velocities measured by serval (Zechmeister et al. 2014); the barycentric velocities, and the systemic velocity γ (Table 1) were additionally taken into the account. The values of shifts from these two sources on the first night are shown in the bottom panel of Fig. 1 (see also Fig. B.1 for the second and third nights). The average S/N of the VIS spectra are 52, 40, and 37 for the first, second, and third night, respectively.

Absorption features of the Earth’s atmosphere were corrected using telluric models generated with Molecfit (Smette et al. 2015; Kausch et al. 2015). In addition, to investigate the telluric emission features superimposed on the Na I D₁ and D₂ doublet, we used sky spectra collected simultaneously with the CARMENES fiber B. Figure 2 illustrates the sky emission spectra, which shows telluric sodium emission lines right in the center of the stellar sodium lines in the average of all spectra.

These emission features are airmass-dependent and can therefore be related to the Earth’s atmospheric sodium (Chapman 1939; Kolb & Elgin 1976; Noll et al. 2012) or they can be contamination by low-pressure sodium street lamps (the lights from the nearby cities of Granada and Almería are visible from Calar Alto and are sometimes reflected by high, thin cirri). The telluric emission on the second observing night was a few times stronger than that on the first night. In addition, there was an instrumental contamination of flux in fiber B due to cross-talk with fiber A, which appeared in the form of a scaled spectrum of fiber A on the spectrum of fiber B. This contamination can be seen in the sky emission spectra in the lower panels of Fig. 2 (the continuum is not straight and is similar in shape to the stellar spectral profile). Therefore, we were not able to simply remove the sky emission features by subtracting each spectrum with the corresponding fiber B spectrum. Instead, we first ignored seven spectra with airmass larger than 1.6, and then masked the emission feature as explained in Sect. 3.2.2. We conclude that the sky emission effect in the first night (after removal of those with airmass larger than 1.6) is negligible, because our final outcomes, before and after masking, are the same.

To normalize the spectra we use a method similar to that used by Khalafinejad et al. (2017). We first computed an averaged spectrum, and then, to obtain the continuum profile, we found the maximum data points in wavelength bins of ~10 Å, ignoring the region around the targeted line (e.g., for Na D lines from ~5887 to 5902 Å and for Hα from ~6563 to 6567 Å). To correct for outliers or any remaining spike, we substituted any point deviating by more than 2σ with the average of its neighboring points. We then linearly interpolated the maximum points to construct the shape of the continuum and divided each individual spectrum by this interpolated curve to flatten the averaged spectrum. After this stage, we again normalized each spectrum by dividing each individual spectrum by the median of the points in the continuum. We performed the analysis for all nights in parallel. However, the second night was highly contaminated by telluric features (see above), and the S/Ns of the second and third nights were lower than that from the first night. We therefore dropped the second and third nights from the final modeling stage, although the results from these nights are shown in Appendix B for comparison.

Fig. 1 Airmass (top) and sum of the radial (stellar) and barycentric velocities (bottom) for each exposure during the first observing night with CARMENES. In this analysis, we discarded seven exposures with airmass higher than 1.6 to reduce the effects of contamination by the telluric features. The dashed lines represent the times of first, second, third, and fourth contact. This figure can be compared with Fig. B.1.

2.2 Long-term photometry with STELLA

The host star WASP-69 was found to show signs of stellar activity by Anderson et al. (2014). The presence of spots in the stellar photosphere causes periodic variations when they rotate in and out of view. The stellar rotation period was estimated to be about 23 d and the photometric variation showed a semi-amplitude of about 1% (Anderson et al. 2014). To characterize the presence of spots on the visible stellar hemisphere during our spectroscopic transit observations, we initiated a long-term photometric monitoring campaign covering the first and second transit observations.

We employed the robotic 1.2 m STELLA telescope (Strassmeier et al. 2004) located at the Observatorio del Teide in Tenerife (Spain), and its wide-field imager WiFSIP (Weber et al. 2012). Observations were taken from May to September 2016 (53 nights) and from June to November 2017 (40 nights) in nightly blocks of three 8 s exposures in Johnson B and three 6 s exposures inJohnson V. The data were reduced and the light curves extracted following the procedure described by Mallonn et al. (2018). In short, the image frames were trimmed, and bias- and flat-field corrected with the standard STELLA pipeline. We then performed aperture photometry using SOURCE EXTRACTOR (Bertin & Arnouts 1996). We tested different aperture shapes (circular versus flexible elliptical) and different aperture sizes to finally employ the aperture minimizing the dispersion in the light curves. We built an artificial comparison star for target differential magnitudes as the flux sum of multiple individual sources. The final choice was again made by a minimization of the light-curve scatter. We verified that the comparison star selection had no significant influence on the photometric variability signal of WASP-69 by manually testing different sets of comparison stars.

In Fig. 3 we show the final light curve after binning the three exposures per filter of an individual observing block. During the second night the star appeared dimmer compared to the first night. In other words, the spot filling factor on the visible stellar hemisphere was smaller on night 1 than on night 2. One could argue that the first night of observation was obtained during a period of maximum brightness of the star, which implies that the spots were at minimum. Unfortunately, we were not able to carry out any photometric observations simultaneous to the third CARMENES observing night in 2020.

Additionally, we re-measured the stellar rotation period. For that, we employed the generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) and the corresponding false-alarm probabilities (FAPs) to search for significant peaks from the available photometric B- and V -band light curves from both 2016 and 2017. We obtained peaks at 24.34 d and 24.95 d above the 0.1% FAP threshold for the 2016 B- and V -band data sets, respectively. We interpreted the ~24 d signal as the stellar rotation period, as it is very close to the 23.07 ± 0.16 d period reported by Anderson et al. (2014). As a reference, periodograms and phase-folded light curves related to this observation are shown in Fig. A.1. In comparison, for the 2017 data sets, the GLS periodograms show peaks at 48.24 d above the 10% FAP level in B band and above the 0.1% FAP level at 11.70 d in V data, which could be considered as the first harmonics of the ~24 d signal. The one-day difference between the measured rotation periods can probably be attributed to differential rotation.

Fig. 2 Left panel, top: average CARMENES stellar spectrum on the first night around the Na I D1 and D2 doublet. Left panel, bottom: average of all sky spectra (uppermost spectrum in black), and the lowest and highest airmass points, respectively (second (blue) and third (red) spectra, respectively). To avoid any offset between the wavelength solutions in fibers A and B, we use the pixel grid instead of wavelength for the x-axis. The profiles are shifted by an arbitrary offset for clarity. Right panels: same as left panels but for the second night.

Fig. 3 Normalized flux of WASP-69 in B (blue) and V (green) filters from our long-term photometric observations with STELLA. The dashed lines mark our two first CARMENES observing nights.

3 Transmission spectroscopy

In this section we describe the data analysis of the high-resolution transmission spectroscopy performed with CARMENES (for all of the three nights) around the strong lines of Na I D₁ and D₂, Hα, Ca II infrared triplet (IRT), and K I λ7699 Å. At the end of this section we discuss the results related to the detected signals.

3.1 Stellar activity and excess light curves

The time evolution of the equivalent widths of some prominent stellar lines (e.g., Ca II H&K, Hα, Ca II IRT) is an indication of stellar chromospheric activity (e.g., Linsky et al. 1979; Notsu et al. 2013; Klocová et al. 2017). To investigate the possible influences of the stellar activity on the exoplanetary signals, we first investigated the temporal evolution of the core of the Hα, Ca II IRT, and Na D lines, and compared them to each other. In this section, we use the term “light curve” to refer to the time-series of the line cores (line indices). For this purpose, we first look for the center of each line by fitting a Gaussian, integrating the flux in a 1.5 Å pass-band centered on the core of the line, and dividing it by the integrated flux of the reference band in the continuum (Khalafinejad et al. 2018). In principle, a Lorenzian profile can better explain the line shape compared to a Gaussian profile; however, due to the relatively large scatter of the data here, neither one has an advantage over the other. We normalized the light curves by dividing the data points by the median of the out-of-transit level. We estimated the uncertainty with the scatter of the out-of-transit data. The time-series obtained by this method are known as “excess light curves” and they have the same pattern as the time evolution of the equivalent width of thelines. Figure 4 shows the corresponding excess light curves for the first night for the Hα, Ca II IRT, and Na D lines individually. Excess light curves for nights 2 and 3 can be seen in Fig. B.2. The star did not show any strong chromospheric flare during any of the nights. The excess light curves for the Ca II IRT together with He D₃ were also shown by Nortmann et al. (2018) for the first and the second nights. Our re-analysis of the Ca II IRT in these nights agrees with their result.

However, the excess light curves display a mild level of stellar activity that reached a maximum value of about 1% in the cores of Hα and the Ca II IRT. On night 1 (Fig. 4), the Na D lines behaved similarly to Hα and the Ca II IRT. On nights 2 and 3, as shown in Fig. B.2, the Na D lines showed larger fluctuations compared to Hα and the Ca II IRT. One possible explanation for these fluctuations is a lower S/N at 5892–5898 Å compared to the other lines at redder wavelengths. In general, the excess light curves do not show any clear atmospheric absorption trend, and so we used other approaches for further investigation of the exoplanetary atmosphere.

Fig. 4 Relative flux time evolution of stellar line cores of Na I D2 (top panel, orange) and D1 (bottom panel, orange) compared with Hα (gray) and the Ca II IRT (blue) on night 1. The vertical dashed lines mark the times of first, second, third, and fourth contact, andthe horizontal dashed line is unity. This figure can be compared with Fig. B.2.

3.2 Transmission spectroscopy at Na I D₁ and D₂

3.2.1 Division approach at Na D

To obtain the transmission spectrum using high-resolution data, a simple method is to divide the average of in-transit spectra by the average of out-of-transit spectra, called “master-out”. However, because the planet is moving around the star, this approach must be corrected by taking into account the radial velocity of the planet. We obtained the transmission spectrum of WASP-69 b using the method described in detail by Khalafinejad et al. (2018). Briefly, we divided each in-transit spectrum by the master-out and then aligned the residual spectra based on the radial velocity of the exoplanet. At the end, we co-added all residual spectra, computed the average, and considered it as the exoplanetary transmission spectrum for each transit observation (Figs. 5 and B.4). In the transmission spectrum of the first night, shown in the middle panel of Fig. 5, the core of the Na I D₂ line shows an absorption at the level of about 2%. However, Na I D₁ does not show a signature of any considerable absorption. For comparison, we also show the stellar spectrum in the first panel, and the transmission spectrum by Casasayas-Barris et al. (2017) from HARPS-N observations in the bottom panel of Fig. 5. The results from the first night of observation with CARMENES agree well with those from the HARPS-N data. The error bars on each data point were obtained by propagating the errors associated with the raw spectra. For the investigation of the physical properties of the atmosphere, we merged the transmission spectrum of the first CARMENES night with that of the HARPS-N data (see Sect. 4). As mentioned in Sect. 2.1, the second night was highly affected by the telluric features and in the third night the S/N was relatively low. Finally, the transmission spectrum of the second night shows large scatter with no sign of exoplanetary absorption. The results from the third night show possible signs of exoplanetary absorption, but with more scatter than the first night. Hence, we exclude these nights from our atmospheric modeling in Sect. 4.2.

3.2.2 Two-dimensional map of residuals at Na D

Another method for visualizing the absorption features of exoplanetary atmospheres is to investigate the matrix of residual spectra obtained from the division of each individual spectrum by the master-out stellar spectrum. In Fig. 6, we show the map of residual spectra at each sodium line for the first night of CARMENES observations (see Fig. B.3 for the second and third nights). The expected exoplanetary signatures are those beneath the dotted diagonal lines. The vertical feature in the background in each panel is the lower S/N region originating from the core of the stellar sodium lines. We did not see a clear trail of exoplanetary Na D absorption in these maps. The lowest pixels could be due to the telluric emission or the variability of stellar emission core and we needed to make sure that the entire signal in the transmission spectrum (obtained through the division approach, see Sect. 3.2.1) was not based on only these few “potential outliers”. Hence, we removed the exposures containing very low pixels on the expected exoplanetary trail and repeated the analysis. We observed no significant changes in the shape of the transmission spectra after this procedure. In addition, we completely masked the narrow band contaminated by the telluric emissions based on Fiber B spectra. The masked region is shown under the shaded area in Fig. 6. Masking the whole region did not change the general shape of thetransmission spectrum or the absorption level; it simply led to larger error bars in the transmission spectrum for the points in the core of the line, because fewer residual spectra formed this region. This confirmed that the skyemission effect (after removal of the seven last spectra with the highest airmass) was within the error bars of the signal. At the end, we did not mask the region to keep the highest possible number of in-transit spectra.

Considering that the 2D map does not reveal any exoplaneary signature, the best way to detect a signal is to co-add the signal in all in-transit exposures, as explained in Sect. 3.2.1.

3.3 Transmission spectrum at Hα and other lines

For all nights, we investigated the exoplanetary Hα line using the division approach and the 2D map of residuals, as for the Na D lines. We found no signatures of Hα except for the third night, which displayed a weak signature (Fig. 7). Both the transmission spectrum and the matrix plots show hints of signatures of Hα absorption for that night. The 2D map shows a relatively uniformly distributed clump of low pixels on the expected exoplanetary trail, which is similar to the He I λ10830 Å absorption feature shown by Nortmann et al. (2018, their Fig. S2). In the transmission spectrum, there is also a signature of absorption at Hα. However, we refrain from claiming a detection because fluctuations at a similar level are visible in the vicinity of this line. In addition, no feature was detected at the Hα line in the other nights, particularly in the first night that had the highest S/N.

With a method similar to the analysis of the Na D and Hα lines, we investigated the Ca II infrared triplet (IRT) and K I λ7699 Å lines. Using the division approach and the 2D representation of the residual spectra in the vicinity of these lines, we observed no sign of absorption for any of the lines.

Fig. 5 Top panel: stellar spectrum (telluric absorption corrected) of WASP-69 around the Na I doublet. Middle panel: planet transmission spectrum of WASP-69 b in the same wavelength range. The gray data points show the division of in-transit by out-of-transit data, corrected for the changing orbital velocity of the planet. The red filled circles are the data binned by 0.25 Å and the red error bars are the standard deviation of each group of binned data points. Bottom panel: planet transmission spectrum obtained with HARPS-N by Casasayas-Barris et al. (2017), shown for comparison. Both CARMENES and HARPS-N data show a similar pattern.

Fig. 6 Two-dimensional map of residuals for both Na I D2 (left) and D1 (right) lines for the first night. The green horizontal dashed lines show the times of first, second, third, and fourth contact and the diagonal dotted line represents the expected trail of the exoplanetary atmospheric absorption. In each panel, the vertical axis is time (similar to the exposure number). Blue and red colors stand for residual fluctuations in the direction of absorption and emission, respectively. The tested masked region is under the shaded area.

3.4 Signal strength at Na D

In order to compare our results with previous studies, we applied a Gaussian toy model to the Na D doublet region of our transmission spectrum to estimate the strength of the absorption signal. We used the sum of two Gaussian profiles centered on the lines of the Na D doublet and applied a Markov chain Monte Carlo (MCMC) procedure using emcee (Foreman-Mackey et al. 2013) to determine the best model parameters and their uncertainties. The details of this procedure including the prior and posterior values are described in Appendix C.

We measured absorption signals of 0.7 ± 0.1% at Na I D₂ and 0.3 ± 0.1% at Na I D₁, both in a passband of 1.5 Å. In other words, we detected a 7σ-level excess absorption at Na I D₂ and a 3σ-level excess absorption at Na I D₁. For estimating the D₂ /D₁ ratio, we took the product of σ and amplitude for each line and obtained a D₂ /D₁ ratio of 2.5 ± 0.7.

In our analysis, the significance of the detection at sodium D₁ is less than for D₂. This is consistent with the results of Casasayas-Barris et al. (2017), who do not detect an absorption at Na D₁. In general, the absorption cross-sectionof the D₂ line is up to two timeslarger than in D₁ (Draine 2011; Gebek & Oza 2020).

In our case, the D₂ /D₁ ratio is about 2.5 ± 0.7, which tends to be closer to the upper limit of what we expect. We further investigate the physical properties of the atmosphere inferred from these lines in Sect. 4.2. The reason for the non-detection of sodium by Deibert et al. (2019) could be related to the fact that the signal at the sodium D₁ line is relatively weak and the stronger signal from the sodium D₂ line is not enough to appear in the cross-correlation method.

We ignored the center-to-limb variation and the Rossiter-McLaughlin effect. Based on Casasayas-Barris et al. (2017), the amplitude of the former is only a small fraction (about one-tenth at 1.5 Å) of the expected absorption depth. In addition, the star is not a fast rotator (v sini = 2.2 ± 0.4 km s⁻¹) and the planetary orbital velocity (~127 km s⁻¹) is much faster than stellar rotation. Therefore, Casasayas-Barris et al. (2017) did not find a significant Rossiter-Mclaughlin effect.

Fig. 7 Left: transmission spectrum around Hα for the third night. The black line is the scaled and shifted stellar spectrum, the gray dots with gray error bars are the residuals of the division, and the red filled circles with red error bars are the data binned by 0.3 Å (here, the error bars are the standard deviation of each group of binned data points). Right: same as Fig. 6 but for Hα.

3.5 Signal significance at Na D

To check the probability of a false positive in our transmission spectra, we applied a bootstrap (or empirical Monte Carlo) analysis similar to those implemented by Redfield et al. (2008) and Seidel et al. (2020). We extracted subsamples from both in-transit and out-of-transit spectra and randomly made three sets of transmission spectra. The first group was the division of a random in-transit spectrum by a random out-of-transit spectrum (“in-out”), the second group a random division of an out-of-transit spectrum by an out-of-transit spectrum (“out-out”), and the third group a random division of an in-transit spectrum by an in-transit spectrum (“in-in”).

Similar to the procedure in Sect. 3.4, we fitted a Gaussian function to each of the division residuals and estimated the absorption signal at theD₂ line (i.e., the area under the Gaussian). For a robust detection, the generated transmission spectra in the second and third groups should showno absorption, meaning that their distribution of the absorption signal must be around zero. However, we expect the absorption signal distribution for the first group to be negative.

The distribution of the absorption signals after 20 000 iterations is shown in Fig. 8. By measuring the mean of these distributions, we measured the peak of the “in-out” absorption signal at the negative value of –0.0030 ± 0.0001, while the other two groups were distributed around zero with the mean value of 0.0002 ± 0.0001 for “out-out” and 0.0005 ± 0.0002 for “in-in”. We considered the false-positive likelihood as the standard deviation of the “out-out” distribution multiplied by the square root of the fraction of out-of-transit spectra of the total number of spectra to account for the sample selection (Redfield et al. 2008; Astudillo-Defru & Rojo 2013; Wyttenbach et al. 2015; Seidel et al. 2020). In this case, our false-positive likelihood was 0.4%, which was still small enough to conclude that we measured a genuine absorption signal from the planet.

4 Atmospheric modeling

We characterized the atmosphere of WASP-69 b using a two-step retrieval approach. First, we performed a low-resolution atmospheric modeling and then we used the output parameters of the best-fit model, including the continuum level information (i.e., haze opacity, atmospheric temperature, and base pressure) as priors for our high-resolution modeling around the Na D features.

4.1 Low-resolution atmospheric modeling

To date, several groups have modeled the WASP-69 b atmosphere in retrieval. Tsiaras et al. (2018) and Fisher & Heng (2018) used the WFC3/HST observations and retrieved the planetary (isothermal) temperature, radius, water abundance, and cloud properties. Their conclusions were similar, showing a radius close to that of Jupiter, a muted water feature, the presence of non-gray aerosols, solar water abundance, and a relatively low atmospheric temperature (492 ± 53 K and 658 K, respectively in each work) consistent with a high albedo. Recently, Murgas et al. (2020) combined these WFC3/HST data with their OSIRIS/GTC optical observations and explored Rayleigh scattering in addition to the atmospheric carbon-to-oxygen ratio (C/O) and metallicity. These authors also investigated the activity of WASP-69 and the presence of stellar spots, testing whether the observed slope seen in the data (the increase in the apparent planet-to-star radius ratio towards the blue wavelengths) could be a consequence of unocculted stellar spots instead of aerosol particles present in the planetary atmosphere (Pont et al. 2013; McCullough et al. 2014). Their analyses showed that the slope seen in the data is more likely caused by hazes present in the atmosphere of WASP-69 b rather than by stellar active regions, revealing a planet with strong Rayleigh scattering, a subsolar C/O ratio, high metallicity, and a super-solar water abundance. Furthermore, they retrieved much higher temperatures than previous studies (1203 K and 1227 K for the cases with and without stellar spots, respectively), which are difficult to reconcile with standard models of planetary atmospheres, even assuming a cloud-free atmosphere (Carone et al. 2021, see also Table 1 in Sect. 1).

Fig. 8 Distributions of the empirical Monte Carlo bootstrap analysis of the core of Na D features on the first night. The “in-in” (blue) and “out-out” (red) distributions are centered around zero (no sodium detection), but the randomized in-out (green) distribution shows a detection. The dashed lines show the center of the Gaussian fit over each distribution.

4.1.1 Modeling setup

To assess the atmospheric properties of WASP-69 b from low-resolution spectroscopy, we performed several retrieval analyses with PyratBay (Python Radiative-transfer in a Bayesian framework, Cubillos & Blecic 2021; Blecic et al., in prep.) on the joint WFC3/HST(Tsiaras et al. 2018) and OSIRIS/GTC (Murgas et al. 2020) datasets. PyratBay is an open-source retrieval framework that generates 1D atmospheric models of planetary temperatures, species abundances, altitude profiles, and cloud coverage, using up-to-date knowledge of the atmospheric physical and chemical processes and alkali, molecular, collision-induced Rayleigh and cloud opacities. The code can work in both forward and retrieval modes, utilizing self-consistency or different parameterization schemes when including atmospheric processes. For this study, we applied basic physical assumptions without exploring the hierarchy of models available in the PyratBay framework. We generated several modeling setups for the purpose of providing the baseline information and inputs for the high-resolution analysis presented in Sect. 4.2. As part of a coordinated effort, Blecic et al. (in prep.) will present detailed forward and retrieval analyses covering different levels of model complexity and discuss their effect on the WASP-69 b spectrum.

In our analysis, we also combined the WFC3/HST (Tsiaras et al. 2018) and OSIRIS/GTC (Murgas et al. 2020) data and explored additional models. We performed two retrievals using two different datasets. Both retrievals included the WFC3/HST data from Tsiaras et al. (2018) (denoted WFC3 in the following text), but the first included only the OSIRIS/GTC data from Murgas et al. (2020) Table A.1 (denoted O1 in the following text), while the second also included the data from their Table A.2 listing the observations around the Na line (denoted as O2). Apart from the continuum information, we also wanted to put constraints on the sodium abundance, which we could then use as inputs for the high-resolution analysis with CARMENES and HARPS-N data (Sect. 4.2). To constrain the sodium volume mixing ratio, we calculated the sodium resonant-line opacities following the formalism from Burrows et al. (2000), using resonance Na D doublet line parameters from the Vienna Atomic LineData Base (Piskunov et al. 1995).

Following the conclusion of Murgas et al. (2020), we assumed that the slope seen in the WFC3/HST and OSIRIS/GTC data comes from the aerosol particles present in the planetary atmosphere and not from unocculted stellar spots. We therefore included a Rayleigh scattering model (Lecavelier Des Etangs et al. 2008) that allowed us to explore and constrain both the enhancement factor (f_Ray) and the power-law index (α_Ray) of the scattering cross-section opacity. We also included a simple cloud model allowing MCMC the option to set a fully opaque cloud deck below a certain pressure (a free parameter of the model, p_top).

Although several molecules have spectral features in the WFC3/HST wavelength bands in the hot-Jupiter temperature range (MacDonald & Madhusudhan 2017), in this analysis we considered only water as the main molecular source of opacity and retrieved its volume mixing ratio (abundance). Our initial tests including all relevant species in the model revealed that additional species have no influence on the WASP-69 b spectra, and that the water-only model is significantly favored, with the probability ratio of e^ΔBIC∕2 ≃ 10¹⁸. This choice of near-infrared absorbers did not affect the continuum at the Na D doublet wavelengths, which is important for our high-resolution analysis.

Water is one of the most abundant species at the WASP-69 b equilibrium temperatures explored in Table 1 (see also Blecic et al. 2016). It also has the most pronounced spectral signatures at the wavelengths of our observations. To include water opacity, we used the ExoMol molecular line-list data from Polyansky et al. (2018). As this database consists of billions of line transitions, we applied the Repack package (Cubillos 2017) to extract only the strongest lines that dominate the opacity spectrum between 300 and 3000 K. This approach causes a difference of less than 1% in dex compared to the original line lists (Cubillos & Blecic 2021), but significantly improves the computational speed. Together with water and sodium opacities, our models also included collision-induced absorption from H₂ -H₂ from Borysow et al. (2001) and Borysow (2002), and H₂ -He from Richard et al. (2012).

In transmission, stellar rays travel through the atmospheric limb probing only the high altitudes (low pressures) of the planetary envelope, usually down to the 0.1 bar level. Additionally, low-resolution observations over a narrow wavelength range such as those used in our analysis do not carry enough information to allow us to further constrain the complex temperature profiles (Rocchetto et al. 2016; Barstow et al. 2013). Thus, when retrieving the atmospheric temperature, we assumed an isothermal temperature profile.

In all retrievals, we also included a planet-to-star radius ratio offset parameter between the WFC3 and OSIRIS observations, δ_WFC3-OSIRIS. Offset values, on the order of those retrieved by Murgas et al. (2020, see their Sect. 5) have a strong and often crucial influence on the retrieved atmospheric properties. There are several possible reasons why such a δ_WFC3-OSIRIS offset could appear between these observations: it could have originated from different system parameters adopted during the WFC3 and OSIRIS data reduction procedures (Alexoudi et al. 2018, 2020), from different instrumental systematic errors that were not fully addressed (e.g., May & Stevenson 2020; Gibson et al. 2012), and/or from stellar activity. Murgas et al. (2020) explored possible values of the δ_WFC3-OSIRIS offset due tothe different system parameters, and calculated the remaining offset assuming that it comes from stellar activity. Our analysis in Sect. 2.2, Fig. 3, and Sect. 3.1 also confirms activity in WASP-69 star.

To address this, we initially allowed δ_WFC3-OSIRIS to be a free parameter in our model. However, none of our retrievals was able to put constraints on this parameter, at the same time causing multiple correlations between other free parameters (we explore these correlations and degeneracies in our following paper Blecic et al., in prep.). Here, we adopted δ_WFC3-OSIRIS = 479 ppm (we shift the WFC3 data upwards), which corresponds to the offset retrieved by Murgas et al. (2020, see their Table 4) for their case without stellar spots. Our decision to choose this value is based on the fact that it is comparable to the variability of the stellar flux of ~2.4% found by Anderson et al. (2014), assuming that such variability occurs in regions not occulted by the planet. We note that by adopting a fixed offset in this analysis, we neglected some uncertainties caused by correlations and degeneracies.

Our model atmosphere had 50 equally spaced layers in log pressure, ranging from 10² bar to 10⁻⁸ bar, and we included all molecular species relevant for the hydrogen-dominated atmospheres with solar C/O ratio (H₂, He, H₂O, CH₄, CO, CO₂, HCN, N₂, NH₃), as well as sodium. We used the TEA code (Blecic et al. 2016) to calculate the chemical equilibrium abundances of all the species at 1000 K and used them for the initial model atmosphere assuming that species abundances remain constant with altitude. In retrieval, we leave H₂O abundance as the only free parameter in our model, while we kept other species abundances at initial values. We set the transit planetary radius at 0.01 bar pressure level and calculated the altitude profile assuming hydrostatic equilibrium. The system parameters used to generate the spectra are listed in Table 1.

We constrained the atmospheric parameters listed in Table 2 assuming wide uniform or log-uniform priors for all parameters. We explored the posterior parameter space with the Snooker differential-evolution MCMC code (DEMC, ter Braak 2006; ter Braak & Vrugt 2008), and obtained between two and six million samples (after discarding the initial 10 000 iterations). This ensured that the Gelman & Rubin (1992) statistics remained at ~1.01 or lower for each free parameter. To determine the best-fit model among those generated on the same dataset, we used the Bayesian information criterion, BIC = χ² + k lnN.

4.1.2 Retrieval results

Figure 9 shows the best-fit spectrum model for the case of WFC3 and O1 data. In this case, we retrieved the planet temperature, radius, water abundance, Rayleigh slope, and strength factor, and the top pressure of an opaque cloud deck. The mean parameter values with 1σ uncertainties are given in Table 2, while the best-fit parameter values are shown in the posterior density plots in the left panel of Fig. 10. In the right panel of Fig. 10, we show the pair-wise correlation plots with posterior densities. Figure 11 shows the best-fit model and posterior density plots for the case of WFC3 and O2 data. In this case, in addition to the parameters listed above, we also retrieved the sodium abundance. Figures 9 and 11 show good fits to the data for both cases.

We note a steep Rayleigh slope from ~ 1.2μm towards shorter wavelengths. The transmission functions confirm that the atmosphere becomes opaque above 10⁻⁴ bar for the near-infrared observations, and above 10⁻⁵ bar for the optical observations. This behaviour is expected, as Rayleigh scattering becomes very pronounced below 1 μm when aerosol particles with sizes smaller than or equal to the wavelength of light cause the atmosphere to become opaque for the incoming radiation. Rayleigh opacity usually scales with λ⁻⁴. Mallonn & Wakeford (2017) performed a comparative study of the scattering slopes of all planets in Sing et al. (2016) together with HAT-P-32 b and GJ 3470 b, and found that the majority of them have slope values log α_Ray > −4 (between –1and –3), with a couple of outliers (HD 189733 b and GJ 3470 b) having log α_Ray < −5. Later studies by Ohno & Kawashima (2020) and May et al. (2020) found that Rayleigh slopes can indeed be steeper than –4. These latter authors discussed the idea that opacity gradients can be naturally generated in hot Jupiters with species other than H, H₂, and He, such as photochemical hazes, when eddy mixing in the planetary atmosphere is very efficient. These findings are in line with those of Pinhas et al. (2019) and Welbanks et al. (2019), who obtained a median spectral index of log α_Ray ≲−5 for different sets of hot Jupiter planets. Both of our analyses (Table 2 and Figs. 10 and 11) reveal a Rayleigh slope on the order of logα_Ray ~−5, steeper than the slope caused by the tiny aerosol particles.

None of our retrievals returned strong constraints on the planetary temperature. The left panel of Fig. 10 shows that the 1σ region spans around 800 K (red section, top left panel). The mean value, however, is close to the planetary equilibrium temperature, assuming albedo A_Bond = 0 and efficientenergy redistribution. The radius of the planet and the Rayleigh slope are most robustly constrained, revealing a planet with a radius very close to that of Jupiter and strong Rayleigh scattering.

The water abundance seems to be slightly bi-modal (solar and super-solar) in the first case (Fig. 10, O1 data), with a tendency towards large values, logH₂O ~−2, and loosely constrained in the second case (Fig. 11, O2 data), spreading over several orders of magnitude between 10⁻⁴ (solar) to 10⁻¹ (super-solar) in volume mixing ratio. WASP-69 b has a Jupiter-like radius, but a very small mass (0.26 M_Jup). Welbanks et al. (2019) showed that there is a metallicity trend of increasing water abundances with decreasing mass going from gas-giant planets to mini-Neptunes. In particular, low-mass planets such as HATP-26 b, WASP-39 b, and WASP-127b show super-solar water abundance on the order of log H₂O ~−2 and smaller. Welbanks et al. (2019) suggested that this trend could be a consequence of the different formation pathways from those usually assumed for hot Jupiter planets.

In the second case (O2 data), we also added sodium as an additional source of opacity and allowed its abundance to be a free parameter of the model. This retrieval returned results akin to those in the first case, with most of the parameters having similar values and confidence ranges (Figs. 9 and 11; see also Table 2). The spectrum, temperatureposteriors, transmission functions, and correlation plots for the second case also looked comparable to the plots from the first case (Figs. 9 and 10). The sodium abundance, of interest here, is loosely constrained to log Na ≃−11. Considering that this abundance is negligible, causing very low sodium opacity, the best-fit model in Fig. 11 does not show any trace of the sodium doublet in the spectrum. Sodium, if present and abundant, would otherwise have a distinct signature at 5895 Å.

Our results (Table 2) agree well with results from previous analyses (Sect. 4.1.1), and in particular with the results of Murgas et al. (2020) who used the same datasets. Our retrieved planetary radius, temperature, and Rayleigh parameters are consistent with theirs. Murgas et al. (2020) also found a high atmospheric metallicity (together with a low C/O ratio compared to solar), causing a super-solar water abundance on the order of our retrieved value.

Fig. 9 Left: best-fit spectrum when the WFC3 and O1 are included. Black points denote the data with uncertainties, while the red points denote the best-fit spectrum integrated over the bandpasses of the observations, which are shown in gray at the bottom of the panel.

Fig. 10 Left: marginal posterior distributions (histograms) when the WFC3 and O1 data are included. The shaded red areas denotethe 68% confidence intervals of the respective distributions. The light green dashed vertical lines denote the retrieved best-fit parameter values. Strong non-linear correlations seen in the right panel of this figure produce that some of the best-fit parameter values fall outside of the red credible regions, and cause gaps in the water confidence interval, reflecting a bi-modal posterior distribution. Right: pair-wise correlation plots with posterior densities.

Fig. 11 Left: best-fit spectrum when the WFC3 and O2 data are included. Black points denote the data with uncertainties, while the red points denote the best-fit spectrum integrated over the bandpasses of the observations, which are shown in gray at the bottom of the panel. Right: marginal posterior distributions (histograms). The shaded gray areas denote the 68% confidence intervals of the respective distributions. The purple dashed vertical lines denote the retrieved best-fit parameter values. The pair-wise correlation plots are similar to those shown in Fig. 10, displaying same nonlinear correlations and degeneracies.

4.2 High-resolution atmospheric modeling

Our analysis in Sects. 3.4 and 3.5 suggested the feasibility of teasing out information from the current high-resolution datasets. Likewise, in Sect. 4.1 we show that the information content in the low-resolution spectrum is statistically significant. In this section, we discuss steps in order to characterize the upper thermosphere of WASP-69 busing atmospheric parameters retrieved from low-resolution datasets as input, and combining them with the CARMENES (first night) and HARPS-N observations.

Normalization of the continuum is a common step in the data reduction of high-resolution spectra of exoplanetary atmospheres. However, such normalization comes at the price of losing continuum information, which in turn causes additional degeneraciesin the retrieved values (Fisher & Heng 2019). A common degeneracy in high-resolution studies occurs between the retrieved abundance of the species, reference pressure, and atmospheric temperature (e.g. Brogi & Line 2019), where different combinations of these parameters could result in similar line profiles. This issue becomes particularly significant when the S/N is relatively low. However, a precise line shape at high S/N could break the degeneracies to some extent (e.g., Benneke & Seager 2012; Heng & Kitzmann 2017; Brogi & Line 2019; Gibson et al. 2020, ; see also Welbanks et al. 2019 for a similar discussion on low-resolution retrieval degeneracies). To break these degeneracies and estimate the continuum at the Na D doublet, we employed the model parameters from low-resolution spectra as estimated in the previous section (see the parameters best-fit values in histograms of Figs. 11 and 10) as priors for the high-resolution atmospheric model.

We used a Python implementation of emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013), called MCKM, to effectively sample the parameter-space. MCKM handles any arbitrary number of model parameters over regular (Cartesian) or irregularly spaced steps and any arbitrary number of data points and their respective covariance matrices. The observational uncertainties can be estimated through a Gaussian process using a variety of kernels. We assumed uninformative priors to initialize 1000 walkers in the MCMC process and iterated until the solution converged.

Table 3 summarizes the model parameters, that is, Na abundance, thermospheric onset pressure, thermosphere temperature, effective turbulent velocity, and their priors. We used the core radiative transfer solver of petitRADTRANS (Mollière et al. 2019) to calculate forward radiative transfer models and included broadening due to turbulent velocity. In order to model the effect of turbulent velocity, we assume a Gaussian velocity distribution, assuming the effective turbulent velocity as the variance of the distribution and zero bulk velocity. We set up the high-resolution model by discretizing the 1D atmosphere to 150 layers, linearly spaced in altitude. We set the lower boundary at 1 bar and the upper boundary at 10⁻⁹ bar to resolve the centroids in grazing geometry appropriate for transmission spectroscopy.

We also clipped the high-resolution data to exclude data points far away from the signal and speed up the retrieval. The clipped data are shown in Fig. 12. Including the discarded points did not change the results and the posteriors were robust against the choice of clipping wavelengths outside the used area.

To explain large D₂ /D₁ ratios (i.e., 2.5 ± 0.7; see Sect. 3.4), it has been hypothesized that D₂ and D₁ may probe different thermosphere regions with different temperatures. Other possibilities are the effect of nonlocal thermodynamic equilibrium (nonLTE) processes (Lind et al. 2011) or an optically thin sodium contribution in the composition of the thermosphere (Gebek & Oza 2020). We examined these hypotheses by fitting our radiative transfer model to the D₂ and D₁ lines separately and by comparing the posterior distributions of the resulting model parameters.

Table 3 also summarizes the posterior distributions of the parameters and Fig. 12 illustrates the best-fitted models to the data as well as our retrieved posterior probability distributions of the model parameters. We find that the retrieved parameter values are statistically consistent, which suggested that both lines may originate from the same region and at similar temperatures. Also, given our results, there was no need for nonLTE processes to explain the observed D₂ /D₁ ratio.

Our retrieved sodium abundances from the D₂ and D₁ lines suggest that the majority of their absorption occurs at optically thin regimes. This scenario is consistent with the conclusion of Gebek & Oza (2020) because our finding supports their main assumption (i.e., both D₂ and D₁ lines probe the same region and arise from the same sodium population) and also supports the optically thin condition for the sodium lines (i.e., our retrieved low abundances). However, it is unclear if the conclusion of these latter authors, and hence our comparison, would hold true if D₂ and D₁ probed different regions of the atmosphere. Moreover, Baranovsky & Tarashchuk (2004) showed that the D₂ /D₁ ratio could be significantly affected by the multiple scattering under optically thick conditions and the D₂ /D₁ ratio could decrease to very small values. Further observations are required to reach a higher S/N and to further constrain the posteriors and address this question.

The posteriors on the thermosphere onset pressure suggest that both D₂ and D₁ lines probe pressures lower than the millibar regime. Our retrieved temperatures also appear to be consistent with 6000 K as inferred for the thermosphere of the ultra-hot Jupiter WASP-189 b (e.g. Yan et al. 2020). However, WASP-69 b is not as highly irradiated as WASP-189 b and their thermospheres might differ in nature. More data are needed to reduce the uncertainties in thermosphere temperatures to ascertain more reliable thermosphere properties. It is worth noting that obtaining constraints on the thermosphere of exoplanets via high-resolution spectroscopy is a topic at the forefront of exoplanet science and perhaps obtaining some constraints on the upper atmosphere for such a relatively cool exoplanet could set the stage for subsequent investigation and the improvement in characterization of exoplanet thermospheres.

We also find that an effective turbulent velocity on the order of 5–10 km s⁻¹ is necessary to explain the line broadening in both the D₂ and D₁ lines. These values are in line with the velocities reported by Lampón et al. (2021), which suggests consistency in retrieved upper thermospheric velocities on irradiated gaseous exoplanets. However, these velocities are slightly higher than the velocities estimated by 3D general circulation models (e.g., Showman et al. 2010; Pierrehumbert & Hammond 2019; Wang & Wordsworth 2020), which might hint at the existence of other broadening mechanisms in addition to the rotation broadening. Nevertheless, all these values are considerably lower than the 40 km s⁻¹ reported by Seidel et al. (2020). Again, higher S/N spectra are required to unambiguously investigate the nature of these broadenings.

Fig. 12 Top: bayesian inference of the upper thermospheric conditions of WASP-69 b. Comparison of the best-fit model spectra of Na I D2 (left panel, shaded green area) and Na I D1 (middle panel, shaded blue area) with discarded and used high-resolution spectroscopic data in gray and orange, respectively, and retrieved temperature profiles based on the posteriors (right panel). The dashed vertical line in red marks the retrieved temperature using low-resolution spectra (Sect. 4.1). Uncertainties are marked by vertical dashed lines at the 16, 50, and 84% quantiles. Bottom: posterior probability distributionof high-resolution model parameters as well as the marginalized distribution for each parameter for Na D2 (shaded green area) and D1 (shaded blue area) in WASP-69 b. The marginalized uncertainties are marked by vertical dashed lines at the 16%, 50%, and 8 4% quantiles. In the density maps, 1σ and 2σ are given as 1 − e−0.5 ~ 0.393 and 1 −e−2.0 ~ 0.865, respectively,as is common for multivariate MCMC results.

5 Summary

In this work, we performed a combined high- with low-resolution transmission spectroscopy analysis to obtain a more complete picture of the properties of the atmosphere of WASP-69 b. Particularly, we aimed to break the degeneracy in abundance and atmospheric temperature, which is present in any ground-based high-resolution-only modeling due to lack of information on the continuum of the transmission spectrum. We presented CARMENES optical and near-infrared spectroscopy during three transits of the low-density hot Saturn WASP-69 b. We searched for planetary atmospheric features from 0.52 to 1.71 μm using the transmission spectroscopy technique, with special emphasis on the Na I D₁ and D₂ doublet. We also added to our analysis the available low-resolution data of this target, as well as the transmission spectra obtained with HARPS-N, and performed atmospheric retrieval on both low- and high-resolution data. We used the fitting parameters of our low-resolution model as priors for our high-resolution modeling.

Due to insufficient data quality and a high level of contamination by sky emission, we disregarded the transit data of the second and third nights in our analysis of the Na D lines. Using the first night data only, we measured Gaussian amplitudes of 3.2 ± 0.3% (D₂) and 1.2 ± 0.3% (D₁), and σ of 0.13 ± 0.01 Å (D₂) and 0.14 ± 0.2 Å (D₁). These values correspond to detection at 7σ and 3σ significance at Na I D₂ and D₁, respectively.Our sodium detection is in agreement with previous measurements (Casasayas-Barris et al. 2017). In addition, although strong signatures of an escaping atmosphere were observed from the upper atmosphere of WASP-69 b through investigation of the He I triplet lines (Nortmann et al. 2018), we found no clear signature of escape in the optical region of the transmission spectrum.

In the analysis of the Hα line, only the third night data showed some hints of this feature. However, further investigations with data of better quality are required to confirm this feature.

Our low-resolution atmospheric analysis revealed a Jupiter-radius planet with strong Rayleigh scattering, solar to super-solar water abundance, and a highly muted sodium feature, in agreement with the findings of previous low-resolution studies by Murgas et al. (2020). The parameters of this model related to the continuum were used in our high-resolution modeling.

In our high-resolution observations, we measured a D₂ /D₁ ratio of 2.5 ± 0.7 and investigated the possible physical conditions that lead to this result. Further, in our high-resolution atmospheric modeling, we fit the D₂ and D₁ lines individually and compared their atmospheric parameter posteriors. We find that the retrieved temperatures (T_th ~ 6000 ± 3000 K) and abundances (log Na ~ 7 ± 2) of these two lines agree with each other to better than 1σ. Therefore, the observed high D₂ /D₁ ratio is consistent with LTE conditions and the non- or weak detection of D₁ is a consequence of obscuration of this line by optical opacity (i.e., haze). We also find that an effective turbulent velocity of 5–10 km s⁻¹ is needed to explain the line broadening in both D₂ and D₁.

Acknowledgements

We thank G. Zhou, J. Seidel and A. Wyttenbach for a fruitful scientific discussion. CARMENES is an instrument at the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l’Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the Deutsche Forschungsgemeinschaft (DFG) through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. We acknowledge financial support from the DFG through priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (KH 472/3-1) and through grant CA 1795/3, NASA through ROSES-2016/Exoplanets Research Program (NNX17AC03G), the Klaus Tschira Stiftung, the European Research Council under the European Union’s Horizon 2020 research and innovation program (694513), the Agencia Estatal de Investigación of the Ministerio de Ciencia, Innovación y Universidades and the ERDF through projects PID2019-109522GB-C5[1:4] and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Instituto de Astrofísica de Canarias (SEV-2015-0548), Instituto de Astrofísica de Andalucía (SEV-2017-0709), and Centro de Astrobiología (MDM-2017-0737), and the Generalitat de Catalunya/CERCA programme.

Appendix A Long-term photometric observation

Periodograms and phase-folded light curves related to the photometric observations in 2016 in B and V bands are shown in Fig. A.1

Fig. A.1 Left panels: GLS periodograms for the WASP-69 data in B and V bands from the 2016 observing season. Horizontal lines indicate false-alarm probability levels FAP = 0.1, 0.01, and 0.001. Right panels: Light curves folded to the best periods, as noted in the plots. A sine wave for the best periods is fitted to the data. The blue and green colors show one complete cycle in B and V filters, respectively.

Appendix B Analysis of the data from the second and third nights

In our atmospheric analysis, we did not consider the data from the second and third nights of observations with CARMENES. As mentioned in Sect. 2.1, these spectra are highly contaminated by the telluric emission features in the sodium region on the second night, and the data quality and S/N are not sufficient in the visible channel for both nights. However, for reference and to present stellar chromospheric activity in the cores of Hα and the Ca II IRT, we show the results of these nights in this section. Figure B.1 shows the airmass and the spectral shifts, Fig. B.2 the time evolution of the cores of the activity indicators Hα and Ca II IRT, Fig. B.3 the residual maps around both sodium D-lines, and finally Fig. B.4 the transmission spectrum around the sodium features on each night. The residual map (Fig. B.3, bottom panels) and the transmission spectrum (Fig. B.4, bottom panel) on the third night show some tentative signatures of Na absorption, albeit at low significance, due to the relatively large scatter of the data points around the sodium lines.

Fig. B.1 Same as Fig. 1 but for nights 2 (top) and 3 (bottom).

Fig. B.2 Same as Fig. 4 but for nights 2 (first and second panels) and 3 (third and fourth panels).

Fig. B.3 Two-dimensional map of residuals (same as Fig. 6) for both the D2 (left panels) and D1 (right panels) lines for nights 2 (top) and 3 (bottom).

Fig. B.4 Same as middle panel of Fig. 5, but for nights 2 (top) and 3 (bottom).

Appendix C Gaussian toy model

To measure the absorption signal at the sodium lines and compare it with previous work, we applied a Gaussian toy model to the transmission spectra around the sodium features. This model is a combination of two Gaussian functions, each centered on one of the sodium lines. As there is a constant wavelength difference between the cores of each Na D line, we only needed one parameter for the positions of the centers of the Gaussian peaks (μ). Here, we selected μ as the center of D₂ line. We employed 200 walkers, with 1000 chains each, where the initial positions were drawn from a Gaussian distribution around our best estimates. We imposed uniform priors to all free parameters. We allowed a burn-in phase of ~50% of the total chain length, beyond which the MCMC converged. The posterior probability distribution was then calculated from the latter 50% of the chain. The priors and values of the best-fit parameters and their uncertainties are shown in Table C.1. The posterior distribution parameters are shown in Fig. C.1.

Fig. C.1 Top panel: Best-fit combined Gaussian model. Bottom panel: posterior distribution of the combined Gaussian model parameters.

The best model is plotted over the data in the top panel of Fig. C.1, and the correlation diagram of the posterior probability distribution is shown in the bottom panel of the same figure. As Fig. C.1 shows, there is a correlation between the width and amplitude parameters of each line. In addition, there is a bi-modality for D₁, related to the high noise level at this line.

References

All Tables

All Figures

Fig. 1 Airmass (top) and sum of the radial (stellar) and barycentric velocities (bottom) for each exposure during the first observing night with CARMENES. In this analysis, we discarded seven exposures with airmass higher than 1.6 to reduce the effects of contamination by the telluric features. The dashed lines represent the times of first, second, third, and fourth contact. This figure can be compared with Fig. B.1.

In the text

Fig. 2 Left panel, top: average CARMENES stellar spectrum on the first night around the Na I D1 and D2 doublet. Left panel, bottom: average of all sky spectra (uppermost spectrum in black), and the lowest and highest airmass points, respectively (second (blue) and third (red) spectra, respectively). To avoid any offset between the wavelength solutions in fibers A and B, we use the pixel grid instead of wavelength for the x-axis. The profiles are shifted by an arbitrary offset for clarity. Right panels: same as left panels but for the second night.

In the text

	Fig. 3 Normalized flux of WASP-69 in B (blue) and V (green) filters from our long-term photometric observations with STELLA. The dashed lines mark our two first CARMENES observing nights.
In the text

	Fig. 4 Relative flux time evolution of stellar line cores of Na I D2 (top panel, orange) and D1 (bottom panel, orange) compared with Hα (gray) and the Ca II IRT (blue) on night 1. The vertical dashed lines mark the times of first, second, third, and fourth contact, andthe horizontal dashed line is unity. This figure can be compared with Fig. B.2.
In the text

Fig. 5 Top panel: stellar spectrum (telluric absorption corrected) of WASP-69 around the Na I doublet. Middle panel: planet transmission spectrum of WASP-69 b in the same wavelength range. The gray data points show the division of in-transit by out-of-transit data, corrected for the changing orbital velocity of the planet. The red filled circles are the data binned by 0.25 Å and the red error bars are the standard deviation of each group of binned data points. Bottom panel: planet transmission spectrum obtained with HARPS-N by Casasayas-Barris et al. (2017), shown for comparison. Both CARMENES and HARPS-N data show a similar pattern.

In the text

Fig. 6 Two-dimensional map of residuals for both Na I D2 (left) and D1 (right) lines for the first night. The green horizontal dashed lines show the times of first, second, third, and fourth contact and the diagonal dotted line represents the expected trail of the exoplanetary atmospheric absorption. In each panel, the vertical axis is time (similar to the exposure number). Blue and red colors stand for residual fluctuations in the direction of absorption and emission, respectively. The tested masked region is under the shaded area.

In the text

	Fig. 7 Left: transmission spectrum around Hα for the third night. The black line is the scaled and shifted stellar spectrum, the gray dots with gray error bars are the residuals of the division, and the red filled circles with red error bars are the data binned by 0.3 Å (here, the error bars are the standard deviation of each group of binned data points). Right: same as Fig. 6 but for Hα.
In the text

	Fig. 8 Distributions of the empirical Monte Carlo bootstrap analysis of the core of Na D features on the first night. The “in-in” (blue) and “out-out” (red) distributions are centered around zero (no sodium detection), but the randomized in-out (green) distribution shows a detection. The dashed lines show the center of the Gaussian fit over each distribution.
In the text

	Fig. 9 Left: best-fit spectrum when the WFC3 and O1 are included. Black points denote the data with uncertainties, while the red points denote the best-fit spectrum integrated over the bandpasses of the observations, which are shown in gray at the bottom of the panel.
In the text

Fig. 10 Left: marginal posterior distributions (histograms) when the WFC3 and O1 data are included. The shaded red areas denotethe 68% confidence intervals of the respective distributions. The light green dashed vertical lines denote the retrieved best-fit parameter values. Strong non-linear correlations seen in the right panel of this figure produce that some of the best-fit parameter values fall outside of the red credible regions, and cause gaps in the water confidence interval, reflecting a bi-modal posterior distribution. Right: pair-wise correlation plots with posterior densities.

In the text

Fig. 11 Left: best-fit spectrum when the WFC3 and O2 data are included. Black points denote the data with uncertainties, while the red points denote the best-fit spectrum integrated over the bandpasses of the observations, which are shown in gray at the bottom of the panel. Right: marginal posterior distributions (histograms). The shaded gray areas denote the 68% confidence intervals of the respective distributions. The purple dashed vertical lines denote the retrieved best-fit parameter values. The pair-wise correlation plots are similar to those shown in Fig. 10, displaying same nonlinear correlations and degeneracies.

In the text

Fig. 12 Top: bayesian inference of the upper thermospheric conditions of WASP-69 b. Comparison of the best-fit model spectra of Na I D2 (left panel, shaded green area) and Na I D1 (middle panel, shaded blue area) with discarded and used high-resolution spectroscopic data in gray and orange, respectively, and retrieved temperature profiles based on the posteriors (right panel). The dashed vertical line in red marks the retrieved temperature using low-resolution spectra (Sect. 4.1). Uncertainties are marked by vertical dashed lines at the 16, 50, and 84% quantiles. Bottom: posterior probability distributionof high-resolution model parameters as well as the marginalized distribution for each parameter for Na D2 (shaded green area) and D1 (shaded blue area) in WASP-69 b. The marginalized uncertainties are marked by vertical dashed lines at the 16%, 50%, and 8 4% quantiles. In the density maps, 1σ and 2σ are given as 1 − e−0.5 ~ 0.393 and 1 −e−2.0 ~ 0.865, respectively,as is common for multivariate MCMC results.

In the text

Fig. A.1 Left panels: GLS periodograms for the WASP-69 data in B and V bands from the 2016 observing season. Horizontal lines indicate false-alarm probability levels FAP = 0.1, 0.01, and 0.001. Right panels: Light curves folded to the best periods, as noted in the plots. A sine wave for the best periods is fitted to the data. The blue and green colors show one complete cycle in B and V filters, respectively.

In the text

	Fig. B.1 Same as Fig. 1 but for nights 2 (top) and 3 (bottom).
In the text

	Fig. B.2 Same as Fig. 4 but for nights 2 (first and second panels) and 3 (third and fourth panels).
In the text

	Fig. B.3 Two-dimensional map of residuals (same as Fig. 6) for both the D2 (left panels) and D1 (right panels) lines for nights 2 (top) and 3 (bottom).
In the text

	Fig. B.4 Same as middle panel of Fig. 5, but for nights 2 (top) and 3 (bottom).
In the text

	Fig. C.1 Top panel: Best-fit combined Gaussian model. Bottom panel: posterior distribution of the combined Gaussian model parameters.
In the text