© ESO 2022

1 Introduction

During a planetary transit, part of the stellar surface is eclipsed by the orbiting planet. This causes a distortion of the spectral lines from the disc-integrated star and a change in their associated radial velocities (RV) that is known as the Rossiter– McLaughlin (RM) effect. Because the RVs partly depend on the stellar rotation, the RM waveform is sensitive to the star–planet alignment, which provides information on the dynamical history of the system. For example, planets in aligned orbits may be indicative of a dynamically gentle planet-disc migration history, whereas misaligned planets may have experienced a more violent migration, such as planet-planet or star–planet scattering via the Kozai–Lidov mechanism (Triaud 2018, and references therein).

If a star experiences differential rotation (DR) between the poles and the equator and/or convection in the stellar photosphere, these variations are naturally encoded in the line-of-sight velocities and are therefore also encoded in the RM signal. For example, our Sun rotates about 20% faster at the equator than at the poles, which equates to a local change in velocity of ~400 m s^–1 across the full range of latitudes. Additionally, the convective envelope of our Sun means that the surface is corrugated: at disc centre, we see the very tops of the granules and the very bottoms of the intergranular lanes, but towards the limb, the tops and bottoms are obstructed from view and we begin to see the granular walls. Furthermore, some granules in the forefront of our line of sight obscure those in the background, and so on. The result are centre-to-limb variations (CLV) in a range of observed properties from continuum intensity to stellar line shape and net position or shift, with net velocity shifts of several hundreds of m s^–1 (see Löhner-Böttcher et al. 2018, and references therein, for more details). As spectrograph precision increases, ignoring these more subtle effects (as is often the case when modelling the RM effect) can inject systematic biases in our star–planet interpretations (Czesla et al. 2015; Cegla et al. 2016b; Casasayas-Barris et al. 2020; Chen et al. 2020). Alternatively, if we account for these effects in RM modelling, we can overcome these biases and learn more about the physics on the stellar surface. We therefore explore the optimal star–planet parameter space here in which transiting planets, in conjunction with the reloaded Rossiter–McLaughlin technique (RRM) (Cegla et al. 2016a), can be used to probe both stellar surface differential rotation and centre-to-limb (net) velocity variations induced by stellar surface convection.

The RMM measures the spatially resolved stellar spectrum behind a planet as it transits, making it ideal for probing variations induced by the stellar surface. Using this approach, we can not only easily account for DR and limb-dependent convective effects, but also do not need to make any assumption about the local line shapes. This latter aspect is thanks to the direct subtraction of the in-transit data from out-of-transit observations (after the spectra have been normalised to a known transit light curve to compensate for the loss of absolute flux due to the Earth’s atmosphere; see Cegla et al. 2016a for more details). Nonetheless, the orientation of the star–planet system also dictates which areas of the star are probed during the planet transit. Hence, some star– planet orientations are more or less favourable to detecting these subtle stellar surface effects.

Naively, it is often assumed that systems in which the planet has a projected obliquity of 90° are the most favourable for detecting stellar DR because the planet has a higher likelihood of transiting more stellar latitudes; there is some evidence to support this idea, as this was indeed the finding of Serrano et al. (2020). However, Serrano et al. (2020) explored a relatively sparse parameter space (with the classical, disc-integrated RM modelling approach), simulating only four projected obliquities, three stellar inclinations, and five impact parameters. We here simulate a representative transiting hot Jupiter (similar to HD 189733 b) at all possible ranges of projected obliquity, stellar inclination, and impact parameter through a RMM analysis. Our goal is to determine the optimal parameter space for detecting and disentangling a solar-like stellar DR and CLV in the net convection-induced line-of-sight velocities.

The foremost goal of our study is to determine in a first instance the star–planet system parameters that allow the easiest and most accurate detection of stellar DR. Additionally, we aim to determine the optimal parameter space in which limb-dependent variations in the net convective velocities may be isolated from the DR.

In Sect. 2, we present the model with which we simulate (and recover) the RRM signal from a transiting planet, including the introduction of stellar surface DR, limb-dependent net convective signal, and photon noise. Our attempt to probe the convection effects is explained further in Sect. 2.2. We discuss our findings in Sect. 3 and conclude in Sect. 4.

2 Simulation and model

To accomplish our goal, we simulated a star–planet system during transit as observed from the Earth. In the optimal (empirical) cases, the target star is bright and the transit is well constrained by an existing light curve. For simplicity, we also assumed a circular orbit. We simulated a typical hot Jupiter because the high planet-to-star ratio of hot Jupiters boosts the signal-to-noise ratio (S/N) for the RRM analysis and their short orbital periods facilitate obtaining multiple transits. Additionally, they are also easier to observe within a single night, allowing for a generous out-of-transit baseline.

