© ESO 2021

1 Introduction

Stellar rotation is a fundamental physical parameter in stellar astrophysics and plays an important role in the formation and evolution of stars (Kraft 1967; Skumanich 1972; Kawaler 1988). Rotation controls stellar magnetism, mixing in the stellar interior, and tidal interactions in close binary systems. Moreover, the relationship between rotation and magnetic activity has important implications for the detectability and characterisation of planets orbiting solar-like stars (Maxted 2016; Collier Cameron 2018). Such stars make up the vast majority of the Kepler exoplanet targets and also allow us to investigate stellar variability due to magnetic activity. In particular, starspots and active regions in the stellar photospheres modulate the stellar flux on the rotation period (Walkowicz & Basri 2013).

The Kepler space mission found several multi-planet systems, one of which is the star orbiting Kepler-30 (KOI-806), which has nearly solar mass and radius, an effective temperature of 5500 ± 55 K, and a rotation period of ~16.0 days as measured using the Lomb–Scargle periodogram (Sanchis-Ojeda et al. 2012; Lanza et al. 2014). These properties make it a very interesting young solar analogue, motivating our choice to investigate its rotation and activity behaviour. Moreover, it is orbited by three planets named Kepler-30b, Kepler-30c, and Kepler-30d with masses of 9.2, 536 (~1.7 Jupiter masses), and 23.7 Earth masses, radii of 3.75, 11.98, and 8.79 Earth radii (Panichi et al. 2018), and orbital periods of 29.3, 60.3, and 143.3 days (cf. Sanchis-Ojeda et al. 2012), respectively. Based on the values of radius and mass, Panichi et al. (2018) estimated that Kepler-30b is a Neptune-like planet rather than a super-Earth, Kepler-30c is a Jovian planet, while Kepler-30d is classified as a Neptune-mass planet, with bulk densities given by ~0.96 g cm⁻³, ~1.71 g cm⁻³, and ~0.19 g cm⁻³, respectively. In this system, the orbits of the planets are on the same plane, almost perpendicular to the stellar spin, a behaviour similar to our Solar System. As a consequence, we observe Kepler-30 almost equator-on. This provides a fundamental constraint to reduce the degeneracies in the spot modelling of the star, an additional motivation for applying our approach (see Sect. 4). Using the spot modelling, Bonomo & Lanza (2012) and Lanza et al. (2019) performed an analysis for the star Kepler-17, a solar-like star younger than the Sun hosting a hot Jupiter that occults starspots during transits. In both studies, the authors analysed the activity and rotation of this star using starspots as tracers, comparing the longitudes of the spots mapped from the out-of-transit photometry with those of the spots occulted during transits to validate their spot-modelling method. The conclusions of these studies are that spot modelling would allow us: (i) to derive active longitudes where spots potentially form; (ii) to estimate the lifetime of active regions; (iii) to determine a lower limit for the stellar differential rotation, and; (iv) to evaluate short-term and long-term activity cycles depending on the length of the available time-series.

Inspired by the peculiar properties of the Kepler-30 system and the availability of almost four years of Kepler high-precision photometry,we decided to invest our efforts in an unprecedented analysis of this system. In particular, the present work presents a novel a joint analysis exploiting spot modelling and multifractal methods, both widely used in the astrophysical context, but never confronted with each other.

Over recent years, several methods have been applied to investigate the (multi)fractal structure of time-series in astrophysical systems (e.g. de Freitas et al. 2016, 2017, 2019a,b; Belete et al. 2018; de Franciscis et al. 2019). The estimation of local fluctuations and long-term dependency of astrophysical time-series is a problem that has recently been studied to understand the effect of long-memory processes caused by stellar rotation (de Freitas et al. 2019a,b). As discussed by de Freitas et al. (2013a), the global Hurst exponent (Hurst 1951) is a powerful (multi)fractal indicator for analysing CoRoT and Kepler time-series, which is mainly used to elucidate the persistence due to rotational modulation in the series themselves. Recently, a novel technique proposed byGu & Zhou (2010), known as the multifractal detrended moving average (MFDMA) method, was applied for the multifractal characterisation of Kepler time-series (de Freitas et al. 2017, 2019a,b). The MFDMA method filters out the local trends of non-stationary series by subtracting the local means. In addition, the MFDMA method can be used to investigate the local fluctuations in time-series using an a priori (fixed) timescale. According to Wang et al. (2014), fixing an a priori scaling range may lead to a crossover falling within the scaling range by mistake, and therefore the results could be biased. A method to avoid mistakes due to fixed scaling ranges is to measure the multifractal properties considering the entire timescales. To this end, a new approach named the multiscale MFDMA method (MFDMAτ) is introduced in this paper (see Sect. 6).

Our paper is structured as follows. Observations and technical details on the Kepler mission are presented in Sect. 2. In Sect. 3, we process and prepare the Kepler-30 time-series as obtained by the pre-search data conditioning (PDC) and simple aperture photometry (SAP) pipelines. In particular, we investigate the impact of the PDC pipeline over stellar variability when compared to the SAP data reduction, giving an overview of the characteristics of the PDC and SAP data using quality flags and discussing the difference between these pipelines. In Sects. 4, 5, and 6, we present our methods: spot modelling, the classical MFDMA method, and a new approach named the multiscale MFDMA method (MFDMAτ), respectively. In Sect. 7, we apply our set of methods to the Kepler-30 data and present the results of spot modelling and then introduce the results of the application of the multifractal analysis. In Sect. 8, we present our conclusions and final remarks.

2 Observations

From 2009 to 2013, the Kepler mission performed 17 observational runs, also referred to as quarters, of ~ 90 days each, composed of long-cadence (data sampling every ~30 min; Jenkins et al. 2010a) and short-cadence (sampling every ~ 1 min) observations (Van Cleve et al. 2010; Thompson et al. 2013). The mission public archive provides two types of data called simple aperture photometry (SAP) and pre-search data conditioning (PDC) time-series. The SAP flux is the sum of all calibrated background-subtracted flux values of the pixels belonging to the target aperture and for which a basic calibration has been performed, while the PDC flux is further corrected to remove instrumental and systematic effects by fitting a linear combination of the so-called co-trending basic vectors (CBV; Van Cleve & Caldwell 2009; Kinemuchi et al. 2012). These vectors describe the common trends observed in the 50% less variable and more mutually correlated targets falling within each CCD of the focal plane (Smith et al. 2012; Stumpe et al. 2014).

Considering our interest in the out-of-transit variations on timescales comparable to the stellar rotation, we make use of the long-cadence data of Kepler-30 publicly available at the MAST archive¹. These cover four years and the mean relative accuracy of each data point is ~ 260 ppm.

3 Kepler-30 preparation and pre-processing

The PDC pipeline was designed to maximise the chance of detecting planetary transits and not to study stellar intrinsic variability on timescales much longer than the transit duration. Therefore, the removal of the instrumental and systematic trends by fitting a linear combination of CBVs often results in overfitting and removal of the intrinsic target variability. The effect is particularly pronounced for the brighter and more intrinsically variable targets and leads to a remarkable attenuation of the components of the variability with timescales longer than 15–20 days (Stumpe et al. 2014; Gilliland et al. 2015). This is particularly relevant for our analysis. On the other hand, the SAP time-series, although presenting some residual instrumental effects and long-term systematic trends, is not affected by this problem. Therefore, we investigate the variability of Kepler-30 by analysing both the SAP and PDC time-series.

Because of the flux offsets between adjacent quarters and the long-term instrumental trends within each quarter (cf. Jenkins et al. 2010b), the time-series corresponding to different quarters were separately treated as follows:

Firstly, in each quarter, residual outliers were removed as described in Sect. 2 of De Medeiros et al. (2013). This procedure resulted in less than 0.2% of the data points being flagged as outliers and discarded. Furthermore, each Kepler time-stamp has a quality flag which alerts the user to systematic or instrumental effects that may affect the quality of the photometric measurement. We discarded all the datapoints with the quality flag SAP_QUALITY different from zero because they could potentially be affected by such systematic errors (Kinemuchi et al. 2012). In this way, the percentage of data points removed was less than 2%.

After applying this procedure, we checked its impact by calculating the cross-correlation coefficient for zero lag (i.e. when the time-series are not shifted) of the series before and after removing the outliers anddatapoints with a non-zero quality flag. To this end, we use the MATLAB function crosscorr². The values of the coefficients found were 0.92 for the PDC data and 0.94 for the SAP one. The gaps in our time-series are short enough to have a negligible impact on our analysis. Therefore, we decided to retain them and to avoid gap filling because artifacts could arise from the necessary interpolation processes.

Secondly, a third-order polynomial was fitted to the data to detrend the time-series separately for each quarter, as suggested by de Freitas et al. (2019b). The planetary transits had a very small impact on the polynomial fits given the small planetary radii and their long orbital periods. Our interest is in investigating the out-of-transit fluctuations on the timescale of stellar rotation or longer, and so the deformation of the transit profiles is not an issue, and they are removed in the next step. Figures 1 and 2 show the results of the two steps described above.

In the third step, the WOTAN package proposed by Hippke et al. (2019) was used to detrend the out-of-transit part of the signal in preparation for the transit removal in the fourth step (see below) as can be seen in the top panels of Figs. 3 and 4 (red lines). To this end, we use the iterative biweight method with a window of 0.75 days as proposed by Hippke et al. (2019). The detrended flux for both pipelines can be seen in the middle panels of Figs. 3 and 4.

In the fourth step, planetary transits were removed according to the parameters given in Table 1 of Sanchis-Ojeda et al. (2012). We also performed a visual inspection and verified that the obtained out-of-transit series has no evident transit signal left after the application of our procedure. The gaps introduced by discarding the transits do not exceed 9 h, which is the duration of the transit of the most distant planet Kepler-30d. These gaps have a negligible effect on our spot modelling because the rotation period of the star is ~ 16 days. In the same way as mentioned above, we did not perform any interpolation process to fill in the gaps because the amount of data removed is negligible. Examples of the final results of our procedure can be seen in the bottom panels of Figs. 3 and 4.

Lastly, the final time-series is obtained by combining the 17 quarters and consists of four years of data. From each quarter series, the median value was subtracted and then all the residual quarter series were stitched together. The final time-series are shown in the upper panels of Figs. 5 and 6, where the median of the errors of the single photometric measurements is 2 × 10⁻⁴ in relative flux units. The number of data points is N = 60007 and N = 60 018 for the PDCand SAP time-series, respectively.

The final PDC and SAP time-series and their Lomb–Scargle periodograms (cf. Lomb 1976; Scargle 1982) with the significance levels are illustrated in Fig. 5. For both time-series, there are two main peaks at 7.8 and 16.8 days. Nevertheless, for SAP time-series, a lower peak at 46.6 days also appears. As mentioned above, the SAP pipeline preserves other variabilities that are attenuated or removed by the PDC pipeline. Consequently, the residual time-series (RTS; i.e., the difference) between theSAP and PDC time-series can reveal the origin of that peak (see upper panel of Fig. 6).

