© ESO 2019

1. Introduction

The interaction of convection and rotation produces differential rotation and magnetic fields in the Sun. On timescales of billions of years, magnetic fields affect the rotation of our star by the angular momentum loss associated with its magnetized stellar wind. A better understanding of these complex processes and their interconnections can be obtained by comparing the Sun with other late-type stars, in particular with younger Sun-like stars that show a higher level of magnetic activity (e.g. Brun et al. 2015; Brun & Browning 2017). Space-borne telescopes, such as CoRoT (Auvergne et al. 2009) or Kepler (Koch et al. 2010), allow us to monitor late-type stars photometrically in the optical passband with a relative accuracy down to 10⁻⁵ − 10⁻⁴ with integration times from minutes to hours for time intervals up to 3–4 years, thus providing unprecedented datasets to study their magnetic activity.

Magnetic fields produced in the stellar interior by hydromagnetic dynamos emerge into the photosphere where they modify the transport of energy and momentum giving rise to cooler and hotter patches called starspots and faculae, respectively (e.g. Gondoin 2008; Strassmeier 2009), that modulate the optical flux integrated over the stellar disc owing to their intrinsic evolution and the rotation of the star. This flux modulation can be modelled to extract information on the locations of the surface brightness inhomogeneities and their evolution. Recent works have explored the possibility offered by starspots as tracers of stellar rotation and differential rotation (e.g. Mosser et al. 2009; Walkowicz et al. 2013; Santos et al. 2017), while a general review of the different approaches to spot modelling can be found in, for example, Lanza (2016) and a comparison between spot models and sunspot group observations is provided in Lanza et al. (2007). Models are better constrained when the inclination of the stellar spin axis to the line of sight is known, such as in the case of stars that have a transiting planet (e.g. Winn et al. 2005; Nutzman et al. 2011).

Among the stars with a transiting hot Jupiter, Kepler-17 is one of the targets with more extended and precise transit observations (Müller et al. 2013) allowing for the spots occulted during the transits to be mapped (Estrela & Valio 2016; Valio et al. 2017) and for the inclination of the stellar spin axis to be constrained (Désert et al. 2011). Moreover, it is a young star of G2V spectral type making it an ideal candidate for solar–stellar connection studies. It has an estimated age of ≲1.8 Gyr, a mean rotation period of ∼12 days, and is accompanied by a planet with a mass of 2.47 ± 0.10 Jupiter masses, a radius of 1.33 ± 0.04 Jupiter radii, and an orbital period P_orb = 1.48571 days (Bonomo et al. 2012).

Bonomo & Lanza (2012) analysed ∼500 days of public Kepler data available at that time; now the availability of the latest Kepler data release with a high-precision photometric time series extending for ∼1500 days calls for a new modelling of this star to study its activity and rotation using starspots as tracers. Moreover, we compare the longitudes of the spots mapped from the out-of-transit photometry with those of the spots occulted during transits, providing independent confirmation of our results and giving constraints on the spot latitudes that cannot be obtained with alternative methods in the case of Kepler-17. In such a way, we investigate the differential rotation of our target, the phenomenology of its active longitudes, and its activity cycles.

The presence of a close-by giant planet affects the properties of Kepler-17, notably its rotation that is used to estimate its age by applying the gyrochronology technique (Barnes 2007, 2010). We look for features associated with star-planet interaction in the photometric time series and investigate the tidal evolution of the system to clarify the difference in the evolution of the stellar rotation with respect to a single star such as our Sun. Moreover, the peculiar evolution of the rotation of Kepler-17 affects the flux of the stellar high-energy radiation experienced by the planet during its lifetime.

2. Observations

The Kepler 95 cm telescope was designed to continuously look at a fixed field in the Cygnus constellation to detect planetary transits. Every three months, the spacecraft is rolled by 90° about its line of sight to keep its solar panels pointing towards the Sun. Each of these periods is called a quarter in Kepler nomenclature. Because of the rotation of the focal plane, each target falls on a different CCD during different quarters, so the observations must be reduced quarter-by-quarter. The usual cadence of Kepler observations is 1765.5 s (long cadence), although for a subset of interesting targets, such as those showing transits, the cadence is reduced to 58.5 s (short cadence). We shall analyse photometry acquired in long cadence because we are interested in the activity and rotation of Kepler-17 both of which have timescales of the order of several days, while short cadence has been used to observe planetary transits and detect spots occulted during transits (Estrela & Valio 2016; Valio et al. 2017).

In the mission archive¹ there are the time series of the electron counts in the individual CCD pixels within a pre-defined area around the image of each target in the focal plane and two light curves. The first is obtained by summing the flux falling within a subset of the pixels included in the target pixel files; a correction for the background flux is also applied. This is called the simple aperture photometry (SAP) light curve. The second time series is obtained by additional processing to remove instrumental and systematic effects and is called the pre-search data conditioning (PDC) light curve (Stumpe et al. 2012, 2014).

We downloaded from the mission archive all the long-cadence SAP time series of the latest data release (Data Release 25) of Kepler-17 covering 14 quarters out of a total of 18. This time series is affected by outliers and systematic instrumental effects that have been corrected in the PDC time series. However, the PDC time series shows a significant reduction of the amplitude of intrinsic stellar variability and sometimes a distortion of the modulations on time scales longer than a few days because it has been designed to detect planetary transits, not to preserve the intrinsic stellar variability (see Gilliland et al. 2015). These effects are particularly relevant for Kepler-17 because it shows a large light curve amplitude of the order of 0.05 mag with typical modulation timescales ≥10 days. For these reasons, we decided not to use the PDC time series in the present analysis and derived two light curves starting from the SAP time series.

The first was derived by means of a procedure called ARC2 introduced by Aigrain et al. (2017) that warrants a better preservation of the intrinsic stellar variability while removing discontinuities, outliers, and instrumental effects by making use of the co-trending basis vectors (CBVs) computed by the PDC pipeline. Co-trending basis vectors describe instrumental effects for each target and are based on the systematic trends observed in targets that are close on the CCD to the given target and are similar in flux. They also take into account instrumental effects by making use of the spacecraft telemetry information. Therefore, they provide the best available description of the systematic trends affecting the time series of a given target. Up to eight CBVs are used by the PDC pipeline to perform its correction often leading to an overcorrection of the instrinsic stellar variability on timescales much longer than those characteristics of planetary transits, and to an injection of additional noise on those timescales. On the other hand, the ARC2 pipeline of Aigrain et al. (2017) was designed to preserve the intrinsic target variability and reduce the injected noise in the correction process as much as possible and is particularly valuable in the case of bright and remarkably variable stars such as Kepler-17 whose light curves are modified by the PDC pipeline.

In our application, we started from the SAP time series considering only the datapoints with a SAP_QUALITY flag equal to zero, that is without any detected problem during their acquisition (Jenkins et al. 2016). The ARC2 pipeline includes a Bayesian method to find the best weights to calculate the linear combination of the CBVs to correct each quarter time series (see Sect. 3.1 in Aigrain et al. 2017) as well as a criterion to select the optimal number of CBVs to be used. The latter is based on balancing the reduction of the normalized light curve amplitude, that comes from removing systematics by adding successive CBVs, against the increase of the short-term noise resulting from the same operation (see Sect. 3.2 of Aigrain et al. 2017, for a quantitative description of the criterion).

We removed the planetary transits from the light curve provided by ARC2 using the ephemeris of Müller et al. (2013), discarding the datapoints before 0.05 and after 0.05 in phase of the first and the fourth contacts, respectively. Each quarter was then normalized to its median value and residual outliers were removed by applying a 3σ clipping to the residuals obtained by subtracting a smoothed version of the light curve obtained with a boxcar filter with a width of 294 min, that is, ten consecutive datapoints. A total of 340 datapoints were flagged as outliers and discarded. The final light curve consists of 40 653 datapoints and ranges from BJD_TDB ² 2454964.512 to 2456423.980 covering a total of 1459.469 days with four main gaps (see Fig. 1, top panel red plot). The median of the error of the data points is 2.54 × 10⁻⁴ in relative flux units.

Fig. 1. Top panel: light curve of Kepler-17 as obtained by applying the de-trending as in Bonomo & Lanza (2012) (in green) and the light curve as provided by the ARC2 pipeline (in red). Each light curve is normalized to the median flux within each quarter after removing the planetary transits. The green datapoints have been shifted upward by 0.0625 for clarity. Bottom panel: difference between the two light curves. The quarters where only the first CBV has been applied to obtain the ARC2 light curve are indicated with blue datapoints, while the quarters where the first two CBVs have been applied are indicated in pink (see the text).

To allow a straightforward comparison with the results previously obtained by Bonomo & Lanza (2012), a second light curve was derived from the SAP time series following their approach. In this case, we choose to consider all the datapoints with a finite flux value, not limiting ourselves to those with a SAP_QUALITY flag equal to zero, because a different flag does not indicate in general that the datapoint is unusable as described in the Kepler Archive Manual Sect. 2.3.1.1. In brief, the procedure used to obtain the analysed light curve was as follows. Data were corrected quarter-by-quarter by discarding the steep variations after the safe modes or data download links; we then removed the points in transits as in the case of the ARC2 light curve. Subsequently, long-term trends of clear instrumental origin were removed by fitting a parabola. Finally, each quarter was normalized to its median value and the 3σ clipping procedure used for the ARC2 light curve was applied to flag and remove outliers. In such a way, a total of 1219 datapoints were discarded. The final light curve consists of 47 376 datapoints ranging from BJD_TDB 2454964.512 to 2456424.001 for a total of 1459.489 days with three main gaps. The number of datapoints is greater than in the ARC2 light curve because there are points covering the first gap in the ARC2 light curve and some others covering some small gaps in other intervals of the light curve (see Fig. 1, top panel green plot) thanks to our less strict selection criterion including datapoints with non-zero SAP_QUALITY flag. The median of the error of the data points is 2.22 × 10⁻⁴ in relative flux units.

