© ESO 2022

1. Introduction

Radial velocities (RVs) and photometry are the most efficient methods used to detect and characterize exoplanetary systems¹. In particular, the first exoplanet detected orbiting a main-sequence star (using the RV method, Mayor & Queloz 1995) pushed forward the development of new state-of-art spectrographs (e.g., Mayor et al. 2003; Pepe et al. 2021). The improvement of instrumental and observational capabilities made it possible, not only to detect planets, but also to characterize them and, in particular, their atmospheres (e.g., Sing et al. 2009; Wyttenbach et al. 2015; Nikolov et al. 2018; Ehrenreich et al. 2020).

During the transit of an exoplanet, the stellar radiation is filtered by the day-night terminator of the exoplanet’s atmosphere. When crossing the stellar disk, the radiation is blocked by opacity sources such as atoms, molecules, and small dust particles, which mimic an increase in the radius of the planet (e.g., Barman 2007; Seager & Deming 2010). By measuring the radius as a function of the wavelength, it is possible to recover the transmission spectrum of the planetary atmosphere (Seager & Sasselov 2000).

Transmission spectroscopy has now been broadly used to detect both narrow- and broadband features (e.g., Kreidberg et al. 2014; Sing et al. 2016; Hoeijmakers et al. 2019, 2020; Allart et al. 2020; Ehrenreich et al. 2020; Casasayas-Barris et al. 2021; Borsa et al. 2021). A notable example is the claimed discovery of the sodium doublet signature in the atmosphere of the exoplanet HD 209458b (Charbonneau et al. 2002) with the NASA/ESA Hubble Space Telescope (HST), which opened a new era in the study of exoplanetary atmospheres². From the ground, the first molecular signature for the same element was detected by Redfield et al. (2008) from the optical transmission spectrum of the exoplanet HD 189733b. Furthermore, other molecules and atomic species such as H₂O, CO, CH₄, CO₂, Li, K, Fe, or TiO have been observed, both from space and the ground, in several exoplanets (Wilkins et al. 2013; Wakeford & Sing 2015; Sing et al. 2016; Tabernero et al. 2021).

In broadband, Rayleigh scattering is responsible for the increase in opacity in the blue visible region of a spectrum, as a result of the interaction of small particles with light. The characteristic slope, which which originates in this mechanism, has been detected in several hot Jupiters (Pont et al. 2008; Wakeford et al. 2013; Nikolov et al. 2018). The Rayleigh and/or Mie scattering mask the presence of features from the deeper atmospheres, which include the broadened alkali K and Na species, particularly in the visible, but also in the near-infrared (e.g., Wakeford et al. 2013; Wakeford & Sing 2015; Sing et al. 2016).

The Rossiter-McLaughlin (RM) effect (Holt 1893; Rossiter 1924; McLaughlin 1924) can be used as a tool to measure changes in the radius of an exoplanet with wavelength, the so-called chromatic Rossiter-McLaughlin (Dreizler et al. 2009). This alternative way of retrieving a transmission spectrum was explored for the first time by Snellen (2004), where the author measured, from ground-based observations, the increased amplitude of the RM effect around Sodium D-lines when compared with the average amplitude for the exoplanet HD 209458b. Later, Di Gloria et al. (2015) applied the same technique to HD 189733b using the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003) data and found the broadband signature of Rayleigh scattering, which has been confirmed by Oshagh et al. (2020).

In this paper, we present CaRM, a novel implementation of the chromatic Rossiter-McLaughlin effect (CRM). The CRM, although low resolution, has the advantage of capturing the broadband transmission spectrum of the exoplanet atmosphere. With high-resolution techniques, the “renormalization” of the photometric continuum wipes out any information it contains in broadband.

CaRM has already been used in Santos et al. (2020) for the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO, Pepe et al. 2010, 2013, 2021) data of HD 209458b. In this work, CaRM is used to analyze HARPS and ESPRESSO data of HD 189733b and WASP-127b, respectively. In Sect. 2 we describe the principles behind the RM effect and how it can be used to measure the planet radius at different wavelengths. We also present our implementation of the CRM effect in detail. In Sect. 3, we explore and assess the effect of considering different RM models and limb-darkening laws. Finally, we detail how we employed CaRM to retrieve the transmission spectrum of the exoplanets HD 189733b and WASP-127b (Sects. 4–6). For the former, we modeled the transmission spectrum to, in addition, address the significance of the broadband feature that is observed.

2. CaRM

2.1. The chromatic Rossiter-McLaughlin effect

The RM effect is an anomaly in the RV curve of a star observed during a transit or eclipse. When the planet or star crosses the stellar disk, it covers regions with different stellar rotational RVs. This causes an unbalance of the integrated RVs for each element of the surface of the star, and a shift in its measured RV. The amplitude (and shape) of the RM effect contains information about the radius ratio between planet and star, the impact parameter, the sky-projected spin-orbit angle, and spin velocity of the star (Triaud 2018):

$$\begin{array}{c}\hfill A_{\mathrm{RM}}\propto \frac{2}{3}D(\lambda )vsin{i}_{\star }\sqrt{1-b^2},\end{array}$$(1)

where v sin i_⋆ is the equatorial velocity of the star projected over the line of sight, D(λ) = (R_p(λ)/R_⋆)², and b the impact parameter which corresponds to the sky-projected distance between the center of the stellar disk and the position of the planetary disk at conjunction. Since the value of D can be a function of wavelength (λ), this equation shows that a measurement of the amplitude of the RM effect can be used to measure the planet radius variations as a function of the wavelength, that is to say to retrieve the transmission spectrum of the exoplanet.

CaRM, for chromatic Rossiter-McLaughlin, is a code developed to measure amplitude variations of the RM anomaly as function of wavelength. For this, the user has to provide the RVs as a function of the wavelength. Then, the data in each wavelength bin, defined by the user, are fitted with an RM model. From the fit, the planet’s radius is computed and, as a result, the broadband transmission spectrum of the exoplanet is returned.

2.2. Data input