In Fig. 6 (bottom panel), we plot the Lomb–Scargle periodogram of the difference time-series (RTS). Only a few points were lost during the process of subtracting the time-series. The Lomb–Scargle periodogram of Fig. 6 clearly shows that the period of 46.6 days, as well as the periods of 23.1 and 33.9 days, are preserved by the SAP pipeline. The Lomb–Scargle technique has a clear disadvantage, because it is limited to analysing inthe frequency space and therefore is not capable of resolving the time-evolution of a given periodicity. In addition, this technique does not allow us to analyse the behaviour of the variability, or of the light fluctuations in particular, over different timescales and amplitudes. To this end, we now consider the spot modelling (cf. Sect. 4) and the MFDMA approaches (cf. Sects. 5 and 6).

Fig. 1 Kepler-30 PDC time-series. Top panel: PDC time-series (black dots) and polynomial third-order fits. Bottom panel: final PDC time-series adjusted by polynomial functions. The signatures of three planets are maintained.

Fig. 2 As in Fig. 1 but for SAP time-series.

Fig. 3 Detrending of Kepler-30 with the iterative biweight filter using a window size of 0.75 days (red line), which allows detection of the transits. Normalised (top), detrended (middle), and residual flux (bottom), shown for only a segment of the total time-series. Black dots indicate the PDC pipeline of Kepler-30 after the first three steps of our procedure. Detrended flux is a combination of noise background and planet transits, indicating the trending of the time-series potentially due to rotational modulation of the star.

Fig. 4 Same as Fig. 3 for SAP pipeline.

Fig. 5 Final PDC (first panel) and SAP (third panel) time-series adjusted according to the steps of Sect. 3. The corresponding Lomb–Scargle periodograms are shown at the bottom. The dashed-dotted blue line marks the 99% confidence level. For SAP data, there is a significant peak to 46.6 days also found in Fig. 6. As can be seen, the PDC time-series has a flat periodogram for periods longer than 20 days confirming that the PDC pipeline tends to remove all the long-term trends.

Fig. 6 Difference between the SAP and PDC time-series, where a peak of 46.6 days can be seen. Second and third peaks of 23.1 and 33.9 days are also seen.

4 Spot modelling of Kepler-30 time-series

The small v sin i ~ 1.94 ± 0.25 km s⁻¹ of Kepler-30(Fabrycky et al. 2012) prevents the application of Doppler imaging techniques to this star (Strassmeier 2009), and therefore the only possibility of mapping its surface is a spot modelling based on the inversion of its light curve.

Spot modelling is a notoriously ill-posed problem because it tries to extract a bi-dimensional map from a one-dimensional time-series. To account for the inherent degeneracies of the solution, we apply a regularisation to our inversion consisting in introducing a priori information to select the unique spot map that maximises the configurational entropy of the spot pattern for a given levelof the chi square of the fit. We choose the configurational entropy to be maximal when the star is unspotted, that is, when the spot map has the least spotted area.

The fundamental assumption of our spot model is that the spot pattern remains stable for a time interval Δt_f during which the flux received from the star is modulated only by stellar rotation. A moderate random starspot evolution on a timescale significantly shorter thanΔt_f is not a problem because it contributes to the noise level and does not affect the spot modelling thanks to the application of the maximum entropy (ME) regularisation (see Appendix C). However, sizeable changes in the spot pattern on a timescale comparable with Δt_f can produce systematic errors in the spot map. An optimal selection of Δt_f is therefore relevant for our modelling and is discussed in Appendix C together with the determination of the other parameters of our model.

A major advantage of modelling a star with transiting planets such as Kepler-30 is the resulting possibility to measure the obliquity of its spin axis projected on the plane of the sky, which allows us to constrain the inclination i of the stellar spin to the line of sight (Sanchis-Ojeda et al. 2012). Given the low projected obliquity, we assume the stellar spin to be normal to the plane of the orbit of the more massive close-by transiting planet Kepler-30c, that is, virtually perpendicular to the line of sight, thus fixing the geometry of the rotating star to be mapped (see Appendix C for the adopted value of the inclination). We note that in the case of a star with an unknown inclination, spot modelling is much more uncertain because of the strong degeneracy between the inclination and the latitudes of the spots.

As we observe Kepler-30 virtually equator-on, spot modelling is unable to provide information on the latitude of the starspots because the duration of their visibility along a stellar rotation is independent of their latitude (we neglect the effect of latitudinal differential rotation, the value of which is unknown a priori). Therefore, we can only map the distribution of the spotted area versus the longitude and the time along successive intervals of duration Δt_f. As a consequence, we consider two-dimensional (2D) spot maps only as an intermediate step in our modelling and derive the distribution of the spotted area versus longitude from them, that is, we compute 1D maps of the filling factor by integrating over the latitude in the 2D maps. Such longitude distributions of the filling factor can be directly related to the light modulation under the assumption that the spot contrast is constant. This greatly reduces the degeneracies associated with the non-uniqueness of light-curve inversion.

We include the effect of solar-like faculae in our model, but their contribution to the flux modulation is found to be very small and can be neglected for a star as active as Kepler-30 where dark spots dominate (cf. Radick et al. 2018, see also Appendix C). More details of our spot modelling and the impact of the modelling parameters on the spot maps are presented in Appendix C.

Recent works have questioned the validity of spot models in their assumption of a few discrete spots on the stellar surface (e.g. Basri & Shah 2020). Our model is based on a continuous distributions of spots and solar-like faculae on the stellar photosphere and has been extensively tested using comparisons to the Sun (Lanza et al. 2007) and to active stars with transiting planets by comparing the spot maps from light-curve inversions with the real distributions of the sunspot groups and the spots detected by occultation during planetary transits, respectively (Silva-Valio & Lanza 2011; Lanza et al. 2019). Thanks to such comparisons, we are confident that our spot modelling is capable of reconstructing the distribution of the spot area versus longitude with a resolution of ≈ 60° in the case of Kepler-30. In consideration of the systematic differences introduced by the Kepler pipeline between the SAP and the PDC time-series, we apply our modelling to both of them and compare the results to look for artifacts produced by those systematic differences.

5 Multifractal description of time-series

5.1 Basic concepts and definitions

The light variations in the time-series of Kepler-30 display a wide range of timescales ranging from tens of minutes, characteristic of photospheric convection, to the rotational modulation due to photospheric brightness inhomogeneities (starspots), up to months or years characteristic of stellar activity cycles. The rotational modulation signal is quasi-periodic and non-stationary owing to theevolution of starspots on the surface of the star. To characterise the fluctuations on the shortest timescales, which are mostly produced by stellar convection, we can use methods of multifractal time-series analysis; we briefly introduce these below. We mainly follow the treatment by Kantelhardt (2015) to whom we refer the reader for a more detailed account because here we limit ourselves to a minimal introduction with the purpose of making this paper accessible to readers with no previous knowledge of the field.

To characterise the wide frequency spectrum of convective fluctuations, it is useful to look for self-similarity properties in the fluctuations themselves. Specifically, we define Δf(t) as the light fluctuation at time t, measured with respect to some constant mean value. A rescaling of the time by an arbitrary factor a, i.e., t → at, may require a re-scaling of the fluctuations by a factor a^H, that is, Δf(t) → a^HΔf(at), in order for the re-scaled fluctuations to follow the same statistical distribution as the initial ones. A time-series that obeys this scaling property for an arbitrary factor a is called a self-affine or self-similar time-series and the exponent H is called the Hurst exponent of the time-series. Several stochastic fluctuations, including those associated with turbulent convection, satisfy this definition.

The Hurst exponent H characterises the kind of self-similarity. The simplest example is the time-series of the mean-square displacement x² (t) of a point performing a random walk with a normally distributed step size. In this case, the time scaling t → at requires a corresponding scaling x → a^1∕2x in order to preserve the statistical distribution of the variable x. In other words, the time-series x(t) exhibits a self-similar appearance with H = 0.5, that is, by expanding the time axis by a factor a and the space axis by a factor a^H, the distribution of the statistical fluctuations of x(t) is the same as the original one. Such a scale self-similarity (or scale invariance) is one of the fundamental properties that define fractals.

Self-affine time-series are persistent in the sense that a large fluctuation is likely followed by another large fluctuation and a small fluctuation by a small one. The persistence holds on all the timescales for which the self-affine scaling holds. In more complex systems, the persistence may change according to the timescale being stronger on some timescales and weaker on others. For example, the weather time-series are usually persistent on timescales of hours or days, but the degree of persistence is higher on shorter timescales. Moreover, the characteristic persistence timescale generally changes with the season, usually being longer during winter and summer and shorter in autumn and spring.

Considering a stochastic self-affine time-series x(t) with the variable sampled at equidistant times t_i, with i = 1, 2, …, N, we define the time-series of the increments Δx_i ≡ x(t_i) − x(t_i−1). In the case of the simple random walk model considered above with H = 0.5, the Δx_i are independent of each other, but when H > 0.5, a positive Δx_i is likely to be followed by a positive increment (persistence), while when H < 0.5, a positive Δx_i is likely to be followed by a negative one (anti-persistence). The degree of persistence can be quantified by the auto-correlation function in the case of a stationary time-series, that is, one with constant mean and standard deviation. If we denote the timescale as n, also indicating a number of datapoints in a uniformly sampled time-series, the autocorrelation is defined as: $$C(n)=\frac{1}{{\sigma }^2}\frac{1}{N-n}{\sum }_{i=1}^{N-n}\Delta x_i\Delta x_{i+n},$$(1)

where σ² is the variance of the time-series of the increments Δx_i.

When the Δx_i are uncorrelated, C(n) = 0 for n > 0. A short-range correlation can be described by an exponentially decaying C(n) ∝ exp(−n∕n_d), where n_d is the decay timescale. For example, this is the case of an autoregressive process defined by the linear recursion Δx_i = cΔx_i+1 + ϵ_i, where c = exp(−1∕n_d) and the ϵ_i are normally distributed random deviates.

In the case of long-term correlations, the asymptotic behaviour of C(n), that is, its scaling for large n, can be described by a power law: $$C(n)\propto n^{-\gamma },$$(2)

where the exponent 0 < γ < 1 characterises the long-term correlations. It can be shown that a long-term correlated, i.e. persistent, behaviour of Δx_i leads to a self-affine scaling of the x(t_i) characterised by a Hurst exponent H = 1 − γ∕2. Therefore, the Hurst exponent can be used to characterise the autocorrelation function of the increments of the time-series.

Ideally, the estimation of the exponent γ can be made by computing the autocorrelation function of the increments or the power spectrum of the original time series x(t_i) (Kantelhardt 2015). In practice, these approaches are doomed to fail because of the noise that dominates the autocorrelation and the power spectrum in the asymptotic regime of large n values to be considered when deriving γ. Therefore, more robust statistical approaches have been developed for such a determination. These are based on the so-called fluctuation functions F_q(n) that measure the scaling of the fluctuation amplitude with the timescale n, where q is a parameter defining the moment of the statistical distribution of the x(t_i) measured bythe function itself (see Appendix A for examples of fluctuation functions). The fluctuation function based on the standard deviation corresponds to q = 2 and plays a special role in the analysis.

