The Hot Neptune WASP-166~b with ESPRESSO I: Refining the Planetary Architecture and Stellar Variability

In this paper, we present high-resolution spectroscopic transit observations from ESPRESSO of the super-Neptune WASP-166~b. In addition to spectroscopic ESPRESSO data, we analyse photometric data from {\sl TESS} of six WASP-166~b transits along with simultaneous NGTS observations of the ESPRESSO runs. These observations were used to fit for the planetary parameters as well as assessing the level of stellar activity (e.g. spot crossings, flares) present during the ESPRESSO observations. We utilise the Reloaded Rossiter McLaughlin (RRM) technique to spatially resolve the stellar surface, characterising the centre-to-limb convection-induced variations, and to refine the star-planet obliquity. We find WASP-166~b has a projected obliquity of $\lambda = -15.52^{+2.85}_{-2.76}$$^{\circ}$ and $v\sin(i) = 4.97 \pm 0.09$~kms$^{-1}$ which is consistent with the literature. We were able to characterise centre-to-limb convective variations as a result of granulation on the surface of the star on the order of a few kms$^{-1}$ for the first time. We modelled the centre-to-limb convective variations using a linear, quadratic and cubic model with the cubic being preferred. In addition, by modelling the differential rotation and centre-to-limb convective variations simultaneously we were able to retrieve a potential anti-solar differential rotational shear ($\alpha \sim$ -0.5) and stellar inclination ($i_*$ either 42.03$^{+9.13}_{-9.60}$$^{\circ}$ or 133.64$^{+8.42}_{-7.98}$$^{\circ}$ if the star is pointing towards or away from us). Finally, we investigate how the shape of the cross-correlation functions change as a function of limb angle and compare our results to magnetohydrodynamic simulations.


INTRODUCTION
Stellar magnetic activity and surface phenomena, such as flares, spots, differential rotation (DR) and centre-to-limb convective velocity variations (CLV), can alter the observed stellar line profiles and induce Doppler shifts which can cause biases in the calculations of planetary companion properties (e.g. Saar & Donahue 1997;Oshagh et al. 2016;Cegla et al. 2016b,a). Furthermore, DR plays a critical role in dynamo processes which are largely responsible for the generation of magnetic fields (e.g. Kitchatinov & Olemskoy 2011;Karak et al. 2020). Therefore, understanding these phenomena is not just important for exoplanet characterisation, but for magnetic activity as a whole.
When a planet transits a host star, a portion of the starlight is blocked in the line of sight and a distortion of the velocities is observed, known as the Rossiter-McLaughlin (RM) effect (see Rossiter (1924); McLaughlin (1924) for original studies and Queloz et al. (2000) for the first exoplanet case). The Reloaded RM (RRM) technique (see Cegla et al. 2016a) isolates the blocked starlight behind the planet to spatially resolve the stellar spectrum. The isolated starlight from the RRM can be used to derive the projected obliquity, , (i.e. the sky-projected angle between the stellar spin axis and planetary orbital plane). If the planet occults multiple latitudes, we can determine the stellar inclination (by disentangling it from sin ). Additionally, if you know the stellar rotation period (P rot ) then you can use this in combination with sin to determine the stellar inclination * . Then we can use this to model for and determine the 3D obliquity, .
Here, is of great importance as it can provide insights into planetary migration/evolution and avoids biases introduced from only having . Additionally, the isolated starlight can be modelled to account for both the stellar rotation behind the planet and any CLV. Sun-like stars which possess a convective envelope have surfaces covered in granules, bubbles of hot plasma which rise to the surface (blueshift), before cooling and falling back into intergranular lanes (redshift). The net convective velocity shift caused by these granules changes as a function of limb angle (i.e. from the centre to the limb of the star) due to line-of-sight changes and the corrugated surface of the star. For example, towards the limb the granular walls become visible (due to geometrical effects) and the bottom and tops of the granules become less visible. Depending on the stellar type, this can cause the net convective blueshift (CB) to be less than compared to disk centre, where the walls of the granules are less visible and dimmer than the top and bottoms. This centre-to-limb CB has been observed on the Sun with RVs changing on the level of 100ms −1 (see Dravins 1982) and on other stars including the K-dwarf HD 189733 (see Czesla et al. 2015;Cegla et al. 2016a). Overall, these velocity shifts can impact the RM effect which is used to determine the projected obliquity (Cegla et al. 2016b). Therefore, ignoring the velocity contributions caused by granulation can bias or skew these measurements and impact our understanding of planet formation and evolution. In addition to the RRM, other techniques utilising transiting exoplanets have been used to probe, detect and map the presence of starspots on the surface of stellar hosts (see e.g. Silva 2003;Wolter et al. 2009;Sanchis-Ojeda & Winn 2011;Morris et al. 2017;Močnik et al. 2017;Zaleski et al. 2020, and references therein).
In our Solar System, all eight planets orbit the Sun in line with its rotation and are nearly aligned, with obliquities ∼7 • ), suggesting all of the planets formed within a protoplanetary disk. Other star-planet systems have been found to be well-aligned and overall ∼ 120 of starplanet systems with measured obliquites are aligned. For example Knudstrup & Albrecht (2022) measure = -2±6 • in the HD 332231 system (a warm giant planet orbiting an F8 star, P orb = 18.7 d). They argue that measurements of obliquities are more meaningful for long period planets as a result of tidal interactions between the star and planet changing the star-planet alignment. Christian et al. (2022) investigate the effects of wide visual binary companions on the formation and evolution of exoplanets. They find that planets to these binary systems tend to be aligned with the binary, demonstrating the alignment is not only due to the formation of the planet. There have been misaligned systems found with obliquities greater than 30 • such as WASP-8 b (Queloz et al. 2010;Bourrier et al. 2017) and WASP-127 b (Allart et al. 2020) in retrograde orbits ( ∼ 180 • ; approximately 11 known), and WASP-79 b (Brown et al. 2017) with a polar orbit ( ∼ 90 • ; approximately 8 known). It has been suggested these misaligned orbits could be a result of Kozai migration (e.g. GJ 436 b Bourrier et al. 2018) where the orbit of a planet is disturbed by a second body causing the argument of pericentre to oscillate around a constant value resulting in an exchange between its eccentricity and inclination, coupled with tidal friction (e.g. Fabrycky & Tremaine 2007;Nagasawa et al. 2008). Additionally, other mechanisms include Secular resonance crossing to explain the polar orbits of Neptune mass planets, magnetic warping which can cause young protoplanetary disks to tilt towards a perpendicular orientation and high eccentricity tidal migration which is often used to explain hot Jupiters (see Albrecht et al. 2021, for more details). Therefore, classifying the obliquities of planetary systems is a vital tool for understanding how planets arrive in their orbits and what their fates will be over time (see the review by Albrecht et al. 2022).
For this study we focus on the WASP-166 b system which is a transiting super-Neptune, within the Neptune desert, with a mass of p = 0.101 ± 0.005 M Jup and radius p = 0.63 ± 0.03 R Jup (Hellier et al. 2019). The star is a bright, = 9.36, F9 main sequence dwarf with an age of 2.1 ± 0.9 Gyr. It was first discovered by WASP-South and later followed up by HARPS and CORALIE with a total of 220 and 41 RVs respectively (see Hellier et al. 2019, for the full analysis). Hellier et al. (2019) analysed RM observations with the HARPS spectrograph and found WASP-166 b to be aligned with the stellar rotation axis with = 3±5 • . In another study, Kunovac & et al. (2022) used the same HARPS observations to re-analyse the system using the RRM technique. They also found the system was aligned with = 1.3 +2.6 −2.0 • . WASP-166 has also been central to various other studies: Bryant et al. (2020) use simultaneous photometric transit observations of WASP-166 b from TESS and NGTS and found the multi-scope precision of NGTS can match TESS. Finally, Seidel et al. (2020) use HARPS and then ESPRESSO (Seidel et al. 2022, same observations as in this paper) observations of WASP-166 b to detect neutral sodium in the planet atmosphere at a 3.4 level and tentative evidence of line broadening. This suggests winds are responsible for pushing sodium further into space causing the bloated nature of this highly irradiated world.
In this paper we apply the RRM technique on newly acquired ESPRESSO observations of the WASP-166 b system to spatially resolve the stellar surface, by subtracting in-transit spectra from out-of-transit spectra, to characterise stellar differential rotation and centre-to-limb convection-induced variations, and to determine the star-planet obliquity. In §2 we detail all of the photometric and spectroscopic observations used in this study. §3 then details the transit analysis of the TESS sectors 8 and 35 and new NGTS photometric lightcurves where the transit parameters of the system are derived. Finally, we discuss the RRM analysis and results in §4 and §5 respectively.

