© ESO, 2017

1. Introduction

1.1. Analysis of the Rossiter-McLaughlin effect

The occultation of a rotating star by a planet removes the profile portion that is emitted by the hidden portion of the stellar photosphere from the apparent stellar line shape. This distortion of the spectral lines, known as the Rossiter-McLaughlin (RM) effect (Holt 1893; Rossiter 1924; McLaughlin 1924) traces the trajectory of the planet across the surface of the stellar disk and is thus sensitive to the velocity of the stellar photosphere and to the alignment between the spins of the planetary orbit and the stellar rotation (namely the obliquity). Measuring the spin-orbit alignment is particularly important because it can feed into theories on planetary migration and evolution. To date, the sky-projected obliquity has been measured for more than 90 planetary systems¹ (Albrecht et al. 2012; Crida & Batygin 2014). Most of these measurements have been carried out using Doppler spectroscopy, as the distortion of the stellar lines and their corresponding cross-correlation function (CCF) induces deviations to the Keplerian radial velocity (RV) of the star during the transit that can be fitted with analytical formulae (e.g., Ohta et al. 2005; Giménez 2006; Hirano et al. 2011; Boué et al. 2013). However, because it only fits this RV anomaly and not the full spectral CCF, this technique (hereafter referred to as the velocimetric analysis of the RM effect) can be prone to significant biases (e.g., Triaud et al. 2009; Hirano et al. 2010; Brown et al. 2017). For example, fitting the RV anomaly with analytical or numerical models that assume the local stellar profile is a constant Gaussian can potentially inject biases of up to 20−30° in the projected obliquity (Cegla et al. 2016b).

In recent years, Doppler tomography has been used to study the alignment and properties of a growing number of planetary systems (e.g., Collier Cameron et al. 2010; Gandolfi et al. 2012; Bourrier et al. 2015). This technique relies on the decomposition of the observed CCF profile into its different components (i.e., the stellar and instrumental profiles, the limb-darkened rotation profile, and the missing starlight signature caused by the planet). As such, it removes some of the biases of the velocimetric analysis and can provide additional constraints on the planet-to-star radius ratio and intrinsic stellar line profile. Yet even Doppler tomography assumes a constant, symmetric line profile across the stellar disk, and typically ignores the effects of differential rotation and velocity shifts due to stellar oscillations and granulation.

This motivated Cegla et al. (2016a) to develop a new technique to analyze and model the RM effect (the Rossiter-McLaughlin effect Reloaded technique, hereafter reloaded RM). By using the planet as a probe, the reloaded RM technique isolates the local CCFs from the regions successively occulted during the transit, with no particular assumptions about the shape of the intrinsic stellar photospheric lines. Many properties can be derived from these local CCFs, including the RV of the stellar surface along the transit chord. In turn, these RVs can be fitted with a simple model to derive the sky-projected obliquity and stellar rotational velocity. Furthermore, because it accounts for differential rotation, the reloaded RM technique can be sensitive to the stellar latitudes occulted by the planet. This can be used to lift the degeneracy between the equatorial velocity and stellar inclination, therefore yielding the value of the true 3D obliquity.

Cegla et al. (2016a) benchmarked the reloaded RM technique on the HD 189733b system and were able to constrain the differential rotation of the host star. However, HD 189733b transits its star at a relatively low latitude (~40°), along a nearly aligned orbit, and the local CCFs measured along the transit chord showed very similar properties. This motivated our search for a planet occulting more diverse regions of its host star, making it even more conducive to the application of the reloaded RM technique.

1.2. The WASP-8 system

Among all systems with known obliquities, we identified WASP-8 as one of the most suitable targets. It is a binary stellar system 87 pc away from Earth, where the A component is a bright G-type star (V = 9.8) hosting two exoplanets, WASP-8b and WASP-8c, and the B component is a faint M dwarf. The two stars have a common high proper motion (111 mas/yr) and their angular separation (4.495′′) sets a minimum separation of 440 au (Evans et al. 2016). WASP-8b was first detected by Queloz et al. (2010) through transit observations with the WASP-south telescope, and then was followed up in radial velocity with the Coralie spectrograph. Queloz et al. (2010) further reported a linear trend in this radial velocity follow-up, which was later revealed by new Keck HIRES measurements as the reflex motion from the outer companion WASP-8c (Knutson et al. 2014; P = 4339$\stackrel{\mathrm{+}\mathrm{850}}{{-390}_{}}$ days, M sin i = 9.5$\stackrel{\mathrm{+}\mathrm{2.7}}{{-1.1}_{}}$ M_jup assuming a circular orbit). WASP-8b is a Jupiter-size planet (1.04 R_jup; Queloz et al. 2010) with a mass of 2.24 M_jup (Knutson et al. 2014). It is on an eccentric (e = 0.3), 8.16 days orbit. Its spectroscopic transit was measured with the HARPS spectrograph (Queloz et al. 2010), revealing the radial velocity anomaly arising from the RM effect caused by a retrograde, highly misaligned orbit. The availability of archival HARPS transit spectra for WASP-8b, their high signal-to-noise ratio (S/N), and the large planet-to-star surface ratio (δ ~ 1.3%) makes this system a target of choice for the application of the reloaded RM technique. From its large impact parameter (b = 0.6) and large sky-projected obliquity (| λ | > 90°) we further expect to probe very different regions of the host star surface, as the transit chord roughly goes from the sky-projected spin axis at ingress to the stellar equator at egress. Finally, with an age of 4 Gy old and a G 6-type, WASP-8 is also a good candidate to search for differential rotation similar to our Sun.