We focused on an HD 189733-like system, as it is one of the brightest transiting hot Jupiter systems known to date and the first application of the RRM, wherein Cegla et al. (2016a) were able to loosely constrain the DR and CLV; α was found to be >0.1 and >0.2 with 99.2% and 91.7% confidence, respectively. Additionally, while the observations agreed with predictions from 3D MHD simulations, the observed local cross-correlation functions (CCF) did not exhibit significant convective variations. The following sky-projected obliquity and the stellar inclination were also recovered: λ ≈ –0.4 ± 0.2° and $i_{*}\approx {92}_{-4}^{+12\circ }$ (and a true 3D obliquity $\psi {\approx }_{-4}^{+12\circ }$

The three main parameters that we considered in our simulations are the stellar inclination i_*, the projected obliquity λ, and the orbital inclination i_p (equivalent to the impact parameter b for a fixed stellar radius and semi-major axis). Because of the way our model is set up, the planet always transits in front of the bottom half of the star from the perspective of an observer on the line of sight. Furthermore, we also studied the influence that the projected equatorial rotational velocity U_eq sin i_* and the visual magnitude V have on our observations, as the RRM signal is amplified by a larger u_eq sin i_* and because the visual magnitude V controls the S/N. Lastly, we briefly explore the impact of a difference in the cadence of observations on the results.

As we follow the same coordinate system as Cegla et al. (2016a), we varied i_* from 0° to 180° and λ from –180° to 180° in steps of 2^° and 4^°, respectively. Between subsequent runs, we also varied the orbital inclination such that the impact parameter varied between 0 to 1 in steps of 0.05. Table 1 lists the system parameters that we held fixed for the HD 189733-like system in our simulations, including the limb-darkening coefficients u₁ and u₂.

Our fitting loop operates as follows: for every orientation on a grid defined by λ, i_*, and b and using the fixed parameters described in Table 1, the stellar velocity behind the planet at each point of the transit was calculated assuming a brightness-weighted solar DR law, following Cegla et al. (2016a). These data served as our reference point, that is, they correspond to the true RVs behind the planet.

To generate realistic velocity uncertainties for the local, occulted regions along the transit chord, we devised a model derived from the relation between surface RVs and the corresponding planet-occulted CCFs (independently of specific observations). The scaling of this formula is based on the local velocities determined from empirical observations of the HD 189733 b transits analysed in Cegla et al. (2016a). The precision of the local RV is dictated by the amount of signal behind the planet, by instrumental sensitivity, the brightness of the star, and the exposure time. The formula is also scaled by a constant offset C, associated with instrumental setup and environmental conditions (that we assumed constant). Hence, our local RV uncertainties e(RV, µ) are defined as $$e\left(\text{RV,μ}\right)=\sqrt{\frac{1}{\text{gain}}}\cdot {\left(\frac{1}{R_{*}}\right)}^2\times {10}^{\frac{V}{25}}\cdot \frac{1}{t_{\exp }}\cdot \frac{LC\left(\text{μ}\right)}{1-LC\left(\text{μ}\right)}\cdot C,$$(1)

where the gain was set to 1 for HARPS and to 6 for ESPRESSO (equivalent to flux’/flux, and it assumes the same exposure time), R_* is the stellar radius, t_exp is the exposure time in seconds, and the light curve (LC) is a function of µ = cos θ, with θ the centre-to-limb angle.

To simulate photon noise, we added Gaussian noise to the local velocities with a standard deviation equal to the value of the uncertainties, scaled from the empirical HARPS observations using Eq. (1). For reference, the local velocity uncertainties in the HD 189733 HARPS observations were close to 50 m s^–1 near disc centre and a few hundred m s^–1 near the limb (Cegla et al. 2016a). To restrain the scope of the study, we did not consider any other sources of noise (instrumental or stellar) in our model. The impact of pressure-mode oscillations for Sun-like stars are likely to be largely mitigated with the optimal exposure time and/or binning (Chaplin et al. 2019; Kunovac Hodžić et al. 2021); however, the effects of temporal convection evolution (i.e. granulation) is likely to increase the difficulty of recovering DR and CLV and will be examined in a forthcoming study. Throughout this analysis, we injected a solar-like DR, with a differential rotation rate α = 0.2. An example of the simulated observations for a given orientation is shown in Fig. 1 by the black markers.

2.1 Considering only differential rotation

In the first part of our study, we restricted our analysis to the sole contribution of DR (without CLV). The simulated data were fit using two different models: one model assumed DR, and the other assumed solid-body rotation (SB). In both models, we held the impact parameter b fixed as we recognised that allowing this parameter to float freely may lead to degeneracy between fitting parameters. Throughout this paper, we place ourselves in idealised cases to better understand the detectability of CLV and DR with system properties (i.e. we assume the brightest systems, with the best light curves available). The DR model has four free parameters: the projected obliquity Å, the stellar inclination i_*, the equatorial velocity u_eq, and the differential rotation coefficient α. When we fit the SB model, we had two free parameters: λ and v_eqsin i_*. In this case, the parameters i_* and α were fixed to 90° and 0, respectively (such that v_eq = v_eqsin i_*). In both models, we used non-linear least squares to fit the data. The fitting algorithm needs a set of starting parameter values, one for each free parameter, which we provided as a random value within a certain range centred on the true value (typically half the possible range for a given parameter). The fitting process allows each free parameter to take any value within the bounds of the parameter space, independently of the value of any other free parameter. Sections 3.1.1 and 3.2.1 provide more detail about the accuracy of the fits.

To compare the DR and SB fits, we first calculated the X² for each. Computing the Bayesian information criterion (BIC) for each fit allowed us to prefer the model with the lower of the two values whilst being conscious of overfitting the data with too many free parameters. A difference ∆BIC >2 is already considered positive, while ∆BIC >6 indicates a strong preference for the smaller BIC model (Kass & Raftery 1995). Hence, in the graphs presented in this paper, any difference in the BIC value higher than 10 is rounded down to that number to facilitate interpretation. Figure 1 shows an example of an RV curve, with the simulated data and the two fits using our DR and SB models.

2.2 Considering a centre-to-limb variation in convection

The second goal of our study is to investigate the impact of CLV on the detection of DR and to determine when these two effects may be differentiated; we use CLV hereafter to refer to the net limb-dependent convection-induced radial velocity shifts. To do this, we injected limb-dependent convective velocities into our previously simulated observations and tried to fit a variation of our two models on this adjusted input.

The injected convective velocities were derived from the magnetically quiet (net 0 G) solar simulations presented in Cegla et al. (2018). Therein, the authors constructed 3D magnetohydro-dynamic (MHD) simulations and synthesised the Fe I 6302 Å absorption line profile. The MHD simulations were inclined at nine centre-to-limb positions in steps of 0.1 µ and corresponded to six snapshots separated in time by 20 minutes each. To emulate exoplanet observations that use CCFs, the limb-dependent velocities were determined by cross-correlating all line profiles to a single profile at disc centre and taking the centroid of these resulting CCFs. Hence, these velocity shifts encapsulate all centre-to-limb variations resulting from the 3D stellar atmosphere; this includes both net shifts from which velocity flows are visible along the line of sight and shifts arising from changes in the line profile shape or asymmetry. As constructed, the velocities are relative to disc centre and thus do not encode the net convective blueshift present at disc centre; this is by design as the disc-integrated net convective blueshift is subtracted off by construction in the RRM (when the systemic and nightly radial velocity offsets are removed). We averaged these velocities over the time series to isolate the net changes across the limb, and fit a fourth-order polynomial to extrapolate to any limb position. The convective velocities were then simply added to the simulated stellar rotational velocities from Sect. 2, denoted hereafter as either DR + CLV or SB + CLV.

Since the MHD simulations are based on a solar twin and only one line profile was synthesised (albeit a typical line, representative of many lines that dominate most CCF template masks), the injected convection model might not be entirely representative of the velocities from a large spectral bandpass or for other stellar types. Moreover, because we used a polynomial fit to a time-average of snapshots, we did not account for temporal variations here. Nonetheless, this investigation provides a starting ground for estimating our ability to recover a limb-dependant signal.

To recover the CLV in the simulated observations, we considered two more free parameters (in addition to those described in Sect. 2): the coefficients c₁ and c₂ of a second-order polynomial that depends on the centre-to-limb angle; the polynomial coefficient c₀ was derived mathematically from c₁ and c₂ (see Cegla et al. 2016a for details). We dropped the third and fourth coefficients in the fitting procedure as we found that with the level of injected photon (Gaussian) noise, we would be overfit-ting the data (this is consistent with empirically reloaded RM analysis in the literature); this simplification also increased the computational speed.

In Fig. 2, we show an example of one star–planet orientation where we included the convective contribution in the reference data. All other parameters being equal, the shape of the curve is clearly radically different from the plot in Fig. 1, where no convective effects are considered. The largest differences can be noted at the limbs, where the convection contributes towards a relative redshift; although the net convective effect at disc centre is a blueshift of a few hundred m s^–1, which dramatically decreases towards the limb, this is a relative effect that is largely subtracted off when we remove the systemic velocities in the empirical observations.

It is important to note that fitting with CLV increases the number of free parameters, which is heavily penalised in the BIC value. We recall that a positive change in the BIC indicates a preference for the more complex model, and because a ∆BIC>2 is already significant, we do not give precise values above ∆BIC = 10 to facilitate visual interpretation of the graphs.

Fig. 1 Simulated local stellar surface velocity below the planet along the transit chord as a function of the phase (bottom) and stellar radius (top) for a given star–planet orientation (λ = 0°; i* = 150°; b = 0.2; V = 10). The best-fit DR model is shown in blue, and the best-fit SB model is plotted in orange; a positive BIC indicates a preference for the more complex DR model.

Fig. 2 Simulated local stellar surface velocity below the planet along the transit chord as a function of the phase (bottom) and stellar radius (top) for a given star–planet orientation (λ = 0°; i* = 150°; b = 0.2; V = 10), wherein both a solar DR and CLV are injected; compare to Fig. 1 to note the CLV impact (overplotted here in light colours), particularly the relative redshifts near the limbs. The corresponding DR + CLV model fit is shown in blue, and the CLV + SB fit is plotted in orange; a positive BIC indicates a preference for the more complex DR + CLV model.

3 Results and discussion

In the following section, we present our findings for the optimal parameter space for detecting and disentangling DR and CLV using the RRM technique. First we present our results when stellar rotation alone was considered. Then we present our analysis when a solar-like DR and CLV are considered. To help determine the optimal parameter space for which a DR fit is preferred to SB (or DR + CLV versus SB + CLV), we display our results in two and three dimensions, using colour-coded heat maps.

The 2D heat maps are retrieved by holding the impact parameter b constant while we varied the stellar inclination i_* (on the vertical axis) and the stellar projected obliquity λ (on the horizontal axis). They represent the model preference over all star–planet orientations considered for a particular value of b and are colour-coded based on the BIC difference between the DR and the SB fit models; a darker shade of green (i.e. more positive ∆BIC) indicates a significant preference of the DR fit (or DR + CLV) over the SB fit (or SB + CLV). On the other hand, light green and white indicate that DR (or DR + CLV) cannot be preferred over SB (or SB + CLV) with confidence. All negative ∆BIC returns in which SB is preferred were set to zero in the plots to more easily identify the regions in which a DR detection is possible. As a reminder, ∆BIC values higher than 10 were set to 10 to facilitate the presentation of the results in the figures.

Comparison between the various 2D heat maps gives a clear representation of the potential of retrieving DR (or DR + CLV) as the impact parameter increases. To provide a convenient way to compare the results for multiple values of b, we generated a 3D representation by stacking the 2D plots over the whole range of impact parameters. To facilitate the interpretation, we projected the 3D structure that emerges onto the corresponding planes. The label “fill” associated with the 3D plots provides a way to quickly appreciate the fraction of the parameter space below the top surface that holds orientations favourable to the detection of DR. Technically, because the parameter space is divided into cells for which we know the preference for DR over SB, it is trivial to average this value over a whole stack of equivalent orientations (varying only the impact parameter). A lighter colour therefore corresponds to the presence of holes in the stack below the top surface, which is defined as the highest impact parameter for which the DR fit is preferred over the SB fit. The colour scale of the projected contours on the two vertical planes corresponds to the fraction of the orthogonal stack in the parameter space for which DR is preferred.

3.1 Stellar rotation alone

We compared DR and SB alone (i.e. without CLV). We recall that the DR model has four free parameters (λ, i_*, v_eq, and α) and the SB model has only two (λ and v_eq sin i_*). Initially, we set the injected v_eq sin i_* and V magnitude according to the values listed in Table 1. We plot a series of 2D heat maps at increasing values of b to probe the DR preference in the parameter space.

As Fig. 3 demonstrates, for low values of b, most of the parameter space allows a clear detection of DR, except in a few areas in which i_* = 90°. Most of these regions of poor DR detection correspond to stellar orientations in which the transiting planet crosses the fewest stellar latitudes. We also note two white streaks in the left side plot at λ ≈ ±90°, which is a sign of a small degeneracy. At λ = ±90°, the planet occults only one longitude value, and at b = 0, this means that the planet crosses a region of zero stellar rotation; hence, at low impact parameters and perpendicular projected alignments, the occulted stellar rotation can fall below the measurement precision. It is also important to note that for low-impact parameters, there can be degeneracies in the observed transit light-curve fit, such as between a/R_⋆, i_p and T 0; hence, a more precise light curve may be needed to retrieve the results empirically. Using ESPRESSO and its higher precision, we note that the streaks of Fig. 3 narrow down and even fill in around i_* ≈ 90° (Fig. 7 and Fig. 8 are useful to visualize the difference in sensitivity between a HARPS and an ESPRESSO setup, with the latter drastically improving detection compared to the former).

The picture changes, however, as the impact parameter increases, and our significant preference for DR decreases rapidly. At high values of b, DR can only be reliably detected in distinct areas of the parameter space, listed below:

• i_* < 90° and 120° < |λ| < 180°

• i_* > 90° and |λ| < 60°

• i_* ≈ 90° and λ ≈ ±90°.

Around these hot-spots (i.e. areas where DR detectability is high), the DR model is preferred over SB up to high values of the impact parameter b. The hot-spots likely arise because they correspond to configurations in the system in which the planet crosses a maximum of latitudes on the star, that is, regions in which the effects of the differential rotation are heightened. When we map the range of latitudes that the planet crosses during the transit, we indeed observe a maximum around λ = ±90^° and i_* = 90^°, regardless of the impact parameter. Furthermore, as the planet always transits in front of the lower half of the star in our model (due to our coordinate system), regardless of the stellar orientation, the range of latitudes crossed is higher when the star is farther away from the pole (i_* > 90°) than when it is pole-on (i_* < 90°), which explains the vertical asymmetry in Fig. 3. This effect is more prevalent at lower impact parameters. At higher impact parameters, the majority of the transits occurs near the limbs of the star, where the S/N is worse. Another major factor is the (assumed) instrumental precision, as an ESPRESSO/VLT-like setup is much better at picking up DR at high values of b than a setup with HARPS-like precision (cf. Fig. 7 and 8). Other factors that affect our ability to retrieve DR at a high impact parameter are the absolute magnitude of the rotational shear α, the projected stellar rotational velocity v_eq sin i_*, and the sampling rate of the observations, for example, as we discuss in Sects. 3.1.2 and 3.1.4.

To give a point of comparison in the subsequent discussion, we present the 3D parameter space for a HD 189733-like system in Fig. 4. This figure corresponds to the stacking of several 2D plots like those presented in Fig. 3 and therefore presents the same features.

Figures 3 and 4 demonstrate that, in theory, a large parameter space is expected in which we may prefer a DR model fit over a SB scenario. However, these figures do not convey the accuracy of the DR fit as there are no constraints on the rate of DR, that is, α. Hence, it is possible that an incorrect DR model could be preferred and therefore inject systematic biases into our understanding of solar-type stars.

Fig. 3 Two-dimensional heat maps illustrating the preferred parameter space (colour-coded by the change in BIC) for detecting stellar DR as a function of projected obliquity and stellar inclination for three different impact parameters (b = 0.1,0.4, and 0.7 from left to right). For a given system, as the impact parameter increases, the parameter space for preferring a DR model over SB decreases and isolated hot-spots appear. These results have been obtained with α = 0.2, and the parameters that were held constant in the three figures are listed in Table 1.

Fig. 4 Three-dimensional heat map of the likelihood of detecting DR, α = 0.2. The parameters held constant are listed in Table 1. The 3D surface represents the highest impact parameter at which DR can be retrieved for the corresponding stellar orientation. The structure of the volume below this surface is represented by way of the projections on the respective planes and according to the colour bar. A darker colour indicates a lower preference for DR in the projected area.

3.1.1 Accuracy of fits

To investigate the accuracy of our DR fits, we highlight in Fig. 5 only the regions in which both a DR fit is preferred and α agrees within 0.05 of the injected value. We find that the outcome is very similar, that is, the aforementioned distribution of hot spots is clearly visible. In particular, we observe that at low impact parameters, our ability to accurately constrain DR is only slightly worse (about λ = 0° and λ = ±90°), but it decreases significantly as the value of b increases. Comparing the right plots in Figs. 3 and 5, we observe a clear deterioration in the sharpness of the boundaries between high- and low-confidence areas (green and white areas). The areas in which DR was easily preferred previously are now more sparsely populated, which indicates that the DR fits are less reliable in these regions. Regarding the accuracy of the fits on λ, i_*, v_eq and α (when such parameters are fitted), our results are remarkably close to the injected value for the most part. Starting with λ, the mean offset between the fitted value and the real value is only 0.10°, with a standard deviation of 4.81°: the small offset means that the fitted results are distributed equally around the injected value, with minimum bias. However, the width of the distribution means that the fit is not perfectly accurate. These values are slightly higher on i_*, with a mean offset between fitted and real values of –0.20° and a standard deviation of 21.85°. While the fitting accuracy on v_eq looks worse with a mean offset of –6.71 km s^–1 and a standard deviation of 39.41 km s^–1, when we reconstruct the value of v_eqsin i_*, we find that the predicted value has a mean of 3.36 km s^–1 (standard deviation: 1.76 km s^–1) that is surprisingly close to the injected value of 3.25 km s^–1 (cf. Table 1). Regarding α, we find a mean offset of –0.09 and a standard deviation of 0.26. As Fig. 5 illustrates, the fitted values on α show a higher dispersion than the injected parameter, especially at a higher impact parameter.

Fig. 5 Two-dimensional heat maps illustrating the preferred parameter space (colour-coded by the change in BIC) for accurately detecting stellar DR as a function of projected obliquity and stellar inclination for three different impact parameters (b = 0.1,0.4, and 0.7 from left to right). Accurate detections are those for which the recovered α agrees with the injected value within 0.05. The behaviour is similar to Fig. 3, where accuracy is ignored, but the cold spots are even larger; regions that change colour between these two figures should be taken with caution in empirical analyses. These results have been obtained with α = 0.2, and the parameters that we held constant in the three figures are listed in Table 1.

Fig. 6 Two-dimensional heat maps illustrating the preferred parameter space (colour-coded by the BIC) for detecting stellar DR as a function of projected obliquity and stellar inclination for two different sampling rates (left: 100 s, right: 300 s). These results have been obtained with b = 0.5, V = 10 and α = 0.2, with the other parameters set to the values listed in Table 1.

3.1.2 Impact of the exposure time

We assumed an exposure time of 100s for our simulated data so far, which is a fairly optimistic assumption. Figure 6 shows the expected difference between this 100 s cadence and that of a 300 s exposure time, which is more in line with expectations for empirical observations of one transit for a bright target.

Increasing the exposure time increases the S/N and should yield better results, but the decreased sampling rate balances the increase in accuracy. Even an exposure time of 600s only leads to marginally better results. Furthermore, a long exposure time runs the risk of not sampling ingress and egress, which is more detrimental to DR detection than the high S/N. Integrating readout delay in the simulation has the effect of decreasing the sampling rate without increasing exposure time, but does not produce any significant difference.

3.1.3 Impact of brightness

Our ability to prefer a DR model is driven by both the precision and sampling rate of our measurements with respect to the absolute magnitude of the DR. Increases in the S/N will naturally lead to an increase in the measured precision on the local velocities. Hence, we explored the potential impact of the V magnitude by scaling our uncertainties following Eq. (1). In particular, we investigated scenarios for V = 10 and V = 12, which values are still considered fairly bright for transiting exoplanet systems, but are substantially dimmer than an HD 189733-like system (V ≈ 7.7).

Figure 7 demonstrates that it may be extremely difficult to measure a solar-like DR on even V = 10 systems using a HARPS-like setup, at least given the cadence herein. Detecting DR in these conditions would be exceedingly hard for any impact parameter b > 0.2. However, if we assume an ESPRESSO/VLT setup, we obtain very similar results as before (i.e. a V = 7.7 host star observed with a HARPS-like setup, see Fig. 8); this illustrates the powerful impact of increasing the precision through a collection of resolving power and increased instrumental precision. Figure 8 is remarkably analogous to Fig. 4. Decreasing the magnitude further and setting V = 12, we now approach the limit of an ESPRESSO/VLT-like setup, as shown in Fig. 9. This last figure is almost identical to Fig. 7, in that detecting DR in such a system seems to be beyond reach for now, at least with one transit.

Fig. 7 Three-dimensional heat map of the likelihood of detecting DR. using HD 189733 system parameters and V = 10, with a HARPS-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.

Fig. 8 Three-dimensional heat map of the likelihood of detecting DR. using HD 189733 system parameters and V = 10, with an ESPRESSO-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.

3.1.4 Impact of the stellar rotation rate

Alternatively, it is also easier to detect DR when the absolute magnitude of the effect is larger. Hence the limits on brightness and instrumental setup or precision can be somewhat mitigated depending on the v_eq sin i_*. of the target. Figures 7–9 have all been obtained using the v_eq sin i_* listed in Table 1, and we demonstrate that while an ESPRESSO/VLT like setup is viable for slowly rotating, dim stars (up until V ≈ 12), the same cannot be said of HARPS. However, if we consider a much faster rotator (i.e. v_eq sin i_* = 15 km s^–1), a HARPS-like setup becomes feasible, as Fig. 10 shows.

It is interesting to see in Fig. 10 that while the crests reach almost as high as they did previously, in the nominal HD 189733 case (cf. Fig. 4), the troughs in the 3D surface are much deeper. This suggests that while DR detection is definitely possible in the case of a dim but fast rotator using a HARPS-like setup, it is much more sensitive to the parameters of the star–planet orientation, and our ability to retrieve DR may be limited to a much reduced set of targets.

Fig. 9 Three-dimensional heat map of the likelihood of detecting DR, using HD 189733 system parameters and V = 12, with an ESPRESSO-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.

Fig. 10 Three-dimensional heat map of the likelihood of detecting DR, using HD 189733 system parameters, veq sin i* = 15 km s–1 and V = 10, with a HARPS-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.

3.2 Stellar rotation and centre-to-limb convection variation

In this second part, we now consider the compounded impact of both DR and CLV. We are interested in determining whether DR + CLV or SB + CLV is preferred and for which parameter space. We also examine whether SB + CLV can be confused for pure DR and pure DR be confused for pure CLV.

Starting with the first question, we evaluated the likelihood of retrieving the injected data (DR + CLV) as opposed to the simpler model (SB + CLV). In this configuration, we recall that the DR + CLV model has six free parameters (λ, i_*, v_eq, and α from the DR model described in Sect. 3.1, plus c₁ and c₂), while the SB + CLV model has four (λ and v_eq sin i_*, plus c₁ and c₂ as well). In the same way as for the case without CLV (cf. Sect. 3.1), Fig. 11 illustrates the preference for DR + CLV over SB + CLV for three values of the impact parameter b.

As before, the darker the shade of green, the more likely it is that we can retrieve the injected data. The main hot-spots to retrieve DR + CLV are:

• i_* > 110° and |λ| < 50°;

• i_* < 50° and |λ| > 100°;

• i_* ≈ 90° and |λ| ≈ 60°;

• i_* ≈ 90° and |λ| ≈ 120°.

Despite similarities with the results from the previous section, we observe some important changes around i_* ≈ 90° and λ ≈ +90°: Whereas we observed two distinct hot-spots in Fig. 3 (without CLV) that appeared at moderate to high impact parameters, now a more complex structure emerges. Instead of the two ellipses around λ ≈ +90° and 60° < i_* < 120°, four arms appear to wrap around the same parameter space, and the centre area (i_* = 90° and λ = +90°) remains empty. This result might be a sign that at this point, DR is being confused with the CLV and that consequently, the less complicated SB + CLV model is the preferred fit to the data points for this specific orientation.

Another major difference is that without CLV, these same hot-spots around i_* ≈ 90° and λ ≈ ±90° can still be observed at high impact parameters (cf. right plot of Fig. 3). However, when CLV is considered, the right plot of Fig. 11 shows that the four arms fade away at a lower impact parameter compared to the other hot-spots that are present in the centre plot of the same figure. Hence, when CLV is present in the injected data, retrieving it along with DR is more difficult the higher the impact parameter.

The 3D heat map (cf. Fig. 12) shows that the density map in the bottom panel is more contrasted than the case without CLV (cf. Fig. 4). While DR + CLV can be reliably detected up to high values of b around particular hot-spots, it is much less favoured between them, even at low impact parameters: there are dead spots in the parameter space when CLV is considered that were not present without CLV. This means that our ability to detect DR is much more dependent on the orientation of the system than before.

Fig. 11 Two-dimensional heat maps of the degree of preference for the DR + CLV model over the SB +CLV model, depending on the impact parameter b (b = 0.1,0.4, and 0.7 from left to right). These results have been obtained with α = 0.2, and the parameters that were held constant across the three figures are listed in Table 1.

3.2.1 Accuracy of fits

Like the previous subsection (cf. Sect. 3.1), we are interested in the accuracy of the fit on α where DR + CLV is preferred to SB + CLV. We compare Fig. 11 and Fig. 13 to determine whether the hot-spots that we previously identified result in a reliable detection of the differential rotation. The area in which both the DR + CLV model is preferred and the difference between the fitted value of a and the injected value is less than 0.05, is less expensive than in the first series of plots (DR-only). The sharpness of the aforementioned hot-spots declines sharply: while the underlying structure is very similar between Fig. 11 and Fig. 13, it becomes more fuzzy in the latter. The conclusion that we drew in the DR-only case seems to apply in this case as well: a particular area might be favourable to retrieve the injected model on average, but given a single point in the parameter space, the fit might not be reliable.

To verify whether the CLV could be mistaken for DR and introduce a bias in the fit on a, we ran several sets of simulations (see Table 1 for the system configuration), injecting SB + CLV (setting the parameter α = 0). In the first case, we compared a SB fit with a DR fit. Our results show that without considering the accuracy of the fit on the value of a, a DR model explains the injected SB + CLV better than the SB model and could therefore point towards a true SB + CLV being confused for DR (if the CLV contribution to the fit were not considered). We suggest that this ambiguity would certainly be lifted if the fitted values of α could be checked against a credible estimate of the true value. However, the challenge is indeed that a sensible estimation of the true value of α is needed to begin with, which is a problematic prospect when constraining the value of this parameter is the primary intention of the study.

Interestingly, despite the inconclusive results of this approach, almost all of the best fits we obtained when trying to fit an injected SB + CLV with a DR model led to a recovered α = ±1. Such extreme results should certainly be interpreted with caution and might point towards the conclusion that trying to fit a pure DR model on SB + CLV will produce inadequate effects. In the second scenario, we take our case with the injected SB + CLV velocities and compare an SB + CLV model fit with a DR model (i.e. SB + CLV vs DR, as opposed to the previous SB vs. DR test), which is interesting because both models have the same number of free parameters. Without imposing any conditions on any parameters, we find that the DR fit cannot explain the data better than the SB + CLV fit. Hence, if the true stellar surface is explained by SB + CLV, it may be possible to prefer a DR model over an SB model (albeit an inaccurate DR), but the best fit will still be from a SB + CLV model. Regardless of fitting an SB, DR, or SB + CLV model to our model star (with SB + CLV injected), the recovered projected obliquity remained robust. From these results, we conclude that while the possibility of CLV being mistaken for DR cannot completely be ruled out, it seems to be very limited. For the accuracy of the fits on λ, i_* and veq, our results are in line with those we obtained previously. We are able to retrieve the injected parameters with good precision: The mean offset between the fitted and the injected value for λ is slightly worse than in Sect. 3.1.1 at –0.48° with a standard deviation of 13.45°. The increased width of the offset distribution shows that the fit is not quite as reliable as in the DR case. For i_* we find a mean offset of 0.16° with a standard deviation of 20.02°. Despite a smaller mean offset on v_eq at –1.5 km s^–1 (standard deviation 26.51 km s^–1) compared to Sect. 3.1.1, the mean reconstructed value for v_eqsin i_* is farther away from the injected value, at 3.57 km s^–1 (standard deviation 3.12 km s^–1) instead of 3.25 km s^–1 (cf. Table 1). Lastly, the fitted α, despite a mean offset to the injected value of –0.05 (standard deviation 0.18), is not as reliable as in the DR-only case, especially at high impact parameters (cf. Fig. 13).

Fig. 12 Three-dimensional heat map of the likelihood of retrieving DR + CLV as compared to SB + CLV. A darker colour indicates a lower preference for DR + CLV in the projected area. These results have been obtained with α = 0.2, and the parameters that were held constant are listed in Table 1.

Fig. 13 DR + CLV preference over SB +CLV depending on b (b = 0.1,0.4, and 0.7 from left to right), considering the proximity of the fitted value of α to the injected value. Accurate detections are those where the recovered α agrees with the injected value within 0.05. These results have been obtained with α = 0.2, and the parameters that were held constant in the three figures are listed in Table 1.

3.2.2 Impact of the sampling rate

When we examined the impact of exposure time or cadence on our ability to retrieve DR + CLV, our observations are largely the same as in the DR only case. We conclude that while the increase in S/N resulting from the longer exposure time helps, the benefits are largely negated by the lower sampling rate. Figure 14 illustrates this point.

Complementary results support the fact that while the two effects we discussed cancel out on a single transit, obtaining multiple transits helps to avoid the negative effect of a lower sampling rate while maintaining the full benefit of the higher S/N.

3.2.3 Impact of brightness and stellar rotation rate

For the influence of brightness and the rotation rate for the preference for the DR + CLV over the SB + CLV fit, our conclusions are once again very similar to the first section. A drop in brightness to V = 10 or V = 12 would prove very challenging, even considering the precision of an ESPRESSO/VLT setup.

Comparison of Fig. 15 with Fig. 10 shows that if CLV is present in the injected data (at the rate injected here, following a solar MHD simulation of the Fe I 6302 Ź), a solar-like DR is unlikely to be retrieved in a practical way using a HARPS-like setup if the star is not observed pole-on or pole-away, even in the case of a fast rotator (v_eq sin i_* = 15 km s^–1 in both figures). The fraction of the parameter space in which the system orientation favours detection is much more limited when CLV is injected in the model.

4 Conclusion

To take full advantage of the improved precision of recent instruments (e.g. ESPRESSO), a more detailed understanding of the stellar surface is important. Throughout this paper, we have presented our results for the optimal parameter space for detecting stellar differential rotation and centre-to-limb convective variations. To do this, we simulated a typical transiting hot Jupiter (similar to HD 189733) and varied several parameters, including the projected obliquity, stellar inclination, and impact parameter.

In the first part of this paper, we examined the likelihood of retrieving an accurate estimate of the differential rotational sheer a depending on system parameters, leaving out any contribution of the CLV in the injected data. We find that DR detection is null for λ = ±90° at low impact parameters (~b ∈ [0,0.2]) for all i_*. The main hot-spots for detection are 120° < |λ| < 180° for i_* < 90°, and |λ| < 60° for i_* > 90° on average, with the planet transiting the lower half of the star in our coordinate system. Additionally, at moderate impact parameters, DR detection is likely at |λ| ≈ 90° for i_* ≈ 90°, which corresponds to the configuration in which the largest number of latitudes are crossed. At higher impact parameters (~b > 0.7), DR detection is increasingly difficult even in the regions described above. We also find that the estimate on a is reliable and only drops in accuracy near the edges of the hot-spots and with increasing values of b.

We also observe that increasing the exposure time does not drastically improve detection, as the increase in S/N is offset by the loss of sampling. Unsurprisingly, multiple transits will alleviate the adverse impact of the lower sampling rate.

Regarding the influence of brightness and rotation rate, we observe that the detection limit for a HARPS-like setup is probably about V = 10 for a moderately slow rotator such as HD 189733, while V = 12 seems to be the limit considering an ESPRESSO-like setup. If we consider faster rotators (e.g. ~ 15 km s^–1), a HARPS-like setup may detect DR beyond V = 10 (if it is photon-noise dominated), although within a more limited parameter space.

In the second part of this paper, we investigated the impact of a CLV component on the injected data. Similarly to the first case, detection still is compromised for λ = ±90° at low impact parameters (for all i_*). However, the emerging structure in the parameter space is different and more complex than previously. The main hot-spots for detecting DR + CLV are |λ| > 100° for i_* < 50° and |λ| < 50° for i_* > 110°, which are visible up until b = 0.7, that is slightly worse than in the first case. Additionally, we observe secondary hot-spots at |λ| < 50° for i_* < 60°, |λ| > 110° for i_* > 120°, as well as |λ| ≈ 60° and |λ| ≈ 120° for i_* ≈ 90°, which subside at lower impact parameter than the others as b increases. Surprisingly, DR + CLV detection is null in the region around |λ| ≈ 90° and i_* | ≈ 90° that we identified in the first part, maybe due to degeneracy between the two components.

When comparing the fitted value on α with the injected parameter, we find that the hot-spots near i_* = 90° are practically non-existent at any impact parameter, constraining DR + CLV detection to the edges of the parameter space. The other conclusions we reached in the first part regarding the influence of brightness, exposure time, and rotation rate are still applicable, with the major difference being that reliable detection of DR + CLV is confined to more well-delimited regions of the parameter space, and to lower impact parameters overall. For these reasons, detection for V > 10 might prove challenging even with an ESPRESSO-like setup, especially if the number of available transits is limited.

Lastly, we established that CLV is unlikely to be mistaken for DR if reasonable constraints can be placed on the value of a. It seems that when CLV and DR might be confused, the DR fit produces largely unreliable results that would not pass scrutiny.

Fig. 14 Two-dimensional heat maps of the parameter space (b = 0.5), with the preference for DR + CLV over SB + CLV colour-coded in green for two different sampling rates (left: 100 s, right: 300 s). These results have been obtained with α = 0.2, and the parameters that were held constant across the three figures are listed in Table 1.

Fig. 15 Three-dimensional heat map of the likelihood of detecting DR + CLV, using HD 189733 system parameters, veq sin i* = 15 km s–1 and V = 10, with a HARPS-like precision. These results have been obtained with α = 0.2, and the parameters that were held constant are listed in Table 1.

Acknowledgements

The authors acknowledge the financial support of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). This research has made use of NASA’s Astrophysics Data System Bibliographic Services. HMC was supported during this project from both an NCCR PlanetS CHEOPS Fellowship and a UKRI Future Leaders Fellowship (MR/S035214/1). VB acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project SPICE DUNE, grant agreement No 947634).

