An ESPRESSO view of HD 189733 system. ⋆ Broadband transmission spectrum, differential rotation, and system architecture

Context. The development of state-of-the-art spectrographs has ushered in a new era in the detection and characterization of exoplanetary systems. The astrophysical community now has the ability to gain detailed insights into the composition of atmospheres of planets outside our solar system. In light of these advancements, several new methods have been developed to probe exoplanetary atmospheres using both broad-band and narrow-band techniques. Aims. Our objective is to utilize the high-resolution and precision capabilities of the ESPRESSO instrument to detect and measure the broad-band transmission spectrum of HD 189733b’s atmosphere. Additionally, we aim to employ an improved Rossiter-McLaughlin model to derive properties related to the velocity fields of the stellar surface and to constrain the orbital architecture. Methods


Introduction
The detection and characterization of exoplanets depend on several techniques that enable us to uncover the subtle signatures left by planets in the signals of their host star.One of the primary methods for detecting exoplanets is radial velocities (RVs), which is used to measure the star's reflex motion around the system's barycenter and thereby detect planetary companions.Precisely measuring RVs is often challenging, and requires stateof-art, high-resolution spectrographs (e.g.Mayor et al. 2003;Pepe et al. 2021) with long-term stability.However, RVs are not limited to detection alone; they are also valuable for the atmospheric characterization of exoplanetary systems, particularly during transits or occultations.
A transiting exoplanet covers stellar regions with varying brightness, spectral content and velocities.This is a result of ⋆ Based on Guaranteed Time Observations (GTO) collected at the European Southern Observatory under ESO programme 1102.C-0744 by the ESPRESSO Consortium.
the limb-darkening effect, the presence of stellar activity such as spots and planes, and the intrinsic rotation of the star.In photometric observations, this results in a decrease in the observed brightness of the star.When measuring RVs, it manifests as an anomaly caused by the planet blocking areas of the stellar surface with different projected velocities.This RV variation, first observed in binary stars, is known as the Rossiter-McLaughlin (RM) effect (Holt 1893;Rossiter 1924;McLaughlin 1924).The RM effect, for transiting exoplanets, was first measured for the HD 209458 system (Queloz et al. 2000) and, more recently, successful attempts have been made to measure the RM anomaly within the solar system, such as the Earth (e.g.Molaro et al. 2015;Yan et al. 2015) and Venus (e.g.Molaro et al. 2012).
The RM effect can be used as a tool to complement the orbital geometry derived from the Keplerian motion outside transit, providing a direct way to measure the projected spin-orbit angle.Some notable applications of these orbital measurements include statistical studies about orbital tilts (e.g.Fabrycky & Winn 2009;Albrecht et al. 2022;Mancini et al. 2022), and simultaneous measurements of the spin-orbit angles in multi-planetary systems (Bourrier et al. 2021).Additionally, since the planet covers different portions of the star along its track during a transit, the RM curve can also be explored to measure, e.g., the existence of differential motions of the stellar surface (e.g Cegla et al. 2016b) or even to probe the spectra of the star behind the planet (Dravins et al. 2021).
Focusing on the exoplanet studies, RVs have become a source for the robust detection of atmospheres and the chemical species that compose them, using a variety of techniques (e.g.Sing et al. 2009;Wyttenbach et al. 2015;Nikolov et al. 2018;Ehrenreich et al. 2020;Tabernero et al. 2021), primarily during transits.When an exoplanet with an atmosphere transits, the radiation from its host star is filtered in the evening and morning terminators, encoding information about the composition and physical properties of the underlying atmospheric processes in the observed stellar spectrum.These properties are a function of atmospheric pressure and temperature.At high altitudes, where the pressure is lower, processes such as atmospheric dynamics, the presence of clouds, hazes, and temperature inversions dominate over the typical chemical reactions timescales (e.g.Moses 2014).If we examine even higher altitudes, the intense radiation fields induce photochemical reactions that determine the atmospheric composition on these layers.In this region, at visible wavelengths, we frequently observe the signature of ionized alkali metals or molecules (e.g.Sedaghati et al. 2021;Azevedo Silva et al. 2022;Seidel et al. 2022).
In this paper, we analyze the HD 189733 system using highresolution ESPRESSO data.The star of this system is a K2V dwarf (Gray et al. 2003) located at a distance of approximately 19.3pc from Earth.It has a V-band magnitude of 7.6 1 and belongs to the group of variable stars known as BY Draconis (Koen et al. 2010).The known planet orbiting this star, HD 189733b, was one of the first hot Jupiters discovered (Bouchy et al. 2005) and has since been extensively studied due to its favorable planet-to-star radius ratio.It was the first exoplanet to have its surface temperature mapped (Knutson et al. 2007), and one of the first (along with HD 209458b) to have its atmosphere measured using spectroscopy with infra-red data from Spitzer (Grillmair et al. 2007).From the visible to the IR, this planet is known for the characteristic decrease in radius with increasing wavelength.This phenomenon is thought to be the caused by the interaction of very small particles (smaller than the wavelength which is being observed) in the upper atmosphere, and it is commonly referred to as Rayleigh scattering (e.g.Pont et al. 2008;Lecavelier Des Etangs 2006).The presence of this effect partially attenuates atomic and molecular signatures of the transmission spectrum.
Our goal in this work is to measure the broad-band transmission spectrum of HD 189733b using the high-resolution capabilities of ESPRESSO.To achieve this, we will use the chromatic Rossiter-McLaughlin (CRM, Di Gloria et al. 2015) method implemented in CaRM (Santos et al. 2020;Cristo et al. 2022).This approach for retrieving transmission spectra was first used by Snellen (2004), who attempted to measure an increase in RM amplitude near the sodium D lines.
In section 2, we provide a summary of the observations and data.Section 3 describes the method used to retrieve the transmission spectrum, along with an explanation pf the relevant effects in the RM model.This is followed by an analysis of the white-light fit in the subsequent section.Finally, in section 5, we 1 The stellar atmospheric parameters and chemical abundances of HD189733 can be found in Sousa et al. (2018).
examine the transmission spectrum to search for the presence of Rayleigh scattering and heavy metal signatures.