In Sect. 2, we describe the spectroscopic transit dataset of WASP-8b and how it is processed to extract the CCFs of the local stellar photosphere. In Sect. 3 we present our reloaded RM model for the WASP-8 system, and the measurements derived for the system architecture and stellar surface properties. Section 4 discusses changes in the local CCF shape across the stellar disk and their interpretation based on simple star-planet simulations. In Sect. 5 we compare different techniques used to analyze the RM effect. We conclude and discuss the significance of our results in Sect. 6. Properties of the WASP-8 system that were used in our analysis as fixed parameters are given in Table 1.

2. Observations and data reduction

2.1. HARPS transit spectra

We used the transit observations of WASP-8b (Queloz et al. 2010) obtained with the HARPS echelle spectrograph mounted on the ESO 3.6 m telescope in La Silla, in Chile. 75 exposures of 300 s duration were secured during the single night of October 4, 2008, covering the full duration of the transit (47 exposures) and providing a good baseline for the unocculted host star both before (7 exposures) and after (20 exposures) the transit. During the measurement of the transit sequence a significant change in telescope focus occurred at orbital phase −0.0014, leading to a better accuracy and lower dispersion on the parameters derived from the reloaded RM analysis (e.g., Fig. 2). Observations were reduced with the version 3.5 of the HARPS Data Reduction Software, yielding spectra with resolution 115 000 covering the region between 380 nm and 690 nm. The reduced spectra were passed through an order-by-order cross-correlation with a G2-type mask function, weighted by the depth of the lines. Since the output CCFs are oversampled with a step of 0.25 km s^-1, for a pixel width of about 0.8 km s^-1, we kept one in four points in all CCFs prior to their analysis. The CCFs are defined in the solar system barycenter rest frame.

Fig. 1 Top panel: local CCFloc profiles (with colors, chosen for viewing ease) and residuals between the out-of-transit CCFDI and the master-out (gray). Bottom panel: residuals between the CCFloc profiles and their Gaussian fits. The dotted black line indicates the star rest velocity.

2.2. The reloaded RM technique

We used the reloaded RM technique developed by Cegla et al. (2016a) to isolate the starlight from the regions occulted by the planet during its transit. Hereafter we distinguish between the local CCFs (CCF_loc), which correspond to the light emitted by the planet-occulted regions, and the disk-integrated CCFs (CCF_DI), which corresponds to the light emitted by the entire stellar disk (minus the planet-occulted regions during the transit). First, the CCF_DI obtained from the HARPS spectra were shifted in velocity space, correcting for the Keplerian motion of the star that was induced by WASP-8b (using the orbital properties in Table 1). We do not account for WASP-8c, as it has a long period (4339 d; Knutson et al. 2014) and induces a negligible reflex motion on the host star during the ~9 h span of the observations (~8 cm s^-1). We then separated the CCF_DI secured outside of the transit and those within, using the ephemeris given in Table 1. The CCF_DI outside of the transit were co-added to build a single “master-out” CCF_DI, whose continuum was normalized to unity. The continuum is defined as regions farther than 10 km s^-1 from the star rest velocity, which corresponds to γ_⋆ in the barycentric frame. Fitting the master-out with a Gaussian profile, we found its centroid at γ_⋆ = −1537.5 ± 5.5 m s^-1. The first spectrum of the night was excluded from our analysis because the residuals between its out-of-transit CCF_DI and the master-out revealed a spurious signature, which was probably caused by twilight illumination. Using the value obtained for γ_⋆, all CCF_DI were then shifted to the stellar rest frame. Since the HARPS observations are not calibrated photometrically, each in-transit CCF_DI has to be continuum scaled to reflect the planetary disk absorption. To do this, we used the photometric light curve fitted to the WASP-south photometry by Queloz et al. (2010). Since the quadratic limb-darkening coefficients were not available in their paper, we retrieved them using the EXOFAST calculator² (Eastman et al. 2013) with Queloz et al. (2010) values for T_eff, log g, and [Fe/H]. Note that we used the Queloz et al. (2010)T_eff value for consistency, but a higher temperature was obtained by Evans et al. (2016). Our light curve, which was calculated with the batman package (Kreidberg 2015), matches well that in Queloz et al. (2010). Finally the local CCF_loc associated with the planet-occulted regions were retrieved by subtracting the scaled in-transit CCF_DI from the master-out (Fig. 1).