The scaling of the fluctuation function with n is characterised by a generalised Hurst exponent or Holder exponent h(q) that depends on the moment of the distribution sampled by F_q(n). To define the fluctuation function F_q(n), we start from the root-mean-square (RMS) fluctuation function F_ν(n), which is defined by $$F_{\nu }(n)={\left\{\frac{1}{n}{\sum }_{i=1}^n{ϵ}_{\nu }^2(i)\right\}}^{\frac{1}{2}},$$(3)

where ϵ_ν(i) is the deviation of the ith datapoint from a suitable moving average introduced to remove the long-term variation in the time-series (see Appendix A for its precise definition). Starting from the RMS fluctuation function, we define the fluctuation function of order q as $$F_q(n)={\left\{\frac{1}{N_n}{\sum }_{\nu =1}^{N_n}{F}_{\nu }^q(n)\right\}}^{1/q}\text{for}q\ne 0$$(4)

and, for q = 0, $$\ln \left[F_0(n)\right]=\frac{1}{N_n}{\sum }_{\nu =1}^{N_n}\ln [F_{\nu }(n)],$$(5)

where N_n is the number of non-overlapping segments fora given segment of size n in the time-series, i.e., N_n ~ N∕n, where N is the total number of datapoints in the series (see Appendix A for the precise definition of N_n). For larger values of n, the fluctuation function scales following a power-law given by $$F_q(n)~n^{h(q)},$$(6)

which allows us to define the generalised Hurst exponent h(q).

The previously defined Hurst exponent corresponds to the value attained for q = 2, that is, H ≡ h(2). By definition, when h(q) is constant, we speak of a fractal (or monofractal) time series, while when h depends on q we have multifractal time series. In the former case, the scaling exponent of the fluctuation self-affine law is independent of the amplitude of the fluctuations, while in the latter case it depends on their amplitude because different values of q correspond to different moments that sample fluctuations of large or small amplitudes in different ways (cf. Eq. (A.5)).

The scaling of the fluctuations with their amplitude can be quantitatively characterised by means of the multifractal spectrum and its characteristic parameters as detailed in Appendix A. Multifractal time-series can be produced, for example, by simple, non-linear, recurrence laws, such as the so-called logistic law x_i+1 = λ_llx_i(1 − x_i), where the coefficient λ_ll falls in a range leading to chaotic behaviour (e.g. Luo & Han 1992).

5.2 Origins of multifractality and the impact of quasi-periodic modulations in the time-series

The multifractal character of a stochastic time-series is manifested by different scaling exponents for small and large fluctuations on the same timescale n. This behaviour can be originated by the simultaneous presence of different long-term correlations and persistences – that depend on the amplitude of the fluctuations – and/or by their statistical distribution that, even in the absence of long-term correlations, deviates strongly from a Gaussian (normal) distribution. In the latter case, the distribution of the fluctuations usually shows fat tails, as in the case of a Cauchy distribution, for example.

It is possible to discriminate between these two sources of multifractality by considering the so-called surrogates of the original time-series. A first kind of surrogate is simply produced by randomly shuffling the datapoints over the set of the N indexes i. In this way, we keep the statistical distribution of the datapoints, but remove all short- and long-term correlations, producing a new time-series with h(q) close to 0.5, that is the Hurst exponent of a random time-series consisting of uncorrelated points. A second kind of surrogate can be obtained by taking the Fourier transform of the original time-series and constructing a modified Fourier transform with the same modulus, but a phase randomly extracted in the [−π, π] interval. The backward transformation of this modified transform into the time domain gives a time-series whose points retain most of the correlations present in the original time-series, but whose distribution is nearly Gaussian, thus effectively eliminating the effect of non-Gaussianity in the original distribution (e.g. Kantelhardt 2015; Wang et al. 2014).

By comparing the h(q) function and the multifractal spectrum of the original time-series with those of the above surrogates, we can extract information on the dominant source of multifractality. Specifically, when the dominant source of multifractality is time correlation, the h(q) function of the randomly shuffled series will be remarkably different from that of the original series, while the h(q) of the phase-randomised series will be similar. On the other hand, when the dominant source of multifractality is a non-Gaussian statistical distribution of the datapoints, the converse will be true.

From a more general point of view, time-series can result from both deterministic and stochastic processes, the former usually being responsible for long-term trends or oscillations. Those trends should be removed to allow the best determination of the properties of the stochastic fluctuations. This is usually done by fitting the trends with polynomials of different degrees in the calculation of the fluctuation functions (e.g. Kantelhardt 2015) or by subtracting a moving average, as in our case (see Appendix A). However, in the case of astrophysical systems, and of magnetically active stars in particular, a clear separation between deterministic long-term trends and short-term stochastic fluctuations is not always possible. For example, the rotational modulation of the flux due to photospheric starspots contains a stochastic component produced by the evolution of the individual starspots often appearing randomly in time and/or longitude (e.g. Lanza et al. 2019, and references therein).

In light of such an impossibility of performing a clear separation between the different components of the photometric time-series of magnetically active stars, the adopted approach has been that of including the rotational modulation and the evolution of the surface brightness inhomogeneities as additional sources of multifractality. In other words, the multifractal analysis itself has been used as a technique to characterise those processes in addition to standard techniques such as the Lomb–Scargle periodogram or autocorrelation to measure stellar rotation periods (e.g. McQuillan et al. 2013) or spot modelling to derive the evolution of the surface features (Lanza 2016). This approach has been extended even to the detection and characterisation of strictly periodic phenomena such as planetary transits in the presence of noise (e.g. Agarwal et al. 2017).

In a sample of magnetically active solar-like stars observed by the CoRoT space telescope, de Freitas et al. (2013a) found a correlation between the Hurst index and their rotation period by assuming that their time-series were monofractal. In a subsequent work on the same sample, de Freitas et al. (2016) established that their time-series are indeed multifractal, that is, the correlations produced by magnetic activity and rotation are better described by assuming a generalised Hurst exponent h(q) and adopting the full set of parameters introduced in Appendix A to characterise the multifractal spectrum, in particular Δα.

Later, de Freitas et al. (2017) analysed a sample of 34 M dwarf stars previously investigated by Mathur et al. (2014) to test the behaviour of the Hurst exponent against the magnetic activity indicator S_ph. The latter index is based on the standard deviation of the flux in the Kepler passband computed over non-overlapping time intervals of k mean rotation periods P_rot after removing the contribution of the photon shot noise to the standard deviation. An average index ${\langle S_{\text{ph}}\rangle }_k$ is then defined as the arithmetic mean of all the individual S_ph indexes determined from each of the kP_rot intervals along the time-series. Therefore, it includes the contributions of all the sources of light variations from the shortest timescales up to kP_rot, including the rotational modulation. On the other hand, the index H = h(2) is defined as the scaling exponent of the RMS fluctuation function F₂(n) ∝ F_ν(n) ~ n^H which is evaluated with the fluctuations obtained after subtracting a moving average from the time-series. This gives a greater weight to the shorter timescales, while reducing the contribution of the rotational modulation and longer timescales. Therefore, with a generally adopted value of k = 5 (Mathur et al. 2014), the correlation between the H index and the ${\langle S_{\text{ph}}\rangle }_k$ index is expected to be low because of such a difference in the respectively dominant timescales. This expectation was confirmed by de Freitas et al. (2017) who found a slight anti-correlation between the two indexes with a Pearson coefficient of − 0.33 when considering the above-mentioned sample of M dwarf stars (see their Fig. 6).

On the other hand, even within a limited range in the rotation period (< 15 days), de Freitas et al. (2017) confirmed the strong correlation between H and P_rot and found the Hurst exponent to be an indicator of magnetic activity. Recently, de Freitas et al. (2019a,b) extended the field of application to time-series showing more than one significant rotational periodicity, which can be considered as an indication of stellar differential rotation. Using a large sample of Kepler active stars, de Freitas et al. (2019a) showed that in the relationship between the Hurst exponent and the relative amplitude of the differential rotation ΔP∕P a strong trend of increasing ΔP∕P toward larger H-index is observed. In addition, the study shows that the correlation is stronger for the most active stars in the sample.

6 MFDMA analysis

Recently, de Freitas et al. (2019a) used a multifractal analysis to investigate the multi-scale behaviour of a set of ~ 8000 active stars that were observed by the Kepler mission. de Freitas et al. (2016, 2017, 2019a,b) showed that the multifractal detrending moving average (MFDMA) algorithm, which was developed by Gu & Zhou (2010)³ and Tang et al. (2015), is a powerful technique that provides valuable information on the fluctuations of a time-series. A complete description of the MFDMA algorithm can be found in Appendix A.

Several authors (e.g. Gieratowski et al. 2012; Wang et al. 2014) have suggested that a single scaling exponent (H = h(2)) is inadequate to describe the internal dynamics of complex signals, and consequently, classical (multi)fractal methods (e.g. DFA and MFDFA) have been modified for an estimation of the temporal dependence of the spectrum of scale exponents. These improved methods areused to better quantify the short- and long-range correlations, which is important in our case because multifractality in our time-series arises mostly because of the correlations (cf. Sect. 4.2 and Table A.1).

Inspired by these studies, we decided to introduce another method to calculate a time-dependent scaling exponent h(q, τ), where τ is a timescale (see below). We call it the multiscale MFDMA method (MFDMAτ) and it allows us to extend the investigation of stellar variability by including a time dependence in the description of the multiscale fluctuations. MFDMAτ is relatively immune to additive noise and non-stationarity, including the non-stationarity due to the occurrence of events of a different dynamics.

To explore the multiscale fluctuations in flux, we introduce the increments Δx(t, τ) = x(t + τ) − x(t), where τ varies from 29.4 min (Kepler cadence) to ~60 days. We apply the MFDMA analysis to the increment time-series Δx(t, τ), thus obtaining a Hurst exponent h(q, τ) that is also a function of the timescale τ, and in this way introduce the MFDMAτ method. The consideration of the increment time-series in the MFDMAτ approach strongly reduces the effects of all the stationary modulations with a period equal to τ or one of its multiples. In other words, this method allows us to isolate the effects of the fluctuations after removing those periodicities. Given the remarkable rotational modulation present in our time-series, this method will allow us to study the properties of the fluctuations after removing the modulation itself with an appropriate selection of the timescale τ.

We explored different values of τ using a plot of x(t + τ) versus x(t) following an approach reminiscent of that used to study multifractal attractors by means of plots of the trajectory of the dynamical system in a bi-dimensional phase space. In Fig. 7, we plot the case corresponding to the optimal value τ = 8.7 days, which shows a repetitive, quasi-stationary pattern organised around a fixed point with a behaviour similar to that of an attractor.

According to Henry et al. (2001), the optimal delay time τ can be determined using the approach of topological properties based on the behaviour of the autocorrelation function. There are various proposals for choosing an optimal delay time using the autocorrelation function. We choose the best value for τ as that corresponding to the maximum of the mean value of the autocorrelation coefficient. In the subplot of Fig. 7, we show that the mean autocorrelation coefficient of the stochastic process x(t), defined as E[x(t + τ)x(t)] (for all t), is maximum for τ = 8.7 days.