The difference between the two light curves is plotted in the bottom panel of Fig. 1. The flux values of the light curve obtained with the approach of Bonomo & Lanza (2012) were linearly interpolated at the epochs of the ARC2 time series to compute the flux difference. Such an interpolation was in order because there were 492 datapoints of the ARC2 light curve that did not have a corresponding datapoint at the same time in the other light curve owing to the different criteria adopted to reject points affected by systematics close to data gaps or steep variations (automated removal in the case of ARC2, manual removal in the other case). The correction computed by the ARC2 pipeline makes use of only the first CBV for eight quarters and of the first two CBVs for the remaining six quarters, providing a better preservation of the intrinsic stellar variability than the PDC light curve. This is confirmed by the small difference with the light curve computed with the approach of Bonomo & Lanza (2012) for which the long-term trends inside each quarter were corrected by a simple parabolic fit and the steep variations detected by eye simply eliminated. This produces a flatter light curve, that is, showing less modulation of its mean level on timescales comparable with a quarter or longer. The parabolic shape of the difference between the two light curves in some quarters (see Fig. 1, lower panel) corresponds to quarters where only the first CBV was used for the correction of the ARC2 light curve.

3. Methods

The spot modelling approach applied in the present study is the same already introduced in Sect. 3 of Bonomo & Lanza (2012) to which we refer the reader for details. In brief, the surface of the star is subdivided into 200 surface elements that contain unperturbed photosphere, dark spots, and solar-like faculae. The specific intensity of the unperturbed photosphere in the Kepler passband is assumed to vary according to a quadratic limb-darkening law:

(1)

where I₀ is the specific intensity at the centre of the disc, μ = cosθ with θ being the angle between the local surface normal and the line of sight, and a_p, b_p, and c_p are the limb-darkening coefficients in the Kepler passband.

The dark spots are assumed to have a fixed contrast c_s ≡ I_spot(μ)/I(μ) in the Kepler passband, where I_spot is the specific intensity in the spotted photosphere. The fraction of a surface element covered by dark spots is given by its filling factor f. The faculae are assumed to have a fixed contrast c_f = 1.115 at the limb that varies linearly with μ becoming unity (no flux perturbation) at the centre of the disc. In this way, they mimic the contrast behaviour of solar photospheric faculae, at least in a rough and average sense. The ratio Q of their area to that of the dark spots is fixed, so that their filling factor is Qf. In our model, Q always appears in combination with c_f in the product c_fQ. Therefore, it is sufficient to vary Q in order to change the contribution of the faculae in our model (cf. Lanza et al. 2007; Lanza 2016).

This model is fitted to a segment of the light curve of duration Δt_f (see Sect. 4) by varying the filling factors of the individual surface elements that can be represented as a 200-element vector f. Therefore, the model has 200 free parameters and suffers from non-uniqueness and instability due to the effect of photometric noise. To select a unique and stable solution, we apply a maximum entropy (ME) regularization by minimizing a functional Z that is a linear combination of the χ² and of a suitable entropy function S:

(2)

where λ >  0 is a Lagrangian multiplier that controls the relative weights given to the χ² minimization and the configuration entropy of the surface map S in the solution. The expression of S is given in Eq. (5) of Bonomo & Lanza (2012); it is maximal when the star is unspotted, that is, when all the elements of the vector f are zero. In other words, the ME criterion selects the solution with the minimum spotted area compatible with a given χ² value of the best fit to the light curve. When the Lagrangian multiplier λ = 0, we obtain the solution corresponding to the minimum χ² that is unstable. By increasing λ, we obtain a unique and stable solution at the price of increasing the value of the χ². An additional effect is that of making the residuals between the model and the light curve biased towards negative values because we reduce the spot filling factors by introducing the entropy term (see Lanza et al. 1998; Lanza 2016, for details).

The information on the latitude of the spots is lacking in our maximum-entropy maps because the inclination of the stellar spin axis is very close to 90° (cf. Sect. 4), which makes the transit time of each feature independent of its latitude. Therefore, we limit ourselves to mapping the distribution of the filling factor versus the longitude.

The optimal value of the Lagrangian multiplier λ is obtained by imposing that the mean μ_reg of the residuals between the regularized model and the light curve verifies the relationship (Bonomo & Lanza 2012; Lanza 2016):

(3)

where σ₀ is the standard deviation of the residuals of the unregularized model, that is, that computed with λ = 0, with N being the number of datapoints in the fitted light curve interval of duration Δt_f.

The optimal value of Δt_f is not known a priori and must be determined with an analysis of the light curve itself because it is related to the lifetimes of the active regions in a given star. We shall adopt a unique value of Δt_f for the entire light curve of Kepler-17 because the ratio Δt_f/P_rot, where P_rot is the stellar rotation period, controls the sensitivity of the spot modelling to active regions located at different longitudes as discussed by Lanza et al. (2007).

The optimal value of the facular-to-spotted area ratio Q is also derived from the light curve best fit. To find the best values of Δt_f and Q, we use a simple spot model consisting of three active regions and a varying background level that was introduced in the case of the Sun (cf. Sect. 3 of Lanza et al. 2003) and previously applied by Bonomo & Lanza (2012) to Kepler-17. It has a much smaller number of free parameters than the maximum-entropy model, allowing us a faster exploration of the Δt_f-Q parameter space to look for the combination that minimizes the χ² of the entire light curve.

In principle, Δt_f and Q are not parameters of the same kind because the former is the duration of the individually fitted time intervals, while the latter is a model parameter. However, Δt_f affects the determination of Q, because Q is constrained by the different facular contrasts between the centre of the disc and the limb. In other words, the determination of Q is affected by the duration of the intervals during which the model active regions are close to the centre or to the limb, which in turn depend on Δt_f. For example, too short a Δt_f does not allow the model active regions to move all along their chords across the stellar disc, meaning that they cannot span the full range between the centre of the disc and the limb and provide the best constraint on Q. In this case, the dependence of the χ² of the entire light curve on Q is characterized by random oscillations. For this reason, we optimize Q for several fixed values of Δt_f considering the range where the total χ² of the three-spot model depends regularly on Q and look for the minimum χ² in the Δt_f-Q space.

4. Stellar parameters

The basic stellar parameters, that is, mass, radius, effective temperature T_eff, and surface gravity log g, are taken from Bonomo et al. (2012) and are the same adopted by Bonomo & Lanza (2012). They do not directly enter into our geometric spot model, except for the computation of the relative difference ϵ_rot between the polar and the equatorial axes of the ellipsoid used to represent the surface of the star. Its value is obtained by a simple Roche model assuming rigid rotation with a period of P_rot = 12.01 days as in Bonomo & Lanza (2012). The gravity darkening effect associated with ϵ_rot ∼ 4.7 × 10⁻⁵ is of the order of a few times 10⁻⁶ mag, and therefore it can be neglected in our model.

The best fit to the transits of Kepler-17b can be used to extract information on the inclination i of the stellar spin axis to the line of sight and on the quadratic limb-darkening coefficients. Müller et al. (2013) provided a refined analysis of the first seven quarters of Kepler short-cadence observations adopting the stellar parameters of Bonomo et al. (2012) and compared the limb-darkening parameters derived from the transit fit with those given by model atmospheres. The precision of the limb-darkening parameters derived from the fitting of the average transit profile is remarkably high because Kepler-17 has a photometry with a high signal-to-noise ratio. Nevertheless, the theoretical linear coefficient is significantly different from that derived from the transit fitting. The discrepancy is more clearly evident in the recent analysis by Maxted (2018) who considered all the short-cadence quarters available in Kepler Data Release 25 and suggested that the responsibility lies with the remarkable activity of Kepler-17 (cf. his Fig. 4). In other words, the active regions on the stellar surface, most of which are not resolved even by methods applied to detect spot occultations during transits, make the limb-darkening profile of Kepler-17 significantly different from that computed from model atmospheres. In our spot modelling, we shall adopt the limb-darkening coefficients as derived from the transit fitting, thus modelling the out-of-transit light curve in a way consistent with that adopted for the occulted spots. However, we shall investigate the impact of this choice by computing additional models with the theoretical limb-darkening coefficients to show that our main results are not critically dependent on those coefficients (cf. Sect. 5 and Appendix A). Specifically, we shall adopt the coefficients as derived by Müller et al. (2013) from a model atmosphere with T_eff = 5787 K, log g = 4.45, and solar metallicity (cf. their Sect. 5.1).

In consideration of his use of the full dataset of the latest Kepler Data Release, we assume the transit fit by Maxted (2018) as the best for our purposes and derive the quadratic limb-darkening coefficients from his model. Specifically, given that Maxted (2018) adopted a non-polynomial limb-darkening law of the form I(μ)/I₀ = 1 − c(1 − μ^α) with free parameters c and α, we fitted our quadratic law in Eq. (1) to 10⁵ realizations of this equation with the values of the coefficients c and α drawn from their a posteriori distributions. The quadratic fit is good, except for a relative deviation of ∼2−3 % close to the limb (μ >  0.8), which can be neglected because the photometric effect of starspots close to the limb is negligible owing to the strong reduction of their projected area by foreshortening. In such a way, we obtain the limb-darkening coefficients in Table 1 as the medians of their a posteriori distributions. We point out that we cannot adopt the same functional form of the limb darkening as in Maxted (2018) because our computer code should have been completely rewritten given that it makes use of the quadratic form of the limb-darkening law to speed up calculations.