References

All Tables

All Figures

	Fig. 1 Simulated local stellar surface velocity below the planet along the transit chord as a function of the phase (bottom) and stellar radius (top) for a given star–planet orientation (λ = 0°; i* = 150°; b = 0.2; V = 10). The best-fit DR model is shown in blue, and the best-fit SB model is plotted in orange; a positive BIC indicates a preference for the more complex DR model.
In the text

Fig. 2 Simulated local stellar surface velocity below the planet along the transit chord as a function of the phase (bottom) and stellar radius (top) for a given star–planet orientation (λ = 0°; i* = 150°; b = 0.2; V = 10), wherein both a solar DR and CLV are injected; compare to Fig. 1 to note the CLV impact (overplotted here in light colours), particularly the relative redshifts near the limbs. The corresponding DR + CLV model fit is shown in blue, and the CLV + SB fit is plotted in orange; a positive BIC indicates a preference for the more complex DR + CLV model.

In the text

Fig. 3 Two-dimensional heat maps illustrating the preferred parameter space (colour-coded by the change in BIC) for detecting stellar DR as a function of projected obliquity and stellar inclination for three different impact parameters (b = 0.1,0.4, and 0.7 from left to right). For a given system, as the impact parameter increases, the parameter space for preferring a DR model over SB decreases and isolated hot-spots appear. These results have been obtained with α = 0.2, and the parameters that were held constant in the three figures are listed in Table 1.