The CCF_loc profiles spectrally and spatially resolve the photosphere of the star along the transit chord. They have a lower continuum flux near the limb of the stellar disk because of the assumed limb darkening and the partial occultation from the planetary disk during ingress and egress (Fig. 1). We analyzed the shape of the CCF_loc and its possible variations across the stellar disk by fitting independent Gaussian profiles to each CCF_loc and deriving their RV centroid, contrast, and full width at half maximum (FWHM). The flux errors assigned to the CCF_loc were derived from the standard deviation in their continuum flux, and the uncertainties on the derived parameters correspond to the 1σ statistical errors from the Levenberg-Marquardt least-squares minimization. As can be seen in Fig. 1, the residuals between the CCF_loc and their fit yield a low dispersion and no particular features, showing that the assumption of a Gaussian profile for the CCF_loc is satisfactory. Hereafter we exclude from our analysis exposures for which the CCF_loc was either not detected or was detected with a contrast that is three times lower than the dispersion of the residuals in the continuum. This concerns the exposures at orbital phases lower than –0.01 (two at ingress) and higher than 0.0095 (three at egress; see the upper panel in Fig. 2).

Fig. 2 Upper panel: map of the CCFloc series, as a function of orbital phase (in abscissa) and radial velocity in the stellar rest frame (in ordinate). Colors indicate flux values. The vertical and horizontal dashed black lines indicate the mid-transit time and stellar rest velocity, respectively. The four vertical dotted black lines show the times of the contacts. The in-transit CCFloc correspond to the average stellar line profiles from the regions occulted by WASP-8b across the stellar disk. Middle panel: RVs of the stellar surface regions occulted by the planet (blue points). The green curve is the RV model obtained for the obliquity and rotational velocity in Queloz et al. (2010). The red curve corresponds to our best fit for the solid-body rotation case. Lower panel: residuals between the measured surface RV and our best-fit shown in the middle panel.

3. Analysis of the stellar surface velocity field

3.1. Method

The RVs of the planet-occulted regions are shown in Fig. 2. These RVs were interpreted using the model developed by Cegla et al. (2016a). This model consists of a single formula that depends on the position of the planet-occulted regions, the sky-projected obliquity (λ), the stellar limb darkening, the stellar inclination (i_∗), and the equatorial rotational velocity (v_eq), differential rotation, and convective velocities of the star. Our definitions for the coordinate system and angle conventions are the same as in Cegla et al. (2016a, see their Fig. 3). The stellar inclination (in 0–180°) is the angle counted positive from the line-of-sight toward the star spin axis, while λ (in −180−180°) is the angle counted positive from the star spin axis toward the orbital plane normal, and projected in the plane of sky. Different scenarios for the reloaded RM model are used in Sect. 3.2, depending on whether or not we account for differential rotation and convective effects. In each case, the theoretical RV of a planet-occulted region is calculated as the brightness-weighted average of the model RVs sampled over a square grid with resolution about 7 × 10^-3 the stellar radius (equal to R_p/31). We excluded contributions from cells that do not lie beneath the planet. Furthermore, because a transiting planet can move significantly during the duration of an exposure, the occulted regions should not be modeled by a simple disk. To account for this blur, we oversampled the chord between the planet position at the beginning and end of each exposure with a resolution R_p/20. The final theoretical RV, which is comparable to the measured RV, is the brightness-weighted average of all RVs from the planet grid at every oversampled position. We found, however, that the blur is negligible during the 300 s long exposures of WASP-8b.

We fit the reloaded RM model parameters using the Metropolis-Hasting Markov chain Monte Carlo (MCMC) algorithm detailed in Bourrier et al. (2015) (see also Tegmark et al. 2004; Ford 2006). We applied an adaptive principal component analysis so that step jumps take place in an uncorrelated space, allowing us to better sample the posterior distributions in the eventuality of non-linear correlations between parameters (Díaz et al. 2014). Jump parameters are described in the following sections depending on the scenario investigated. Except when specified, uniform priors were used on all parameters. We analyzed the system with multiple chains, started at random points in the parameter space, and with a jump size adjusted to get an acceptance rate of ~25%. We checked that all chains converge to the same solution, before thinning them using the maximum correlation length of all parameters. Finally, we merged the thinned chains so that the posterior distributions contain a sufficient number of independent samples. We set the best-fit values for the model parameters to the medians of the posterior probability distributions and evaluated their 1σ uncertainties by taking limits at 34.15% on either side of the median. We compared the results from the different reloaded RM models by calculating the Bayesian information criterion for the best fit (we fit 42 RV points).

3.2. Stellar rotation and obliquity results

Before the application of the reloaded RM model, information can be retrieved from a visual analysis of the local RV series (Fig. 2). The measured velocities decrease with orbital phase, going from positive to negative values, which shows that WASP-8b is on a retrograde orbit crossing first the redshifted regions of the stellar disk (moving away from the observer) and then its blueshifted regions. Furthermore, the strong asymmetry with respect to the phase and velocity origin indicates a significantly misaligned orbit with more time spent occulting the blueshifted regions.