The inclination of the stellar spin axis is assumed equal to that of the orbital plane as derived from the transit fitting as in Bonomo & Lanza (2012). Adopting the parameters of the model of Maxted (2018), we obtain the median value listed in Table 1, which is also compatible with that of Müller et al. (2013) and with the estimate coming from the rotation spectral line broadening νsini, P_rot and the estimated stellar radius (Bonomo & Lanza 2012). The mean rotation period P_rot is adopted as in Bonomo & Lanza (2012) for a straightforward comparison with their results.

The contrast of the dark spots was assumed to be equal to that of sunspot groups, that is c_s = 0.677, by Bonomo & Lanza (2012). However, the recent work by Valio et al. (2017) provides a direct measure based on the modelling of starspots occulted during transits. Their modelling gives the mean value c_s = 0.55 ± 0.17 adopted in the present analysis. However, in view of the large differences in c_s among the different spots modelled by Valio et al. (2017), in Sect. 5 and in Appendix A, we explore the effects of the variation of the spot contrast on our results by computing regularized models with the two extreme values c_s = 0.38 and 0.72.

The facular-to-spotted area ratio Q was derived using the simplified three-spot model to fit the entire light curve of Kepler-17 by selecting the value that gives the minimum total χ². The plot of the ratio between the χ² and its minimum versus Q is shown in Fig. 2 for the ARC2 light curve and in Fig. 3 for the light curve de-trended as in Bonomo & Lanza (2012) assuming Δt_f = 8.733 days. The ratio of the χ² to its minimum is statistically distributed according to the Fischer-Snedecor statistics F as (cf. Lampton et al. 1976, Sect. VI):

(4)

Fig. 2. Ratio of the χ2 of the best fit to the entire ARC2 light curve to its minimum value vs. the parameter Q, i.e. the ratio of the area of the faculae to that of the dark spots in active regions. The horizontal dashed line indicates the 95% confidence level for determining the interval of acceptable Q values.

Fig. 3. As in Fig. 2 but for the best fit to the light curve obtained with the de-trending of Bonomo & Lanza (2012).

where p is the total number of free parameters in the model and N_T the total number of datapoints in the fitted time series. Equation (4) allows us to estimate the confidence interval of the parameter Q that depends on the maximum value of the ratio corresponding to a given confidence limit (Lampton et al. 1976).

The best fit to the ARC2 light curve shows some oscillations of the total χ² as a function of Q, probably related to some residual suppression of the facular modulation by the ARC2 pipeline. Using the original PDC light curve, the oscillations dominate the plot and prevent the determination of the optimal Q value because there is no clearly determined global minimum. This happens because solar-like faculae produce a photometric signal only when they are close to the limb, that is, only in limited intervals of the rotational modulation produced by active regions (Lanza 2016). Such a tiny signal is easily suppressed by the PDC de-trending leading to insufficient information to constrain the value of Q in the PDC light curve. On the other hand, the ARC2 pipeline, although still based on the use of the CBVs derived by the PDC pipeline, applies only the first one or a linear combination of the first two CBVs in an attempt to preserve the intrinsic stellar variability on all the accessible timescales. In such a way, most of the facular signal is preserved allowing us to constrain the value of Q. The light curve de-trended with the method by Bonomo & Lanza (2012), that removes only long-term trends by a parabolic fit, gives an even cleaner result. From both the light curves, we derive an optimal value of Q = 2.4 that shall be adopted for our analysis. It is different from the value Q = 1.6 found by Bonomo & Lanza (2012) as a consequence of the smaller value of the spot contrast c_s used in the present modelling as we verified by running another minimization with the previous c_s = 0.667. Given the impact of the value of Q on our results, in Sect. 5 and Appendix A, we also consider spot models with Q = 1.0 and Q = 4.0, which are well beyond the 95% confidence regions in Figs. 2 and 3, to explore the effect of varying Q on our models.

The duration of the individual segments of the light curves fitted with the three-spot model has been kept at Δt_f = 8.733 days in all of the above models, that is, the same value given in Bonomo & Lanza (2012). When we increase Δt_f by only 2.5%, the total χ² found by minimising with respect to Q becomes significantly worse, increasing beyond the value corresponding to the 95% confidence region as computed with the statistics in Eq. (4). This suggests that the previous Δt_f value is still the optimal one. Together with the choice of the same rotation period P_rot, this has the advantage of yielding spot maps directly comparable with those obtained by Bonomo & Lanza (2012). The remarkable sensitivity of the χ² to Δt_f is caused by the best fits of some individual time intervals that become significantly worse when adopting a longer Δt_f, likely as a consequence not only of the short spot lifetimes, but also of the difficulty experienced by the numerical optimization routine in reaching the same minimum χ² when the number of datapoints is increased.

5. Results

5.1. Light-curve models

The best fit of the ARC2 light curve without regularization (λ = 0) was computed with the parameters in Table 1. The minimum and maximum of the residuals are −0.00110 and 0.00160, respectively, with an arithmetic mean of 1.751 × 10⁻⁷ in relative flux units. The best fit to the distribution of the residuals with a Gaussian has a mean of −3.428 × 10⁻⁶ and a standard deviation σ_ARC2 0 = 2.333 × 10⁻⁴. This value is close to the photometric accuracy of the datapoints indicating that the spot model is generally able to fit the light curve down to the level of the photon shot noise (cf. Sect. 2).

The composite regularized best fit obtained with the ARC2 light curve is plotted in Fig. 4. The light curve and the best fit have been normalized to the maximum of the light curve. The intervals plotted in the three panels of Fig. 4 have different durations because long gaps have been excluded from our plots with the exception of that in the bottom plot that was left to avoid a fourth plot with too short an interval. The minimum and maximum residuals are −0.00163 and 0.00175, respectively, with a mean of −1.325 × 10⁻⁵ in relative flux units. The distribution of the residuals is plotted in Fig. 5. The best fit with a Gaussian has a mean of −1.412 × 10⁻⁵ and a standard deviation of 2.689 × 10⁻⁴. The negative mean of the residuals is a consequence of the regularization that tends to reduce the spotted area as much as possible leading to a best fit systematically higher in flux than the data points. The convergence criterion is verified within 7% in all the cases, with more than 75% of the individual intervals verifying it within 2%.

Fig. 4. Top panel: light curve detrended with the ARC2 pipeline and fitted with our composite maximum-entropy regularized spot model and the parameters listed in Table 1. The observed flux, normalized to its maximum value, is plotted vs. time (filled dots) and the best fit is superposed (solid line). Lower panels: residual of the regularized best fit vs. time (filled dots).

Fig. 5. Histogram of the distribution of the residuals of the regularized best fit to the light curve in Fig. 4 (solid red histogram) and its best fit with a Gaussian (solid green line). The vertical dotted line marks the zero value.

The best fit of the light curve obtained with the de-trending approach of Bonomo & Lanza (2012) and without regularization (λ = 0) was computed with the parameters listed in Table 1. The residuals of the best fit range between −0.00230 and 0.00241 with a mean of 3.491 × 10⁻⁷ in relative flux units. The largest residuals are found close to data gaps and are probably associated with the residual systematic trends before and after the gaps that the approach by Bonomo & Lanza (2012) is not capable of correcting as efficiently as the ARC2 pipeline. The best fit with a Gaussian to the distribution of the residuals has a mean of −2.631 × 10⁻⁶ and a standard deviation σ_BL 0 = 2.504 × 10⁻⁴ in relative flux units.

The composite regularized best fit of the light curve obtained with the approach of Bonomo & Lanza (2012) is plotted in Fig. 6. The light curve and the best fit have been normalized to the maximum flux of the light curve. The minimum and the maximum of the residuals are −0.00603 e 0.00260, respectively, while the mean of the residuals is −1.326 × 10⁻⁵ in relative flux units. The distribution of the residuals is plotted in Fig. 7 together with a Gaussian best fit with a mean of −1.415 × 10⁻⁵ and a standard deviation of 3.067 × 10⁻⁴. While the mean is comparable with that of the fit to the ARC2 light curve, the standard deviation is about 10% larger as indicated by the greater fraction of relatively larger residuals. The convergence criterion is verified with a maximum deviation of 6% for all the individually fitted intervals of duration Δt_f = 8.733 days, with more than 75% of the intervals verifying the criterion within 3%.

Fig. 6. As in Fig. 4 but for the light curve de-trended with the approach of Bonomo & Lanza (2012).

Fig. 7. As in Fig. 5 but for the light curve de-trended as in Bonomo & Lanza (2012) and plotted in Fig. 6.

In contrast with the distribution of the residuals of the light curve fitted by Bonomo & Lanza (2012), which showed a positive tail in excess of the Gaussian best fit, the present distribution is well fitted by a Gaussian and has a symmetric shape. This difference is likely due to the improved correction of the residual systematic trends in the latest release of Kepler data.

It is interesting to explore possible periodicities in the time series of the residuals of the light-curve best fits. To this aim, we analyse the residuals of the unregularized best fit because they have no systematic bias introduced by the regularization. In Fig. 8, we plot the Generalized Lomb-Scargle periodogram (GLS; see Zechmeister & Kürster 2009) of the residuals of the light curve de-trended with the approach of Bonomo & Lanza (2012), chosen because it has more datapoints and less gaps than the ARC2 light curve. Given that the removal of the transits may affect the power spectrum, we checked the spectral window of the time series finding low side lobes at frequencies of ∼2 and ∼3.25 day⁻¹ that are away from the orbital frequency of 0.672 day⁻¹, showing that the transit removal does not significantly affect our periodogram. The relative duration of the discarded intervals is in fact 12.7% of each orbital period, which explains why the impact on the spectral window is limited when applying the GLS periodogram, specifically designed to treat time series with gaps.

Fig. 8. Generalized Lomb-Scargle periodogram of the residuals of the unregularized best fit to the light curve de-trended with the approach of Bonomo & Lanza (2012). The vertical red dashed line indicates the orbital period of Kepler-17b, while the red dotted line indicates the synodic period (see the text).