Observations
Two transits of HD 189733b were observed with ESPRESSO during the nights of August 11 and 31, 2021, as part of the guaranteed time of observation (GTO) under the program 1104.C-0350(T).ESPRESSO is a high-resolution fiber-fed spectrograph that covers the visible range from roughly 380 to 788 nm, distributed over 170 slices (2 slices correspond to one spectral échelle order).They were carried out at UT1 using the HR21 observing mode, with a 1" fiber, a spatial binning factor of two, and with a resolving power of approximately 140 000.
The observations were performed using two fibers, Fiber A pointed to the target, and Fiber B targeted at the sky.The integration time for each observation was set to 300s during both nights.For the first and second nights, a total of 41 and 43 data points were obtained, respectively, resulting in a total effective observation time of 3h 25m and 3h 35m.This resulted in an average SNR (signal-to-noise ratio) of 156 and 167, around 580 nm, and an average uncertainty for the RV measurements of approximately 32 cm s −1 on both nights.A summary of the observations is provided in Tab. 1, and the plots with the RV measurements, their uncertainties, as well as SNR and airmass variation, can be found in Fig. 1.
The data was reduced using the ESPRESSO data reduction pipeline, DRS, version 2.3.1 and 3.0.0.We proceeded with the RVs reduced with 2.3.1 since the latest version is seemingly more prone to jitter.The Cross-Correlation functions were computed using a K2 mask.The RVs were derived by fitting a Gaussian function to the sky-subtracted CCFs for each slice.One spectrum from the end of the second night was removed due to low SNR.

Simultaneous EulerCam photometry
We observed two full transits of HD 189733b with the Euler-CAM photometer (Lendl et al. 2012) at the 1.2m Euler-Swiss telescope located at La Silla observatory.The observations were carried out on 10 August 2021 and 30 August 2021 in the Gunn r ′ filter with an exposure time of 30 s and 10 s, respectively.The EulerCam data were reduced using the standard procedure of bias subtraction and flat-field correction.The transit lightcurves were obtained using differential aperture photometry, with a careful selection of reference stars and apertures that minimize the final light curve RMS.To account for correlated noise that affects the photometric data due to observational, instrumental and stellar trends, we used a combination of polynomials in several variables (time, stellar FWHM, airmass, coordinate shifts and sky background).The system parameters were obtained using CONAN (Lendl et al. 2017(Lendl et al. , 2020)), a Markov Chain Monte Carlo (MCMC) framework, by fitting for R P /R * , b, T 14 , P and T 0 , assuming wide uniform priors.The quadratic coefficients and their uncertainties for the photometric filter were calculated with the LDCU2 routine (Deline et al. 2022) and allowed them to vary in the fit with Gaussian priors.We also took into account additional white noise by adding a jitter term for each light curve.2.0 Fig. 1.Radial velocities 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, seeing evolution during the night, the signal-to-noise ratio at around 580 nm and airmass as a function of the number of days since the first epoch.The red dot represents the observation that was removed, from the second-night set, due to low SNR.

Method
The development of high-precision spectrographs and sophisticated data has enabled the measurement of RVs with unprecedented precision.Consequently, we can now measure the RM effect with higher detail.However, achieving this level of precision comes with challenges.At the sub-meter per second precision, there are second-order effects in RVs that must be modeled to avoid bias in estimation of orbital parameters.
In this section, we describe how we take advantage of CaRM modularity to retrieve the broad-band transmission spectra.To model the RM effect, we used a version of SOAP (e.g.Boisse et al. 2012;Oshagh et al. 2013;Dumusque et al. 2014;Akinsanmi et al. 2018;Zhao & Dumusque 2023) similar to the described in Serrano et al. (2020).This is an addition to the already available to use models proposed by Boué et al. (2013) and Ohta et al. (2005) (ARoME and RMcL implemented by Czesla et al. 2019).As we will discuss later, this approach allows us to address some approximations made in the previous models that simplify greatly the description of the stellar surface.