3.2.1. Solid-body scenario

To further analyze the system, we first applied the reloaded RM model assuming a solid-body (SB) rotation for the star. In that case, the absence of differential rotation leads to a degeneracy between v_eq and sin i_∗, which prevents us from determining the exact stellar latitudes transited by the planet. The two jump parameters of the MCMC are thus the sky-projected obliquity λ and the projected rotational velocity v_eq sin i_∗. Their posterior probability distributions are shown in Fig. A.1, along with their marginalized 1D distributions that are well represented by Gaussians and allowed us to derive tight constraints on the inferred best-fit values (Table 2). The RV model obtained for these values provides a very good fit to the data (Fig. 2), and the corresponding system architecture is shown in Fig. 3. We note that the residual at phase near −0.005 is significantly larger than its uncertainty. It is possible that the planet occulted a stellar spot during this exposure. However the planet would have had to be grazing the spot to yield such an isolated variation in phase, and the amplitude of the variation seems too large to be caused by a spot (~800 m s^-1). Furthermore, the exposure does not show similar deviations to the contrast and FWHM time series (Fig. 4). Thus, the outlying RV value could rather be a statistical deviation and we did not exclude the corresponding exposure from our analysis.

Fig. 3 Upper panel: projection of WASP-8 in the plane of sky for the best fit in the differential-rotation scenario. The stellar disk is colored as a function of its RV field with thin black lines showing iso-RV curves equally spaced. The stellar spin axis is inclined and shown as a black arrow going from the south pole to the north pole (dashed when within the star). The normal to the orbital plane is shown as a green arrow, and the orbital trajectory is shown as a green curve. In the solid-body scenario, the spin axis orientation is degenerate and iso-RV curves would be parallel to the projected stellar spin. Lower panel: surface RVs of the regions occulted by the planet, as a function of orbital phase (increasing from right to left to match the retrograde motion of the planet). The measurements (blue points) are well reproduced in both the SB (solid red curve) and DR (dashed red curve) scenarios

Fig. 4 Contrast (upper panel) and FWHM (lower panel) of the local CCFloc from the planet-occulted regions as a function of orbital phase. A significant decrease in contrast occurs along the transit chord, while the FHWM remains stable overall. The solid horizontal lines correspond to the contrast and FWHM of the disk-integrated CCFDI measured outside of the planet transit. Vertical dotted lines indicate the planet contacts. The gray line in the upper panel shows our best fit for the local contrast, assuming it varies linearly from the pole to the equator of the star.

We then accounted for center-to-limb convective variations, including in the model a polynomial law in μ = cos(θ) (θ is the center-to-limb angle) for the convective velocities (CV). It is not possible to detect the constant component of the CV because it is degenerate with the true stellar rest velocity in the γ_⋆ measurement of the master-out CCF_DI, which was subtracted from the data. The model brightness-weighted CV integrated over the stellar disk must therefore equal zero, which implies that the zeroth order coefficient of the polynomial law (c₀) can only be determined by the higher order coefficients, which are added as jump parameters to the MCMC. We found that the difference in BIC does not favour the inclusion of a linear CV law in our model (Table 2), with its c₀ and c₁ coefficients not significantly different from zero. Nonetheless, the best-fit model corresponds to a decrease in CV from center-to-limb by about 200 m s^-1, which is of the same order as solar variations (on the 100 s of m s^-1 level; Dravins 1982). The inclusion of a quadratic CV law in the model yielded similar BIC value as the linear law (Table 2), and we thus conclude that the data are most consistent with a constant CV across the stellar disk.

3.2.2. Differential rotation scenario

In a second step we considered a model with differential rotation (DR), using a solar-like DR law Ω(θ_lat) = Ω_eq(1−α sin²(θ_lat)), where Ω is the rotation rate (equal to Ω_eq at the equator), α is the relative differential rotation rate, and θ_lat is the stellar latitude. Constraining the DR is one way to determine the exact regions of the star transited by the planet, and hence disentangle the true 3D spin-orbit geometry from projection effects. Indeed, in this case v_eq and sin(i_∗) can be treated as independent parameters, as they influence the model stellar surface velocity differently (see Eq. (8) in Cegla et al. 2016a). The jump parameters of the MCMC are thus α, λ, v_eq, and cos(i_∗). We did not allow α to go beyond 1 for physical reasons. We used cos(i_∗) with a uniform prior in [−1; 1] so that i_∗ can be bound directly to the range 0–180° (Cegla et al. 2016a). We also used sin(i_∗) as a prior to reflect the geometric probability of observing a given stellar inclination, but removing this prior does not significantly change the discussion below. The results of the MCMC, shown in Fig. A.2, yield asymmetric posterior probability distributions (PD) for all model parameters except λ. Therefore, in addition to the median and its associated 68.3% confidence interval (Sect. 3.1), we also report in Table 2 the highest density interval (HDI), which contains 68.3% of the distribution mass such that no point outside the interval has a higher density than any point within it. The HDI are equivalent to the range of values that encompass 68.3% of the distribution on each side of the median when the distribution is Gaussian.