Fig. 7 Time-delayed PDC time-series x(t + τ) with the optimal delay τ = 8.7 days vs. the original time series x(t). The subplot shows that τ = 8.7 days corresponds to the maximum value of the autocorrelation function of the time-series.

7 Results and discussions

We present the results of our spot modelling first and then introduce the results of the application of the multifractal analysis to the original photometric time-series of Kepler-30 as well as to the residuals of the spot models. In this way, we can apply the results of the spot modelling to elucidate those of the multifractal analysis, discussing the characteristic rotational and activity timescales found in the data.

7.1 Best fit of the Kepler-30 light curves

In Fig. 8, we plot the PDC light curve together with the unregularised best fit obtained with our model and the corresponding residuals. The unregularised best fit provides the minimum chi square and effectively reproduces the spot light modulation giving residuals where it has been filtered out.

The distribution of the residuals is shown in Fig. 9 and is almost Gaussian with a mean μ_res very close to zero (μ_res = −2.18 × 10⁻⁷) and a standard deviation σ₀ = 3.78 × 10⁻⁴ in relative flux units, slightly larger than the photon shot noise. This indicates the presence of other processes that contribute to the random noise level, probably associated with surface convection and magnetic fields affecting the stellar flux.

A GLS periodogram of the residuals is shown in Fig. 10. Our spot model strongly reduces any variability with periods longer than about 6 days leaving the variability below 3–4 days almost intact. The maximum of the GLS is reached at a period of 2.565 days. The peak is about two times higher than the noise level as indicated by the other nearby peaks, and therefore is probably not very significant – usually a peak reaches a significance of 99 percent when its height is at least four times the noise level –, but we shall discuss its possible origin in Sect. 7.3.1.

The unregularised best fit of the SAP light curve is very similar to that of the PDC light curve except for a few larger residuals probably associated with some small uncorrected instrumental effects. By fitting the residual distribution with a Gaussian, we find a slightly larger mean μ_reg = −1.737 × 10⁻⁶, but a standard deviation of σ₀ = 3.746 × 10⁻⁴ in relative flux units which is virtually identical to that of the PDC residual distribution. Also the periodogram of the residuals is very similar and is not shown here. The residuals of the PDC and SAP light-curve fits will be used to study the multifractal properties of the stellar flux variability without the additional complication introduced by the light modulation produced by starspots. Analysing the residuals of the unregularised best fits allows us to avoid possible systematic effects introduced by the a priori assumptions adopted for the ME regularisation that increase the standard deviation and make the mean of the residuals systematically negative (see Appendix C).

Fig. 8 Upper panels: PDC light curve of Kepler-30 (solid black dots) and the unregularised best fit obtained with our spot model (red solid line). Lower panels: residuals of the best fit.

Fig. 9 Distribution of the residuals of the unregularised best fit to the PDC light curve in Fig. 8.

Fig. 10 GLS periodogram of the residuals of the unregularised best fit to the PDC light curve in Fig. 8.

7.2 Spot maps and variation of the spotted area

The spot map obtained from the ME-regularised best fit of the PDC light curve is plotted in Fig. 11, where we plot the spot filling factor versus the longitude and the time. The adopted longitude reference frame rotates with the mean rotation period of the star of 16.0 days and the longitude increases in the direction of stellar rotation. Therefore, spots rotating with the mean rotation period of 16 days stay at a constant longitude on this plot. On the other hand, spots rotating faster than the mean rotation period show a longitude that increases with time, while those rotating slower show a longitude that decreases with time. We see that the evolution of the starspots of Kepler-30 is rather fast with sizeable variations of the filling factor over timescales of 20–30 days (the time resolution of our spot modelling is Δt_f = 11.963 days). We introduce a modified Julian date as MJD = HJD − 2454900. During the first half of the modelled time series, that is, up to MJD ~ 800, spots cluster around two main active longitudes at ~0° and ~200°, while a smaller and intermittent active longitude is present in between them at ~ 100°. Starting from MJD ~300, the main active longitudes migrate towards higher and lower longitudes, respectively, probably as a consequence of the decay of the spots rotating with the mean rotation period and the formation of new spots that rotate faster and slower, respectively. The intermediate active longitude is approached by the main longitude initially at ~ 200° and looses its identity. ForMJD ≳800, the active longitudes are less well defined and the pattern is characterised by individual spots with lifetimes ranging from ≈ 30 to ≈ 200 days that rotate at different rates. The spot map obtained with the regularised best fit to the SAP light curve is similar to that displayed in Fig. 11 and is plotted in Appendix C.

The migration of the different spots in longitude can be used to estimate the surface differential rotation if we assume that the migration is caused by the location of those spots at different latitudes on a differentially rotating star. Given the lack of information on the spot latitudes, we can only derive a lower limit to the surface shear. The estimate should be regarded with great caution because the apparent migration can be the result of the starspot evolution rather than of the differential rotation. In the case of Kepler-30, the evolution of individual starspots occurs on timescales shorter than or comparable to those of the migration of the active longitudes, and therefore this source of systematic error is certainly present and casts a further uncertainty on our estimate.

Looking at the spot map in Fig. 11, the best suitable feature with which to estimate the differential rotation is the active longitude, which starts its migration at MJD ~ 400 at a longitude of ~ 30° and reaches a longitude of ≈100° at MJD~550. Taking into account an uncertainty of at least ±40° in the longitude shift between those two dates owing to the intrinsic width and spreading of the active longitude in time, we estimate a relative surface shear of ΔΩ∕Ω ~ 0.020 ± 0.012, where Ω is the angular velocity of rotation, the mean value of which corresponds to the mean rotation period of 16 days. A similar estimate based on the spot map obtained from the SAP time-series (cf. Fig. C.1) gives a comparable result. In any case, we stress that this estimate is particularly uncertain, especially because of the fast spot evolution seen inside the active region with individual spot lifetimes of tens of days, which is comparable with the time resolution of our spot modelling. Other active longitudes show a less clearly defined migration pattern and are more affected by the intrinsic evolution of the spots, and therefore they are not considered for our estimate.

The total spotted area, obtained by summing up the filling factor over all the longitudes, is plotted in Fig. 12 against time. Only the values obtained from time intervals Δt_f without significant data gaps have been plotted in Fig. 12, in order to avoid systematic errors associated with the ME regularisation that tends to decrease the spotted area in the time intervals containing gaps. To find if the distribution of the datapoints within each interval Δt_f is uniform, we subdivide it into five equal subintervals and count the number q_i of datapoints in each subinterval (i = 1,.., 5). If the ratio (max{q_i}− min{q_i})∕max{q_i}≤ 0.2, the corresponding interval Δt_f is considered as not affected by gaps and its area value is considered in our analysis. Despite the fact that we disregard the area values obtained from intervals with significant gaps, we still see some points with deviations up to ± 10 percent from the mean value in Fig. 12. Such relatively large deviations are due to flux discontinuities at the epochs when successive intervals are stitched together after a repointing of the Kepler telescope or to some small deviations from the convergence criterion adopted to fix the level of regularisation (see Appendix C).

A GLS periodogram of the spotted area shows a peak at 33.8 days with a false-alarm probability (FAP) of 0.29 according to the analytical formula by Zechmeister & Kürster (2009). The corresponding best-fitting sinusoid is plotted in Fig. 12 with thearea values affected by the systematic effects mentioned above clearly deviating from it. Repeating the same analysis with the spot model of the SAP light curve (see Appendix C), which is not affected by the systematic effects introduced by the PDC pipeline, we find the same periodicity with a period of 33.9 days and a FAP = 0.024, confirming the presence of a short-term modulation in the total spotted area of Kepler-30, reminiscent of the short-term cycles found in CoRoT-2 (Lanza et al. 2009; Zaqarashvili et al. 2021).

Such a short-term cycle of ~34 days is relatively close to the synodic period P_syn = 35.2 days of the planet Kepler-30b as computed with a rotation period P_rot = 16.0 days and an orbital period P_orb = 29.334 days from $P_{\text{syn}}^{-1}=P_{\text{rot}}^{-1}-P_{\text{orb}}^{-1}$. This period also appears in the periodogram of Fig. 6 and is close to the period of the modulation of the spotted area. In addition, the 23.1-day period (see Fig. 6) is close to the synodic period of Kepler-30c. However, an interpretation in terms of tidal or magnetic star-planet interactions (cf. Lanza 2012, 2013) is problematic because of the large star–planet separations that range from ~ 42 to ~ 121 stellar radii for the three planets⁴. More precisely, both the tidal torque on the star and the energy density of its coronal magnetic field decrease with distance as ${(a/R_{\text{s}})}^{-6}$, where a is the orbital separation and R_s the radius of the star.

Fig. 11 Distribution of the spot filling factor vs. the longitude and time for the ME regularised spot model of the PDC light curve of Kepler-30. The minimum of the filling factor corresponds to dark blue regions, while the maximum is rendered in orange (see the colour scale close to the lower right corner of the plot). We note that the longitude scale of the horizontal axis is extended beyond the [0°, 360°] interval to better follow the migration of the starspots.

Fig. 12 Total spotted area on Kepler-30 as derived from the spot modelling of the PDC light curve versus the time (green dots). A sinusoid with a period of 33.784 days corresponding to the maximum of the GLS periodogram is overplotted (red solid line).

Fig. 13 Log-log plot of the fluctuation functions Fq(n) (black circle for q = 2) calculated for the final PDC time series presented in Fig. 5. The red curves correspond to q between −5 and 5 in steps of 0.2. Vertical dashed lines mark three domains of the fitting windows for the small n between 29.4min and 8.7 days, for middle n between 8.7 and 16.8 days and, for the large n greater than 16.8 days. The green dashed line gives the average slope H = h(2) = 0.95 for timescales shorter than 8.7 days.

7.3 Multiscale multifractal analysis

7.3.1 Fluctuation functions

In Figs. 13 and 14, we plot the fluctuation functions F_q (n) against the timescale n for different values of q for the PDC and SAP time-series, respectively. In this logarithmic plot, the scaling relation given by Eq. (6) becomes a straight line in the intervals of n where h(q) is constant. A value of n separating two consecutive intervals where h(q) is constant for a given q is called a crossover. We see a crossover in the plots for the values of q > 0 at 8.7 days as indicated bythe vertical dashed line. The crossover corresponds to the optimal delay time as determined with the procedure described in Sect. 6 (cf. Fig. 7). The change in slope is made clear by the green dashed line that gives an average slope H = h(2) = 0.95 for timescales shorter than8.7 days. The small decrease in the fluctuation functions following the crossover and extending up to the rotation period of 16.8 days is barely significant and is due to the variation in amplitude produced by the use of a moving average to approximate what is actually a nearly sinusoidal modulation in the computation of the fluctuation function itself (see Appendix A). Once the rotation period is reached, the amplitude of the fluctuation function saturates at a nearly constant value, except for some small oscillations produced by the sinusoidal modulation of the signal (cf. Appendix B). In other words, the slope h beyond the crossover is almost zero because the variability due to the rotational modulation dominates over the stochastic fluctuations. We reiterate that the fluctuation function with q = 2 is based on the variance of the fluctuations in the time-series.