Modeling the RM effect
In a rotationally symmetric star without activity, the integrated projected velocity fields towards the observer cancel out exactly the portion associated with the stellar surface rotating away.However, when a planet transits, the planetary disk blocks light from the host star and the stellar surface behind it.Consequently, it becomes possible to measure the integrated velocity imbalance resulting from the unblocked portion of the star.The measured RV amplitude is primarily a function of the area being blocked (planet radius and atmospheric height), the rotation velocity of the star, and the impact parameter (Triaud 2018).
Over the years, several attempts have been made to model the classical RM effect with increasing accuracy.Most of these approaches try to express the RM anomaly though analytical formulations, which are constrained by the simplifying assumptions and symmetry conditions.These assumptions are primarily made to make the integration time to obtain the RV profile becomes practical, but it is also important to note that second degree phenomena create solutions that are not analytically exact.
Stellar inclination relative to the sky plane u i Limb-darkening coefficients (T e f f , σ T e f f ) (K) Effective temperature and uncertainty (log(g), σ log(g) ) Surface gravity and uncertainty (z, σ z ) (dex) Metallicity and uncertainty Notes.The first set corresponds to orbital and planetary properties.The second constitutes the set of stellar parameters, from which P rot , V CB , α B , α C , i ⋆ and u i can be free parameters of the RM modeling.
One such limiting assumption is that the underlying "quiet star"3 CCF remains constant across the stellar disk and can only be subject to displacements resulting from the projected rotational velocity (in longitude).As such, no latitudinal variations such as those associated with the differential rotation can be accounted for.
To overcome this problem, in this paper, we use an alternative formulation to measure the impact in RVs of a transiting exoplanet using the SOAP code.In short, the code simulates the star as a 2D disk with a grid with a user-defined resolution.We adopt a grid of 600 × 600 which strikes balance between speed and accuracy.For each point on the grid, the code computes the velocity shift to be applied to the "quiet-star" CCF, photometrically scaled to match the limb-darkening at the grid position.This CCF, the default from SOAP, was obtained by cross correlating an observation of the Fourier Transform Spectrograph (FTS) with a G2 HARPS template.The RV measurement results from the Gaussian fit to the sum of all CCFs on the grid.We exploit this numerical point-by-point RV shift computation to include the effect of center-to-limb variations (CLVs) induced by the convective blueshift (CB) and differential rotation.This changes and additional updates will be presented on a forthcoming publication of SOAP.

Convective Blueshift
The convective blueshift (Jewell 1896) is an RV signal that arises from the granular nature of stars with an external convective layer.On average, the hot and bright rising plasma in the granules contributes more significantly to the integrated radial velocities compared to the darker and cooler gas that sinks in the inter-granular spaces.The first degree changes induced by the radial component of the CB are the result from the projection effect along the disk and the photometric weight of the limb darkening.Similar to the rotational velocity fields, a transiting exoplanet introduces an unbalance in the perceived projected velocities, due to the CB, which is superimposed on the RM effect.
To incorporate the effects of the CLVs induced by the CB, we implemented in SOAP the first-order approximation of the effect described by Shporer & Brown (2011).We note, however, that more sophisticated approaches exist and can potentially better model the CB effect.However, these approaches involve additional physics (and the corresponding assumptions) based on magnetohydrodynamical simulations of the stellar surfaces (Cegla et al. 2016a).We consider that the first-degree polynomial description of the CLVs is a good compromise between model complexity and capturing the bulk of the CB effect.Nonetheless, this model has significant limitations, as it neglects the presence of meridional flows, variations in the CB strength and line shape along the disk for different chemical species (e.g.Liebing et al. 2021) and differential rotation.We expect this approximation to be reasonable for slow rotating solar-type stars, where the FWHM is greater than the rotation speed.Despite, even within this group, significant deviations from the profiles given by MHD models are possible (e.g.Cegla et al. 2016a).
To compute the CB-induced CCF shifts, we consider an initially limb-darkened perfect sphere.In each cell SOAP grid, we assume a constant CB velocity perpendicular to the surface.The magnitude at each point is radially symmetric and results from the projected component by the angle between the normal to the stellar surface and the line of sight (LOS) θ: V CB cos θ where V CB is the local CB velocity.In the literature, the measurement of the convective blueshift velocity is often represented by the solarscale factor S (e.g.Liebing et al. 2021).This quantity represents the ratio between the CB amplitude of the star and that of the Sun, where the solar value corresponds to −350 m s −1 .

Differential rotation
Differential rotation was initially detected on the Sun through the observation of variations in the migration rate of spots at different latitudes.Spectroscopically, it manifests as a latitudedependent shift of the spectral lines, resulting in measurable RV shifts (e.g.Livingston 1969).
Solar-like star, are known to possess strong magnetic fields that give rise to activity that is manifested on the surface (Hale 1908).In alpha-omega dynamos, such as the Sun, magnetic fields are generated by the dynamo that is located within the convective envelope (Parker 1955).The dynamo mechanism itself relies on the presence of turbulent velocities and differential rotation, both of which are observed in the Sun.
To account for the latitudinal effect of the differential rotation, we incorporate the formulation presented in, for example, Reiners (2003) or Gray (2005): where Ω(l) represents the angular velocity component at the latitude l, Ω eq denotes the angular rotation velocity at the equator, and α B represents the differential rotation coefficient.For a given star with differential rotation, the impact on the RM profile is not only greater in amplitude for higher coefficients but also more pronounced (and less degenerate with the other parameters) if the planet's transit path traverses a wider range of stellar latitudes (Roguet-Kern et al. 2022).
Table 3. Set of priors for the white-light fit with CaRM for HD 189733b.
Parameter Prior V sys † (km s −1 ) G(−2.18, 0.1) R p /R ⋆ U(0.14, 0.17) G(−0.85, 0.32) P rot (days) G(11.953, 0.1) Notes.In the table, the symbol G represents a normal distribution, where the first value corresponds to the mean and the second the standard deviation.Similarly, the symbol U represents a uniform distribution, with the respective lower and upper boundaries.
(*) We present the prior for σ W in units of m s −1 instead of the logarithmic, as it is more intuitive to perceive the range in velocity units. (†) The fitting process is performed independently for each dataset.