The input data can be given in the form of standard CCF files, as delivered by their respective pipelines for ESPRESSO and HARPS, or text files (one for each bin). For ESPRESSO, the orders are duplicated thanks to the use of a pupil slicer delivering two simultaneous spectra (or two slices) for each spectral order (Pepe et al. 2021).

For the automated retrieval of the RVs for HARPS and ESPRESSO, the code takes as input the ranges of orders or slices, which define the bins, as well as the data paths and the DRS fits’ files extensions for the CCFs. For each bin, it performs a flux weighted (by the inverse variance) sum of the corresponding CCFs. The flux uncertainties are computed with the standard error for the weighted mean for each point³. The RV corresponds to the mean of the Gaussian fit of the resulting CCFs associated to each wavelength bin. For ESPRESSO, the error is estimated by taking the gradient of the CCF flux (in the RVs’ axis) and summing the inverse square of the former weighted by the CCF errors. For HARPS, there is no information about the uncertainty of the flux in each CCF point. In this case, we employed the Bouchy et al. (2001) method directly to each spectral order (which is read from the S2D files). Both methods are equivalent and are used to measure the contribution of the photon noise to the uncertainty in RV measurements.

For the text files, in the table format, the first line defines the wavelength range and the second provides the identification of each data column. The table itself is expected to contain the observation times (BJD), the RVs, and respective uncertainties (km s⁻¹). An example is provided together with the code⁴.

2.3. The models

The RVs of a star hosting planets are composed of the joint velocity of the system, the Keplerian motion, and, during the transit, the RM effect. In CaRM, the Keplerian component is modeled by a sinusoidal function (i.e., in this version only a circular orbit is considered) with a semi-amplitude K, an orbital period P, and a shift of the mid-transit time, in terms of orbital phase, Δϕ₀. A constant velocity, V_sys, is added to account for the systemic velocity. For the “RM” component, we employed two distinct models: ARoME (Boué et al. 2013) and Ohta et al. (2005) (implemented in PyAstronomy, Czesla et al. 2019).

The measurement of the Doppler shift of spectra can be made using, for example, the iodine cell technique (Butler et al. 1996) or the CCF technique that uses the cross-correlation between the observed spectrum and a mask that indicates the location and strength of the known stellar lines (Baranne et al. 1996). ARoME is a code developed to accommodate these observation techniques. The version we included in CaRM is an adaptation with Cython of the original library⁵. As inputs, the code uses the planet to star radius ratio (R_p/R_⋆), the mid-transit mean-anomaly that is derived from the phase assuming an argument of periastron of π/2, the semi-major axis of the orbit (a/R_⋆), the orbital inclination (i_p), the projected angle between the spin-axis and the orbital plane (λ), the projected rotational velocity, and the stellar macro-turbulence parameter ζ_t. This last parameter combined with the stellar rotation velocity is added to the width of the Gaussian fit to the stellar surface CCF behind the planet (β₀) plus the width that results from the Gaussian fit to the out-of-transit CCFs (σ₀).

The model available in the PyAstronomy package (Czesla et al. 2019) uses a geometrical approach to describe the RM anomaly based on the description of Ohta et al. (2005) for RVs derived from template matching. The projected velocity toward the observer is modeled as the sum of the projected proper motion and the rotational contribution of the surface. The model uses the equatorial stellar angular velocity (Ω), which we computed using the rotational velocity of the star (V_rot) with a fixed stellar radius, and the inclination of the stellar rotation axis (i_⋆).

Both RM models require one to provide a set of limb-darkening coefficients. For ARoME, it accepts both quadratic or nonlinear limb-darkening coefficients as input. In contrast, PyAstronomy only works with a linear law. The impact of using different limb-darkening laws is explored in Sect. 3. We use the limb-darkening toolkit LDTk⁶ (Parviainen & Aigrain 2015), which is a package that automates the computation of the limb-darkening profiles using specific intensity spectra generated from the PHOENIX spectral library (Husser et al. 2013). This code allows the fit of limb-darkening profiles in specific, user provided, wavelength bins. To retrieve the coefficients, LDTk uses, in addition to the wavelength interval, the effective temperature of the star T_eff, the logarithm of the stellar surface gravity log(g), and the stellar metallicity [Fe/H]. We adopted an approach similar to the photometric transmission spectroscopy, where the limb-darkening profile is calculated integrating over the full range of the wavelength bin. A word of caution is needed, though, since the RM effect is a consequence of the shape variation of the stellar lines in a given wavelength range. Wavelength trends on the strength of individual spectral lines would change the limb-darkening profiles, which can potentially mimic chromatic variations.

CaRM also offers the possibility to model the activity-induced RV signal using two approaches: adding a linear phase-dependent slope to the data and adding a GP model. The former takes into account long-term trends due to stellar rotation and the effect of stellar activity in the out-of-transit slope, while the latter can be used to try to model in-transit variability (e.g., occulted spots or granulation).

For the linear slope, if the transit window is provided, the code can identify the in-transit and out-of-transit measurements. This is fitted simultaneously with the RM curve.

For the GP analysis, CaRM calls the celeriteRealTerm kernel (Foreman-Mackey et al. 2017). This requires the definition of the log of amplitude and timescale, log(a) and log(τ), of the residual RV signal. The GP kernel takes the following functional form:

$$\begin{array}{c}\hfill k(\varphi )=a{e}^{-\varphi \tau },\end{array}$$(2)

where ϕ is the orbital phase.

2.4. Implementation

CaRM is developed in a blend of Python3 and Cython3 with common standard packages such as Numpy and Scipy. The RM models are chosen by the user employing an input file. This file consists of a set of instructions to fit both the white light (i.e., using the full passband of the instruments) and the chromatic data. There we define the stellar and orbital properties combined with the selection of the RM model. As input, we can select the preferred limb-darkening law and provide stellar parameters such as the radius, the effective temperature, log-surface gravity, and metallicity (Table 1 and the documentation available on CaRM’s page). Furthermore, in the same file, we provide the inputs and the initial estimates of the model parameters, or if each parameter must be fitted or not with a given prior. The range of the priors and the possibility of some parameters to be independently fit of each night can be enabled. The code assumes that if the prior is not declared, the corresponding parameters are fixed with the value provided by the user input values.