OBSERVATIONS
For the analysis of the WASP-166 system we use both photometric data from TESS and NGTS as well as spectroscopic data from ESPRESSO. Here we give details of the observations and data reduction in more detail, where a summary of the observations can be found in Table 1

Photometric Data
To isolate the stellar spectrum behind the planet, in-transit spectra are directly subtracted from out-of-transit spectra. This is made possible by scaling the CCFs using photometric lightcurves. As a result, for WASP-166 b we utilise six TESS and two NGTS photometric transits to determine the transit, planet properties, and to scale the ESPRESSO CCFs.

TESS Photometry
WASP-166 was observed by the Transiting Exoplanet Survey Satellite (TESS: Ricker et al. 2014) in 2-min cadence capturing a total of seven transits during Sector 8 and Sector 35 between the 2 to 27 February 2019 and 9 February to 6 March 2021, respectively. We accessed the 2-min lightcurves produced by the TESS Science Processing Operations Centre (SPOC) pipeline (Jenkins et al. 2016) and use the PDCSAP_FLUX time series for our analysis. We used the transit ephemerides from Hellier et al. (2019) to mask out the in-transit data points and then ran a Lomb-Scargle analysis to search for and determine a stellar rotation period. Our analysis yielded no significant periodic signals. From a visual inspection of the light curve it is clear that there is some low level out-of-transit variability (see Figure B1); however, it is unclear whether this is astrophysical or instrumental in nature. To prevent this variability affecting our transit fit, we flatten the light curve using a spline fit, again masking out the in-transit data points to avoid the spline affecting the transit shapes. The methodology used may underestimate the errors on the planetary properties but this does not affect the RRM analyis in this work. The TESS photometry and detrending spline are shown in Figure B1.