CaRM
CaRM4 is a semi-automatic code written in Python, designed to extract the broad-band transmission spectra of exoplanets using the chromatic RM method (Cristo et al. 2022).The code takes as input CCF files from HARPS and ESPRESSO (which are products from the default DRS) or user-specified text files containing the order/slice ranges, the observation times, radial velocities, and their uncertainties (organized by columns).
To initialize the code, the user modifies a setup file, which includes the path to the folders containing the reduced data (the CCF files) that correspond to the observations.CaRM scans the folders for CCF files with a user-defined suffix.It then converts the data to the specifies text file structure, employing lists with the spectral format of each instrument.At this step, the RVs are computed by performing a Gaussian fit to the weighted sum, by the variance, of the CCFs as defined in the range list.The uncertainties are computed assuming a photon-noise-limited observations, following the method described in Bouchy et al. (2001) but adapted to directly measure them from the CCF.
The input file allows the user to choose which model use to perform RM RV anomaly fit.Stellar parameters such as effective temperature, stellar surface gravity, and stellar metallicity are not directly fed into the models but are used instead to fit the limbdarkening coefficients usingLDtK (Parviainen & Aigrain 2015).
The code performs the model fitting using wither MCMC implementation emcee (Foreman-Mackey et al. 2013) or a dynamic nested sampler Dynesty (e.g.Speagle 2020; Koposov et al. 2022).Priors are defined as a dictionary, associating each parameter to be fitted with its prior distribution.Joint fits of specific parameters can be performed if multiple data sets are provided.The same prior definition scheme is employed for the subsequent fits performed on different wavelength ranges.The code assumes that the first element of the range list corresponds to the white light, and it uses this data to compute the colorindependent parameters with the highest SNR.The transmission spectrum is then constructed, after each chromatic RM variation is performed, with the measurements of the planet-to-star ratio as a function of the wavelength.

The white light fit
The chromatic Rossiter-McLaughlin effect relies on accurately fitting of the RM anomaly as a function of wavelength.Therefore, we begin by fitting the white light RVs (which correspond to the full bandpass of ESPRESSO) to constrain the colorindependent parameters (Tab.2) with the highest SNR the data can offer.

Priors and assumptions
For the white-light fit, we assumed a circular orbit that is expressed by a Keplerian with semi-amplitude K and a systemic velocity V sys .We allow V sys to vary as a free parameter with a uniform prior, as it can be influenced by nightly offsets, tel- luric contamination, and stellar activity, which may alter its measured value 5 .The out-of-transit slope, associated with the Keplerian orbit of the planet, and often parameterized by its semiamplitude K, is known to be affected by stellar activity (Boldt et al. 2020).Notably, K is not given directly in SOAP but as a function of the planetary mass.Given that we expect the planetary mass to remain constant during the transit, the Gaussian prior only accommodates the impact of stellar activity.For the mid-transit time, a shift term is introduced with a Gaussian prior.The mean value of the prior was set to zero, and a standard deviation was chosen to be compatible with the timescale of the transit.The mid-transit shifts can be originated from the uncertainties in the determination of the mid-transit time, which are amplified by the number of orbits since the epoch used as reference.Although we anticipate this value to be minimal due to extensive studies of this planet, other sources of mid-transit shifts, such as small-scale transit time variations (TTVs) or small orbital eccentricities, cannot be ruled out.We adopt uniform priors for α B , i ⋆ , V CB , and σ W since these parameters are poorly constrained in the literature or unknown for this target or these particular observations.Gaussian priors based on literature values were used for the remaining parameters, with the standard deviation matching the reported uncertainties.The width of the P rot prior was increased since it is used to compute the rotational velocity of the star, which presents some variability in the literature (model 1 or M1).It is important to note that in our analysis, we selected a source from the literature that obtained the rotation period value using an alternative method based on the evolution of stellar spots.This choice was made because it is possible that fitting the RM models alone may not be sufficient to fully resolve the degeneracies between the parameters in this well-aligned planetary system.By incorporating additional constraints from spot evolution studies, we aim to improve the accuracy and reliability of our results.In addition, we also tested a broader prior for the rotation period to access rotation periods that significantly deviate from the value that we used as reference before (model 2 or M2).Table 3 summarizes and completes the description of the set of free parameters and the respective priors we adopt.