In the text

Fig. 4 Three-dimensional heat map of the likelihood of detecting DR, α = 0.2. The parameters held constant are listed in Table 1. The 3D surface represents the highest impact parameter at which DR can be retrieved for the corresponding stellar orientation. The structure of the volume below this surface is represented by way of the projections on the respective planes and according to the colour bar. A darker colour indicates a lower preference for DR in the projected area.

In the text

Fig. 5 Two-dimensional heat maps illustrating the preferred parameter space (colour-coded by the change in BIC) for accurately detecting stellar DR as a function of projected obliquity and stellar inclination for three different impact parameters (b = 0.1,0.4, and 0.7 from left to right). Accurate detections are those for which the recovered α agrees with the injected value within 0.05. The behaviour is similar to Fig. 3, where accuracy is ignored, but the cold spots are even larger; regions that change colour between these two figures should be taken with caution in empirical analyses. These results have been obtained with α = 0.2, and the parameters that we held constant in the three figures are listed in Table 1.

In the text

	Fig. 6 Two-dimensional heat maps illustrating the preferred parameter space (colour-coded by the BIC) for detecting stellar DR as a function of projected obliquity and stellar inclination for two different sampling rates (left: 100 s, right: 300 s). These results have been obtained with b = 0.5, V = 10 and α = 0.2, with the other parameters set to the values listed in Table 1.