We compared our results with the measurements of | α_Kep | obtained for Kepler G-type stars by Balona & Abedigamba (2016); because they use photometric methods, only the absolute value of α_Kep could be inferred. These authors derived an analytical formula to approximate | α_Kep |, which matches well their measurements for stars in the range 6000−6600 K (see their Eq. (2) and Fig. 11). Using this formula, we found that WASP-8 (T_eff = 6100 K; Ω_eq between 0.24 and 0.41 rad d^-1, as derived from our HDI for v_eq) would be expected to show relative differential rotation rates between 0.19 and 0.28. It is thus noteworthy that the PD of α and v_eq in our DR best-fit scenario (Fig. A.2) show a common peak at about 0.25 and 2 km s^-1, respectively.

Nonetheless a wide range of values for α and v_eq are allowed by the data. This is linked to the high stellar inclination obtained in the DR scenario, as one can see in Fig. 3. The combination of the stellar north pole pointing away from the observer and a large sky-projected obliquity means the planet always occults regions that are either close to the projected stellar equator (where differential rotation is lowest) or are near the stellar spin axis (where RVs are low, and were measured with larger errors during ingress). As a result, RVs of the planet-occulted regions in the DR scenario are very similar to that of the SB scenario, even when α is large (see lower panel in Fig. 3). Accordingly, the PD of v_eqsini_⋆ derived from the equatorial velocity and stellar inclination in the DR scenario yields similar results as in the SB scenario (Table 2). This is also the case for λ because variations of the obliquity directly affect the stellar longitudes (and thus also the iso-RV curves) crossed by the planet, which cannot be compensated for by variations of the other parameters.

Even though the best fit for the DR scenario yields a slightly lower χ² than in the SB scenario, the increase in BIC penalizes its additional free parameters (Table 2). The same conclusion was reached using the Akaike information criterion. Therefore, we cannot conclude with certainty to the differential rotation of WASP-8 and we consider the solid-body model as the most plausible. The iso-RV curves deviate more significantly from the SB case near the limbs of the star (Fig. 3) and, thus, measurements of the surface RV at higher S/N during ingress and egress are crucial in confirming differential rotation.

4. Spatial variations of the stellar spectra

4.1. Analysis

We analyzed the variations in the shape of the local CCF_loc along the transit chord, by examining the FWHM and contrast of their Gaussian fits (Sect. 2.2). We found that the contrast decreases almost linearly during the transit, from about 75% at ingress down to about 40% at egress (Fig. 4). We checked that varying the fixed system properties (Table 1) by up to three times their uncertainty did not remove this trend. Performing the reloaded RM analysis with the master-out computed from CCF_DI taken either only before the transit or only after the transit, also did not remove the trend. Furthermore, these different master-out CCFs were found to be consistent between themselves (i.e., with differences in RV centroid, contrast, and FWHM lower than 1σ). The FWHM remains roughly constant while the planet is fully in front of the stellar disk, but increases by about 1 km s^-1 at ingress and then again at egress, although the noise on the local FWHM is too large to confirm these variations. Also, the fully in-transit CCF_loc have a lower FWHM than the master-out CCF_DI (7.253 ± 0.02 km s^-1; Fig. 4). This is likely due to the rotational broadening in the disk-integrated CCF.

4.2. Interpretation

Instrumental effects, even the partial defocusing of the telescope during the first half of the transit (Sect. 2), are unlikely to explain such a continuous and large variation in the measured contrast. The reloaded RM technique is sensitive to the light curve used to normalize the CFFs. In particular, a larger transit depth can decrease the contrast of the extracted CCF_loc. However, we found that to obtain a constant local contrast during the transit of WASP-8b, the light curve would have to be abnormally asymmetric. Indeed, by comparison with the Mandel & Agol light curve (Mandel & Agol 2002) used by Queloz et al. (2010), the required light curve would have to be more than 1% deeper during ingress and the first part of the transit and up to 0.5% shallower during the last part of the transit. Such variations are excluded by the photometric measurements obtained by Queloz et al. (2010). Furthermore, these variations spread over the whole transit would likely arise from wide bands of spots and plages that strongly vary in contrast with stellar latitude, which would be inconsistent with the assumption of a constant local contrast from the planet-occulted regions. Finally, even the required asymmetric light curve would not change significantly the RV time series, except during ingress and egress where the deviations to the Mandel & Agol light curve would be the largest. Thus, we find it more likely that the measured contrast variation traces variations in the depth of the lines emitted by the stellar photosphere that could arise, for example, from the effect of the stellar magnetic field. Using theoretical line profiles generated in the frame of a previous study (Cegla et al. 2016a, albeit with a G2V spectral type in the MHD simulations), we investigated whether changes in the magnetic field strength could yield the observed contrast variations. However, even for the Fe I-617.3 nm line that is particularly sensitive to the magnetic field, an increase in field strength from 20 to 500 Gauss only changed the contrast by a few percent, which is much lower than the measured 35%.