Results
We run CaRM selecting Dynesty to perform the fit using a nested sampling approach.The number of live points is set to 500.The convergence is evaluated by monitoring the estimate of the logarithmic evidence and terminating the run when the variation is below 1%.This yielded 33 056 posterior samples for model M1 and 38 336 for model M2 during the static sampling phase.The corner plots illustrating the posterior distributions of the fitted parameters can be found in The fit for M2 is similar to M1 but they show higher dispersion (80 cm s −1 ).Based on this, we decide to adopt M1 for further study.It is important to note, however, that achieving lower residuals at the expense of over-fitting is a potential concern.To address this, we have verified that the dispersion of the residuals for M1 is still more than twice the expected photon noise.
There is no significant observable correlation between the residuals of the first and second nights.It exists, however, points that deviate significantly from the average RVs, which could be due to occultation events of active regions on the stellar surface, such as spots or plages.Notably, a larger radial velocity variation can be observed on the second night just after ingress, but there Notes.The uncertainties represent the 68.27% confidence interval around the median value. (*) We provide here the prior for σ W , in m s −1 , and not the logarithmic since it is more intuitive to perceive the range in velocity units ( †) Independently fit for each data set.
is no corresponding variation in the first night.The signal of the RV variation suggests a possible spot-crossing event.
To investigate the origin of the significant deviation in the residuals of the first night, we utilized simultaneous ESPRESSO and EulerCam observations.We found a good agreement both in magnitude and phase, b by modeling the residual radial velocities using SOAP with an occultation of a single spot at a stellar longitude of −30 • and latitude of 0 • .We assumed spot size of 0.1% and a temperature contrast of 660K.This produces a negative RV variation with an amplitude of approximately 3 m s −1 and a corresponding positive flux excess of 1300 ppm.For the photometric measurements, we rebinned the flux observations to approximately match the RVs at the same phase, as shown in Figure 2. The photometric measurements are compatible with the spot occultation scenario.However, since the effect in the flux is low when compared with the average error bars, we cannot rule out other possibilities.
Figure 2 displays the best-fit model obtained from combining the data of the two nights, as well as the best fit to the observed data.The joint model fit of the nights results in a RV scatter of 67 cm s −1 , which is approximately twice the predicted photon noise level.The error bar that includes the jitter in quadrature closely matches this value.By comparing our retrieval with these previous studies, we demonstrate the significance of incorporating models that account for center-to-limb variations (CLVs) and differential rotation to avoid potential biases (Cegla et al. 2016a;Cegla et al. 2016b).Additionally, the presence of strongly correlated signals in the residuals has been highlighted in previous works (e.g., Triaud et al. 2009;Cegla et al. 2016b;Cristo et al. 2022).

The stellar surface and true spin-orbit angle
The stellar surface has a significant impact on the shape of the RM curve, as we discussed before.Therefore, it is essential to analyze how effectively we can retrieve the fundamental parameters that describe using RV observations modeled with SOAP.
Figure 3 illustrates the posterior distribution diagrams for several keys parameters that are used to model the stellar surface.In particular, the equatorial rotation period of the star, the stellar rotation axis inclination, the differential velocity ratio, the local convective blueshift amplitude velocity, and the spin-orbit angle for M1 (right plot).
Our analysis yields a rotation period P rot ≈ 11.45 ± 0.09 days.Assuming a spherical star, this can be translated converted in linear rotation velocity dividing the equatorial perimeter by the rotation period.We estimate V eq = 3.38 ± 0.06 km s −1 .Since HD 189733b's orbit is aligned the planet crosses only a small number of stellar latitudes, which can be approximated reasonably well by the average latitude.To compare our results with other studies, we compute the average velocities of the stellar surface behind the planet, using the planetary and stellar inclinations, as well as the differential rotation of the stellar surface at the average latitude crossed by the planet.We approximate the latitudes crossed by the planet using the impact parameter, following a similar approach to Cegla et al. (2016b), with the expression V eq sin i ⋆ (1 − α(a cos i P ) 2 ).
Measurements of the projected rotation velocity for HD 189733 vary significantly in the literature.For instance, stellar activity photometric modelling yields 2.97 ± 0.22 km s −1 (Winn et al. 2006), while line broadening analysis by Bouchy et al. (2005) suggests 3.5 ± 1 km s −1 .Our approach is similar to the one adopted by Triaud et al. (2009) 6 where V eq sin(i ⋆ ) is computed from an RM model.In Triaud's paper, however, their best solution for the rotational velocity (3.316 +0.017 −0.068 km s −1 ) produces clear wave-like residuals and they don't consider that the star may have differential rotation or CB.The authors adjust the RM model to account for the observed trends, which results in a lower estimate V eq sin(i ⋆ ) = 3.05.Also, Cegla et al. (2016b) applied an analytic model of the Doppler shifts behind the planet, which includes the effect of CLVs and differential rotation, and fitted it to the residual CCFs.They find V eq = 4.50 +0.51  −0.49km s −1 and V eq sin(i ⋆ ) ≈ 3.3 km s −1 .Our solution for the sky-projected spin-orbit angle (λ) of HD 189733b is consistent with previous literature.We estimate λ = −1.00+0.22  −0.23 , which is in agreement with values ranging from −1.4 ± 1.1 • reported by Winn et al. (2006) to −0.35 ± 0.25 • reported by Campante et al. (2016).These results suggest a higher statistical probability for the stellar axis to be close to the orbital inclination.Statistical analyses of spin-orbit misalignment, such as those conducted by Campante et al. (2016) or Fabrycky & Winn (2009) provide valuable insights into the distribution of spin-orbit angles and their implications.Our finding of λ ≈ −1.00 supports the notion that the stellar axis is more likely to align closely with the orbital inclination, based on these statistical analyses.
In our model, we fit the stellar rotation axis angle.This is an important measurement to constrain models of planetary formation and evolution.For example, the angle between the Sun's rotation axis and the ecliptic plane is 7.15 • (e.g.Beck & Giles 2005) which represents only a small deviation from it.For the HD 189733 system, we retrieve a stellar axial rotation tilt of 71.87 +6.91 •  −5.55 • .This result is in line with the prediction from Henry & Winn (2007) which computes i ⋆ > 54 • with 95% confidence, with a most probable value of 65 • .An additional source for measurements of the stellar inclination axis can be found in Cegla et al. (2016b).Comparing the results for the model similar to ours (one parameter CB and differential rotation), they find 92.0 +11.0 • −3.8 • which is significantly different from our prediction as  2010) and Cegla et al. (2016b).Additionally, we compare these measurements with the differential rotation coefficients obtained from the analysis of photometric data of 24,124 Kepler stars (Reinhold & Gizon 2015).We represent these measurements as black dots.For the period, we selected the minimum one provided by the authors since it corresponds to the equatorial period for stars exhibiting solar-like differential rotation.
their posterior distribution clearly prefers higher values for the angle.
We additionally derive the differential velocity ratio α B ≈ 0.29 +0.12 −0.13 with 93.4% confidence for α B > 0.05.This means that, at the poles, the star rotates with a velocity 2.58 ± 0.35 km s −1 which corresponds to a rotational period of 14.9 ± 2.0 days.In Cegla et al. (2016b) we can find both one and two parameters (using a similar parametrization as can be found on Eq. 1) models to describe the latitudinal change of the stellar surface velocity.The authors find α B > 0.1 with 99.1% confidence and amplitude that ranges between 0.3 and 0.86 with one parameter (derived from HARPS data).A more precise determination is available on Fares et al. (2010), where it is derived α B ≈ 0.278±0.093from polarimetry.We compare these measurements with our solution and the observed distribution of differential rotation with the stellar rotational period in Fig. 4 and find a good agreement with the literature values.
HD 189733 is a K dwarf and, as such, it is expected to have a granular surface like the Sun but with a lower contrast between the granule centers and the inter-granular lanes.Using the formulation for the CLVs produced by the CB, the local convective blueshift must be sub-solar (or S factor lower than one) (Liebing et al. 2021).We are able to measure V CB ≈ −211 +69 −61 m s −1 which represents a scale factor S ≈ 0.60 +0.20  −0.18 which is sub-solar as expected.Despite the value, in RV amplitude for the CLVs created by the CB, being lower than the differential rotation, it contributes to a deformation symmetrical about the mid-transit to the RM profile (Fig. D.1).This characteristic allows a confident determination of the value for the CB itself, with V CB < −50 m s −1 with a level of 99.8%.We compare our result with Liebing et al. (2021) in Fig. 4. We depict the Doppler maps generated by the combination of the effects described in this section, as well the contribution of each of them, in