A call of the code starts to build the “results” file in which the parameters provided by the user are stored (see Fig. 1). The times of observation for each epoch are converted into orbital phases using the input values for the period and the mid-transit time. The passbands and stellar parameters are used to compute the limb-darkening coefficients, for the specified law, with LDTk. For the combined white light and chromatic analysis, the first data set that is run is for the white-light RVs. After this fitting, taking advantage of the higher S/N, it is possible to use the computed parameters to constrain the chromatic run (assuming they are wavelength independent). The best-fit model values can be obtained, in an automated way and by user definition, as the posterior maximum probability or median solutions. The prior definition, the chains (which are picked from a randomly drawn distribution from the priors), and log likelihood are provided to the Markov chain Monte Carlo (MCMC) implementation emcee (Foreman-Mackey et al. 2013). Each evaluation of the model is computationally expensive; however, with the parallelization potentialities, the sampling is accelerated, reducing the overall time of computation.

Fig. 1. Processes flow scheme of CaRM. The top gray area represents the user inputs and the bottom is for the last process execution in a run. The RVs can be directly read from text files or, alternatively, computed from the CCFs files for HARPS and ESPRESSO.

The MCMC runs in two distinct phases. The burn-in is used to sample the posterior distributions and narrow the range of values of each parameter. After this, the final position of the walkers is saved and used as the initial position for the production phase. The fit of the chromatic data has a similar process. The major distinction is the use of the updated parameters, fixing the constant parameters from the white-light RV fitting. After the code runs for all the passbands, the information is stored for analysis.

3. The limb-darkening impact

Limb darkening is a radial effect produced by the observation of atmospheric stellar layers at different temperatures. From the perspective of the observer, the effective optical depth decreases with increasing radius due to lower gas density and a shorter line-of-sight distance through the star.

Espinoza & Jordán (2016) studied the impact of the different limb-darkening laws (Appendix A) in transit light curves. They conclude that, comparing the limb darkening recovered for a number of Kepler targets, there are biases introduced that are particularly relevant for noncentral transits. This can affect the determination of the stellar inclination, semi-major axis (a/R_⋆), and the planet radius (R_p/R_⋆). The proposed solution to partly solve these biases is to let the limb-darkening coefficients free with high polynomial order limb-darkening laws.

To verify if the choice of the limb-darkening law can produce false transmission spectroscopy features, we assessed the impact in RVs, assuming that the intrinsic law that models the stellar limb-darkening profile of the simulated system is nonlinear. For this, we fit the limb-darkening laws (whose profiles are compared in Fig. A.1) with the LDTk package for a simulated typical system with a hot Jupiter (parameters in Table B.1). The impact of the laws in the RM was computed with the average temperature range (i.e., from minimum to maximum) of the FGK stars. We then compared the Rossiter-McLaughlin anomalies, computed with ARoME, to quantify the differences in RVs between linear and quadratic laws with the underlying nonlinear law (Fig. 2).

Fig. 2. Stellar temperature impact, in RVs, in the limb-darkening laws. Top row: RM curve (computed with ARoME), with increasing spin-orbit misalignment, for three temperatures in the temperature range of the FGK stars. Middle and bottom row: residuals’ variation (in tenths of cm s−1) between a linear and a quadratic law. The RV scale is constant to emphasize the variation of the residuals’ amplitude with stellar temperature.

The RVs’ comparisons show the same global behavior, with the quadratic law presenting the lowest residuals overall. In addition, there is an increasing trend in the amplitude of the linear residuals, whilst there is a slight decrease for the quadratic law. ESPRESSO, for example, has achieved a precision of < 30 cm s⁻¹ (Suárez Mascareño et al. 2020) on multiple observations, despite the precision being able to be lower for single measurements (depending, e.g., on the magnitude or spectral type). For the simulated hot-Jupiter system, the peak difference (< 25 cm s⁻¹) is below the precision of most state-of-the-art spectrographs. As such, we would not be able to measure differences that originate from the choice of the limb-darkening law and which could produce some chromatic variability.

4. Observations and data reduction

To apply CaRM to the observations of HD 189733b and WASP-127b, we computed the RVs using the CCF method. The CCFs from each observation were computed by the HARPS and ESPRESSO data reduction software packages (DRS version 3.4 and 2.2.8, respectively), and they are provided per order and per slice (ESPRESSO).

4.1. HD 189733b

The HD 189733b data consists of three sets of observations taken on 8 September 2006, as well as 20 July and 29 August 2007 with HARPS (Fig. 3). This corresponds to the same data analyzed in Di Gloria et al. (2015). HARPS is a fiber-fed spectrograph installed on the 3.6 m telescope Cassegrain in La Silla. It covers a spectral range between 378 nm and 691 nm, distributed over 71 echelle orders. The observations were acquired under the programs 072.C-0488(E), 079.C-0828(A) (PI: Mayor), and 079.C-0127(A) (PI: Lecavelier des Étangs). The exposure time was set to 900 s for the first night and 300 s for the following one. The detector was set for the fast readout and low gain mode, which offers a lower overhead between exposures (23 s when compared with 180 s for the high gain mode) at the expense of increased readout noise.

Fig. 3. RVs of HD 189733b retrieved from the CCF header for the white light. The error bars are not visible since they are smaller than the dimension of the markers. In the same column, the associated S/N near 580 nm and the airmass is a function of the number of days since the first epoch. The average signal-to-noise is 178, 102, and 116, respectively.

The CCFs were produced cross-correlating with a K0 spectral template binary mask, closely matching the spectral type of the host star, HD 189733 (K2 V from Gray et al. 2003). This choice is is different from that that adopted by Di Gloria et al. (2015), which used the reduction available in the ESO archive. That reduction used a different mask, G2, for the last two nights.

4.2. WASP-127b