We find that our spot modelling accounts for the light modulations with timescales longer than ∼3.5 days as indicated by the almost complete disappearance of power in the periodogram for periods longer than that limit. There is a peak at the orbital period with a false-alarm probability of ∼10⁻⁸ as estimated with the analytic formula of Zechmeister & Kürster (2009). This suggests a variability with the orbital phase that can be attributed to light reflection and secondary eclipses as discussed by Bonomo et al. (2012). However, we are not in the position to characterize such effects because we are analysing only long-cadence data the time sampling of which is too coarse to give precise information on these phenomena, in particular on the secondary eclipses. We also indicate in Fig. 8 the synodic period P_syn calculated as , where we adopted a mean rotation period P_rot = 12.01 days. Variability with the synodic period could be an indication of magnetic star-planet interactions (Lanza 2008, 2012). In the present case, the peaks close to P_syn are not particularly prominent suggesting that there is no detectable effect of the planet on stellar variability in the present dataset. Similar results are obtained from the analysis of the residuals of the regularized best fits, but they are not shown here.

5.2. Longitude distribution of active regions, comparison with occulted spots, and differential rotation

In Fig. 9 we plot the distribution of the filling factor of the starspots f versus the longitude and the time for the regularized spot maps obtained with the ARC2 light curve. The origin of the longitude is at the meridian pointing toward the Earth on BJD_TDB 2454964.512 and the longitude is increasing in the same direction as the stellar rotation and the orbital motion of the planet. The reference frame is rotating with the star with a period P_rot = 12.01 days as in Bonomo & Lanza (2012). In the following, we measure the time t starting from a reference epoch, that is, we introduce a time t′ = t − 2454900.0 measured in BJD_TDB.

Fig. 9. Distribution of the spot filling factor vs. the longitude and time as derived by our maximum-entropy spot model of the long-cadence ARC2 light curve. The maximum of the filling factor is indicated by the yellow-orange colour, while the minimum by dark blue (see colour scale in the lower right corner). We note that the longitude scale is repeated beyond the [0° ,360°] interval to better follow the migration of the spot features. White circles mark the longitude and time of the spots occulted by the planet during transits as detected by modelling transit profiles observed in short cadence (Valio et al. 2017). Their size is proportional to their flux deficit D as defined in Sect. 5.2 (see the text). Data gaps without enough observations to compute a spot map along the individual Δtf intervals are indicated by a black bands. We note that short-cadence transit data are available during the first gap of the light curve extracted by the ARC2 pipeline.

Spots occulted by the planet during transits have been modelled by Valio et al. (2017) by analysing the short-cadence Kepler photometry whose time series began after the first transit detection, that is, approximately after t′ ∼ 290 days. The belt covered by the planet during its transits extends from −7.0° to 8.4° in latitude, assuming that the transit chord covers the northern hemisphere of the star, given that the impact parameter is b = 0.012 ± 0.010 and the radius of the planet 0.1335 ± 0.0001 stellar radii according to Maxted (2018). A similar result is obtained with the transit model by Valio et al. (2017), who took into account the systematic effects of the spots on the transit profile obtaining a slightly larger radius than that of Maxted, that is, 0.138 ± 0.001 and an impact parameter b = 0.10 ± 0.01 stellar radii, leading to an occulted belt between −3.6° and 13.8° in latitude, again assuming a transit chord on the northern hemisphere. In Fig. 9, the occulted spots are marked as white open circles the size of which is proportional to their flux deficit defined as , where r_s is the radius of the spot as derived by the duration of its occultation and I_spot/I the ratio of its specific intensity to the unperturbed intensity as derived from the height of the photometric anomaly produced by the spot itself (the “bump” along the residual transit profile).

The map of the distribution of the spot filling factor as derived from the light curve de-trended with the approach of Bonomo & Lanza (2012) is displayed in Fig. 10 together with the spots occulted by the planet during its transits. The ARC2 pipeline discarded the datapoints of an interval centred around t′ ∼ 400 days that were instead retained by the simpler de-trending algorithm applied by Bonomo & Lanza (2012), therefore the corresponding map shows only three major gaps.

Fig. 10. As in Fig. 9 but for the spot modelling of the light curve de-trended with the method of Bonomo & Lanza (2012).

The two maps in Figs. 9 and 10 are remarkably similar. This is due to the fact that the maps are based on best fits to short individual time intervals of Δt_f = 8.733 days, along which the light variations are comparable in the two time series, thus the different trends of the light curves on longer timescales do not strongly affect the spot distributions. However, we see that the map based on the ARC2 light curve shows a slight preference for a greater filling factor, in particular for values closer to the maximum, as indicated by the more extended red and orange areas in Fig. 9. In Appendix B and Fig. B.1, we show an enlargement of Fig. 10 to better show the migration of the spot pattern and the association between spots mapped from the out-of-transit light curve and spots occulted during transits.

Considering Figs. 9 and 10, during the first time interval (t′ ≲ 550), two main active longitudes are apparent, one beginning at ≈50° longitude and slowly migrating to ≈0°, the other beginning at ≈200° and staying approximately stationary in the adopted reference frame. Another active longitude appears in between, characterized by a remarkable intermittency. Individual starspots have a duration ranging from ≈10 to ≈40−50 days, which is remarkably shorter than the duration of the active longitudes (cf. Fig. B.1). This pattern is very similar to that mapped by Bonomo & Lanza (2012) using the Kepler long-cadence data available at that time (cf. their Fig. 4).

To measure the association between the longitude distribution of the spots as derived from the out-of-transit light curves and the distribution as mapped from the occultations during transits, we define the cross-correlation coefficient ρ_cc as:

(5)

where s_oot is the out-of-transit longitude distribution of the spot filling factor f and s_occ the distribution of D of the occulted spots, both mapped onto N_L = 20 equal longitude bins of 18°; the overbar indicates the mean value; and ℓ ∈ [ − 10, 10] is the lag index, the longitude lag being given by Δλ = 18° × ℓ. The distributions s_occ and s_oot are treated as circular datasets, that is, they repeat themselves beyond an interval of 360°. Because Δt_f <  P_rot, we consider the mean of two consecutive out-of-transit spot distributions and smooth them to a resolution of 54° to derive s_oot. Similarly, to have a complete longitude coverage along the chords occulted by the planet, we average the distributions of the occulted spots along four consecutive transits to compute s_occ. The association between the two distributions is measured by the cross-correlation coefficient at zero lag, that is, ρ_cc(0), which is plotted against time in Fig. 11 for our two light curves. The correlation is zero for t′ ≲ 290 because no transits were observed. We note that ρ_cc(0) is equal to the Pearson linear correlation coefficient r between s_oot and s_occ as defined in Sect. 14.5 of Press et al. (2007). Estimating the significance of the correlation, that is, the probability of obtaining the given r or a larger one in the case of a chance association, is difficult because the statistical distributions of the correlated variables are in general not known. A good alternative is to resort to the Spearman or rank-correlation coefficient r_S for which an analytic evaluation is possible (cf. Press et al. 2007, Sect. 14.6). Therefore, we can use r_S to compute the significance of ρ_cc(0).

Fig. 11. Cross-correlation coefficient at zero lag (see text) between the distributions of the starspots as obtained from the ME models of the out-of-transit light curves and the spots occulted during transits as mapped by Valio et al. (2017). The considered ME spot longitude distributions are those obtained from the light curve of Bonomo & Lanza (green diamonds) and the ARC2 light curve (orange triangles).

We can apply Eq. (5) also to evaluate the migration of the spot pattern between two consecutive out-of-transit spot longitude distributions. Specifically, the longitude lag Δλ that maximizes the cross-correlation ρ_cc between consecutive out-of-transit spot distributions can be used to quantify the longitude migration of the spot pattern occurring between them. This migration is assumed to be produced by the differential rotation when the most prominent spots are not rotating with the period of 12.01 days assumed for the reference frame. In Fig. 12, the migration rate obtained from the derivative of Δλ as a function of time is plotted against time itself.

Fig. 12. Migration rate between consecutive spot pattern distributions as derived from the Bonomo & Lanza light curve (green diamonds) and the ARC2 light curve (orange triangles). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).

Considering Figs. 9 and 10, we see that several spots occulted during transits are concentrated around the active longitudes found by the out-of-transit spot modelling with the coincidence being better for the two more persistent longitudes and during some time intervals, in particular when ρ_cc(0) ≳ 0.25. This coincidence provides independent confirmation of our spot modelling approach and indicates that some of the occulted spots and the active longitudes mapped by the out-of-transit light curve are at a similar latitude. The migration of the trails of the occulted spots towards positive longitudes indicates that they are rotating slightly faster than the main active longitude with a period of ∼11.90 ± 0.04 days. We note that the longitude resolution of the maps of the occulted spots can be as small as a few degrees, while that of the spots mapped from the out-of-transit light curve reaches only ≈50° in the best cases, thus accounting for the lack of a complete coincidence between the two maps (cf. Lanza et al. 2007; Silva-Valio & Lanza 2011) and the relatively low significance of the correlation that ranges from ∼0.1 to ∼0.4 for ρ_cc(0) ∼ 0.25 as derived from the analytical method introduced above.