Chromatic RM fit
Similarly to the white light, we computed the RVs that results from the sum of CCFs that are defined in the specific slice intervals (see Tab. F.1).These RVs reflect, in practice, the impact of the planetary transit on the stellar spectrum, which results in an RM profile that captures the sum of the planetary radius and scale height for the particular wavelength interval.CaRM fits these chromatic Rm s using a priori information from the white light.We fix the wavelength-independent parameters such as the differential rotation velocity ratio, rotational period, or planetary inclination such that biases are not introduced in the chromatic radius determination (the RVs here have a substantially lower SNR).We fit similarly the systemic velocity of the star for each night with a uniform prior to account for a potential wavelength-dependent offset due to stellar activity.For the same reason stated for the white light, in addition to the chromatic variations, we let the planetary mass change around the literature value with comprehensive 10% width.We also give a uniform prior to the convective blueshift since it is expected to be a function of the wavelength (Cegla et al. 2016b), as it results from the contribution of different spectral lines that are formed at different depths in the photosphere (and, as such, with different velocity contributions).To the jitter amplitude, it is attributed a broader prior, when compared with the white-light, since each bin contains less RV information.For the RM fit binby-bin approach used here, we assume that the stellar spectrum shape doesn't change as a function of the distance to the disk center (e.g.Dravins et al. 2021).There may be additional parameters that change as a function of the wavelength that we are not considering.The set of priors for the chromatic priors are summarized in Tab.2022) in green.Shifts were applied to the literature results to match, approximately, the radius level @550nm.
Table 5. Set of priors for the chromatic fit.The elements of the table as the same meaning as Tab. 3.