For both the PDC and SAP time-series, we find a crossover at the timescale corresponding to the first harmonic of the rotation period, that is ~8.7 days. The predominance of the modulation at the first harmonic is seen in the light curves of active stars when two active longitudes on opposite hemispheres are responsible for most of the flux modulation, a behaviour that is confirmed by our spot modelling (see Sect. 7.2). From the perspective of stochastic process analysis, the behaviour we see is reminiscent of an attractor (see Fig. 7) with a phase trajectory that revolves around a period-one fixed point. As a result, the logarithms of the fluctuation functions F_q (n) show a linear regime until ~8.7 days after which they become saturated and begin weakly oscillating. There is a slight difference in the level of oscillation when the time-series PDC and SAP are compared as can be seen in Figs. 13 and 14. Based on the results in Sect. 7.1, the presence of a few larger residuals found in the SAP time-series probably explains theplateau within the range from 8.7 to 16.8 days, because the noise is responsible for dampening the oscillation.

In Fig. 15, we plot the fluctuation functions for the difference time series RTS. In this case, we see a crossover at a longer timescale of ~23.5 days. We can associate it with the characteristic evolutionary timescales of active regions in Kepler-30 because the rotational modulation was remarkably reduced by subtracting the two time-series from each other (cf. Sect. 7.2). We remind that the intrinsic variability on timescales longer than 15–20 days was preserved in the SAP time-series, while it was strongly reduced in the PDC time-series owing to the de-trending applied by the Kepler pipeline (cf. Sect. 3). Therefore, the long-term intrinsic variability stands out in the difference time-series.

The fluctuation functions of the time-series of the residuals of the unregularised spot modelling of the PDC light curve are plotted in Fig. 16. We see two crossovers at ~ 1 and ~ 35 days, where the Hurst exponent H = h(2) changes from 0.73 to 0.10 and then to 0.37, respectively. Therefore, the fluctuations with timescales shorter than ~1 day are characterised by persistence, while those on longer timescales show a value of H characteristic of antipersistent time-series. The timescale of ~1 day is characterised by a remarkable increase in the power level in the GLS power spectrum of the residuals (cf. Fig. 10) and the multifractal analysis reveals a change in the persistence property of the fluctuations on only that timescale. This corresponds to the characteristic turnover time of supergranular convective cells in the Sun and Sun-like stars (cf. Meunier et al. 2015). The other crossover at ~35 days coincides with the possible short-term cycle in the area of the starspots (cf. Sect. 7.2) and the change in the slope of the fluctuation functions may indicate changes in the properties of convection on that timescale of stellar activity. The fluctuations functions of the spot model residuals of the SAP light curve are very similar and shall not be discussed here. In both the cases, we do not find any particular feature of the fluctuation functions at the peak period of 2.565 days in the GLS periodogram of the residuals, suggesting that it is not significant and not associated with a change in the behaviour of the stellar light fluctuations.

We based our conclusions on the fluctuation functions F_q (n) with positive q where the effect of large-scale fluctuations is amplified with respect to that of small-scale fluctuations in the computation of the moments of order q > 0. Conversely, small-scale fluctuations will have the most important impact on the fluctuation functions F_q (n) when q < 0 because the effect of large-scale fluctuations is attenuated by a negative exponent in the calculation of the momenta (cf. Eq. (A.5)). In all the above plots of the fluctuation functions F_q (n) with q < 0, we see that the exponent of the scaling law, that is, the slope of the plot, changes frequently and abruptly. This indicates that the asymptotic regime where F_q(n) ∝ n^h(q) is not fully reached in those plots, which is likely an effect of a few very small fluctuations that dominate the functions.

Fig. 14 As in Fig. 13 but for final SAP time-series.

Fig. 15 As in Fig. 13 but for final RTS time-series. In this case, only two domains are considered and separated by the vertical dashed line at 23.5 days.

Fig. 16 Fluctuation functions of the residuals of the unregularised spot model of the PDC light curve of Kepler-30.

Fig. 17 Hurst surface h(q, τ) calculated for final PDC time-series. The (h, q) plane corresponds to h(q) calculated with the standard MFDMA method. The colour bar indicates values of h(q, τ), where for q < 0 the higher values are found and for q > 0 the lower ones.

Fig. 18 As in Fig. 17 but for SAP data.

7.3.2 Hurst exponent surfaces

Now, we investigate the effects of a periodic trend using the Hurst surface h(q, τ). With this, we can study the generalised Hurst exponent h behaviour not only as afunction of the q-moment (standard MFDMA method) but also as a function of the timescale τ (MFDMAτ). The Hurst surfaces illustrated in Figs. 17, 18, and 19 calculated for PDC, SAP, and RTS data show abundant features that may be hidden by the standard MFDMA proposed by Gu & Zhou (2010).

For a fixed scale n, when q changes from − 5 to 5, there are downward trends for all h(q, τ) surfaces, except in the domain where n < 8.7 (n < 23.5) days and q < 0 for PDC and SAP (RTS) series. In particular, we can recover the results from the standard MFDMA method, calculating the scaling exponent for each q at the whole scale range n. In this way, varying q, when τ changes from small to large values, all curves of h(q > 0, τ) show a rapid rise to their peak values, then become more or less constant, showing a plateau. However, for h(q < 0, τ), there are oscillations that increase with the increase of |q|, in addition to a stronger jump when τ is small and q < 0 for PDC and SAP and less pronounced for RTS. The presence of bumps at negative q values in the Hurst surfaces in Figs. 17, 18, and 19 may be the result of the amplified impact of small fluctuations that in turn depend on several effects such as the evolution of active regions and the differential rotation.

The main feature of the h(t, τ) plots is the remarkable decrease in the value of h from ≈1 for q ~ −5 to ~ 0.5 for q > 0. This suggests that the small-scale fluctuations, the effect of which on h is amplified for q < 0, have a strong degree of persistence, while the large-scale fluctuations, which dominate h for q > 0 are more similar to an uncorrelated random process. The results based on the portion of the plots with q < 0 should be taken with some caution because of the sizeable fluctuations of h, but these plots suggest a different physical origin for small- and large-scale fluctuations in the Kepler-30 light curve, a conclusion that can be useful to constrain models of stellar microvariability. Crossovers are not clearly evident in the h(t, τ) plots of the SAP and PDC time-series, except for that at τ = 8.7 days already discussed above. This probably happens because the variability is dominated by the rotational modulation. Therefore, the SAP and PDC h(t, τ) surfaces show only a tiny hint of possible slope changes for τ between ≈20 and ≈35 days, which could be associated with the evolution of active regions. On the other hand, the h(t, τ) plot of the RTS time-series shows a clear crossover at ~23 days that extends over an interval of q and is likely indicating the typical active regions evolution timescale as estimated from the spot modelling.

Following the discussion in Sect. 5.2 and in de Freitas et al. (2017), we computed the ${\langle S_{\text{ph}}\rangle }_{k=5}$ activity index for Kepler-30 finding 4712 and 5209 parts per million for the PDC and SAP pipelines, respectively. The larger value found with the SAP data is a consequence of the reduction of the stellar variability in the PDC time-series for timescales longer than ~ 15 days. The Hurst exponent surface h(2, τ) shows a steep increase up to τ = 8.7 days followed by an almost constant plateau, indicating that most of the variability in the time-series of Kepler-30 occurs for timescales shorter thanhalf the rotation period of the star which corresponds to only 10% of the time interval sampled by the ${\langle S_{\text{ph}}\rangle }_{k=5}$ index. Therefore, such an index provides a very limited description of the light variability in the case of Kepler-30, while the Hurst exponent surfaces allow a detailed description of the dependence of the variability on the timescale. The short-term spot cycle of ~ 34 days of Kepler-30 has a period shorter that the five rotations adopted to compute the ${\langle S_{\text{ph}}\rangle }_{k=5}$ index, therefore this index may not be appropriate to sample the activity timescales characteristic of Kepler-30.

Fig. 19 As in Fig. 17 but for RTS data. In this case, random walk-type fluctuations (H > 1) are stronger than PDC and SAP data, indicating more apparent slowly evolving fluctuations in regime q < 0.

7.3.3 Effect on surrogate time-series

As mentioned in Sect. 5.2, multifractality can have two sources: long-term correlations or a fat-tailed probability distribution of the fluctuations. In the case of Kepler time-series, it was already shown that the first source of multifractality is by far the dominant one (de Freitas et al. 2017). We confirm this result in the specific case of Kepler-30 by calculating h(q, τ) for the shuffled and phase-randomised surrogates of the three time-series considered above.

Figure 20 shows the average results of 200 realisations of the shuffled and phase-randomised surrogates of PDC (upper panels), SAP (middle panels), and RTS (bottom panels). In fact, as can be seen in Fig. 20, the shuffling procedure destroys the correlations, that is, the Hurst surface h(q, τ) is flat (⟨h(q, τ) ≈ 0.5⟩, the value of h corresponding to white noise). However, the Hurst surface of phase-randomised series vary slightly, both with the order q and timescale τ (⟨h(q, τ) ≈ 0.82⟩). These findings suggest that the multifractality of rotational modulation is due to both long-range correlation and non-linearity in turn due to fat-tailed probability distributions, but long-range correlations are the main source of multifractality, which is consistent with the results of standard MFDMA.

8 Summary and conclusions

In this study, we analyse the rotation and evolution of photospheric active regions of the moderately young Sun-like star Kepler-30 accompanied by a three-planet system, using both its PDC and SAP time-series by means of an MFDMA-based multifractality analysis approach in a standard and a new multiscale version. In the latter case, the PDC and SAP time-series, as well as the difference RTS data, are examined considering the generalised dependence of the local Hurst exponent on the timescale τ by means of the surface h(q, τ). We also consider the impact of the rotational modulation on the characterisation of the multifractal properties of the light fluctuations, as already investigated by de Freitas et al. (2013b, 2016, 2017, 2019a,b), and show that such an impact depends on the timescale. Furthermore, we applied a maximum entropy spot modelling to extract information on the longitude, the area variation, and the evolutionary timescales of the active regions responsible for the rotational modulation of the stellar flux. Such an approach reveals that the characteristic timescales of stellar activity in Kepler-30 are significantly shorter than the five rotations adopted to compute the ${\langle S_{\text{ph}}\rangle }_{k=5}$ index as defined by Mathur et al. (2014), and therefore such an index is of limited use to characterise the activity of our target.

Our main conclusions can be summarised as follows:

• (i)

  The fluctuation functions F_q(n) show that the multifractal properties of the Kepler-30 time-series have a relationship with the range of scale n, indicating the limitation of the standard MFDMA method which uses a fixed timescale. Then, we systematically investigate the dynamic behaviours of the three time-series, PDC, SAP, and RTS, by applying a new approach here named MFDMAτ;