The best correspondence between the out-of-transit and the occulted spot distributions is found for t′ ≈ 300, 900, and 1100 as indicated in Fig. 11 by ρ_cc(0) ≳ 0.35, which corresponds to an analytical significance better than 0.1. At those times, we see active longitudes in Figs. 9 and 10, while the migration rate of the out-of-transit spot distributions is close to zero or fluctuates between zero and +1.5 deg day⁻¹ in Fig. 12, indicating that the out-of-transit light curve is dominated by low-latitude spots, mostly occulted during transits. On the other hand, when ρ_cc(0) ≲ 0.2, active longitudes are not apparent in the adopted reference frame and we see mostly negative migration rates as, for example, for 750 ≲ t′ ≲ 850 or 1400 ≲ t′ ≲ 1500 (cf. Figs. 12 and B.1). The minimum migration rate is of −(3–4) deg day⁻¹, if we consider couples of consecutive similar measurements and exclude more extreme, isolated values which have smaller correlation coefficients as indicated by the smaller sizes of their plotted symbols. This minimum migration rate corresponds to a rotation period of 13.35−13.85 days that gives a relative differential rotation of ΔP_rot/P_rot ≃ 0.14 ± 0.05 considering an error of ±20° in the measurement of the lag between successive distributions and basing the determination on two consecutive measurements. The differential rotation is solar-like, that is, the lower latitudes rotate faster than the higher latitudes as we deduce from the occultation of the faster rotating spots during transits.

The regular increase of the migration rate as observed in Fig. 12 during 750 ≲ t′ ≲ 850 or 1400 ≲ t′ ≲ 1500 can be interpreted as a consequence of the evolution of the high-latitude spots that produce the negative values of the rate itself. When those spots first appear, they dominate the correlation between successive spot distributions producing the largest negative deviations from the angular velocity corresponding to the rotation period of the reference frame of 12.01 days. Their subsequent decay makes the spots at lower latitudes increasingly more relevant in the cross-correlation until the migration rate crosses the zero value and becomes positive when the latter completely dominate the cross-correlations. From the duration of such phases of ∼100 days, we deduce that high-latitude spots have a maximum lifetime of about three months. Alternatively, individual high-latitude spots may have a shorter lifetime (cf. Fig. B.1), while it is the latitude of the activity belt that steadily migrates towards the equator as we observe in the solar cycle. In this case, we can trace three cycles in Fig. 12 separated by ∼400 and ∼600 days, respectively. These periods agree with the period in the total area of the occulted spots as found by Estrela & Valio (2016) who reported a modulation with a period of 490 ± 100 days.

In Fig. 12, we see some very fast changes in the migration rate of the starspots. This is similar to the observations of the solar rotation period as derived from disc-integrated tracers such as the chromospheric Ca II H&K lines. In those time series, the onset of a new cycle is marked by an abrupt increase of the period as activity disappears at low latitudes and re-appears at higher latitudes (e.g. Donahue & Keil 1995; Hempelmann & Donahue 1997).

We explore the effect of varying the parameters of our spot modelling in Appendix A finding that decreases in the contrast of the spots c_s = I_spot/I or in the facular-to-spotted area ratio Q have the highest impact on our measurement of the differential rotation, reducing its amplitude by approximately a factor of two for the extreme values of those parameters.

5.3. Variation of the spotted area

Our spot modelling of the out-of-transit light curve allows us to determine the variation of the total spotted area versus time by integrating the filling factor over the longitude. The error is estimated from the photometric accuracy of the data points. The presence of gaps inside each individually fitted interval of duration Δt_f affects the total area because the maximum entropy regularization drives the solution towards the minimum spotted area compatible with the data, thus reducing the filling factor at the longitudes that are in view during the gaps in the light curves.

To reduce the impact of this effect on the variation of the total spotted area, we measured the presence of significant gaps along each interval Δt_f. We divided each interval into five equal subintervals and counted the number of data points into each subinterval n_i, with i = 1, …, 5 numbering the subinterval. A measure δ of the inhomogeneous distribution of the data points along the interval Δt_f is defined as δ ≡ [max(n_i) − min(n_i)]/max(n_i). In the case of the area values obtained from the light curve with the de-trending of Bonomo & Lanza (2012), the intervals with δ >  0.2 are discarded giving a total of 83 area measurements unaffected by the gaps over a total of 135 intervals. We note that as the intervals with δ >  0.2 have on average ∼4% less datapoints, they show a similar systematic decrease of the spot coverage values. With a mean spot coverage of ∼0.068, this amounts to a systematic difference of ∼2.7 × 10⁻³ that is comparable with the amplitude of the modulation we detect in the spot coverage itself (see below). Therefore, we choose δ = 0.2 as our acceptance threshold to avoid systematic errors comparable with the amplitude we intend to measure.

For the spot coverage obtained by the best fits to the ARC2 light curve, we had to relax the acceptance criterion and discarded only the intervals with δ >  0.25 to obtain a comparable number of coverage measurements, that is, 70 acceptable values over a total of 133 time intervals. We note that the ARC2 light curve has less data points and more gaps than the light curve de-trended with the approach of Bonomo & Lanza (2012). We therefore base our analysis mainly on the latter light curve.

The plot of the total spotted area versus the time for this light curve is shown in Fig. 13 together with the best-fitting sinusoid with a period of 47.906 days as obtained by the GLS periodogram. The area values put in phase with that period are shown in Fig. 14. The false-alarm probability (FAP) estimated with the analytic formula proposed by Zechmeister & Kürster (2009) is 0.0166. By performing 10 000 shuffling of the area values of the time series, we estimate a FAP of 0.05, which is not significantly different from the analytic estimate. The plot in Fig. 13 shows that the 48 day oscillation is particularly evident for t′ < 600 days, possibly around t′ ∼ 800 days, and in the latest part of the time series, that is for t′ > 1100−1200 days, although with a varying amplitude.

Fig. 13. Total spotted area as derived from the ME best fits to the light curve de-trended with the method of Bonomo & Lanza (2012) vs. time (green filled circles). The error bars have an amplitude of 3σ, where σ is the standard deviation as derived from the photometric accuracy of the datapoints. Values for the intervals with δ >  0.2 have been excluded. The best fitting sinusoid with a period of 47.906 days is superposed to the time series (red solid line).

Fig. 14. As in Fig. 13, but with the area values put in phase with the period of 47.906 days. The best-fitting sinusoid with a period of 47.906 days is superposed to better show the oscillation of the spotted area (red solid line).

To trace this varying periodicity, we apply a Morlet wavelet with the same parameters as in Bonomo & Lanza (2012). The amplitude of the wavelet versus period and time is plotted in Fig. 15. The relative maxima of the power are concentrated around a period of ∼50 days, although there is sometimes power at periods around 30 days. The gaps in the time series affect the Morlet wavelet and can account for the secondary maxima at different periods. A comparison of the wavelet map in Fig. 15 with that in Fig. 7 of Bonomo & Lanza (2012) shows the same overall structure, although the new map extends for a longer time interval and suggests a re-appearance of the periodicity at ∼50 days close to the end of the time series.

Fig. 15. Amplitude of the Morlet wavelet of the total spotted area variation in Fig. 13 vs. period and time. The amplitude was normalized to its maximum value. Different colours indicate different relative amplitudes from the maximum (orange) to the minimum (dark blue) as indicated in the colour scale in the right lower corner.

We explore the impact of the variation of our model parameters on the spot coverage in Appendix A considering different values of the limb-darkening coefficients, spot contrast c_s = I_spot/I, and facular-to-spotted area ratio Q. We find that the ∼48 day periodicity is retrieved in all cases, although its FAP becomes larger for non-optimal values of the Q parameter or the largest values of the spot contrast c_s.

The variation in the total spotted area as derived from the intervals of the ARC2 light curve with δ ≤ 0.25 is plotted in Fig. 16 together with the GLS best-fitting sinusoid. This sinusoid has a period of 48.202 days, but its analytic FAP is 0.1667, while the FAP from 10 000 shuffling is 0.40, likely as a consequence of the lower number of data points and larger fluctuations from one interval to the next. The corresponding Morlet wavelet map is plotted in Fig. 17 and its overall aspect is similar to that of the map in Fig. 15, reinforcing the case for a periodicity in the spotted area of ∼50 days during the first half of the time series. In the second half, the wavelet power is split among several different periodicities, likely as a consequence of the gaps in the time series and a long-term modulation with a period of several hundred days appearing for t′ ≳ 800 days. Such a modulation is not observed in the area time series obtained from the light curve de-trended with the approach of Bonomo & Lanza (2012) that fits and removes a parabolic trend within each quarter, thus filtering out the variations on timescales comparable with the quarter duration of ∼90 days (cf. the photometric time series in Figs. 1, 4, and 6).

Fig. 16. As in Fig. 13, but for the area values derived by the ME best fit to the ARC2 light curve. Intervals with δ >  0.25 have been discarded. The period of the GLS best fitting sinusoid is 48.202 days. We highlight the different scale on the y-axis.

Fig. 17. As in Fig. 15, but for the time series in Fig. 16.

Unfortunately, Kepler data are not useful to search for long-term changes of the mean light level of Kepler-17 to confirm the activity cycle of ∼400−600 days suggested by the different regimes of spot longitude migration (cf. Sect. 5.2). Some hint of a long-term variation of the spotted area may be apparent in Fig. 16, but the lack of a complete correction for the systematic variations from one quarter to the next hampers our attempts to confirm this result. We note that the variation in the spotted area plotted in Figs. 13 and 16 refers to the spots that are unevenly distributed in longitude, because the amplitude of the rotational modulation is insensitive to uniformly distributed spots. In other words, if the cycle of ∼400−600 days is associated with a variation of the area of spots almost uniformly distributed in longitude, it can go undetected in those plots and only a long-term variation of the mean light level would reveal its presence. With a chromospheric index (Bonomo et al. 2012, Sect. 4.2.2), Kepler-17 is at the boundary separating very active stars with predominantly non-axisymmetric spot distributions and active longitudes from less active rotators with an almost uniform distribution of spots in longitude. This may account for the complex behaviour of Kepler-17, probably indicating phenomena common to both kinds of stars (Lehtinen et al. 2016).

5.4. Star–planet interaction: tides