Transmission spectrum retrieval
The transmission spectrum derived from ESPRESSO data, Fig. 5, shows a steep decrease in planetary radii as a function of increasing wavelength.In the wavelength range where they coincide, the transmission spectrum derived in this paper is globally consistent with the spectrum obtained from HARPS data with the same technique (Cristo et al. 2022).Our retrieval, however, suggests an enhanced slope when compared with the spectrophotometric spectrum obtained with the observations from Hubble (STIS).The observed slope is classically attributed to the scattering of particles with a size smaller than the incoming light's wavelength.In literature (e.g.Lecavelier Des Etangs et al. 2008), this radius-wavelength slope is often parametrized as: where α corresponds to the scattering slope and H to the atmospheric scale-height.When α = −4 this kind of interaction is often called as Rayleigh scattering.There are several exoplanets that exhibit 'Super-Rayleigh' (or α<-4) slopes, such as the hot Jupiter WASP-19b with α ∼ −35 (Sedaghati et al. 2017) or the super-Neptune HATS-8b with α ≃ −26 ± 5 (May et al. 2019).Studies of the planet population of Sing et al. (2016) show that, in general, planets have an enhanced slope (e.g.We perform the first analysis of our retrieved transmission spectrum by fitting a linear model to the transmission planet radii as a function of the logarithmic wavelength, using Eq. 2. The result of the fit and confidence bands can be observed in Fig. 5.We measure a decrease in planet radius of ∆R p /R ⋆ = −0.0082± 0.0013 along the wavelength range of ESPRESSO, which corresponds to a scattering slope of −31 ± 5 assuming a scale height of 190 km (Kasper et al. 2018).The slope obtained through this analysis is much higher than values we can find in literature up to date.The reasons behind this excessive slope are partially addressed in Oshagh et al. (2020), where the authors obtain a lower radius variation computing the RVs with a custom mask combined with the effect of stellar activity.Using only custom CCFs they derive α ≃ 9.77 ± 2.72.In this work, we chose to use the default ESPRESSO masks and try to evaluate the effect of stellar activity afterwards.It is not yet certain what the weight of these effects is and whether ESPRESSO masks suffer from the same problems.
To interpret the transmission spectrum, obtained with the CRM technique, we additionally use the forward modeling capabilities of PLATON (Zhang et al. 2019(Zhang et al. , 2020)).We use PLATON assuming an isothermal atmosphere with equilibrium chemistry, divided by default into 500 layers.The physics of the atmosphere is governed by the hydrostatic equation, setting the planetary radius reference pressure to 1 bar.In each layer, the code includes the contribution from gas absorption, collisional absorption, and Rayleigh scattering.The contribution of the gas absorption is computed by solving the radiative transfer equation.The presence of different species changes the absorption coefficient of the different layers, changing the opacity as a result.The Rayleigh scattering is controlled by a parametric law, where it is possible to control the wavelength dependence and the slope strength.In alternative, PLATON supports Mie scattering, which models the contribution of hazes to the transmission spectrum.This model is a function of the imaginary refractive index of the particles, particle size, geometric standard deviation, and fractional scale height.Furthermore, it is possible to select between the already For active stars, the presence of unocculted spots and plages is known to change the perceived planet-to-star radius ratio, (McCullough et al. 2014;Rackham et al. 2018) which not only impacts photometric measurements but also RV measurements during transits (Oshagh et al. 2014;Boldt et al. 2020).In addition, the inflation does not only affect the white-light radius but can also introduce spurious trends on the transmission spectrum that can mimic a true signal from the planet.To account for this, PLATON introduces a model that is controlled by the fraction of the surface that is plagued with activity and the temperature contrast between these regions and the solar surface.One more additional and important source of attenuation of the spectral features are clouds.With PLATON it is possible to define the pressure of the cloud deck.Below it, the atmosphere is fully opaque and as a consequence, there is no contribution from these layers.
We performed several retrievals to understand what is the main mechanism behind the observed slope in the transmission spectrum.The uncertainty was propagated, using a Gaussian prior, for the stellar radius, planetary mass, planet radius, planet temperature, C/O ratio, and atmospheric metallicity.For each of the runs, we computed two distinct inference criteria.ln(z) which corresponds to the natural logarithm of the evidence, that is used to compare directly distinct models, and the χ 2 which is a weighted (by the variance) measure of the distance from the model to the data points, and it can be used to tell us how well the models reproduce the data.
We started by constructing a flat model that includes all species and setting the scattering slope to zero.The uncertainties are propagated for the stellar radius, planetary mass (M p ), planetary radius, and temperature from the literature using Gaussian priors.In addition, we also propagate in the same way for the log-metallicity (log(Z)), C/O ratio, and temperature of the atmo-sphere to the values derived in Zhang et al. (2020).In order to try to reproduce the observed slope in the transmission spectrum, we used the parametric law for the Rayleigh scattering with a uniform prior for the scattering slope and a comprehensive range for the scattering strength k α .The priors are summarized in Tab. 6. Notes.The reference parameters are from Zhang et al. (2020), and the planet-to-star radius ratio results from the white-light fit from the ESPRESSO data.In the last column, the fit median values and uncertainties.
Fig. 6 presents our results for the models with and without scattering.For the flat model (R0) we obtain a log evidence ln(B R0 ) = 154.06± 0.10 and a χ 2 = 63.76.We retrieve the bestfit model with the parametric scattering (R1) with ln(B R1 ) = 160.07± 0.11 and χ 2 = 33.68.The Bayes factor between this model and the flat one is 407 ± 61, which represents very strong evidence against R0 (Kass & Raftery 1995).The model seems also to better reproduce the global trend of the data, especially in the bluest range where it is able to fit the small plateau that is observed.However, it is unable to explain the radius measure-ments centered at 468.975 nm and 680.295 nm which are of by over 2-sigma.This model produces a decrease in radius over the entire range of wavelengths of ∆R p /R ⋆ = −0.0046± 0.0008, which is significantly different from the value we derived before.Increasing the number of species in the model seems to explain, in part, the observed slope and produces a significantly lower α ≃ 22.4 ± 7.