The observations of WASP-127b were acquired on the 24 February and 17 March 2020 under the GTO program ID 1102.C-0744 (PI: F. Pepe). In total, there are 148 RV measurements with 300 s of integration. The spectrograph was set to singleHR21 mode with a spectroscopic HR (which corresponds to the highest prescribed RV precision of the instrument), one UT, slow reading, and 2 × 1 binning. The RVs were extracted using a G8 mask as provided by the ESPRESSO DRS. This same data were analyzed by Allart et al. (2020) to study the atmosphere of WASP-127b in high resolution.

The inspection of the data shows hints of a low amplitude RM effect. The after transit RVs of the first night reveal an increased dispersion and average error when compared with the rest of the data. This is accompanied by a decrease in the signal-to-noise (S/N) and it is caused by an increase in the airmass (see Fig. 4).

Fig. 4. RVs retrieved from the CCF header for the white light and respective error bars for each night for WASP-127b. In the same column, the associated S/N near 580 nm and the airmass is a function of the number of days since the first epoch.

5. The transmission spectrum of HD 189733b

5.1. The white-light RM fit

As presented in Santos et al. (2020) for HD 209458b, the first step to retrieve the transmission spectrum is to fit the white-light data with the RM model. This is performed to fix wavelength-independent parameters using the higher S/N data. We used the ARoME model and adopted, as a reference, most parameters from Triaud et al. (2009), which is the same source from the non model-specific parameters used by Di Gloria et al. (2015) (see Table E.1). We also used the updated stellar parameters, with the exception of the stellar radius, derived in Sousa et al. (2008) and listed in SWEET-CAT⁷ (Santos et al. 2013) which compiles sets of stellar atmospheric parameters derived, when possible, with a homogeneous methodology for planet host stars.

The initial value for the mean systematic velocity of the star V_sys (−1.9 km s⁻¹) was estimated by observation of the data, and it is only a rough estimate. In Fig. 3 we see that, taking the average position of the mid-transit as a reference, this value is quite variable. This variation is expected between observations as a consequence of the nightly zero point, which can have an origin in the different calibrations performed at the beginning of each night for different observations, or as a result of stellar variability⁸. We chose to give a comprehensive uniform prior to account for these variations.

The spin-orbit angle prior was fixed to the Triaud et al. (2009) estimate using the classic RM approach, $-0.85{}_{-0.32}^{+0.28}°$ (using also the same dataset). Other estimates from literature, using different data sets or techniques, such as −1.4 ± 1.1° from Winn et al. (2006), −0.35 ± 0.3° and −0.67 ± 0.3° (independently estimated from two transits) from Collier Cameron et al. (2010) employing the line profile tomography technique, or −0.4 ± 0.2° from Cegla et al. (2016) who used the reloaded RM technique, are compatible with the one adopted and they suggest an aligned orbit. The projected rotational velocity of the star (v sin i_⋆) was first estimated by Bouchy et al. (2005) using the ELODIE spectrograph, from the following line broadening method: 3.5 ± 1.0 km s⁻¹. Cegla et al. (2016) recovered the rotational profile of HD 189733 and excluded the rigid body rotation with a high degree of confidence. They discuss that considering the rigid body rotation in these cases biases the velocities’ (v sin i_⋆) computation toward a lower value, and that it can produce wave-like residuals such as the ones observed by Triaud et al. (2009). We chose to fit the velocity with a Gaussian prior centered on the value from Triaud et al. (2009) after the correction (minimizing the RM residuals), as the RM models used by CaRM do not include stellar differential rotation.

The value adopted for σ₀, that is, the width of the best Gaussian fit, results from the averaged measurement of the out-of-transit CCF widths obtained from a Gaussian fit. We estimate a similar value as Di Gloria et al. (2015), but it is important to note that there is a variability of 0.03 km s⁻¹ between the average out-of-transit of the night with the highest and lowest value. This should not greatly impact the fit; nevertheless, the RM amplitude does not have a strong dependence on this parameter and the variation represents a small percentage of the value we used (Santos et al. 2020). The width of the nonrotating star β₀ depends on measurements of the CCF width of a sample of stars and incorporates the instrumental profile. We used the calibration done for HARPS data, similar to the one presented in Santos et al. (2002).

The white light was set to run with 50 chains with an initial burn-in phase with 1500 steps and 3000 for the production, using as an initial guess the parameters from Table E.1 and the priors from Table 2. The individual fits of the three observations and residuals after subtraction of the best fit model are presented in Fig. D.1. By inspection, there is some structure in the residuals of each night that turns out to be more evident when the residual RVs are combined (Fig. 5). The HARPS data of HD 189733b have been extensively analyzed in the literature (e.g., Triaud et al. 2009; Di Gloria et al. 2015; Cegla et al. 2016), and similar residual structures have been found. Some hypotheses that were proposed to explain this are differential rotation, spot crossing, or convective blueshift (CB).

Fig. 5. HD 189733 RM fit. Top: best fit model (solid red) of the combined observations of HD 189733b with HARPS. The gray area corresponds to a random sample of the posterior about the highest likelihood. The data points present two error bars, the green corresponding to the original RVs’ uncertainties computation from the CCF method and the black after quadratically adding the white-noise level fitted to the data. Bottom: residuals after subtraction of the data from the best-fit. The lighter yellow represents the ingress and egress and the darker the region where the planet is fully in front of the stellar disk. At the top of the residual’s plot, σres represents the average value of the RVs’ residual dispersion, σW is the estimate of the jitter velocity, and σ is the average uncertainty in the residual RVs.

The RM models we include assume a rigid body rotation for the star, which is equivalent to fix the stellar surface rotational velocity. It is known that stars possess differential rotation, in particular the Sun, which is a function of the latitude. Cegla et al. (2016) addressed this issue with the reloaded RM, which does not assume any specific line profile (and then velocity profile) to the occulted areas by the planet. They find that with differential rotation in the model, the wave-like features vanish. This star has also been reported to be active (Bouchy et al. 2005) and multiple events of spot crossing have been detected in photometric transits (Pont et al. 2013). In the context of the RM, similar spot-crossing events (phase aligned) seem unlikely with observations that are taken far apart in time. Furthermore, the effect would have to be present during the complete duration of the transit, as the residuals have the same time span. Finally, CB (e.g., Shporer & Brown 2011), which is caused by the disk-integrated differences in brightness from the hotter upward-flowing material and the downward-flowing one, is a hypothesis that also needs to be considered. In ARoME, the effect of the CB is said to be negligible, as it only adds a constant offset to the RVs and does not modify the RM signal. Despite this, Albrecht et al. (2012) shows that although symmetric, it contributes to distortions in amplitude of the RM signal as a result of distortions of the line profiles from the center to the stellar limb (Dravins et al. 2021).