We did not find evidence of a light modulation associated with a possible star–planet interaction in the residuals of our spot models (cf. Sect. 5.1). However, some effect of the planet on the star is expected because of its mass of ∼2.5 Jupiter masses and its proximity. Tides raised on the star by the planet are an example of such an interaction. Thanks to our determination of stellar rotation, we can derive information on the tidal dissipation inside the G2V star Kepler-17 that is useful to model the evolution of the rotation itself. This can be applied to evaluate the activity level of the star in the past, thus providing information for models of planetary evolution and evaporation (e.g. Murray-Clay et al. 2009), as well as the confidence of stellar age based on gyrochronology.

The orbital angular momentum of the planet is approximately four times the stellar spin angular momentum, while the total angular momentum of the system is only ∼0.6 of that required to reach a synchronous final state, in the hypothesis that the total angular momentum of the system is conserved (Hut 1980). However, the stellar magnetized wind produces a steady loss of angular momentum from the system that accelerates the shrinking of the planetary orbit until the planet is engulfed by the star (cf. Damiani & Lanza 2015).

We investigate the evolution of the stellar rotation and the orbital semimajor axis by applying the simple model of Lanza & Mathis (2016) that includes the wind braking of the stellar rotation using a Skumanich-type law. The star is assumed to rotate rigidly and the strength of the tidal interaction is parameterized by the stellar modified tidal quality factor Q′ (Zahn 2008). A stronger interaction implies a faster dissipation of the kinetic energy of the tides and is parameterized by a smaller value of Q′. We include in the model the evolution of the radius of the star calculated by means of the EZweb interface³ because the tidal torque is proportional to (R/a)⁶, where R is the stellar radius and a the orbit semimajor axis (Zahn 2008). We consider a model for a main-sequence star of mass 1.095 M_⊙ and metal abundance Z = 0.03 because it has a radius of 1.06 R_⊙ at the estimated age of Kepler-17, that is 1.8 Gyr (cf. Bonomo et al. 2012). This value of the mass is different from that derived by Désert et al. (2011) and Bonomo et al. (2012) by fitting different stellar evolution models to the position of the star in the mean density–effective temperature diagram, but it is still within ∼1σ from their mass estimates. We prefer to adopt a stellar evolution model that fits the radius at the putative age of 1.8 Gyr rather than the estimated mass because the radius evolution has a much stronger impact on the tidal evolution of the system, while the ratio of the stellar mass to the planetary mass stays fixed in our model. We assume a circular orbit because tides inside the planet damp any initial eccentricity on timescales of 10−100 Myr, which is much shorter than the age of the star. The obliquity of the planetary orbit is assumed to be zero following the discussion in Sect. 7.2.1 of Désert et al. (2011).

The strength of the tidal interaction in star-planet systems is unknown and there are theoretical reasons to believe that it depends on the ratio between the tidal frequency⁴ and the rotation frequency of the star (Ogilvie & Lin 2007). Observational estimates have been performed only with statistical methods that do not provide information on individual systems, but only an indication of the mean values of Q′ in different regimes (e.g. Bonomo et al. 2017; Collier Cameron & Jardine 2018). Therefore, we model the evolution assuming two constant values for Q′, that is 10⁷ and 10⁸, with a preference for the latter from a theoretical point of view. Specifically, the tidal frequency and the rotation frequency in Kepler-17 are always sufficiently separated so as to avoid the excitation of inertial waves, that is, , thus leading to a weaker tidal interaction between the star and the planet (cf. Ogilvie & Lin 2007). Conversely, in the regime, the excitation of inertial waves that are strongly dissipated in wave attractors would lead to a decrease of Q′ by two to three orders of magnitude (cf. Rieutord et al. 2001; Ogilvie & Lin 2007; Goodman & Lackner 2009).

In Fig. 18 we plot the evolution of the stellar spin and orbit semimajor axis together with the evolution of the stellar radius adopted to compute the first two quantities. We also plot the evolution of the stellar rotation without any tidal torque, that is, under the action of the Skumanich-type wind braking only. We assume that our model can be applied for ages later than ∼0.6 Gyr because younger ages may still show the effects of the initial conditions and of an incomplete internal core–envelope coupling. The stellar age is assumed to be 1.8 Gyr and the presently measured mean rotation period of 12 days is imposed at that age to all our models.

Fig. 18. Top panel: mean rotation period of Kepler-17 vs. time for our tidal evolution model computed with Q′ = 107 (dotted line) and Q′ = 108 (solid line). The evolution of the stellar spin without any tidal torque and assuming a period of 12.01 days at an age of 1.8 Gyr is plotted as the dot-dashed line. The vertical dashed line indicates the adopted age of the star. Middle panel: as in the top panel, but for the orbit semimajor axis. Lower panel: radius of the star adopted to compute the tidal evolution vs. time.

The spin evolution is dominated by magnetic braking up to about ∼2 Gyr for Q′ = 10⁸, while it deviates remarkably from the Skumanich law for Q′ = 10⁷ all along the evolution. This is because the angular momentum exchange due to tides dominates over the loss of angular momentum by the wind braking. This scenario has led to a remarkably longer rotation period of the star in the past than in the case with Q′ = 10⁸ because the planet spun up the star through its tidal interaction, while its orbit decayed. However, even for Q′ = 10⁸, the slope of the braking law is reduced by the tidal interaction with tides that counteract magnetic braking leading to a less steep variation of the rotation period. This suggests that the age estimated by means of standard gyrochronology, which ranges between 1.0 and 1.4 Gyr (Barnes 2010; Barnes et al. 2016), is not accurate for Kepler-17 because of the tidal spin-up induced by its massive and close-by planet. Moreover, the planet could also affect the efficiency of the stellar wind (cf. Cohen et al. 2010; Lanza 2010). Therefore, the present results support the adoption of a modified gyrochronology relationship to evaluate the age of Kepler-17 (cf. Lanza 2010) and provide a rotation evolution scenario in agreement with the age of 1.8 Gyr as estimated by Bonomo et al. (2012).

The expected survival time of the planet is estimated to be ≈0.35 Gyr in the case of the stronger tidal interaction and ≈2.4 Gyr in the case of the weaker interaction. The orbital decay is mainly ruled by the increase of the stellar radius along its main-sequence evolution owing to the remarkable dependence of the tidal torque on R/a.

6. Discussion and conclusions

We have analysed the activity of Kepler-17 using two different approaches to de-trend the systematics present in Kepler time series. We confirm that the PDC pipeline introduces an overcorrection of the subtle light modulations produced by solar-like faculae in late-type stars making PDC timeseries not recommendable to produce spot maps by light-curve inversions (see also Aigrain et al. 2017). Conversely, the de-trending by the ARC2 pipeline suffers from many fewer problems and provides results that are comparable to those derived by the simpler approach by Bonomo & Lanza (2012), at least in the case of this target with high-S/N data and a high level of activity. Nevertheless, the evaluation of the total spotted area from the ARC2 time series can still be affected by some residual trends on timescales ranging from ∼10 to ∼90 days, probably due to the use of a limited number of CBVs to correct the light curve.

Our results provide an extended comparison of the maps obtained from the in-transit and out-of-transit light modulations, that is from two independent datasets and methods. The good correspondence found during certain time intervals indicates that the adopted maximum-entropy approach is capable of reconstructing the overall starspot distribution in longitude and time and gives support to the existence of active longitudes. A similar comparison was made by Silva-Valio & Lanza (2011) validating the spot models for CoRoT-2, but this was based only on ∼150 days of data, while the present comparison is more extended and shows the effects of a likely activity cycle of ∼400−600 days.

The active longitude around ∼200° is remarkable because it lasts for at least 1400 days, although its level of activity is continuously changing as indicated by the varying spot filling factor. Similar long-lived active longitudes are commonly observed on very active rotators such as young solar-type stars (Lehtinen et al. 2016) or the subgiant members of the close active binaries II Peg (Rodonò et al. 2000) or HR 1099 (Lanza et al. 2006). Space-borne photometry has revealed an active-longitude phenomenology similar to that of Kepler-17 in, for example, CoRoT-2 (Lanza et al. 2009) or CoRoT-6 (Lanza et al. 2011) or in the M dwarf GJ 1243 (Davenport et al. 2015). Active longitudes have also been proposed for the Sun, although they are not as evident as in more active and rapidly rotating solar-like stars (e.g. Usoskin et al. 2007). Recent hydromagnetic dynamo models have provided some insight into the physical mechanisms that could produce such active longitudes (Weber et al. 2013).

We found some indication of an activity cycle of ∼400−600 days from the latitudinal migration of the starspots in Kepler-17. Estrela & Valio (2016) found evidence of a similar periodicity in the area of the spots occulted during transits. This cycle could be similar to the solar 11-year cycle, although of remarkably shorter duration. Short activity cycles have recently been found in ι Horologii, an F8V star that hosts a giant planet on a 300 day orbit and shows a cycle of ∼585 days (Sanz-Forcada et al. 2013); and in the young (∼1 Gyr) G1.5V star HD 30495 that has a rotation period of ∼11 days and shows two chromospheric cycles, one of 620 ± 150 days and another of 12 ± 3 years. The short-term modulation is intermittent and does not appear to be related to the longer-term cycle (Egeland et al. 2015). The case of HD 30495, which is similar to Kepler-17 in effective temperature and rotation period, suggests that such short cycles may be characteristic of young Sun-like stars. Other possible examples are HD 76151, a G3V star with a rotation period of ∼15 days and a cycle of ∼920 ± 10 days; and HD 190406, a G1V star with a rotation period of ∼14 days and two cycles of ∼950 ± 10 days and ∼17 years (see Baliunas et al. 1995, 1996).