• (ii)

  The Hurst surfaces reveal that for negative q values there are remarkable fluctuations in the local Hurst exponent h(q, τ) and significant differences for different values of q. Conversely, the positive q values showthat the Hurst surfaces become flat starting from a minimum period (~ 8.7 days) that corresponds to the first harmonic of the rotation period as found by the Lomb–Scargle periodogram in the case of the SAP and PDC time-series, while it corresponds to the evolutionary timescale of active regions (~ 23.5 days) for the RTS time-series, in which the rotational modulation has been almost completely removed. The analysis of the residuals of the spot modelling shows a crossover at ~1 day that coincides with the characteristic turnover time of the supergranules (cf. Sect. 7.3.1), and another at ~ 35 days that corresponds to an increase in the power level of the light fluctuations and a possible short-term cycle in the total area of the starspots (see below), respectively;

• (iii)

  The multifractality of the Kepler-30 time-series is principally due to the long-range correlations with a minor contribution from a broad non-Gaussian probability density distribution of the fluctuations. This result is found by comparing the original time-series with their shuffled and phase-randomised surrogates. Our h(t, τ) plots also suggest that the small-scale fluctuations that dominate the function F_q for q < 0 have a remarkable persistence, while the large-scale fluctuations dominating for q > 0 have a random and uncorrelated behaviour similar to that of white noise, once the effect of the rotational modulation has been removed;

• (iv)

  A maximum entropy spot modelling shows that the photospheric features of Kepler-30 evolve on timescales ranging from 10 to 20 days for individual active regions up to a few hundred days for the longer lived active longitudes. Their migration can be used to estimate a lower limit for the relative surface shear of ΔΩ∕Ω ~ 0.02 ± 0.01, which should be taken with great caution owing to the rapid evolution of the individual starspots, which can mimic the effects attributed to differential rotation. A short-term cycle of about ~ 34 days in the total area of the starspots could be present and its timescale compares well with that found in the difference RTS time-series as well as with a crossover (a slope change) in the fluctuation functions of the residuals of the spot modelling;

• (v)

  Finally, we note a coincidental proximity of some activity timescales with the synodic periods of the two closest planets. The short-term spot cycle of ~34 days is close to the synodic period of 35.2 days of the planet Kepler-30b when a mean rotation period of 16 days is adopted for Kepler-30. From the Lomb–Scargle periodogram of the RTS time-series, a 23.1-day period emerges that is close to the synodic period of Kepler-30c of 22 days. This period also appears in the multifractal analysis as a crossover of the fluctuation functions and is likely associated with the characteristic evolutionary timescales of active regions in Kepler-30 as indicated by the spot modelling. For the most distant planet Kepler-30d, we did not identify any similar coincidence. Nevertheless, the large separations of planets b and c suggest that great caution should be taken when attributing such coincidences to a possible star–planet magnetic interaction (see Sect. 7.2 for more detail).

We point out the relevant timescales and the persistence characteristics of the light fluctuations of Kepler-30, an active Sun-like star. This information can provide constraints for the models of the stellar flux variations based on the effects of magnetic fields and surface convection, thus promising to contribute to the refinement of those models in future investigations.

As a perspective, our multifractal methods are particularly interesting for analysing ongoing TESS and future PLATO high-precision photometric time-series. They represent a useful complement to regularised spot models, and can be used for an investigation of the activity phenomena on timescales comparable with the stellar rotation period or longer, while multi-fractal methods can cover the full range of timescales and characterise the persistence of the fluctuations produced by magnetoconvection on timescales remarkably shorter than the rotation period, which are not accessible to spot modelling. On the other hand, we illustrate how the two methodologies provide comparable results in the range of timescales they have in common. In this case, spot modelling allows us to gain more physical insight into the origin of the characteristic timescales pointed out by the multifractal analysis, at least when they can be associated with starspot evolution or activity cycles. To reduce the degeneracies of spot modelling, it is advisable to target stars whose spin axis inclination can be constrained through the occultations of spots during planetary transits such as in the case of Kepler-30 (Sanchis-Ojeda et al. 2012) or by means of asteroseismology (e.g. Ballot et al. 2006).

Fig. 20 Hurst surface h(q, τ) calculated for the shuffled (left side) and phase-randomised (right side) versions of time-series: PDC (top), SAP (middle), and RTS (bottom). Presented are results averaged over 200 realisations of the surrogate series.

Acknowledgements

The authors are grateful to an anonymous referee for a careful reading of the manuscript and several valuable suggestions that helped them to improve their work. DBdeF acknowledges financial support from the Brazilian agency CNPq-PQ2 (grant no. 311578/2018-7). Research activities of STELLAR TEAM of Federal University of Ceará are supported by continuous grants from the Brazilian agency CNPq. AFL gratefully acknowledges support from the INAF Mainstream Project entitled “Stellar evolution and asteroseismology in the context of PLATO space mission” coordinated by Dr. Santi Cassisi. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. All data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). The authorswould like to dedicate this paper to all the victims of the COVID-19 pandemic around the world.

Appendix A The multifractal background

In the present appendix, we present the steps of the MFDMA algorithm according to Gu & Zhou (2010):

Step 1.

Construct the sequence of cumulative sums of time-series x(t) over time t = 1, 2, 3, …, N, assuming the data points are evenly spaced: $$y(t)={\sum }_{i=1}^{t}x(i),t=1,2,3,\dots ,N,$$(A.1)

where N is the total number of data points in the time-series and the index i is the time index t.

Step 2.

Calculate the moving average function of Eq. (A.1) in a moving window: $$\tilde{y}(t)=\frac{1}{n}{\sum }_{k=0}^{\lceil n-1\rceil }y(t-k),$$(A.2)

where n is the window size, and $\lceil (x)\rceil$ is the smallest integer that is not smaller than argument (x)⁵;

Step 3.

Detrend the series by removing the moving average function, ỹ (i), and obtain the residual sequence, ϵ(i), through: $$ϵ(i)=y(i)-\tilde{y}(i),$$(A.3)

where n ≤ i ≤ N. We divide the residual series ϵ(i) into N_n = int[N∕n − 1] non-overlapping segments of the same size n, where int[x] is the largest integer that is not larger than x. Each segment can be denoted by ϵ_ν so that ϵ_ν(i) = ϵ(l + i) where 1 ≤ i ≤ n and l = (ν − 1)n;

Step 4.

Calculate the root-mean-square (RMS) fluctuation function, F_ν (n), for a segment of size n: $$F_{\nu }(n)={\left\{\frac{1}{n}{\sum }_{i=1}^n{ϵ}_{\nu }^2(i)\right\}}^{\frac{1}{2}}.$$(A.4)

Step 5.

Generate the fluctuation function F_q(n) of the qth order: $$F_q(n)={\left\{\frac{1}{N_n}{\sum }_{\nu =1}^{N_n}{F}_{\nu }^q(n)\right\}}^{\frac{1}{q}},$$(A.5)

for all q≠0, where the qth-order function is the statistical moment (e.g. for q = 2, we have the variance), while for q = 0, $$\ln \left[F_0(n)\right]=\frac{1}{N_n}{\sum }_{\nu =1}^{N_n}\ln [F_{\nu }(n)].$$(A.6)

The scaling behaviour of F_q(n) follows the relationship $$F_q(n)~n^{h(q)},$$(A.7)

where h(q) denotes the Holder exponent or generalised Hurst exponent. Each value of q yields a slope h and, in particular, q = 2 gives the classical Hurst exponent.

Step 6.

Knowing h(q), the multifractal scaling exponent, τ(q), can be computed: $$\tau (q)=qh(q)-1.$$(A.8)

Finally, the singularity strength function, α(q), and the multifractal spectrum, f(α), are obtained via a Legendre transform: $$\alpha (q)=\frac{\text{d}\tau (q)}{\text{d}q}$$(A.9)

and $$f(\alpha )=q\alpha -\tau (q).$$(A.10)

By definition, for a monofractal signal, h is the same for all values of q and is equal to α giving a single-value spectrum f(α) = f(h) = 1, while for a multifractal signal, h(q) is a function of q, and the multifractal spectrum is generally parabolic (see Fig. 2 from de Freitas et al. 2017).

We use the following model parameters to compute the multifractal spectrum, as recommended by Gu & Zhou (2010): N = 30; q ∈ [−5, 5] with a step size of 0.2; the lower bound of segment size n, which is denoted n_min and set to 10; and the upper bound of segment size n, which is denoted n_max and is given by N∕10. An estimateof the standard deviation of the Hurst exponent is provided by Kantelhardt (2015) and is < 0.03 when we consider time-series with more than 10⁴ datapoints, as in our case.

A.1 Multifractal indicators

We use a set of four multifractal indicators that can be extracted from the quantities defined in Eqs. (A.8), (A.9) and (A.10), in particular from the multifractal spectrum f(α) as illustrated in Fig. 2 of de Freitas et al. (2017). In this paper, we refer to this same figure to describe the shape of the multifractal spectrum. Here, we show a list of the main indicators:

Parameter α₀

The α₀ parameter is the value of α corresponding the maximum of the multifractal spectrum f(α). It delivers information about the structure of the process producing the fluctuations, with a high value indicating that it is less correlated and has a fine structure (Krzyszczak et al. 2018). In addition, this parameter is strongly affected by signal variability. This is evidenced when one investigates the different sources of multifractality that are present in an astrophysical signal (see Sect. 5.2).

Singularity parameter Δf_min(α)

The parameter Δf_min(α) characterises the broadness, which is defined as the difference f(α_max) − f(α_min) of the singularity spectrum. This difference provides an estimate of the spread in changes in fractal patterns. Consequently, if Δf_min(α) is positive, the left-hand side of the spectrum is less deep, while a negative value indicates that this side is deeper. On the other hand, when Δf_min(α) is null, the depths of the tails are the same on both sides. As quoted by Ihlen (2012), larger fluctuations (q > 0) imply that the singularities are stronger, whereas the smaller fluctuations (q < 0) indicate that the singularities are weaker (cf. Tanna & Pathak 2014). In other words, multifractal spectra that have a longer right tail than the left one reveal that the structure of the time-series is more regular and less dominated by extreme (maximal) values and therefore the parameter Δf_min(α) can be a useful way to estimate the impact of noise in the periodic signal.

Degree of asymmetry (A)

This index is defined as the ratio: $$A=\frac{{\alpha }_{\text{max}}-{\alpha }_0}{{\alpha }_0-{\alpha }_{\text{min}}},$$(A.11)

where α₀ is the value of α where f(α) is maximal, while α_min and α_max are the minimum and maximum values of the singularity exponent α as defined by Eq. (A.9), respectively. The value of A indicates one of three possible skewness of the singularity spectrum: right-skewed (A > 1), left-skewed (0 < A < 1), or symmetric (A = 1).

Degree of multifractality (Δα)

This index represents the measure of the interval: $$\Delta \alpha ={\alpha }_{\text{max}}-{\alpha }_{\text{min}},$$(A.12)