The high obliquity of the system suggests that WASP-8b crosses many stellar latitudes along the transit chord (Fig. 3). Following an empirical approach, we thus considered a simple model where the contrast of the CCF_loc decreases linearly from the poles to the equator of the star, and the stellar inclination is allowed to vary such that v_eqsini_⋆ = 1.9 km s^-1 (following the SB scenario; Table 2). We then simulated the transit of WASP-8b and calculated the contrast of the average local profiles from planet-occulted regions. We used χ² statistics to compare the measured and simulated contrast series, varying the polar and equatorial contrasts, and the stellar inclination, in our model. We also included in the fit the comparison between the contrasts of the out-of-transit model CCF_DI and the observed master-out CCF_DI (Fig. 4). The best fit was obtained for a polar contrast of 100%, an equatorial contrast of 52%, and a stellar inclination i_⋆ = 112°. The corresponding theoretical contrast along the transit chord is shown in Fig. 4. Interestingly, the inclination we derive falls within the range obtained from the independent analysis of the surface RV (Sect. 3.2.2). In that configuration, WASP-8b would cross the equator of the star at near the center of the transit. Incidentally, this is also the time when the local contrast becomes similar to the contrast of the master-out CCF_DI (0.5423 ± 9 × 10^-4, Fig. 4), which might indicate that the equatorial bands of the stars contribute the most to the shape of the disk-integrated CCF.

We note that a polar contrast of 100% is likely to be unphysical. However this value is a projection of our gradient model at the pole, and the variation in contrast effectively observed and fitted actually goes from about 75% to 50% during the first part of the transit. Also, our gradient model does not explain the continuous decrease of the measured contrast during the last part of the transit. To yield a linear decrease in contrast along the full transit chord the simulated star would have to point toward the observer, so that the planet always occult the same stellar hemisphere. Yet this is inconsistent with the stellar inclination derived from the RV analysis, which corresponds to the star pointing away from the observer. Furthermore, with the star pointing toward us we found that the out-of-transit contrast was overestimated because the star is then mostly visible through its polar regions, which have the deepest line contrast. In conclusion, the measured contrast can be explained during most of the transit if the star points away from the observer and its local contrast decreases linearly at high stellar latitudes. The situation is more complex at the end of the transit, when the planet occults the equatorial bands of the star where the local CCFs are both shallower and broader (Fig. 4). Further investigations, beyond the scope of this paper, will be required to investigate whether Zeeman splitting or photospheric temperature changes could explain these variations.

5. Pitfalls of the velocimetric and tomographic RM analysis

Neglecting chromatic effects, the occultation of the stellar disk by the transiting planet decreases the flux in every part of the disk-integrated CCF_DI by the same percentage. In terms of absolute flux, the decrease is thus smaller at the bottom of the stellar lines than in their continuum, resulting in an apparent (positive) bump in the in-transit CCF_DI that is maximized at the velocity of the planet-occulted region. When the planet transits the stellar hemisphere that is rotating toward the observer, the occulted regions are blueshifted and the bump is located in the blue wing of the CCF_DI. The fit to this distorted CCF_DI with a Gaussian profile thus yields a RV centroid that is redshifted with respect to the Keplerian RV of the star (blueshifted when the planet transits the stellar hemisphere rotating away from the observer). The amplitude of this anomalous shift depends on the amplitude of the bump, i.e., on both the planet-to-star radius ratio and the contrast and width of the intrinsic local CCF behind the planet.

Both the velocimetric and tomographic analysis of the RM effect assume a constant shape for the local CCF_loc across the stellar disk (e.g., Ohta et al. 2005; Collier Cameron et al. 2010), yet the case of WASP-8b show that the local contrast in particular can vary significantly along the transit chord (Sect. 4). An increase in the CCF_loc contrast increases the amplitude of the planet bump in the CCF_DI, resulting in a larger RV deviation to the Keplerian curve (and vice versa for a contrast decrease). In a velocimetric analysis of the RM effect, the planet thus appears to be occulting regions farther away from the projected stellar spin axis than it is in reality. In the case of WASP-8b, the planet transits the redshifted half of the stellar disk for a short time before occulting its blueshifted half (Sect. 3.2), which results in a small blueshifted RV anomaly at the beginning of the transit followed by a larger redshifted RV anomaly, as detected by Queloz et al. (2010). Yet, because the contrast of the CCF_loc decreases steadily during the transit (Fig. 4), we argue that the planet is closer to the spin axis at the beginning of the transit and farther away toward the end of the transit than derived from the velocimetric analysis of Queloz et al. (2010). If correct, this would have led Queloz et al. (2010) to underestimate the absolute value of the sky-projected obliquity, which is in agreement with our measurement of λ = −143.0${\stackrel{\mathrm{+}\mathrm{1.6}}{{-1.5}_{}}}^{\mathrm{◦}}$, larger than the Queloz et al. (2010) value (−123.0${\stackrel{\mathrm{+}\mathrm{3.4}}{{-4.4}_{}}}^{\mathrm{◦}}$) by more than 4σ. We also argue that the variation in local contrast biased their RM-derived estimation of v_eqsini_⋆ (1.59$\stackrel{\mathrm{+}\mathrm{0.08}}{{-0.09}_{}}$ km s^-1) by ~3σ. However we note that our value for v_eqsini_⋆ = 1.90 ± 0.05 km s^-1 is consistent with their spectroscopic value (2.0 ± 0.6 km s^-1), which is derived from the fit of the stellar Fe I lines in the HARPS spectra. In conclusion, variations in the shape of the local CCF_loc cannot be detected in the traditional velocimetric analysis of the RM effect, despite clearly leading to strong biases in the derived properties.