In addition to this possible activity cycle in Kepler-17, we find marginal evidence for an oscillation of the total spotted area with a period of ∼48 days, previously reported by Bonomo & Lanza (2012). The mean total area of the occulted spots found by Valio et al. (2017) is ∼6 ± 4%, agreeing with our mean spotted area. However, these latter authors do not find evidence of the 48 day periodicity, probably because the large variations in the area of the individual spots hamper its detection, making only the periodicity at the mean rotation period clearly apparent (cf. Fig. 6 in Valio et al. 2017). The 48 day spotted area modulation is not associated with a migration of the main latitude of spot formation and is reminiscent of the so-called Rieger cycles in the Sun (e.g. Oliver et al. 1998; Zaqarashvili et al. 2010; Gurgenashvili et al. 2017). These cycles have been attributed to Rossby-type waves propagating in the solar interior that modulate the toroidal magnetic field responsible for the formation of the spots (cf. Zaqarashvili et al. 2010; Gurgenashvili et al. 2016; Zaqarashvili 2018). Similar cycles have been observed in, for example, CoRoT-2 (Lanza et al. 2009), and in some young late-type stars with an age between 4 and 95 Myr investigated by Distefano et al. (2017). Short-term cycles, possibly of Rieger type, have also been investigated using CoRoT (Ferreira Lopes et al. 2015) and Kepler (Arkhypov et al. 2015) time series.

The amplitude of the latitudinal differential rotation in Kepler-17 derived from our spot modelling is only a lower limit because we do not know the latitudes of the spots rotating with different periods. Moreover, the rotation of the overall spot pattern can be different from that of individual starspots as suggested in the case of CoRoT-2 by Fröhlich et al. (2009) and in our case by Fig. B.1. Therefore, our estimate of ∼14 ± 5% relative amplitude of the differential rotation should be taken with some caution. Moreover, this amplitude is reduced to about 8 ± 5% for some values of our spot-modelling parameters (cf. Appendix A). Nevertheless, the solar-like character of the differential rotation, that is, the faster rotation of the equator with respect to the higher latitudes, is well established thanks to the comparison with the spots occulted during the transits that are certainly located at low latitudes. Using the measurements of the rotation periods of those spots and assuming a solar-like differential rotation profile, Valio et al. (2017) estimated a relative pole-equator angular velocity difference of ΔΩ/Ω ≃ 8.0 ± 0.9%, close to the lower limit of our determination.

The range of rotation periods derived from the modulation of the chromospheric flux in late-type stars by Donahue et al. (1996) suggests ΔP_rot/P_rot​​​​ ∼ ​​​​0.13 for Kepler-17, although HD 190406, which has similar spectral type and mean rotation period, shows ΔP_rot/P_rot = 0.21. The large statistical sample considered by Reinhold & Gizon (2015) shows a relative amplitude up to 0.1−0.2 for G-type stars with the mean rotation period of Kepler-17, based on the analysis of the photometric time series of Kepler targets. Recent theoretical models by Brun et al. (2017) predict a solar-like differential rotation for Kepler-17, that is, with the equator rotating faster than the poles. Its fluid Rossby number (see Brun et al. 2017, Eq. (33) and Fig. 22) is R_of ∼ 0.6 giving an expected relative amplitude of the differential rotation between the equator and 60° latitude of ΔΩ/Ω ≃ 0.2. Therefore, we conclude that the amplitude of the latitudinal differential rotation of Kepler-17 as derived from our analysis is in agreement with both observations and theoretical models for stars of similar spectral type and mean rotation period.

In the Kepler-17 system, the tidal interaction between the planet and the star is likely to be relevant. We find that it is capable of modifying the evolution of the stellar rotation by counteracting the braking by the stellar wind. Even if we assume a weak tidal coupling (Q′ ∼ 10⁸), which is favoured in our model, the stellar spin up is significant and makes it impossible to derive a precise age for the star by means of gyrochronology (Barnes 2007, 2010). This is probably the case of several stars hosting massive close-by planets as discussed by Ferraz-Mello et al. (2015) and Damiani & Lanza (2015) for example. Nevertheless, an estimated age of ∼1.8 Gyr (Bonomo et al. 2012) is in agreement with a simple model of the tidal evolution of the system. The same model predicts a rotation period for the star of ∼7.5 days at an age of 0.6 Gyr and an almost constant orbit semimajor axis for the planet over the 0.6−1.8 Gyr time interval. This information can be used to evaluate the evolution of the high-energy stellar radiation flux (e.g. Sanz-Forcada et al. 2011), which affects the evaporation of the planet and is controlled by the rotation rate of the star itself (e.g. Schmitt 2010).

Acknowledgments

The authors are grateful to an anonymous referee for a careful reading of their manuscript and valuable comments that improved the presentation of their results. AFL is grateful for interesting discussions with the colleagues of the team on “Rossby waves in Astrophysics” led by Prof. T. Zaqarashvili and supported by the International Space Science Institute in Bern, Switzerland. AFL and ASB acknowledge support from INAF/Frontiera through the Progetti Premiali funding scheme of Italian Ministry of Education, University, and Research. YN gratefully acknowledges support for the Programa de Doutorado Sanduíche no Exterior (PDSE) from the Brasilian Federal Agency for Support and Evaluation of Graduate Education (CAPES) during his stage at INAF-Catania Astrophysical Observatory (Process no. 88881.134871/2016-01).

References

Appendix A: Effects of changing model parameters on the spot distributions

We explore the effects of varying our model parameters on our main results. Specifically, we assume theoretical limb-darkening coefficients in place of those derived from the fitting of the transits by Maxted (2018); or change the contrast of the spots c_s; or vary the facular-to-spotted area ratio Q. In this investigation, we change one parameter at a time to isolate its effects on our results. We present results for the light curve detrended according to the method of Bonomo & Lanza because it has less gaps than the ARC2 light curve and because the results are very similar. We focus on the migration rate of the spots and on the variation in their total coverage that affect our measurements of the differential rotation and of the activity level, respectively.

A.1. Varying the limb-darkening coefficients

In Fig. A.1, we plot the migration rate versus time as derived from the ME spot modelling with the limb-darkening coefficients obtained from the fitting of the transits according to Maxted (2018) (green diamonds) or from model atmospheres (orange triangles; see Müller et al. 2013, and discussion in Sect. 4). The migration rates of the spot distributions are closely comparable, if we consider only the values between −4 and +2 deg day⁻¹. Therefore, the amplitude and the sign of the differential rotation are not affected, if we adopt the theoretical limb-darkening coefficients.

Fig. A.1. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with the limb-darkening coefficients as derived by fitting the transits (green diamonds) or from model atmospheres (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).

The total spotted area shows a systematic difference of about −2.8% with oscillations not exceeding 0.4%, which do not affect our conclusions concerning possible activity cycles in Kepler-17 (cf. Fig. A.2). Specifically, the GLS periodogram of the modulation of the spotted area has its maximum at a period of 47.916 days, very close to that of 47.906 days obtained with our reference model; also, the FAP is similar.

Fig. A.2. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with the limb-darkening coefficients as derived by fitting the transits (green diamonds) or from model atmospheres (orange triangles). Lower panel: relative difference between the values of the area obtained with the two sets of limb-darkening coefficients.

A.2. Varying the spot contrast

We explore the effect of varying the spot contrast between the extreme values measured by Valio et al. (2017) by modelling spot occultations, that is, c_s ≡ I_spot/I = 0.38 and 0.72; for comparison, in the case of sunspot groups, c_s = 0.67. In Fig. A.3, we see that the spot migration rate has a less negative minimum value of about −2 deg day⁻¹ for c_s = 0.38, but a slightly greater positive maximum value, leading to a smaller amplitude of the relative differential rotation, that is, ΔP_rot/P_rot ∼ 0.08 ± 0.05 when considering a period of 11.90 days for the rotation of the occulted spots (cf. Sect. 5.2). On the other hand, we find ΔP_rot/P_rot ∼ 0.12 ± 0.05 when we consider c_s = 0.72 because a minimum migration rate of about −3 deg day⁻¹ is measured during several cross-correlations (cf. Fig. A.5).

Fig. A.3. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with a spot contrast cs = 0.55 (green diamonds) or cs = 0.38 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).

The total spot coverage changes in a systematic way showing a smaller area and a smaller amplitude of its modulation when the spots are darker, that is for c_s = 0.38 (cf. Figs. A.4 and A.6). However, the relative variations in the individual values of the coverage of ± 2−3% do not affect the main peak of the GLS periodogram that is at 47.916 and 47.962 days for c_s = 0.38 and 0.72, respectively. Only the value of the FAP as given by the formula of Zechmeister & Kürster (2009) is increased to 0.105 in the case of c_s = 0.72.

Fig. A.4. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with cs = 0.55 (green diamonds) or cs = 0.38 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different spot contrasts.

Fig. A.5. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with a spot contrast cs = 0.55 (green diamonds) or cs = 0.72 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).

Fig. A.6. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with cs = 0.55 (green diamonds) or cs = 0.72 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different spot contrasts.

A.3. Varying the facular-to-spotted area ratio

Finally, we explore the effects of the variation in the facular-to-spotted area ratio Q between the extreme values 1.0 and 4.0 that are well beyond the 95% joint confidence interval of this parameter as derived in Sect. 4 by the analysis of the ARC2 and of the Bonomo & Lanza light curves. The spot migration rate shows minimum negative values of about −2 deg day⁻¹ for Q = 1.0 (cf. Fig. A.7), while for Q = 4.0, the minimum is about −3 deg day⁻¹ (cf. Fig. A.9), corresponding to relative amplitudes of the differential rotation ΔP_rot/P_rot = 0.08 ± 0.05 and 0.12 ± 0.05, respectively. The variation of Q produces systematic changes in the longitudes of the model active regions because the relative contributions of dark spots and bright faculae depend on their positions with respect to the centre of the stellar disc at a given rotation phase (see Lanza et al. 2007).

Fig. A.7. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 1.0 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).