where α_max and α_min are defined above. A low value of Δα indicates that the time-series is close to fractal with the multifractal strength being higher when Δα increases (de Freitas & De Medeiros 2009; de Freitas et al. 2017). As mentioned by de Freitas et al. (2016), larger values of Δα denote more complex fluctuations, whereas smaller values indicate that the spectrum tends towards the monofractal limit. According to Makowiec & Fuliński (2010), if Δα is less than 0.05 a monofractal behaviour of the spectrum should be assumed because of the intrinsic precision in deriving such a parameter from a statistics based on a number of datapoints that is in any case limited.

A.2 Results based on the standard multifractal approach

First, we calculate the MFDMA fluctuation functions F_q(n) as a function of window size n (in days) for the three time-series PDC, SAP, RTS. The scale parameter n is varied from 10 to N∕10, and the exponent q is varied from − 5 to 5 in steps of 0.2. The scaling pattern of F_q(n) of original data (red lines) for q = 2 is shown in the top left panels of Figs. A.1, A.2 and A.3. We repeat the analysis for a set of 200 randomly shuffled series as well as for 200 phase-randomised series (blue and green lines, respectively).

For PDC, SAP, and RTS from Kepler-30, the standard multifractal characterisation was performed following the procedures described above and also for the series obtained from surrogates of the original data sets. For all of the time-series, the above four multifractal indicators together with the Hurst exponent H were obtained, and their values are shown in Table 1 where the parameters associated with the original time-series are denoted with an (O) superscript, while those associated with the shuffled and the phase-randomised series with an (S) or (P) superscripts, respectively.

Fig. A.1 Multifractal analysis for PDC time-series following steps 5 and 6 presented in Sect. 4. Top panel: original (in red), shuffled (green), and phase-randomised (blue) data are based on the procedure mentioned in Sect. 4. Left top: multifractal fluctuation function Fq (n) obtained from MFDMA method for only q = 2, indicated henceforth as a big circle. Right top: q-order Hurst exponent (h(q)) as a functionof q-parameter. This panel shows the truncation originated from the leveling of the h(q) for positive q. Left bottom: comparison of the multifractal scaling exponent τ(q) of three data. In this panel, it is possible to identify a crossover in q ~−1. Right bottom: multifractal spectrum f(α) of three time-series, respectively.

Fig. A.2 As in Fig. A.1 but for SAP time-series.

Fig. A.3 As in Fig. A.1 but for RTS data.

Fig. A.4 Deviations of scaling exponent for (top left) PDC, (top right) SAP, and (bottom) RTS data.

For all the series and q-domains, the multifractality due to correlation is stronger than that due to the non-Gaussianity of the distribution of the fluctuations with fat tails as indicated by ‘1’ and the positive values of $\Delta f_{\text{min}}^O(\alpha )$. This value does not imply that there are only long-term correlations in fluctuations, but it implies that non-linearities due to fat-tailed distributions are very weak. As a consequence, a time-independent structure with h(q) ≈ 0.5 is shown in the shuffled time-series (see Figs. A.1, A.2 and A.3). After shuffling, all time-series exhibit a smaller degree of multifractality than the original one. This result can be emphasised by the multifractal spectra (green curves from right top panels). In general, this analysis is not conclusive for explaining the source of this behaviour. In this case, it is necessary to verify the behaviour of τ(q). As illustrated by Fig. A.4, there is a strong dependence of h(q)− h_shuf(q) on q, which is clearest for q < 0, where a deeper right tail occurs in multifractal spectra (see right bottom panels from Figs. A.1, A.2 and A.3). In contrast, the dependence on q can be neglected for the phase-randomised data, except for q < 0 in PDC andRTS data as illustrated by Fig. A.4. Thereby, the two sources of multifractalities appear in q < 0, the domain where the rotational modulation occurs. Basically, the correlations and non-linearities are negligible for q > 0, where the strong fluctuations due to the noise (correlated or not) appear.

Because the presence of rotational modulation and noise can affect the multifractal indicators, we decided to investigate the changes of the multifractal spectrum and compare the indicators calculated from the original series with those obtained from the surrogate series. Consequently, it is possible to find the source(s) that affect the values of α₀, H, Δα and A following the behaviour of their calculated values from the surrogate series.

Firstly, the relation between the values of α₀ in the three time-series show a clear difference among them. It is interesting to note that SAP data have larger ${\alpha }_0^{\text{O}}$ values than the PDC time-series, which means that in PDC data the fluctuations governing the rotational modulation are more correlated and have a less fine structure (a structure more regular in appearance) compared to the fluctuations of the SAP data. For RTS data, the fluctuations are even less correlated and have a more fine structure (a structure less regular in appearance) compared to the processes governing the PDC and SAP data. Secondly, the values of H for the surrogate data indicate that shuffling the data destroys the correlations, and therefore H tends to 0.5, whereas the phase-randomised data recover a value very close to that found for the original series.

In contrast, the degree of multifractality Δα changes more broadly. Because this parameter is connected to the richness of the data structure, we highlight that the Δα for the original time-series of RTS indicates such data may promote some values of the fluctuations, making the signal structure richer. Nevertheless, there is an important detail. As $\Delta f_{\text{min}}^{\text{O}}(\alpha )$ is the smallest value, the broadness of Δα is a mixing between strong and weak fluctuations, and therefore the noise has an important effect over the data, as can be emphasised by the small values of Δα_S. However, following the criterion of Makowiec & Fuliński (2010), the shuffling process does not reduce the series to a monofractal.

It is possible to observe a dissimilarity for the asymmetry parameter A_O among the time-series for both the original and the surrogate data. As shown in Figs. A.1, A.2 and A.3, and Table A.1, it is interesting that generally the spectra are rather right-skewed (which suggests that fine structures are more frequent) or tend to be symmetrical in shape. However, the extreme events become stronger for the RTS data, as the left side of the spectrum of this series is deeper than those from PDC and SAP data.

In conclusion, it is interesting to compare the parameters of the multifractal spectra between the two time-series, namely PDC and SAP. This allows us to address the question of whether the MFDMA method can be applied as an indicator of the changes in the dynamics of fluctuations and therefore of processes occurring inthe stellar atmosphere from which the stellar flux originates. It can be observed that the values of α₀ change only slightly and are greater in the SAP data. On the other hand, the degree of multifractality is more developed in the PDC data, whereas the asymmetry for SAP time-series is slightly more positive. Some differences between both time-series can also be seen when analysing the absolute differences of Hurst exponents for the original and shuffled data or original and surrogates (see Fig. A.4). Even though we do not see a change between both time-series, the source of multifractality due to the contribution of long-range correlations can be seen to be dominant.

Appendix B Fluctuation functions for a sinusoidal signal

We consider a noiseless, purely sinusoidal signal with constant amplitude and phase to study how the fluctuation functions computed with the MFDMA behave under the effect of a sinusoidal modulation. Suyal et al. (2009) investigated a similar case using a re-scaled range analysis, the so-called R∕S method (e.g. Kantelhardt 2015), to evaluate the fluctuations function, but did not discuss the case of the MFDMA method. In our case, we see that the F_q (n) fluctuation functions, plotted as a function of n (see Fig. B.1), show a linear trend followed by an oscillating phase, a behaviour similar to that shown in Figs. 13, 14, and 15. Other factors can lead to deviations from this ideal behaviour, such as noise and/or the superposition of other low-amplitude periodicities. The combination of these contributions should make the simulated time-series closer to the observed ones.

Fig. B.1 Log-log plot of fluctuation functions Fq(n) vs. n for a sinusoidal time-series (shown in the subplot) with a period of ~6.3 days.

Appendix C Spot modelling

To map the photosphere of Kepler-30, we subdivide it into 200 surface elements of side 18° × 18° and assumethat in each element there are dark spots with a filling factor f, solar-like faculae with a filling factor Qf, and unperturbed photosphere with a fillingfactor 1 − (Q + 1)f, where Q is the ratio of the area of the faculae to that of the spots in each element. For simplicity, Q is assumed constant in our model.

The contrast of the dark spots is defined as c_s ≡ I_spot∕I_u, where I_spot is the brightness of the spotted photosphere and I_u that of the unperturbed photosphere; c_s is assumed to be constant. In our case, c_s = 0.85 in the Kepler passband has been derived from the occulted spots producing characteristic bumps along the photometric profiles of the transits (Sanchis-Ojeda et al. 2012).

The facular contrast is assumed to be zero at the centre of the stellar disc and maximum at the limb as we observe in the Sun, that is, c_f = c_{f 0}(1 − μ), where c_f0 = 1.115 is the contrast at the limb and μ ≡ cosψ, where ψ is the angle between the normal to the surface at a given point and the line of sight. We note that the effect of the faculae is parametrised by the product Qc_f0, and thus their contrast and their area ratio are not independent parameters (cf. Lanza 2016). Therefore, we keep c_{f 0} constant and very Q (see below).

The limb darkening of the unperturbed photosphere is expressed by means of a quadratic law: $$I_{\text{u}}(\mu )\propto a_{\text{p}}+b_{\text{p}}\mu +c_{\text{p}}{\mu }^2,$$(C.1)

where a_p, b_p, and c_p are derived from the limb-darkening coefficients given by Sanchis-Ojeda et al. (2012).

Our model assumes that the distribution of the filling factor f over the surface of the star is fixed while fitting the light curve. Because spots are evolving, this means that we cannot apply our model to fit the entire photometric time-series, but we must fit individual intervals of duration Δt_f over which the hypothesis of a fixed spot pattern is satisfied. In stars with a slowly evolving spot pattern, the best choice is to take Δt_f equal to the mean rotation period. This will give a uniform sampling of all the longitudes along each individual rotation (Lanza et al. 2007).

Unfortunately, this is not the case for Kepler-30 because its starspots evolve quite rapidly, thus imposing a shorter time interval to adequately fit its light modulation. Following the method applied in previous modelling of other CoRoT or Kepler targets (e.g. Lanza et al. 2019), we derive the best values of Δt_f and Q by applying a simplified spot model assuming only three discrete spots plus a uniform background (Lanza et al. 2003).

The value of Δt_f is derived by considering that in general the chi square of the best fit to the entire light curve decreases with the decrease of Δt_f. However, the reduction of the chi square becomes small when Δt_f becomes equal to or shorter than the typical timescale of spot evolution. Therefore, we progressively decrease Δt_f until the decrease of the chi square is no longer significant according to a statistical test based on the Fisher-Snedecor statistics and in this way we select the optimal Δt_f. We note that for each trial value ofΔt_f, we compute models with different values of Q spanning a large interval to select the minimum of the chi square. This is made possible by the relatively short CPU time required to compute the best fit with the three-spot model in comparison to the full model assuming a continuous spot distribution.

The value of Q minimising the chi square for the optimal value of Δt_f is adopted for our spot modelling. In the case of Kepler-30, we find that the optimal Δt_f = 11.963 days, that is, ~ 0.75 of the mean rotation period, while Q = 0.5, that is much smaller than the value adopted for the spot modelling of the Sun as a star by Lanza et al. (2007), where Q_⊙ = 9.0. This is in agreement with the finding that dark spots dominate the optical light modulation of stars significantly more active than the Sun as is the case for Kepler-30 (cf. Radick et al. 2018).