Since Doppler tomography compares the profiles of the CCF_DI, rather than their RV centroid, it should be less impacted than a traditional velocimetric analysis regarding the determination of the planet trajectory and v_eqsini_⋆. This is because the signature of the planet is directly included in the tomographic model of the stellar line fitted to the in-transit CCF_DI. Therefore, variations in the contrast and FWHM of the local CCF_loc should not significantly affect the localization of the planet signature. However the usual assumption of a constant shape for the CCF_loc would likely bias both the estimation of the planet-to-star radius ratio and the width of the intrinsic stellar plus instrumental profile.

Alternatively, the reloaded RM technique allows for a simple and efficient way to isolate the CCF_loc and assess its variations. The residuals between the CCF_loc of WASP-8 and their fits do not show significant dispersion or spurious features (Fig. 1), giving us confidence that the shape of the CCF_loc is well approximated by a Gaussian profile along the transit chord. As such, a variation in the CCF_loc contrast or FWHM does not change the RV centroid and, therefore, does not affect our retrieval of the stellar surface RVs.

6. Conclusion

Analysis of the RM effect during the transit of an exoplanet allows us to measure the dynamical properties of the stellar photosphere and to constrain the system architecture through the measurement of the spin-orbit alignment. While the velocimetric and tomographic analysis of the RM effect yield the sky-projected obliquity of a system, the reloaded RM method developed by Cegla et al. (2016a) can potentially constrain the differential rotation of the host star. In that case, breaking the v_eqsini_⋆ degeneracy also yields the measurement of the stellar inclination and true 3D obliquity. Even when the reloaded RM technique is unable to constrain differential rotation, it avoids some of the biases that can affect the more traditional techniques. Indeed, our application of the reloaded RM technique to the WASP-8 system revealed near-linear variations of up to 35% in the contrast of the local line profile along the transit chord. Those variations could not have been detected in the velocimetric analysis of the RM effect performed by Queloz et al. (2010) and led to a significant misestimation of the sky-projected obliquity and projected stellar rotational velocity of the system. Such biases in the measured obliquity, caused by changes in the local CCFs shape, were predicted by Cegla et al. (2016b). Since the reloaded RM approach relies on the analysis of the RV centroid of the local CCF_loc, it is not sensitive to variations in their FWHM and contrast as long as they keeps a symmetric shape. This is the case for WASP-8, for which the local line profile is well approximated by a Gaussian, and we thus refine the properties of the system to λ = −143.0${\stackrel{\mathrm{+}\mathrm{1.6}}{{-1.5}_{}}}^{\mathrm{◦}}$ and v_eqsini_⋆ = 1.90 ± 0.05 km s^-1. Our study does not change the conclusions from Queloz et al. (2010) that WASP-8b is on a strongly misaligned, retrograde orbit, nor does it change the subsequent interpretations in the literature.

We did not find significant convective velocity variations from the center to the limb of WASP-8. We also found no definitive evidence for differential rotation, although the posterior distribution for α in this scenario peaks at about 0.25, which is consistent with the measurements obtained by Balona & Abedigamba (2016) for Kepler stars similar to WASP-8. The constraints obtained on α partially lift the degeneracy between v_eq and i_∗, suggesting possible values in the range 1.8−3.1 km s^-1 and 101−158°. The distributions derived for i_∗ and λ would also yield the true 3D obliquity (see calculation in Cegla et al. 2016a) in the range 106−129°. Even though measurements obtained with higher accuracy during the ingress and egress of WASP-8b will be necessary to assess the differential rotation of its host star, we emphasize that the posterior distribution for the sky-projected obliquity is well defined and yields the same results whether differential rotation is accounted for or not.