Conclusions
In this paper, we used the chromatic Rossiter-McLaughlin effect applied to high-resolution ESPRESSO data to retrieve the transmission spectrum of HD 189733b.
We started by fitting the white-light RVs with a model composed of a Keplerian component and the effect of a transiting exoplanet on the integrated CCF.
We tentatively detect the presence of differential rotation with a confidence of 93.4% and derive an equatorial rotation period of 11.45 ± 0.09 days and a polar period of 14.9 ± 2 days.Additionally, using a first-degree surface brightness (CB) model, we derive a convective blueshift scale factor of S ≈ 0.60 +0.20  −0.18 .We also test a broader range for the rotation period, which increases the confidence for the differential rotation ratio being larger than 0.05 to 99.6%.The median value for the convective blueshift is compatible in both scenarios.
The presence of differential rotation breaks the symmetry of the stellar RV field, which further allows the determination of the stellar tilt.We find i ⋆ ≈ 71.87 +6.91 • −5.55 • and a projected spinorbit angle of λ ≈ −1.00 +0.22 • −0.23 • .In turn, we compute the true 3D spin-orbit angle (ψ ≈ 13.6 ± 6.9 • ) and note that the planet seems to be well aligned.
We analyze the transmission spectrum by first fitting a simple Rayleigh scattering model to the data.We find a greatly enhanced Rayleigh scattering slope, often called Super-Rayleigh, of −31±5 and a radius decrease of ∆R p /R ⋆ = −0.0082±0.0013.Using the forward retrieval software PLATON, with all the atomic and molecular species available, we estimate a significantly lower slope (α ≃ −22.4 ± 7) and radius variation: ∆R p /R ⋆ = −0.0046±0.0008.We reproduce the plateau observed in the blue range (378.05− 443.26 nm), which cannot be explained by scattering alone.The enhanced radius-wavelength slope has been explored in the literature, with some authors being able to emulate the decrease using models of the stellar surface with unocculted cold spots (e.g., McCullough et al. 2014;Kasper et al. 2018;Rackham et al. 2019) or occultation of plages (Boldt et al. 2020).However the origin of the stronger slope observed in groundbased observations requires further investigation.It is plausible that STIS observations may be less sensitive to stellar activity, as they primarily rely on variations in flux contrast on the stellar surface during transits.Additionally, spectroscopic observations can exhibit line profile deformations induced by stellar activity, potentially varying with wavelength.This characteristic could be specific to active stars since a similar approach was used in (Santos et al. 2020) to retrieve the broadband transmission spectrum of HD 209458b did not result in any detectable strong slope in the blue wavelengths.To understand the origin of these differences and determine the more reliable method, it is crucial to conduct high precision simultaneous observations in the future.

Fig. 2 .
Fig. 2. Fit and residuals of the white light data.Top: The best-fit model (solid red line) obtained from the combined data of the nights observed with ESPRESSO.Bottom: Residuals after subtracting the model.The light orange areas represent the phases where the planet is fully inside the stellar disk, while the lighter orange regions correspond to the ingress and egress phases.The data points are represented with two error bars.The green error bar is computed from the CCF, assuming a photon-noise limited observation.The black error bars is obtained by adding in quadrature the value of the green error bars and the jitter amplitude.At the top of the bottom figure, two quantities are presented.The first is the average dispersion of the residual RVs (σ res ) in units of m s −1 .The second is the average value of the black error bars, also given in m s −1 .

Fig. 3 .
Fig. 3. Posterior distribution diagrams for M1 and M2.The corner plots depict the posterior distribution for the equatorial rotation period, stellar axis inclination, differential rotation coefficient, local CB velocity and projected spin-orbit angle.The black contour lines represent (from center to outwards) the confidence intervals enclosing 68.27% 95.45% and 99.73% of the accepted samples.The histograms display the parameter posterior distributions, with the darker dashed line indicating the median value lighter lines delimiting the 1-sigma interval.
Fig. A.1 and Fig. B.1.The fits of the individual nights with the best-fit model and corresponding residuals after subtracting it are in Fig. C.1.

Fig. 4 .
Fig. 4. Estimates of the solar-relative CB and differential rotation ratio, along with a comparison to existing literature.Left: The solar-relative scale factor as a function of the effective temperature for a sample of stars ranging from M to F spectral types (Liebing et al. 2021).The red point and error bars represent the prediction from our best-fit model.Right: Relative differential rotation coefficient as a function of the rotation period of the stellar surface at the equator.We compare our retrieval with the results from Fares et al. (2010) andCegla et al. (2016b).Additionally, we compare these measurements with the differential rotation coefficients obtained from the analysis of photometric data of 24,124 Kepler stars(Reinhold & Gizon 2015).We represent these measurements as black dots.For the period, we selected the minimum one provided by the authors since it corresponds to the equatorial period for stars exhibiting solar-like differential rotation.
5 and the fit for each wavelength bin is represented in Fig. D.2.

Fig. 5 .
Fig. 5. Broad-band transmission spectrum of HD 189733b.Left: Linear fit (black dashed) to the measured planet-to-star radius ratio as a function of the logarithmic wavelength.The red bands represent the one to three σ confidence levels.Right: Comparison of our results with Pont et al. (2013) in blue and Cristo et al. (2022) in green.Shifts were applied to the literature results to match, approximately, the radius level @550nm.

Fig. 6 .
Fig. 6.Fit of the transmission spectrum for HD 189733b.From top to bottom: Transmission spectrum of HD 189733b modeled with PLATON with a flat model (red colors) and a model transmission spectrum with a variable Rayleigh-like slope.Residuals after subtraction of the best-fitting model to the flat model and the Rayleigh-like parameterizations.

Table 1 .
Observation summary of HD 189733b observations.

Table 2 .
CaRM parameters using the SOAP model.

Table 6 .
Set of priors for the transmission spectrum fit (R1) with PLATON.