The median of the posterior distribution and uncertainties of the fitted parameters are in Table 3. The K parameter, with a contribution from the Keplerian amplitude and variability, differs the most for the third night of observations. During this night, the beginning of the transit was not observed and roughly half of the data points correspond to out-of-transit measurements (Fig. D.1). It can be seen in the residuals, in Fig. 5, that structures with an amplitude up to ∼2.5 m s⁻¹ do not correspond to in-transit events. The higher level of activity also affects the jitter amplitude, which is the highest from the three nights. The median value of β₀ is in line with the value reported by Di Gloria et al. (2015). For the projected rotational velocity, the median value of the posterior seems to tend to a higher value, standing between the value reported by Triaud et al. (2009) and Bouchy et al. (2005).

5.2. Transmission spectrum

The hot Jupiter HD 189733b has been extensively studied since its discovery (Bouchy et al. 2005). The bluer side of the optical transmission spectrum of HD 189733b was reported to present a high radius-wavelength slope (Pont et al. 2008; Sing et al. 2008; Di Gloria et al. 2015) which is not often seen in the transmission spectra of exoplanets of this type, a possible signature of Rayleigh scattering. Furthermore, there is no evidence for the detection of wings of the sodium doublet, while their cores are detected (e.g., Wyttenbach et al. 2015; Seidel et al. 2020a). With the low amplitude water feature, the presence of an extended atmosphere dominated by submicron haze particles seems to dominate the spectra from visible to infrared wavelengths (Sing et al. 2008, 2009).

The chromatic Rossiter-McLaughlin method was applied twice to this target in the literature. The work of Di Gloria et al. (2015) found, with the same HARPS data analyzed here, the same strong wavelength radius slope consistent with the observations of Pont et al. (2008). The value of the slope from HARPS data using the RVs obtained through the CCF method is larger when compared with the HST observations. Later, the study of Oshagh et al. (2020) shows that it is possible to obtain a slope closer to the one present on the transmission spectrum derived from HST data, using RVs obtained via template matching (Butler et al. 1996; Anglada-Escudé & Butler 2012). In the same paper, the hypothesis is raised that the standard templates from the HARPS DRS can contain blended lines that contribute to changes in the retrieved RVs, contaminating the transmission spectrum. Alternatively, Oshagh et al. (2014) and Boldt et al. (2020) show that stellar activity, such as spots of faculae, can mimic broadband features such as the slopes from Rayleigh scattering mechanisms. In particular, the activity seems to affect the RVs’ observations producing variable slopes. Nevertheless, the probability of the activity giving an origin to a coherent (in time) strong feature is low, in particular with data from several transits.

To retrieve the transmission spectrum of HD 189733b, we used the standard RVs, computed with the CCF method with the default DRS masks (K0). As described in Sect. 2, we can use two methods that help to mitigate the impact of stellar activity in the retrieved transmission spectrum. In Santos et al. (2020), the slope K was set as a free parameter (for each night). Di Gloria et al. (2015) and Oshagh et al. (2020) fitted the out-of-transit points of each night separately and removed a linear trend before fitting the RM profile. In this context, K must be interpreted not purely as the Keplerian amplitude, but a combination of it with the time-dependent linear slope.

The in-transit observation can also be affected by activity. Potential occultations of stellar spots or faculae are observed in the residuals, with a profile similar to the RM anomaly. This is seen in Fig. 5 and discussed in Sect. 5.1. Di Gloria et al. (2015) tried to mitigate the effect of these residuals by subtracting the residuals that result from the white-light fit from each chromatic bin. With CaRM, we chose not to perform this as it does not account for the chromatic variations of the phenomenon that gives origin to it. Furthermore, the contribution of the convective blueshift is differential along the stellar disk, causing center-to-limb variations (Cegla et al. 2016; Reiners et al. 2016). The line profile is deeper and thinner at the disk center and broader and shallower at the limb. The subtraction of the white-light line profiles from the chromatics creates residuals of variable width and depth (Dravins et al. 2021).

We performed three runs with CaRM to assess the impact on the transmission spectrum of different setups:

• Gaussian processes with a linear slope ξ (run1);

• fit the Keplerian semi-amplitude K (run2);

• fit the linear slope (run3).

The parameters in these runs were fitted for each observation independently, and the transmission spectra are compared in Fig. 6. We preferred run2 (chromatic fit priors in Table 4) since it minimizes the number of parameters in the model and, at the same time, can account for the chromatic variations during transit. The first runs present problems when fitting the out-of-transit slope in particular for the third night, as there is no available data before ingress. Despite this, the comparison between retrievals show a good agreement at the 1σ level (Fig. 6 and Table C.1 for the preferred solution).

Fig. 6. HD 189733b transmission spectra. Top left: comparison of the transmission spectra retrieved with CaRM, with ARoME, for different setups (see text). The bins centered near 545 nm are aligned for better comparison. After, a comparison between our result and Pont et al. (2013) observed with Hubble (STIS), Di Gloria et al. (2015) with HARPS, and Oshagh et al. (2020) which combines HARPS and CARMENES data. Shifts where applied to compare with the various literature sources.

To check the validity of fixing parameters from the white light to the chromatic runs, in particular the ones associated with the ARoME model, we decided to measure the average out-of-transit CCF width as a function of the wavelength. We found a positive slope from blue to red wavelengths, with an average value which corresponds to the white-light σ₀ (Table E.1). To check if this variation affected the transmission spectrum, we provided as input to the code the measured σ₀ instead of fixing it to the white-light value. Since β₀ is correlated with σ₀, these variations also affect its value. We then quadratically added the difference between the bin CCF width and the white-light value. There were no significant changes, within the error bars, on the final transmission spectrum. In addition to this, we also checked if the transmission spectrum could be the result of wavelength-dependent BIS variations. For this, we subtracted the white-light BIS from the value for each bin and compared it with the general shape of the transmission spectrum. The BIS value corresponds to a few percent of the CCF widths; however, we did not find a correlation with the radius variations.