In conclusion, we caution against using a velocimetric analysis to interpret the RM effect of transiting planets, as it is sensitive to variations in the shape of the local stellar CCF_loc that cannot be detected except through a direct study of the spectra. Doppler tomography should be less affected by such variations, although they could still lead to biases in the estimation of the planet-to-star radius ratio and the width of the stellar line profile. It is possible that the properties of many planetary systems analyzed in the literature values have been biased by variations in the shape of the local CCF_loc. The reloaded RM-derived stellar surface velocities can still be biased by variations in the symmetry of the stellar line profile or deviations to the light curve used to normalize the CCF_DI (caused, e.g., by spots). Nonetheless, this method allows for a simple and efficient way to isolate the local stellar CCF_loc, assess their variability, and correclty interpret their RV centroid even in case of variations in their contrast or FWHM. Furthermore, because the reloaded RM technique extracts the average stellar properties behind the occulting planet surface, small Earth-size planets that will be detected by missions like CHEOPS, TESS, and PLATO will act as finer probes of the stellar surface in terms of spatial resolution. Although the lower flux retrieved from behind the surface of small planets will make it more challenging, analyzing the RM effect with this technique will thus have a strong potential when coupled with observations from future generations spectrographs like ESPRESSO, SPIROU, or NIRPS.

Acknowledgments

We warmly thank the referee for their appreciative reading of our paper. This work has been carried out in the frame of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). V.B. and A.W. acknowledges the financial support of the SNSF.

References

Appendix A: Additional posterior probability distributions

Fig. A.1 Correlation diagrams for the probability distributions of the solid-body rotation model parameters. Green and blue lines show the 1σ and 2σ simultaneous 2D confidence regions that contain 39.3% and 86.5% of the accepted steps, respectively. One-dimensional histograms correspond to the distributions projected on the space of each line parameter. The red line and white point show median values.

Fig. A.2 Correlation diagrams for the probability distributions of the differential rotation model parameters. Green and blue lines show the 1σ and 2σ simultaneous 2D confidence regions that contain 39.3% and 86.5% of the accepted steps, respectively. One-dimensional histograms correspond to the distributions projected on the space of each line parameter. The red line and white point show median values.

All Tables

All Figures

	Fig. 1 Top panel: local CCFloc profiles (with colors, chosen for viewing ease) and residuals between the out-of-transit CCFDI and the master-out (gray). Bottom panel: residuals between the CCFloc profiles and their Gaussian fits. The dotted black line indicates the star rest velocity.
In the text

Fig. 2 Upper panel: map of the CCFloc series, as a function of orbital phase (in abscissa) and radial velocity in the stellar rest frame (in ordinate). Colors indicate flux values. The vertical and horizontal dashed black lines indicate the mid-transit time and stellar rest velocity, respectively. The four vertical dotted black lines show the times of the contacts. The in-transit CCFloc correspond to the average stellar line profiles from the regions occulted by WASP-8b across the stellar disk. Middle panel: RVs of the stellar surface regions occulted by the planet (blue points). The green curve is the RV model obtained for the obliquity and rotational velocity in Queloz et al. (2010). The red curve corresponds to our best fit for the solid-body rotation case. Lower panel: residuals between the measured surface RV and our best-fit shown in the middle panel.

In the text

Fig. 3 Upper panel: projection of WASP-8 in the plane of sky for the best fit in the differential-rotation scenario. The stellar disk is colored as a function of its RV field with thin black lines showing iso-RV curves equally spaced. The stellar spin axis is inclined and shown as a black arrow going from the south pole to the north pole (dashed when within the star). The normal to the orbital plane is shown as a green arrow, and the orbital trajectory is shown as a green curve. In the solid-body scenario, the spin axis orientation is degenerate and iso-RV curves would be parallel to the projected stellar spin. Lower panel: surface RVs of the regions occulted by the planet, as a function of orbital phase (increasing from right to left to match the retrograde motion of the planet). The measurements (blue points) are well reproduced in both the SB (solid red curve) and DR (dashed red curve) scenarios

In the text

Fig. 4 Contrast (upper panel) and FWHM (lower panel) of the local CCFloc from the planet-occulted regions as a function of orbital phase. A significant decrease in contrast occurs along the transit chord, while the FHWM remains stable overall. The solid horizontal lines correspond to the contrast and FWHM of the disk-integrated CCFDI measured outside of the planet transit. Vertical dotted lines indicate the planet contacts. The gray line in the upper panel shows our best fit for the local contrast, assuming it varies linearly from the pole to the equator of the star.

In the text

Fig. A.1 Correlation diagrams for the probability distributions of the solid-body rotation model parameters. Green and blue lines show the 1σ and 2σ simultaneous 2D confidence regions that contain 39.3% and 86.5% of the accepted steps, respectively. One-dimensional histograms correspond to the distributions projected on the space of each line parameter. The red line and white point show median values.

In the text

Fig. A.2 Correlation diagrams for the probability distributions of the differential rotation model parameters. Green and blue lines show the 1σ and 2σ simultaneous 2D confidence regions that contain 39.3% and 86.5% of the accepted steps, respectively. One-dimensional histograms correspond to the distributions projected on the space of each line parameter. The red line and white point show median values.

In the text