The inclination of the stellar spin to the line of sight is derived by assuming that the stellar spin is normal to the orbital plane of the planet Kepler-30c that is the largest body in the system after the star. Nevertheless, the orbits of the other planets are virtually coplanar, and so this is equivalent to assuming that the stellar spin is perpendicular to the mean plane of the planetary orbits. We note that the possibility of constraining the inclination of the stellar spin reduces the degeneracies of our spot modelling and is made possible by the presence of transiting planets around the star.

The parameters adopted in our spot modelling are listed in Table C.1. The mass and the radius of the star together with its rotation period are used to compute the surface gravity as a function of the latitude to account for the effects of the gravity darkening. They are of the order of 10⁻⁶ in relative flux units and do not affect our solution.

The effects of the uncertainties in the model parameters on the spot modelling have been discussed in detail by Lanza et al. (2009, 2019). They are negligible for our application because we are mainly interested in the typical timescales of starspot evolution rather than on the absolute values of the spotted area or of the surface differential rotation.

By minimising the chi square of the fit to the light curve with our continuous spot model, we can find a unique spot map for Kepler-30. However, such a map is unstable in the sense that a small change in the data produces a large change in the map because most of the map details come from fitting the photometric noise and the evolution of small spots along each Δt_f time interval. Regularisation can provide a unique map by reducing the fitting of the noise and of the small fluctuations and therefore the consequent instability. This is achieved by imposing a priori assumptions on the spot map that are coded into a regularisation functional S(f) that depends on the vector f, the components of which are the filling factors of the 200 surface elements of our map (see Eq. (3.25) in Lanza 2016, for the expression of S).

Fig. C.1 Distribution of the spot filling factor vs. the longitude and time for the ME regularised spot model of the SAP light curve of Kepler-30. The minimum of the filling factor corresponds to dark blue regions, while the maximum is rendered in orange (see the colour scale close to the lower right corner of the plot). We note that the longitude scale of the horizontal axis is extended beyond the [0°, 360°] interval to better follow the migration of the starspots.

In the case of a regularised best fit, instead of minimising the chi square χ² between theobserved and the modelled fluxes, we minimise a linear combination of the χ² and the regularising functional S, as: $$Z(f)={\chi }^2(f)-{\lambda }_{\text{ME}}S(f),$$(C.2)

where λ_ME > 0 is a Lagrangian multiplier that rules the level of regularisation. When λ_ME is zero, we get the unregularised model that has the minimum χ² and residuals that are symmetrically distributed around a zero mean μ_res. When we increase λ_ME, we increase the χ² above its minimum and the mean of the residuals is no longer zero because the ME regularisation drives the map towards an unspotted photosphere, thus making the residuals systematically negative. In other words, the functional S is designed to be maximal when the star is unspotted, so that, by introducing the regularisation, we drive the map towards that of an unspotted star.

Fig. C.2 Total spotted area on Kepler-30 as derived from the spot modelling of the SAP light curve vs. the time (green dots). A sinusoid with a period of 33.878 days corresponding to the maximum of the GLS periodogram is overplotted (red solid line).

The criterion to fix the Lagrangian multiplier is crucial to decide when the process of regularisation has to be stopped to prevent the χ² from becoming too large and the fit unacceptable. As in the case of similar stars with a photon shot noise of the single datapoint of the order of ~ 1% of the amplitude of the flux modulation produced by starspots, we iteratively increase λ_ME until the absolute value of the mean of the residuals |μ_res| becomes equal to the standard error ϵ_st of the data points in the fitted interval of duration Δt_f. The standard error is computed from the standard deviation σ₀ of the residuals of the unregularised best fit and the number of data points M in each interval Δt_f as ${ϵ}_{\text{st}}={\sigma }_0/\sqrt{M}$. For each interval Δt_f, we determine λ_ME by enforcing the equality |μ_res| = ϵ_st within 5% both for the fit of the PDC and the SAP light curves.

The unregularised best fit of the SAP light curve obtained with our model is very similar to that of the PDC time-series and is not shown here. The regularised spot map is shown in Fig. C.1 and is similar to that obtained with the PDC light curve. The SAP light curve shows a range of amplitudes of the light modulation larger than the PDC light curve because the Kepler pipeline tends to reduce the variability at timescales longer than 15–20 days. As a consequence, the SAP spot map shows a wider distribution of the filling factor than the PDC map, but the location of the spots and their evolution are very similar to those shownin the PDC map.

The total spotted area as derived by the ME regularised spot model of the SAP light curve is plotted in Fig. C.2 as a function of the time. As in the case of the PDC light curve, only the intervals Δt_f with data points of sufficiently uniform distribution are considered to avoid systematic effects produced by the ME regularisation (see Sect. 7.2). A modulation with a period of ~ 33.9 days is apparent as indicated by the sinusoid superposed to the plot.

References

All Tables

All Figures

	Fig. 1 Kepler-30 PDC time-series. Top panel: PDC time-series (black dots) and polynomial third-order fits. Bottom panel: final PDC time-series adjusted by polynomial functions. The signatures of three planets are maintained.
In the text

	Fig. 2 As in Fig. 1 but for SAP time-series.
In the text

Fig. 3 Detrending of Kepler-30 with the iterative biweight filter using a window size of 0.75 days (red line), which allows detection of the transits. Normalised (top), detrended (middle), and residual flux (bottom), shown for only a segment of the total time-series. Black dots indicate the PDC pipeline of Kepler-30 after the first three steps of our procedure. Detrended flux is a combination of noise background and planet transits, indicating the trending of the time-series potentially due to rotational modulation of the star.

In the text

	Fig. 4 Same as Fig. 3 for SAP pipeline.
In the text

Fig. 5 Final PDC (first panel) and SAP (third panel) time-series adjusted according to the steps of Sect. 3. The corresponding Lomb–Scargle periodograms are shown at the bottom. The dashed-dotted blue line marks the 99% confidence level. For SAP data, there is a significant peak to 46.6 days also found in Fig. 6. As can be seen, the PDC time-series has a flat periodogram for periods longer than 20 days confirming that the PDC pipeline tends to remove all the long-term trends.

In the text

	Fig. 6 Difference between the SAP and PDC time-series, where a peak of 46.6 days can be seen. Second and third peaks of 23.1 and 33.9 days are also seen.
In the text

	Fig. 7 Time-delayed PDC time-series x(t + τ) with the optimal delay τ = 8.7 days vs. the original time series x(t). The subplot shows that τ = 8.7 days corresponds to the maximum value of the autocorrelation function of the time-series.
In the text

	Fig. 8 Upper panels: PDC light curve of Kepler-30 (solid black dots) and the unregularised best fit obtained with our spot model (red solid line). Lower panels: residuals of the best fit.
In the text

	Fig. 9 Distribution of the residuals of the unregularised best fit to the PDC light curve in Fig. 8.
In the text

	Fig. 10 GLS periodogram of the residuals of the unregularised best fit to the PDC light curve in Fig. 8.
In the text

Fig. 11 Distribution of the spot filling factor vs. the longitude and time for the ME regularised spot model of the PDC light curve of Kepler-30. The minimum of the filling factor corresponds to dark blue regions, while the maximum is rendered in orange (see the colour scale close to the lower right corner of the plot). We note that the longitude scale of the horizontal axis is extended beyond the [0°, 360°] interval to better follow the migration of the starspots.

In the text

	Fig. 12 Total spotted area on Kepler-30 as derived from the spot modelling of the PDC light curve versus the time (green dots). A sinusoid with a period of 33.784 days corresponding to the maximum of the GLS periodogram is overplotted (red solid line).
In the text

Fig. 13 Log-log plot of the fluctuation functions Fq(n) (black circle for q = 2) calculated for the final PDC time series presented in Fig. 5. The red curves correspond to q between −5 and 5 in steps of 0.2. Vertical dashed lines mark three domains of the fitting windows for the small n between 29.4min and 8.7 days, for middle n between 8.7 and 16.8 days and, for the large n greater than 16.8 days. The green dashed line gives the average slope H = h(2) = 0.95 for timescales shorter than 8.7 days.

In the text

	Fig. 14 As in Fig. 13 but for final SAP time-series.
In the text

	Fig. 15 As in Fig. 13 but for final RTS time-series. In this case, only two domains are considered and separated by the vertical dashed line at 23.5 days.
In the text

	Fig. 16 Fluctuation functions of the residuals of the unregularised spot model of the PDC light curve of Kepler-30.
In the text

	Fig. 17 Hurst surface h(q, τ) calculated for final PDC time-series. The (h, q) plane corresponds to h(q) calculated with the standard MFDMA method. The colour bar indicates values of h(q, τ), where for q < 0 the higher values are found and for q > 0 the lower ones.
In the text

	Fig. 18 As in Fig. 17 but for SAP data.
In the text

	Fig. 19 As in Fig. 17 but for RTS data. In this case, random walk-type fluctuations (H > 1) are stronger than PDC and SAP data, indicating more apparent slowly evolving fluctuations in regime q < 0.
In the text

	Fig. 20 Hurst surface h(q, τ) calculated for the shuffled (left side) and phase-randomised (right side) versions of time-series: PDC (top), SAP (middle), and RTS (bottom). Presented are results averaged over 200 realisations of the surrogate series.
In the text

Fig. A.1 Multifractal analysis for PDC time-series following steps 5 and 6 presented in Sect. 4. Top panel: original (in red), shuffled (green), and phase-randomised (blue) data are based on the procedure mentioned in Sect. 4. Left top: multifractal fluctuation function Fq (n) obtained from MFDMA method for only q = 2, indicated henceforth as a big circle. Right top: q-order Hurst exponent (h(q)) as a functionof q-parameter. This panel shows the truncation originated from the leveling of the h(q) for positive q. Left bottom: comparison of the multifractal scaling exponent τ(q) of three data. In this panel, it is possible to identify a crossover in q ~−1. Right bottom: multifractal spectrum f(α) of three time-series, respectively.

In the text

	Fig. A.2 As in Fig. A.1 but for SAP time-series.
In the text

	Fig. A.3 As in Fig. A.1 but for RTS data.
In the text

	Fig. A.4 Deviations of scaling exponent for (top left) PDC, (top right) SAP, and (bottom) RTS data.
In the text

	Fig. B.1 Log-log plot of fluctuation functions Fq(n) vs. n for a sinusoidal time-series (shown in the subplot) with a period of ~6.3 days.
In the text

Fig. C.1 Distribution of the spot filling factor vs. the longitude and time for the ME regularised spot model of the SAP light curve of Kepler-30. The minimum of the filling factor corresponds to dark blue regions, while the maximum is rendered in orange (see the colour scale close to the lower right corner of the plot). We note that the longitude scale of the horizontal axis is extended beyond the [0°, 360°] interval to better follow the migration of the starspots.

In the text

	Fig. C.2 Total spotted area on Kepler-30 as derived from the spot modelling of the SAP light curve vs. the time (green dots). A sinusoid with a period of 33.878 days corresponding to the maximum of the GLS periodogram is overplotted (red solid line).
In the text