In the text

	Fig. 7 Three-dimensional heat map of the likelihood of detecting DR. using HD 189733 system parameters and V = 10, with a HARPS-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.
In the text

	Fig. 8 Three-dimensional heat map of the likelihood of detecting DR. using HD 189733 system parameters and V = 10, with an ESPRESSO-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.
In the text

	Fig. 9 Three-dimensional heat map of the likelihood of detecting DR, using HD 189733 system parameters and V = 12, with an ESPRESSO-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.
In the text

	Fig. 10 Three-dimensional heat map of the likelihood of detecting DR, using HD 189733 system parameters, veq sin i* = 15 km s–1 and V = 10, with a HARPS-like precision. These results have been obtained with α = 0.2. The parameters that were held constant are listed in Table 1.
In the text

	Fig. 11 Two-dimensional heat maps of the degree of preference for the DR + CLV model over the SB +CLV model, depending on the impact parameter b (b = 0.1,0.4, and 0.7 from left to right). These results have been obtained with α = 0.2, and the parameters that were held constant across the three figures are listed in Table 1.
In the text

	Fig. 12 Three-dimensional heat map of the likelihood of retrieving DR + CLV as compared to SB + CLV. A darker colour indicates a lower preference for DR + CLV in the projected area. These results have been obtained with α = 0.2, and the parameters that were held constant are listed in Table 1.
In the text

Fig. 13 DR + CLV preference over SB +CLV depending on b (b = 0.1,0.4, and 0.7 from left to right), considering the proximity of the fitted value of α to the injected value. Accurate detections are those where the recovered α agrees with the injected value within 0.05. These results have been obtained with α = 0.2, and the parameters that were held constant in the three figures are listed in Table 1.

In the text

	Fig. 14 Two-dimensional heat maps of the parameter space (b = 0.5), with the preference for DR + CLV over SB + CLV colour-coded in green for two different sampling rates (left: 100 s, right: 300 s). These results have been obtained with α = 0.2, and the parameters that were held constant across the three figures are listed in Table 1.
In the text

	Fig. 15 Three-dimensional heat map of the likelihood of detecting DR + CLV, using HD 189733 system parameters, veq sin i* = 15 km s–1 and V = 10, with a HARPS-like precision. These results have been obtained with α = 0.2, and the parameters that were held constant are listed in Table 1.
In the text