The total spot coverage is affected by Q because a larger spot area is required to counterbalance the effect of the larger faculae and reproduce the amplitude of the observed light modulation when Q is increased (cf. Figs. A.8 and A.10). In addition to this systematic variation, there are also fluctuations of relative amplitude of about ±2% with respect to the reference case with Q = 2.4 that do not affect our results on a possible short-term activity cycle. Specifically, we find the maximum of the GLS periodogram of the area time series at periods of 47.955 and 47.962 days for Q = 1 and Q = 4, respectively. Only the analytic FAP is increased to 0.06 and 0.12, respectively, likely as a consequence of the non-optimal values of Q adopted in those models.

Fig. A.8. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 1.0 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different values of Q.

Fig. A.9. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 4.0 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).

Fig. A.10. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 4.0 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different values of Q.

Appendix B: A detailed view of the filling factor map

In Fig. B.1, we show an enlargement of Fig. 10. It shows the changes occurring between t′ ≃ 750 and t′ ≃ 850 days when the spot pattern displayed an overall backward migration produced by a slower rotation, that is, a rotation period longer than 12.01 days, the period of the reference frame adopted to plot the spot map. This backward migration was clearly detected by cross-correlating successive distributions of the filling factors as obtained from the ME models of the out-of-transit light curve (see Fig. 12). However, we see in Fig. B.1 that individual longitudes show different migration rates also outside 750 ≲ t′ ≲ 850 days indicating that they are produced by spots at different latitudes. Therefore, the migration rate given by the cross-correlation is an average over the whole longitudinal distributions. For t′ ≲ 850 days, the correspondence between the spots as mapped by the out-of-transit light curve and those mapped from transit occultations is poor (the cross-correlation at zero lag ρ_cc(0) ≲ 0.15 in Fig. 11) suggesting that most of the former are located outside the occulted belt. Individual spots are short-lived with typical lifetimes of a few tens of days, while active longitudes where spots form and decay are long-lived – see the case of the active longitude around 200° − 250° that disappears at t′ ≈ 770 and re-appears for t′ ≳ 850 days (see also Sect. 5.2).

Fig. B.1. Enlargement of Fig. 10 showing the distribution of the spot filling factor (see the colour scale at the bottom right) vs. longitude and time as obtained by the ME modelling of the light curve de-trended as in Bonomo & Lanza (2012). The spots detected during transits by Valio et al. (2017) are overplotted as white circles the radius of which is proportional to their flux deficit D as defined in Sect. 5.2. The longitude scale goes beyond the interval [0° ,360°] to help us following the migration of the spots.

All Tables

All Figures

Fig. 1. Top panel: light curve of Kepler-17 as obtained by applying the de-trending as in Bonomo & Lanza (2012) (in green) and the light curve as provided by the ARC2 pipeline (in red). Each light curve is normalized to the median flux within each quarter after removing the planetary transits. The green datapoints have been shifted upward by 0.0625 for clarity. Bottom panel: difference between the two light curves. The quarters where only the first CBV has been applied to obtain the ARC2 light curve are indicated with blue datapoints, while the quarters where the first two CBVs have been applied are indicated in pink (see the text).

In the text

	Fig. 2. Ratio of the χ2 of the best fit to the entire ARC2 light curve to its minimum value vs. the parameter Q, i.e. the ratio of the area of the faculae to that of the dark spots in active regions. The horizontal dashed line indicates the 95% confidence level for determining the interval of acceptable Q values.
In the text

	Fig. 3. As in Fig. 2 but for the best fit to the light curve obtained with the de-trending of Bonomo & Lanza (2012).
In the text

	Fig. 4. Top panel: light curve detrended with the ARC2 pipeline and fitted with our composite maximum-entropy regularized spot model and the parameters listed in Table 1. The observed flux, normalized to its maximum value, is plotted vs. time (filled dots) and the best fit is superposed (solid line). Lower panels: residual of the regularized best fit vs. time (filled dots).
In the text

	Fig. 5. Histogram of the distribution of the residuals of the regularized best fit to the light curve in Fig. 4 (solid red histogram) and its best fit with a Gaussian (solid green line). The vertical dotted line marks the zero value.
In the text

	Fig. 6. As in Fig. 4 but for the light curve de-trended with the approach of Bonomo & Lanza (2012).
In the text

	Fig. 7. As in Fig. 5 but for the light curve de-trended as in Bonomo & Lanza (2012) and plotted in Fig. 6.
In the text

	Fig. 8. Generalized Lomb-Scargle periodogram of the residuals of the unregularized best fit to the light curve de-trended with the approach of Bonomo & Lanza (2012). The vertical red dashed line indicates the orbital period of Kepler-17b, while the red dotted line indicates the synodic period (see the text).
In the text

Fig. 9. Distribution of the spot filling factor vs. the longitude and time as derived by our maximum-entropy spot model of the long-cadence ARC2 light curve. The maximum of the filling factor is indicated by the yellow-orange colour, while the minimum by dark blue (see colour scale in the lower right corner). We note that the longitude scale is repeated beyond the [0° ,360°] interval to better follow the migration of the spot features. White circles mark the longitude and time of the spots occulted by the planet during transits as detected by modelling transit profiles observed in short cadence (Valio et al. 2017). Their size is proportional to their flux deficit D as defined in Sect. 5.2 (see the text). Data gaps without enough observations to compute a spot map along the individual Δtf intervals are indicated by a black bands. We note that short-cadence transit data are available during the first gap of the light curve extracted by the ARC2 pipeline.

In the text

	Fig. 10. As in Fig. 9 but for the spot modelling of the light curve de-trended with the method of Bonomo & Lanza (2012).
In the text

Fig. 11. Cross-correlation coefficient at zero lag (see text) between the distributions of the starspots as obtained from the ME models of the out-of-transit light curves and the spots occulted during transits as mapped by Valio et al. (2017). The considered ME spot longitude distributions are those obtained from the light curve of Bonomo & Lanza (green diamonds) and the ARC2 light curve (orange triangles).

In the text

	Fig. 12. Migration rate between consecutive spot pattern distributions as derived from the Bonomo & Lanza light curve (green diamonds) and the ARC2 light curve (orange triangles). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).
In the text

Fig. 13. Total spotted area as derived from the ME best fits to the light curve de-trended with the method of Bonomo & Lanza (2012) vs. time (green filled circles). The error bars have an amplitude of 3σ, where σ is the standard deviation as derived from the photometric accuracy of the datapoints. Values for the intervals with δ >  0.2 have been excluded. The best fitting sinusoid with a period of 47.906 days is superposed to the time series (red solid line).

In the text

	Fig. 14. As in Fig. 13, but with the area values put in phase with the period of 47.906 days. The best-fitting sinusoid with a period of 47.906 days is superposed to better show the oscillation of the spotted area (red solid line).
In the text

	Fig. 15. Amplitude of the Morlet wavelet of the total spotted area variation in Fig. 13 vs. period and time. The amplitude was normalized to its maximum value. Different colours indicate different relative amplitudes from the maximum (orange) to the minimum (dark blue) as indicated in the colour scale in the right lower corner.
In the text

	Fig. 16. As in Fig. 13, but for the area values derived by the ME best fit to the ARC2 light curve. Intervals with δ >  0.25 have been discarded. The period of the GLS best fitting sinusoid is 48.202 days. We highlight the different scale on the y-axis.
In the text

	Fig. 17. As in Fig. 15, but for the time series in Fig. 16.
In the text

Fig. 18. Top panel: mean rotation period of Kepler-17 vs. time for our tidal evolution model computed with Q′ = 107 (dotted line) and Q′ = 108 (solid line). The evolution of the stellar spin without any tidal torque and assuming a period of 12.01 days at an age of 1.8 Gyr is plotted as the dot-dashed line. The vertical dashed line indicates the adopted age of the star. Middle panel: as in the top panel, but for the orbit semimajor axis. Lower panel: radius of the star adopted to compute the tidal evolution vs. time.

In the text

	Fig. A.1. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with the limb-darkening coefficients as derived by fitting the transits (green diamonds) or from model atmospheres (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).
In the text

	Fig. A.2. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with the limb-darkening coefficients as derived by fitting the transits (green diamonds) or from model atmospheres (orange triangles). Lower panel: relative difference between the values of the area obtained with the two sets of limb-darkening coefficients.
In the text

	Fig. A.3. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with a spot contrast cs = 0.55 (green diamonds) or cs = 0.38 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).
In the text

	Fig. A.4. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with cs = 0.55 (green diamonds) or cs = 0.38 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different spot contrasts.
In the text

	Fig. A.5. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with a spot contrast cs = 0.55 (green diamonds) or cs = 0.72 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).
In the text

	Fig. A.6. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with cs = 0.55 (green diamonds) or cs = 0.72 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different spot contrasts.
In the text

	Fig. A.7. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 1.0 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).
In the text

	Fig. A.8. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 1.0 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different values of Q.
In the text

	Fig. A.9. Migration rate between consecutive spot pattern distributions as derived from the ME models of the Bonomo & Lanza light curve with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 4.0 (orange triangles; see our Sect. 4). The size of the symbols is proportional to the cross-correlation coefficient ρcc (cf. Eq. (5)).
In the text

	Fig. A.10. Upper panel: total coverage factor As of the starspots as derived from the ME models of the light curve of Bonomo & Lanza with facular-to-spotted area ratio Q = 2.4 (green diamonds) or Q = 4.0 (orange triangles). Lower panel: relative difference between the values of the area obtained with the two different values of Q.
In the text

Fig. B.1. Enlargement of Fig. 10 showing the distribution of the spot filling factor (see the colour scale at the bottom right) vs. longitude and time as obtained by the ME modelling of the light curve de-trended as in Bonomo & Lanza (2012). The spots detected during transits by Valio et al. (2017) are overplotted as white circles the radius of which is proportional to their flux deficit D as defined in Sect. 5.2. The longitude scale goes beyond the interval [0° ,360°] to help us following the migration of the spots.

In the text