5.2.1. Transmission spectrum modeling with PLATON

To analyze the preferred transmission spectrum, we used PLATON (Zhang et al. 2019, 2020). This tool is a Python package that computes transmission spectra and allows one to fit the transmission models to retrieve atmospheric characteristics using observational data (Zhang et al. 2020).

As a template to run the retrieval, we used the publicly available code example⁹. We considered an isothermal atmospheric profile, which constitutes an adequate approximation (e.g., Lecavelier Des Etangs et al. 2008; Heng & Kitzmann 2017) for the atmospheric layers probed by transmission spectroscopy in the visible, with equilibrium chemistry. The stellar and planetary radius were set with a normal distribution, with the standard deviations and mean set to the values provided in Table E.1, taking a broader radius prior (3σ) to account for shifts in the average radius of the transmission spectrum. We used the value reported by Triaud et al. (2009) for the planet mass (𝒢(1.138, 0.025)) and the posterior values reported in Zhang et al. (2020) for the other parameters of HD 189733b. We used the priors 𝒢(1089, 120) for the temperature of the atmosphere, 𝒢(1.08, 0.23) for the logarithm of the metallicity, and 𝒢(0.66, 0.09) for the C/O ratio. We increased the range of the uniform prior for the logarithm of the scattering factor and maintained the error multiple, 𝒰(−2, 6) and 𝒰(0.5, 5), respectively. The stellar temperature is fixed in the run, with the value also reported in Table E.1.

We tested the retrieval comparing the evidence of a flat model (no radius variation with wavelength) with one containing Rayleigh scattering, fixing the scattering slope parameter to 4 (Fig. 7), and a combination of Rayleigh scattering and Na. The Bayesian evidence was measured for each scenario using the built-in nested sampler to fit the transmission models. The data suggest the model with Rayleigh scattering is preferred over the flat model (with a logarithmic Bayesian evidence of ln(z) = 32.4 ± 0.1 versus ln(z) = 31.7 ± 0.1). The model considering sodium does not present any major difference with the one considering only Rayleigh, and the retrieval outputs an intermediate logarithmic evidence (ln(z) = 32.2 ± 0.1). The transmission spectrum we obtained presents a good agreement when compared with Di Gloria et al. (2015). Nevertheless, the evidence is not sufficient (with a Bayes factor of Δln(z) = 0.7 ± 0.1), due to the low resolution, to confirm the presence of Rayleigh scattering following the Kass & Raftery (1995) criteria (which requires Δln(z) > 3).

Fig. 7. Atmospheric models’ fit of the preferred transmission spectrum with a flat profile (left) and adding Rayleigh scattering (right).

6. The transmission spectrum of WASP-127b

6.1. The white-light fit of WASP-127b

Similarly to the analysis of HD 189733b, we started fitting the combined observation in the full wavelength range of ESPRESSO with the ARoME RM model. The V_sys encompasses a comprehensive interval around the observed average. The R_p/R_⋆ prior is a Gaussian distribution centered in the value provided by Seidel et al. (2020b) and the standard deviation corresponds to the 1σ uncertainty. The spin-orbit misalignment angle, in contrast with the previous study, has no published value using the classic RM analysis. It has been discussed in Cegla et al. (2016) that the reloaded RM analysis tends to give distinct values for the angle, which is attributed to center-to-limb convective variations. For a slow rotating star, such as WASP-127, we found that the residuals are dominated by these variations (producing amplitudes ∼cm s⁻¹ to ∼1 m s⁻¹) and uncertainties in the obliquity from 10° ∼20°. It was chosen, as a result, to let the spin-orbit angle vary using a normal distribution prior with the mean from Allart et al. (2020) and a standard deviation of 20° to account for the differences in the methods (priors summarized in Table 5). The MCMC run parameters were set in the same way as for HD 189733b for all bins. The combined data with the best fit model, as well as the residuals after subtraction from one to the other, are provided in Fig. 8.

Fig. 8. Combined data fit and residuals. The elements in the plots have the same meaning as in Fig. 5.

The posterior distribution of the fitted parameters for the white light is found in Table 6. The median systemic velocities for the observed nights present very similar values. The host star is, in fact, very similar to the Sun, but, with 9.7 ± 1.0 Gyr and a larger radius (Allart et al. 2020), it is at the end of its main-sequence phase and evolving toward the branch of the subgiants. Despite this, the star seems to present a low variability from these observations taken 21 days apart. The spin-orbit angle value, $-150.6{}_{-8.7}^{+8.1}°$, has a difference of approximately 20° when compared with the value from Allart et al. (2020) (Table E.1). This is at the top of the estimated interval by Cegla et al. (2016) for the difference between the values obtained from the classical RM when compared to the reloaded RM. Despite this, the angle solutions are compatible at the 68% level.

6.2. The transmission spectrum of WASP-127b

WASP-127b is an exoplanet that orbits a bright star in a short orbital period (Table E.1). In terms of proximity to the star and radius, it can be classified as a hot Jupiter; however, in terms of mass, it resembles more a member of the Neptune population. With 0.18 Jupiter masses (Lam et al. 2017) and just over a Jupiter radius, it is one of the lowest density planets known. Its high insolation and low mass places it in the so-called Neptune desert (Lecavelier Des Etangs 2006; Mazeh et al. 2016). Lam et al. (2017) estimate a large-scale height for this exoplanet atmosphere (∼2350 km), which makes it a good target for transmission spectroscopy. Despite this, the rotational velocity of the star is low. The amplitude of the Rossiter-Mclaughlin effect is just ∼3 m s⁻¹, which is an order of magnitude lower than the amplitude of the RM detected on HD 189733b.