NGTS Photometry
The Next Generation Transit Survey (NGTS; Wheatley et al. 2018) is a photometric facility consisting of twelve independently operated robotic telescopes at ESO's Paranal Observatory in Chile. This is the same location as ESPRESSO, therefore, both instruments experience the same weather conditions. Each NGTS telescope has a 20 cm diameter aperture and observes using a custom NGTS filter (520 -890 nm).
We observed WASP-166 using NGTS simultaneously with the ESPRESSO transit observations on the nights of 31st December 2020 and 18th February 2021 (see Figure 1). For both NGTS observations, we utilised six NGTS telescopes simultaneously to observe the transit event in the multi-telescope observing mode described in Smith et al. (2020) and Bryant et al. (2020). For both nights, an exposure time of 10 seconds was used and the star was observed at airmass < 2. On the night of 31st December 2020 a total of 9654 images were taken across the six telescopes, and 13845 were taken on the night of 18th February 2021.
The NGTS images were reduced using a version of the standard NGTS photometry pipeline, described in Wheatley et al. (2018), which has been adapted for targeted single star observations. In short, standard aperture photometry routines are performed on the images using the SEP Python library (Bertin & Arnouts 1996;Barbary 2016). For the WASP-166 observations presented in this work, we used circular apertures with a radius of five pixels (25 ). During this reduction comparison stars which are similar in magnitude and CCD position to the target star are automatically identified using the Gaia DR2 (Gaia Collaboration et al. 2016Collaboration et al. , 2018 and parameters found in the TESS input catalogue (v8; Stassun et al. 2019). Each comparison star selected is also checked to ensure it is isolated from other stars. Faint nearby stars within the aperture can contaminate the photometry and make the observed transit appear shallower than it is in reality. We queried the Gaia DR2 (Gaia Collaboration et al. 2018) and found four faint stars ( = 18.61 − −19.87 mag) within 25 of WASP-166. Based on these neighbours we estimate a maximum possible contamination factor of just 0.049%. Therefore, we expect the contamination from these faint stars to not have a significant impact on our results and so, we treat the NGTS data as undiluted throughout the analysis in this paper.

Spectroscopic Data
Two transits of WASP-166 b were observed on the nights of 31st December 2020 (run A) and 18th February 2021 (run B) using the ESPRESSO (Pepe et al. 2014(Pepe et al. , 2021 spectrograph (380 -788 nm) mounted on the Very Large Telescope (VLT) at the ESO Paranal Observatory in Chile (ID: 106.21EM, PI: H.M. Cegla). The ESPRESSO observations were carried out using UT1 for the first night and UT4 for the second under good conditions, with airmass varying between 1.0 -2.5" and 1.0 -1.7" for each run A and run B, respectively, in the high resolution mode using 2 × 1 binning. Exposure times were fixed at 100 s for each night to reach a SNR near 50 at 550 nm (to be photon noise dominated) and to ensure a good temporal cadence, with a 45 s dead-time per observation due to read out. Each run covered a duration of 6h 40m and 6h 15m of uninterrupted sequences covering the full transit duration and includes 2 -3 h of pre and post baseline. A summary of the ESPRESSO observations can be found in Table 1.
The spectra were reduced with version 2.2.8 of the ESPRESSO data reduction software (DRS), using a F9 binary mask to create cross-correlation functions (CCFs) which we use for our analysis. Additionally, the DRS also outputs the contrast (i.e. depth), full width at half max (FWHM: i.e width) and centroid (i.e. radial velocity) of each CCF. In run A the SNR can be seen to increase over the duration of the night which correlates with the decreasing airmass as observing conditions improve. However, in run B the SNR decreases towards the end of the night as a result of passing thin clouds. Despite this, the FWHM and contrast remain steady during both runs and are dispersed around the mean. Overall, the median integrated radial velocity uncertainties for run A and run B are 1.57 ms −1 and 1.77 ms −1 , respectively.

PHOTOMETRIC TRANSIT ANALYSIS
When computing the RRM we spatially resolve the stellar surface by subtracting in-transit spectra from a master out-of-transit spectra. To do this, the spectra are scaled using a modelled transit lightcurve to remove any variations caused by the Earth's atmosphere. This enables a direct subtraction between the spectra and avoids any assumptions on the local stellar line profile shapes. As a result, we need good knowledge of the transit parameters of the system, which we obtain from the spline detrended TESS photometry as the TESS data are more precise for this analysis. We identified that one of the TESS transits falls close to the start of the second orbit of Sector 8. TESS photometry at the start of orbits often suffers from increased systematic trends. Therefore, to avoid these trends affecting our analysis through this transit, we exclude it from the analysis. We fit the TESS photometry (600 -1040 nm) using the MCMC method from the package (Foreman-Mackey et al. 2013) and the transit model (Kreidberg 2015). We take as free parameters: the time of transit centre ( ), the orbital period (P orb ), the planet-to-star radius ratio ( p / * ), the scaled semi-major axis (a/ * ), the orbital inclination ( p ), the limb darkening coefficients ( 1 and 2 ), and a constant flux offset ( 0 ). For the limb darkening, we use a quadratic law and calculate values for the coefficients using LDTK (Parviainen & Aigrain 2015) and the stellar parameters from Hellier et al. (2019). We calculate values of 1, LDTK = 0.3955 ± 0.0023 and 2, LDTK = 0.1523 ± 0.0038, which we adopt as Gaussian priors for the fitting. For the remaining parameters we adopt the Hellier et al. (2019) values as initial positions for the walkers and uninformative priors that ensure a physically realistic solution. We run the MCMC using 20 walkers for 10000 steps per walker, following a burn in processes of 3000 steps. The transit model obtained from this fitting is displayed in Figure 2 and the transit parameters are presented in Table 2.
The NGTS transits (520 -890 nm) are then individually fit to investigate whether the transits are significantly early or late compared to the prediction from the TESS model. This is a vital step for computing the RRM technique as we need the most precise ephemeris possible. For these NGTS fits, we produce a transit model using the TESS fit parameters, and take , p / * , and the limb darkening coefficients as the only free parameters. We compute new limb darkening coefficients of 1, LDTK = 0.4531 ± 0.0026 and  Table 2.
2, LDTK = 0.1492 ± 0.0042 for the NGTS band pass, which we use as Gaussian priors. We also impose a Gaussian prior on p / * of 0.052 ± 0.001 based on the results from the fit to the TESS data.
During the NGTS fitting, we detrend the six independent telescope light curves for each night. We do this by fitting a linear out-of-transit model with time to each light curve simultaneously with the transit model. From these fits, we measure NGTS transit times of C,NGTS−1 = 2459215.739534 ± 0.0009405 and C,NGTS−2 = 2459264.729337 ± 0.000633. These transit times are consistent with the TESS predictions at the levels of 0.83 and 1.47 . Therefore, we do not have any significant indication of timing variations for the two transit events observed with ESPRESSO.
However, slight inaccuracies in the ephemeris at the time of the ESPRESSO observations can bias the RRM analysis (Kunovac Hodžić et al. 2021). Therefore, we fit all the photometric data together to derive the most accurate ephemeris for WASP-166 b. We use the flattened TESS and the two detrended six telescope NGTS observations, as well as the nine telescope NGTS transit observation of WASP-166 b from Bryant et al. (2020). For this final fit, we use the parameters derived from the TESS analysis, and allow and P orb to vary. We find values of = 2458524.40869 ± 0.00030 and P orb = 5.443542± 0.000003 days. The new ephemeris (determined from both NGTS and TESS transits) along with the planet-to-star radius ratio ( p / * ), the scaled semi-major axis (a/ * ), the orbital inclination, ( p ) and the limb darkening coefficients ( 1 and 2 ) (from the TESS alone transit fitting) are detailed in Table 2 as the final planetary properties used in the RRM analysis.
In Section 2.1.2, we mentioned that we expected any dilution of the transit depth due to contamination of the NGTS aperture to not have a significant effect on our results. From this analysis, we derive a transit depth of 3137.6 ± 23.8 ppm. Assuming the contamination factor of 0.049% obtained in Section 2.1.2, we calculate an "undiluted" depth of 3139.1 ± 23.8 ppm. The difference between these two depths is just 1.55 ± 33.62 ppm, clearly confirming that any possible effect had by contamination from the faint stars within the NGTS aperture on the results of this study will be insignificant.
The differences in bandpass between ESPRESSO, NGTS and TESS means there will be slight changes in the limb darkening and p / * . Therefore, we tested the impact these differences would have on the local RVs (see Figure 4) using the LDTK python package to determine limb darkening parameters in the ESPRESSO bandpass and using the p / * determined in the NGTS bandpass from Bryant et al. (2020). Overall, we found there was no significant change in the local RVs and any change was within 1-sigma. Additionally, we checked the impact of any inaccuracies in the ephemeris on the local RVs by using the upper and lower limits of the errors on T 0 and P rot from Table 2; we found this did not have a significant change in the RRM analysis. Finally, the simultaneous NGTS transits are used to ensure there has been no contamination with stellar activity during our ESPRESSO observations.

THE RELOADED ROSSITER MCLAUGHLIN ANALYSIS
To isolate the starlight behind the planet during its transit, we use the reloaded RM (RRM) technique. From here on we will use the term local CCF (CCF loc ) to refer to the occulted light emitted behind the planet and the term disk-integrated CCF (CCF DI ) to refer to the light emitted by the entire stellar disk. Here we will discuss the RRM technique at a high level, therefore, we refer the reader to Cegla et al. (2016a) for further details.

Isolating the Planet Occulted Starlight
Firstly, the CCFs DI from the ESPRESSO observations are shifted and re-binned in velocity space to correct for the Keplerian motions of the star induced by WASP-166 b, using the orbital properties in Table 2. Using the ephemeris and transit duration from Table 2, we separated the CCFs DI into those in-transit and those out-of-transit. We then created a single master-out CCF DI for each night by summing all out-of-transit CCFs DI and normalising the continuum to unity. By fitting the master-out CCF DI for each night with a Gaussian profile we can determine the systemic velocity, , as the centroid, see Table 1. We define the continuum as regions which lie outside the region of ± 10 kms −1 , a combination of the eq sin * and FWHM of the CCFs loc , from the stars rest velocity (i.e. has been subtracted). Before any further analysis, all CCFs DI were shifted to the stellar rest frame by subtracting for each corresponding night from the x-axis (i.e. there was no re-binning). As the ESPRESSO observations are not calibrated photometrically, we have to continuum scale each CCF DI to reflect the change in flux absorbed by the planetary disk. To do this, we normalise each CCF DI by their individual continuum value and then scale them using a quadratic limb darkened transit model using the fitted TESS parameters in Table 2. Finally, we can determine the CCFs loc which represent the CCF of the occulted starlight behind the planet by subtracting the now scaled in-transit CCFs DI from the master-out CCF DI , see Figure 3.

Determining and Modelling the Stellar Velocity of the Planet Occulted Starlight
The CCF loc (both binned and un-binned) in Figure 3 represent the starlight behind the planet along the transit chord, meaning they are both spectrally and spatially resolved. There is a change in continuum flux for CCFs loc which are near the limb of the stellar disk due to limb darkening and partial occultation of the stellar disk by the planet during ingress and egress. To determine the stellar velocity of the occulted starlight, we fit Gaussian profiles using curve_fit from the Python Scipy package (Virtanen et al. 2020) to each of the CCFs loc . There are a total of four Gaussian parameters in our fit including the offset (i.e. continuum), amplitude, centroid, and width. Initial values for these parameters for curve_fit were determined as continuum minus the minimum point of flux for the amplitude, the minimum flux point in velocity space between ± 10 kms −1 for the centroid, the continuum for the offset, and the variance of the data for the width. The continuum was determined as any points which lie outside ± 10 kms −1 and then lie between ± 0.5 of the mean. We applied constraints to both the centroid and offset during the fitting where the centroid was to be between ± 10 kms −1 (a combination of the sin and FWHM of the CCFs loc ) and the offset was allowed to vary between ± 0.5 of the the mean of points which lie outside ± 10 kms −1 . The flux errors assigned to each CCF loc were propagated from the errors on each CCF DI as determined from the version 2.2.8 of the ESPRESSO DRS. The local RVs of the planet occulted starlight are shown in Figure 4, plotted as a function of both phase and stellar disk position behind the planet. Cool main sequence stars show solar-like oscillations, known as pressure modes (p-modes), which are caused by sound waves travelling through the stellar interior. P-mode oscillations are present within the local RVs of WASP-166, therefore, we look to binning them out, see Section 5.2.We removed CCFs with < 0.25 from our analysis, as profiles close to the limb were very noisy and when comparing the depth to the noise the signal was not significant to enable a Gaussian fit, see Figure 3 where they are shown as dashed lines. This resulted in 3 CCFs being removed from the binned 10 minute data analysis and 11 from the un-binned observations.
We fit the local RVs in Figure 4 using the model by Cegla et al. (2016a) (hereafter referred to as the RRM model). This RRM model consists of one key formula (see equation 10 of Cegla et al. (2016b)) which depends on the position of the transiting planet centre with respect to the stellar disk, the sky-projected obliquity ( ), the stellar inclination ( * ), the equatorial rotational velocity ( ), the differential rotational shear ( ), quadratic stellar limb darkening (u 1 and u 2 ), and centre-to-limb convective variations of the star ( conv , see equation 11). We use the same coordinate system and angle conventions as used in Figure 3 of Cegla et al. (2016a). The stellar inclination is the angle from the line-of-sight towards the stellar spin axis (between 0 -180 • ) and is the sky-projected angle between the stellar spin axis and planetary orbital plane(between -180 -180 • ). For WASP-166, we use different scenarios for the RRM model depending on whether or not we account for differential rotation and centre-to-limb convective variations. For every ESPRESSO observation, we calculate the brightness weighted rotational velocity behind the planet following Cegla et al. (2016b) and using the quadratic limb darkening coefficients from the Table 2.
We sample the RRM model parameters using  Table 3 with uncertainties based on the 16th, 50th, and 84th percentiles.

RESULTS
A visual inspection of the local RVs in Figure 4 can provide an indication of the alignment of the WASP-166 system. Overall, the measured local RVs increase with orbital phase from negative to positive as the planet transits the stellar disk. Additionally, there is symmetry within the velocities with respect to phase, suggesting the system is in an aligned orbit with equal time spent crossing the blue and red-shifted regions. In two previous studies measurements for the projected obliquity, , and eq sin * were obtained using HARPS data by Hellier et al. (2019) from spectroscopic line broadening and an MCMC fit to the RM effect ( = 3 ± 5 • and eq sin * = 5.1 ± 0.3 kms −1 ) and RRM modelling by Kunovac & et al. (2022)  • and eq sin * = 4.82 +0.23 −0.21 kms −1 ). All of the results for our model fitting discussed in the following sections can be found in Table 3.

Solid Body Alone Model
To analyse the local RVs further, we first applied the RRM model assuming a solid body (SB) rotation for the star as it has the least number of free parameters. As a result, this leads to a degeneracy between eq and sin * , preventing us from calculating the stellar latitudes occulted by the planet. The two RRM model parameters are, therefore, the sky-projected obliquity, , and the projected rotational velocity, eq sin * . The probability distributions, from the MCMC model fitting, along with the 1D distributions are displayed as a  . The top panel shows the local RVs for the un-binned observations (grey) and for ten minute binned data (colour), determined from the local CCFs of the regions occulted by the planet as a function of phase. The ten minute binned data points are colour coded by the stellar disk position behind the planet in units of brightness weighted (where = cos ). The best-fit model for solid body rotation (SB) for the ten minute binned data (red line) is shown, along with the SB plus centre-to-limb linear (blue), SB plus centre-to-limb quadratic (green) and SB plus centre-to-limb cubic (orange) models. The bottom panel is the residuals (local RVs -model) for all models with colours corresponding to the top panel model lines, with a horizontal line at 0 to guide the eye. Points close to the limb can suffer from systematics (originating from the ephemeris, limb darkening etc.), therefore, we reran the model fitting removing the last point (Phase = 0.012) and found our results were within 1 with the BIC favouring the SB plus centre-to-limb cubic model. corner plot in Appendix D1. We find, for the combined fit to both runs eq sin * = 4.89 ± 0.08 kms −1 and = −4.49 +1.74 −1.73 • . These values are consistent with measurements obtained using HARPS data by Hellier et al. (2019) and Kunovac & et al. (2022). The derived projected obliquity informs us that the system is in an aligned orbit as suspected from a visual inspection of the local RVs. The SB model fit to the data is shown in Figure 4. We used an Anderson Darling test to check if the residuals in Figure 4 were normally distributed and found they were not indicating there is remaining correlated red noise, potentially as a result of magneto convection (i.e. granulation/super granulation) and or residual p-modes.
In addition to fitting both runs simultaneously, we also fit them individually to check their consistency. For run A we derive = −9.60 +2.36 −2.37 • and eq sin * = 4.89 ± 0.11 kms −1 and for run B = 1.21 +2.50 −2.51 • and eq sin * = 4.95 ± 0.11 kms −1 . Overall, eq sin * for both runs is consistent to better than 1 and the precision in is the same. For the combined fit, falls between the two values for both runs individually. The difference in between the two runs is primarily driven by outlier observations at each limb with larger uncertainties.
We run a SB model fit including an uncorrelated white noise term ( ) to check for any unaccounted for noise present in the data as a result of stellar surface variability, see Table 3. From the model, could potentially be picking up on p-modes oscillations which are present within the local RVs. As a result of this, we look to binning out these potential p-modes in Section 5.2.

Solid Body plus Centre-to-Limb Convective Variations Model
We are also interested in how the net convective blueshift (CB) varies across the stellar disk. The net convective velocity shift caused by granules changes as a function of limb angle (i.e. from the centre to the limb of the star) due to line-of-sight (LOS) changes. Since WASP-166 is an F-type main sequence star (T eff = 6050 K) we would also expect to observe the same phenomenon. To model the

Notes:
For all SB models * and are fixed under the assumption of rigid body rotation and the eq column corresponds to eq sin * . For these models we are unable to determine the 3D obliquity, . The BIC of each model was calculated using 2 , therefore, due to different sample sizes it is not possible to compare models between data sets. As a result, the reduced chi-squared ( 2 ) has been added to make these comparisons possible. For clarity, CLV1, CLV2 and CLV3 correspond to centre-to-limb linear, quadratic and cubic respectively. The first four corner plots for the un-binned observation MCMC runs are in an appendix and the remaining are available as supplementary material online.
centre-to-limb convective variations (CLV) we fit the local RVs for CLV and SB at the same time. This was done by using and eq sin * for SB rotation and adding a linear, quadratic or cubic polynomial as a function of limb angle, see Table 3 and Figure 5. According to the BIC of the models, for the un-binned data, the SB plus quadratic CLV is the best fit to the data, resulting from a better fit to the data at disk centre. There is a marginal difference between the SB alone and SB plus linear CLV, with the BIC telling us the difference in fit is not sufficient enough to justify the extra free parameters. Additionally, there is a marginal difference between the SB plus quadratic and the SB plus cubic, where the chi-squared ( 2 ) is slightly lower for the cubic fit but the BIC is higher indicating the extra free parameters are not justified. We would also like to point out that despite WASP-166 b being an aligned system (from the un-binned data), the local RVs at ingress and egress behave differently, see Figure 4. This is unexpected as we would expect to see the same behaviour as the planet transits the limb at ingress and egress. At present we do not have an explanation for this.
As mentioned previously we had potential p-mode oscillations in our local RVs. Radial velocity shifts caused by p-modes can be large enough to mask low mass, long and short period planetary signals. In our case, these p-modes are problematic as they introduce more noise and make it difficult to detect CLV. In Chaplin et al. (2019) they show how fine-tuning exposure times to certain stellar parameters can help to effectively average these oscillations. We used their Python code ChaplinFilter which outputs the ideal exposure times to use; (i) to suppress the total observed amplitude to a 0.1 ms −1 level (29.8 mins) and (ii) corresponding to the Earth-analogue reflex amplitude of the star (21.4 mins). However, our RV precision for the ESPRESSO data is on a 1 ms −1 level, therefore, to get to this we bin our RV errors to an exposure time of ∼10 mins. As a result of this, we binned the local RVs in time to increase our exposure time and increase our SNR, effectively averaging out the p-modes, by 10 mins. Additionally, we also binned the data to 20 mins which is between the ChaplinFilter results and our ten minute bins as anymore binning would reduce our spatial resolution significantly. We refit all SB models to the 10 and 20 min binned local RVs, finding the uncorrelated white noise term ( ) drops to zero, supporting our initial suspicions of p-modes in the un-binned observations.
For the 10 min binned data the best fit model to the data is the SB plus cubic CLV. This can be seen in Figure 5, where a cubic fit characterising the residual velocities (i.e. local RVs -SB model) is shown. The BIC for the quadratic CLV fit is the highest and, therefore, the worst model for the data. This can be seen clearly in Figure 5 where the quadratic fit does not capture the shape of the data, however, the linear fit does capture this better in comparison. For the 20 min binned data, the SB alone is the best fit to the data with the lowest BIC value, however, the SB plus CLV cubic fit has a lower reduced chi-squared but the BIC is higher indicating the addition of free parameters is not justified. This is because this level of binning reduces the sampling of observations at the limb making it difficult to pick out any structure of the CLV. All the SB models for 10 and 20 min binned observations produce a consistent eq sin * of approximately 4.80 kms −1 and within 2 of each other. When comparing all models from the un-binned observations, 10 min and 20 min binned observations, the reduced chi-squared ( 2 , see Table 3) indicates the SB plus cubic CLV model for the 10 min binned observations to be the best fit. All of the 1D distributions for each of the MCMC runs mentioned in this section can be found as supplementary material online. Therefore, the CLV is characterised by the cubic fit to have have a velocity of ∼ -1 to -2 kms −1 at the limb (see Figure 5). In Cegla et al. (2016b) they found for an aligned hot Jupiter with a four day orbit around a Sun-like star, if you neglect contributions of ∼5 kms −1 in CLV there can be uncertainties in the project obliquity of ∼ 10 -20 • . For WASP-166, when comparing between the SB model only and the SB plus cubic CLV of the binned 10 minute observations there is a difference of 9.6 • which is in line with what Cegla et al. (2016b) predicted.

Differential Rotation Models
In the next scenario we apply the RRM model to the un-binned local RVs assuming differential rotation (DR) for the star. If DR is present, we can determine the latitudes on the star which are transited by the planet through disentangling eq sin * and determine the 3D obliquity. Therefore, the MCMC parameters are , , eq and * and the results can be found in Appendix D5. Unfortunately, we have bimodal distributions present in both and * , informing us that there remains a degeneracy within them. This could be due to the spectroscopic transits not being precise enough to separate between if the star is pointing away or towards us. We also ran the MCMC fitting for DR fixing * < 90 • (away) and * > 90 • (towards) to get an estimate on * and finding these models had a higher BIC than the SB models. In Hellier et al. (2019) they estimate a stellar rotation period of 12.1 d from variations in the HARPS RVs and combining the eq sin * fitted to the RM effect with the fitted stellar radius (this was not done using photometry). From this we can disentangle * from eq sin * by using the rotation period and our eq sin * = 4.89 kms −1 , yielding * ∼73.0 • . We can then use this to predict the stellar latitudes transited by the planet, using P rot and eq sin * , which gives a ∼9 • change in latitude.
In Roguet-Kern et al. (2022) they investigate the optimal parameter space to use the RRM technique to detect DR and CLV on a HD 189733-like system (i.e. a hot Jupiter in a circular orbit around a K-dwarf). To do this they use simulations to explore all possible ranges of , * and impact factor, , producing maps of optimal regions. We place WASP-166 with = −4.49 • (from SB model, un-binned), * = 88.0 • (from Hellier et al. 2019, estimation) and = 0.213 (from Table 2) on the heat maps in their Figure 3. According to their maps we find the chance of detecting DR on WASP-166 is slim given these parameters with a change in BIC of <2. We ran the DR plus CLV linear, quadratic and cubic models on the 10 minute binned observations to see if accounting for CLV made a difference when trying to pull out DR, see Table 3. In all of the models there was a bimodal distribution present in * , suggesting we cannot discern if the star is pointing away or towards us. Furthermore, the BIC between the SB plus cubic CLV and DR plus cubic CLV has a difference of ∼1, which isn't a significant enough difference for the DR model to be preferred. Additionally, the differences in 2 and BIC are minimal between the away and towards models meaning we cannot discern between them, see Section C for full details. Therefore, the best model to fit the data is the SB plus cubic CLV.

Variations within the CCF profiles
In section 5.2 we modelled the local RVs from fitting the local CCFs. However, in this section we take a closer inspection of the CCFs themselves and how their shape changes as a function of limb angle. We computed an analysis into the equivalent width (EW) and FWHM (as they represent a measure of changes in the width and height of the profile) of our CCF profiles, using the 10 min binned observations (due to the increase in SNR) from the previous section. In Figure 6, we show the EW and FWHM plots for the 10 min binned data calculated from the Gaussian fitting method explained in Section 4.2. For EW we observe a trend which increases towards the limb and the same for FWHM. We fit both a linear and quadratic relationship to each and derive the R 2 and p-value as a measure of the goodness of fit. The R 2 measures the degree to which the data is explained by the model, where a higher value indicates a better fit. The p-value then indicates if there is enough evidence that the model explains the data better than a null model (i.e. if the p-value is very small, we can reject the null hypothesis). For the EW, R 2 = 0.04 and p-value = 0.18 for a linear fit as a function of and R 2 = 0.12 and p-value = 0.09 for a quadratic fit, where a quadratic is the better fit to the data. Similarly, for the FWHM, R 2 = 0.24 and p-value = 0.001 for a linear fit and R 2 = 0.25 and p-value = 0.004 for a quadratic fit. In this case, the p-value is lower for the linear model despite the R 2 ). The grey data points represent the 10 min binned observations, the red line is a two degree polynomial fit to the data, the green line is a linear fit to the data and the dashed black line is the mean of the data.
being lower in comparison to the quadratic model. Both models are a good fit to the data, however, we will select the model with the lower p-value which has fewer free parameters as we do not need to over complicate the model. Therefore, we argue that the linear model is the best fit to the FWHM data. In Beeck et al. (2013) they use 3D HD simulations to model the convection in a range of cool main sequence stars, then use these models to synthesise various Fe I line profiles. They found similar trends in EW and FWHM for their models of F-(T eff ∼ 6850 K) and G-type (T eff ∼ 5800 K) stars. They attribute the increase in EW towards the limb to be a result of the Fe I lines being stronger due to the increasing temperature of the lower photosphere with respect to optical depth (at the limb our line-of-sight means we see more of the granular walls). For FWHM, an increase towards the limb was observed for all stars in their models due to the increasing component of horizontal velocity flows compared to line-of-sight velocity causing a stronger Doppler broadening. For WASP-166, we observe the same behaviour as Beeck et al. (2013). They synthesise and analyse the Fe I line at 617.3 nm only, however, we are looking at the CCF which is an average of all absorption lines within the CCF template mask including a large number of Fe lines. In a more recent study, Dravins et al. (2018) investigate the centre to limb changes of line widths in HD 189733, a cool K-dwarf, finding the lack of any corresponding observed signature/trend. However, they go on to model the changes in line width for F, G and K-type stars as a function of limb angle (see their Fig 11) finding a change in the region of ∼3 kms −1 . In another study (Dravins et al. 2017), they compute a similar analysis on the hotter G-type star HD-209458 ( eff = 6071 K) finding a change in velocity of between 1 -2 kms −1 in line width towards the limb (see their Figure 17), which is comparable to our measurements in Figure 6. Overall, our finding using observations solidifies that of Beeck et al. (2013), Dravins et al. (2018) and Dravins et al. (2017) who used simulations.
In addition to looking at the EW and FWHM we can also isolate different parts of the CCF to fit individually. For this analysis we focus on fitting the CCF core and CCF wings as each originate from a different depth in the photosphere. We expect asymmetries within the CCFs as we move from disk centre to the limb caused by the net velocities of the granules on the surface of the star. We define the CCF wing to be 50% of the line depth including the continuum . The local RVs for the Gaussian fits to the CCF core (blue) and wing (orange) for the ten minute binned observations. and the CCF core to be the remaining 50% of the line. Due to the sampling of the CCF this was needed in order to have enough data points for a sufficient Gaussian profile to be fitted. The Gaussian fitting was done using the same method described in Section 4.2 with no constraint placed on the continuum of the line core fit. In Figure 7, we show the local RVs for the CCF core and CCF wing of the binned 10 minute observations where a split can be seen towards the limb of the stars between the two at the end of the transit. This shows there are asymmetries within the CCF with the wings appearing more redshifted towards the limb at the end of the transit than the core. We also fitted the local RVs of the wing and core using the MCMC method described in Section 4.2 and the results are in Table 3. The asymmetries in the line produced a rather high for the wing, whereas for the core remained consistent with an aligned system. We investigated the residual RVs (subtracting the local RVs from SB model) however, nothing of significance was found.

CONCLUSIONS
We have conducted an in depth study into the stellar surface variability of WASP-166 while also investigating the impacts this has on the projected obliquity. To do this, we utilised photometric transits of WASP-166 b, two observed with NGTS (simultaneous to ESPRESSO) and a further six observed with TESS, to fit for the transit planetary system properties (see Table 2). We then computed the Reloaded Rossiter McLaughlin (RRM) analysis using two full ESPRESSO transits to determine the local velocity of the star behind the planet, testing these RVs for various models to include solid body (SB) rotation, centre-to-limb convective variations (CLV), and differential rotation (DR), see Table 3. From a SB model fit of the un-binned ESPRESSO observations we found WASP-166 has a projected obliquity of = −4.49 +1.74 −1.73 • and eq sin * = 4.89 ± 0.08 kms −1 which is consistent with WASP-166 b being an aligned system. However, there was evidence of p-modes within our dataset so, we binned the data by 10 minutes to achieve an RV precision of ∼1 ms −1 (this is the expected precision from the exposure time calculator for ESPRESSO). This method effectively averages out the p-modes which can introduce more noise and make it difficult to detect CLV. Fitting the 10 minute binned observations with a SB RRM model yielded = −5.93 +1.99 −2.01 • and eq sin * = 4.77±0.09 kms −1 which remains consistent with our un-binned analysis and WASP-166 b being an aligned system. This is an update on previous studies of WASP-166 b (Hellier et al. 2019;Kunovac & et al. 2022) using HARPS observations where we are in agreement with their values, however, in this study the ESPRESSO observations allow us to be more precise.
In addition to fitting for SB we also looked at models with CLV. The net convective velocity from the centre of the star out to the limb changes due to the convective cells being viewed at different angles from changes in LOS. Therefore, we alter the RRM model to account for this by fitting a linear, quadratic and cubic relation as a function of limb angle. In the 10 minute binned dataset the SB plus cubic CLV model was the best fit to the data with an improvement in the BIC of twenty (see Table 3 and Figure 5 for the plotted trends). The ΔBIC (difference between a particular model and the best model with the lowest BIC) can be used as an argument against the other model. A ΔBIC between 2 and 6 means there is a good argument in favour of the best model and between 6 and 10 the evidence for the best model and against the weaker model is strong. We used an Anderson Darling test to check if the residuals for the un-binned observations in Figure 4 were normally distributed and found they were not indicating there is remaining correlated red noise as a result of magneto convection (i.e. granulation/super granulation). Therefore, we would like to caution the reader on the interpolation of the BIC of the models. For the 20 minute binned observations, the SB alone is the best fit to the data due to the level of binning reducing the sampling of observations at the limb, making it difficult to pick out any structure of the CLV. Overall, ignoring the velocity contributions induced by granulation can skew the measurements of the projected obliquity and, therefore, impact out understanding of planet formation and evolution of the system. This can be seen in our results within this paper where the projected obliquity changes by 9.6 • for 10 minute observations when CLV are accounted for. However, in the case of WASP-166 this change still means the system is aligned and as a result will not change our view of planet formation and evolution. For the un-binned observations = −4.49 +1.74 −1.73 and in 10 minute binned observations = −15.52 +2.85 −2.76 , eq sin * remains consistent within 2 of the SB model results. This is the first time where a tentative CLV has been detected on WASP-166 using the RRM technique.
We attempted to fit a DR model to the un-binned data, however, due to the number of stellar latitudes transited by WASP-166 b and our RV precision we were unable to constrain any DR for this data set. Additionally, we were unable to disentangle * from eq sin * as we could not determine a rotation period for WASP-166 from the available photometric data. From our predictions of the local RVs and keeping = −4.49 • , to pull out DR WASP-166 would need to posses an * between 110 • and 160 • or < 40 • (from Roguet-Kern et al. 2022), where more stellar latitudes would be crossed by the planet and the local RV variation induced by the DR would be greater than the uncertainties. However, we determine a CLV between ∼ -1 --2 kms −1 , therefore, by including this in our modelling we were able to fit and pull out a potential DR shear of ∼ -0.5 (anti-solar) and * either 42.03 +9.13 −9.60 • or 133.64 +8.42 −7.98 • (i.e. either pointing away from or towards us). This further strengthens the need to consider and fit for CLV as it can also impact the detection and characterisation of DR as it is may be masked out by the CLV. Finally, we also inspected how the shape of the CCF changes as a function of limb angle. We found that the FWHM and equivalent width (EW) of the local CCFs for the 10 minute binned data show trends which increase towards the limb of the star by ∼2 kms −1 . We modelled this with a linear and quadratic fit, using the R 2 value to determine the best fit was a quadratic for EW and linear for FWHM. This analysis matches the findings of both Beeck et al. (2013) and Dravins et al. (2017Dravins et al. ( , 2018 who found similar results in simulated line profiles from state-of-the-art 3D HD simulations. Therefore, these types of observations allow us to further validate these simulations for main sequence stars other than the Sun. Finally, we looked at the asymmetry in the local CCFs along the transit chord by fitting for the core and wing separately, finding the wing was more redshifted than the core. This has been seen for the Sun using HD simulations of the Fe I line in Cegla et al. (2018) (see Figure 11) where there is a C-shape within the line bisector which becomes more prominent towards the solar disk. Near the limb we lose the downward bend in the wing which is believed to be the contribution from the down flowing intergranular lanes.
Overall, for WASP-166 we have refined the system parameters and investigated the stellar surface variability, comparing this to and validating 3D HD simulations. We find a centre-to-limb convective variation detection which has an impact on deriving the projected obliquity. Furthermore, by including the centre-to-limb contribution and modelling it together with differential rotation were were able to put limits on the differential rotational shear and stellar inclination. acknowledges support by FONDECYT grant 1201371 and from the ANID BASAL projects ACE210002 and FB210003.

DATA AVAILABILITY
The TESS data are available from the NASA MAST portal and the ESO ESPRESSO data are public from the ESO data archive. NGTS data is available for the two simultaneous transits as supplementary material online with this paper. The third NGTS transit is available to download online from Bryant et al. (2020).