The first complete study of WASP-127b found a cloud-free atmosphere for this puffy planet with a strong Rayleigh-scattering signature and features compatible with TiO and VO (Palle et al. 2017). The tentative signature of sodium is also claimed and was later confirmed at a 5σ level by Chen et al. (2018). Allart et al. (2020) also found signatures of sodium that are compatible both with Chen et al. (2018) and Spake et al. (2021), but the water feature amplitude analysis at 1.3 μm seems to indicate, at least partially, the presence of clouds. This is compatible with the upper limits established by Seidel et al. (2020b) using simultaneous photometry, but with the less precise HARPS data. The two observations of WASP-127b with ESPRESSO have both in-transit and out-of-transit data, which allows for a better determination of the activity-induced slope. We chose then to fix the RV semi-amplitude and fit V_sys, R_p/R_⋆, ξ, and log(σ_W), with the priors from Table 7.

We decided not to use GPs with this data set, as there are no clear structures seen in the residuals (Fig. D.2). The analysis of the transmission spectrum obtained with CaRM (Fig. 9, Table C.2) does not reveal broadband features, in contrast to what has been found for the previous target. At low resolution, it is compatible with the conclusions of Palle et al. (2017), Chen et al. (2018), and Spake et al. (2021). In fact, the transmission spectrum presented in Palle et al. (2017) has a stronger slope in bluer regions that we cannot probe because it is outside the wavelength range of ESPRESSO. Given the flatness of the retrieved spectrum and large error bars, we did not try to model it with PLATON.

Fig. 9. WASP-127b transmission spectrum retrieved with CaRM, with ARoME, compared with other sources in the literature (Palle et al. 2017 with NOT/ALFOSC; Chen et al. 2018 with NOT/ALFOSC and GTC/OSIRIS; and Spake et al. 2021 with HST/STIS and TESS). The gray bars represent the wavelength intervals where the CCFs were masked to mitigate atmospheric contamination. The horizontal bars, in our retrieval, correspond to the binning in wavelength.

7. Conclusions

CaRM is a semi-automatic code that uses the chromatic Rossiter-McLaughlin effect to retrieve the transmission spectra of exoplanets. In this paper, we described the code and studied the variables to be considered to use this method.

We used the code on HD 189733b with HARPS data and found similar results when compared with previous authors using the same method. The white-light RM fit showed that the model is not fully capable of reproducing the observations. The residuals obtained after the subtraction of the best-fit model show an RM-like pattern that has been previously attributed to spot-crossing events, differential stellar rotation, or CB (Triaud et al. 2009; Di Gloria et al. 2015; Cegla et al. 2016; Dravins et al. 2021).

With PLATON, the planetary transmission spectrum was analyzed, and we found that a model containing Rayleigh scattering is preferred over a flat one. We do not have enough evidence with this method to confirm the presence of the haze. Yet, Sing et al. (2016) were able to confirm it using Hubble and Spitzer Space Telescopes observations.

For WASP-127b, the white-light fit shows a rather stable system in terms of the systemic velocity. The spin-orbit misalignment found with the classical RM analysis shows a difference of ∼20° with the reloaded RM solution. The transmission spectrum, at the resolution CaRM is able to retrieve, is rather unremarkable and with no identifiable features. Despite this, when compared with other literature sources, it is more consistent with a flatter profile (Spake et al. 2021) than with a strong Rayleigh scattering in the bluest part of the spectrum (Palle et al. 2017; Chen et al. 2018).

In this work, we show that CaRM is capable of reproducing both past analysis (with HD 189733b) as well as performing new ones with ESPRESSO. However, it is a less sensitive method when compared to broadband transmission spectroscopy from space. The precision on the transmission spectrum depends, in a first degree, on the ratio between the Rossiter-McLaughlin amplitude and the uncertainties on the RV measurements. The uncertainties can be improved by incorporating multiple transits, and better RM models are fundamental to disentangle between the planet and stellar signatures. Taking this in consideration, CaRM proves to be a valuable tool to retrieve and study the transmission spectra of exoplanets using state-of-the-art spectrographs such as ESPRESSO and future HIRES at ELT.

Acknowledgments

The authors acknowledge the ESPRESSO project team for its effort and dedication in building the ESPRESSO instrument. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FIS-AST/32113/2017 and POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 and POCI-01-0145-FEDER-028953; PTDC/FIS-AST/28987/2017 and POCI-01-0145-FEDER-028987. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through Fundação para a Ciência e Tecnologia (FCT). The INAF authors acknowledge financial support of the Italian Ministry of Education, University, and Research with PRIN 201278X4FL and the “Progetti Premiali” funding scheme. N.J.N. acknowledges, additionally, support from FCT through the project CERN/FIS-PAR/0037/2019. J.H.C.M. is supported in the form of a work contract funded by Fundação para a Ciência e Tecnologia (FCT) with the reference DL 57/2016/CP1364/CT0007; and also supported from FCT through national funds and by FEDER-Fundo Europeu de Desenvolvimento Regional through COMPETE2020-Programa Operacional Competitividade e Internacionalização for these grants UIDB/04434/2020 and UIDP/04434/2020, PTDC/FIS-AST/32113/2017 and POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 and POCI-01-0145-FEDER-028953, PTDC/FIS-AST/29942/2017. R.A. is a Trottier Postdoctoral Fellow and acknowledges support from the Trottier Family Foundation. This work was supported in part through a grant from FRQNT. This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. A.S.M., J.I.G.H., C.A.P. and R.R. acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00, and from the Government of the Canary Islands project ProID2020010129. J.I.G.H. also acknowledges financial support from the Spanish MICINN under 2013 Ramón y Cajal program RYC-2013-14875. This 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). This project has received 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, project SCORE, grant agreement No. 851555). V.A. acknowledges the support from FCT through Investigador FCT contract nr. IF/00650/2015/CP1273/CT0001.

References

Appendix A: The limb-darkening laws

There are several parametrizations of limb darkening which can be more adequate depending on the physical properties of the star. Some of those parametrizations, which can be used in CaRM, are the linear (eq. A.1), the quadratic (eq. A.2), and the nonlinear (eq. A.3):

$$\begin{array}{cc}& \frac{I(\mu )}{I(\mu =1)}=1-u_1(1-\mu );\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& \frac{I(\mu )}{I(\mu =1)}=1-u_1(1-\mu )-u_2{(1-\mu )}^2;\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& \frac{I(\mu )}{I(\mu =1)}=1-{\displaystyle \sum _{n=1}^4{u}_n(1-{\mu }^{n/2}).}\hfill \end{array}$$(A.3)

The normalized intensities show the largest difference for μ = 0 for the three temperatures, Fig. A.1, with μ representing the cosine of the angle between the observer line of view and the normal to the stellar surface. We used the cumulative function of the difference between models to show, not only the localized differences expressed by larger slopes, but also the weight they have on the total difference.

Fig. A.1. Stellar temperature impact in the limb-darkening laws. Top row: Comparison between the normalized intensity for a linear and quadratic, and nonlinear limb-darkening laws (top). Bottom row: Cumulative percentage difference for the linear and quadratic laws and the nonlinear.

For the mock system we used, the cumulative difference between the quadratic and nonlinear laws is no larger than 4%, and it increases with increasing stellar temperature. This is similar to the trend that is observed in RVs and shows that a correct modeling of the limb-darkening law is important for future instruments with higher RV precision.

Appendix B: Simulated system parameters

Appendix C: Transmission spectra tables

Appendix D: Individual nights fits

Fig. D.1. Data fit and residuals. The elements in the plots have the same meaning as in Fig. 5.

Fig. D.2. Data fit and residuals. The elements in the plots have the same meaning as in Fig. 5.

Appendix E: Reference values table

Appendix F: Corner plot of HD 189733b

Fig. F.1. Corner plot of the joint white-light fit for the three HARPS observations of HD 189733b. The values and error bars represent the median and 1σ uncertainty for the posterior distributions. The number under the identification of each variable, when present, represents an individual fit of that particular parameter.

Appendix G: Corner plot of WASP-127b

Fig. G.1. Corner plot of the joint white-light fit for the two ESPRESSO observations of WASP-127b.

All Tables

All Figures

	Fig. 1. Processes flow scheme of CaRM. The top gray area represents the user inputs and the bottom is for the last process execution in a run. The RVs can be directly read from text files or, alternatively, computed from the CCFs files for HARPS and ESPRESSO.
In the text

Fig. 2. Stellar temperature impact, in RVs, in the limb-darkening laws. Top row: RM curve (computed with ARoME), with increasing spin-orbit misalignment, for three temperatures in the temperature range of the FGK stars. Middle and bottom row: residuals’ variation (in tenths of cm s−1) between a linear and a quadratic law. The RV scale is constant to emphasize the variation of the residuals’ amplitude with stellar temperature.

In the text

	Fig. 3. RVs of HD 189733b retrieved from the CCF header for the white light. The error bars are not visible since they are smaller than the dimension of the markers. In the same column, the associated S/N near 580 nm and the airmass is a function of the number of days since the first epoch. The average signal-to-noise is 178, 102, and 116, respectively.
In the text

	Fig. 4. RVs retrieved from the CCF header for the white light and respective error bars for each night for WASP-127b. In the same column, the associated S/N near 580 nm and the airmass is a function of the number of days since the first epoch.
In the text

Fig. 5. HD 189733 RM fit. Top: best fit model (solid red) of the combined observations of HD 189733b with HARPS. The gray area corresponds to a random sample of the posterior about the highest likelihood. The data points present two error bars, the green corresponding to the original RVs’ uncertainties computation from the CCF method and the black after quadratically adding the white-noise level fitted to the data. Bottom: residuals after subtraction of the data from the best-fit. The lighter yellow represents the ingress and egress and the darker the region where the planet is fully in front of the stellar disk. At the top of the residual’s plot, σres represents the average value of the RVs’ residual dispersion, σW is the estimate of the jitter velocity, and σ is the average uncertainty in the residual RVs.

In the text

Fig. 6. HD 189733b transmission spectra. Top left: comparison of the transmission spectra retrieved with CaRM, with ARoME, for different setups (see text). The bins centered near 545 nm are aligned for better comparison. After, a comparison between our result and Pont et al. (2013) observed with Hubble (STIS), Di Gloria et al. (2015) with HARPS, and Oshagh et al. (2020) which combines HARPS and CARMENES data. Shifts where applied to compare with the various literature sources.

In the text

	Fig. 7. Atmospheric models’ fit of the preferred transmission spectrum with a flat profile (left) and adding Rayleigh scattering (right).
In the text

	Fig. 8. Combined data fit and residuals. The elements in the plots have the same meaning as in Fig. 5.
In the text

Fig. 9. WASP-127b transmission spectrum retrieved with CaRM, with ARoME, compared with other sources in the literature (Palle et al. 2017 with NOT/ALFOSC; Chen et al. 2018 with NOT/ALFOSC and GTC/OSIRIS; and Spake et al. 2021 with HST/STIS and TESS). The gray bars represent the wavelength intervals where the CCFs were masked to mitigate atmospheric contamination. The horizontal bars, in our retrieval, correspond to the binning in wavelength.

In the text

	Fig. A.1. Stellar temperature impact in the limb-darkening laws. Top row: Comparison between the normalized intensity for a linear and quadratic, and nonlinear limb-darkening laws (top). Bottom row: Cumulative percentage difference for the linear and quadratic laws and the nonlinear.
In the text

	Fig. D.1. Data fit and residuals. The elements in the plots have the same meaning as in Fig. 5.
In the text

	Fig. D.2. Data fit and residuals. The elements in the plots have the same meaning as in Fig. 5.
In the text

	Fig. F.1. Corner plot of the joint white-light fit for the three HARPS observations of HD 189733b. The values and error bars represent the median and 1σ uncertainty for the posterior distributions. The number under the identification of each variable, when present, represents an individual fit of that particular parameter.
In the text

	Fig. G.1. Corner plot of the joint white-light fit for the two ESPRESSO observations of WASP-127b.
In the text