Issue 
A&A
Volume 650, June 2021



Article Number  A40  
Number of page(s)  20  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202140287  
Published online  03 June 2021 
Multiscale behaviour of stellar activity and rotation of the planet host Kepler30
^{1}
INAFOsservatorio Astrofisico di Catania,
Via S. Sofia 78,
95123,
Catania,
Italy
^{2}
Departamento de Física, Universidade Federal do Ceará,
Caixa Postal 6030, Campus do Pici,
60455900
Fortaleza,
Ceará,
Brazil
email: danielbrito@fisica.ufc.br
^{3}
Faculdade de Física â Instituto de Ciências Exatas, Universidade Federal do Sul e Sudeste do Pará,
Marabá,
PA 68505080,
Brazil
Received:
5
January
2021
Accepted:
26
March
2021
Context. The Kepler30 system consists of a G dwarf star with a rotation period of ~16 days and three planets orbiting almost coplanar with periods ranging from 29 to 143 days. Kepler30 is a unique target with which to study stellar activity and rotation in a young solarlike star accompanied by a compact planetary system.
Aims. We use about 4 yr of highprecision photometry collected by the Kepler mission to investigate the fluctuations caused by photospheric convection, stellar rotation, and starspot evolution as a function of timescale. Our main goal is to apply methods for the analysis of timeseries to find the timescales of the phenomena that affect the light variations. We correlate those timescales with periodicities in the star and the planetary system.
Methods. We model the flux rotational modulation induced by active regions using spot modelling and apply the Multifractal Detrending Moving Average algorithm in standard and multiscale versions to analyse the behaviour of variability and light fluctuations that can be associated with stellar convection and the evolution of magnetic fields on timescales ranging from less than 1 day up to about 35 days. The light fluctuations produced by stellar activity can be described by the multifractal Hurst index that provides a measure of their persistence.
Results. The spot modelling indicates a lower limit to the relative surface differential rotation of ΔΩ∕Ω ~ 0.02 ± 0.01 and suggests a shortterm cyclic variation in the starspot area with a period of ~34 days, which is close to the synodic period of 35.2 days of the planet Kepler30b. By subtracting the two timeseries of the simple aperture photometry and presearch data conditioning Kepler pipelines, we reduce the rotational modulation and find a 23.1day period close to the synodic period of Kepler30c. This period also appears in the multifractal analysis as a crossover of the fluctuation functions associated with the characteristic evolutionary timescales of the active regions in Kepler30 as confirmed by spot modelling. These procedures and methods may be greatly useful for analysing current TESS and future PLATO data.
Key words: stars: activity / stars: rotation / techniques: photometric / stars: individual: Kepler30 (KOI806)
© 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 solarlike 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 multiplanet systems, one of which is the star orbiting Kepler30 (KOI806), 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 (SanchisOjeda 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 Kepler30b, Kepler30c, and Kepler30d 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. SanchisOjeda et al. 2012), respectively. Based on the values of radius and mass, Panichi et al. (2018) estimated that Kepler30b is a Neptunelike planet rather than a superEarth, Kepler30c is a Jovian planet, while Kepler30d is classified as a Neptunemass planet, with bulk densities given by ~0.96 g cm^{−3}, ~1.71 g cm^{−3}, and ~0.19 g cm^{−3}, 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 Kepler30 almost equatoron. 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 Kepler17, a solarlike 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 outoftransit photometry with those of the spots occulted during transits to validate their spotmodelling 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 shortterm and longterm activity cycles depending on the length of the available timeseries.
Inspired by the peculiar properties of the Kepler30 system and the availability of almost four years of Kepler highprecision 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 timeseries 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 longterm dependency of astrophysical timeseries is a problem that has recently been studied to understand the effect of longmemory 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 timeseries, 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 timeseries (de Freitas et al. 2017, 2019a,b). The MFDMA method filters out the local trends of nonstationary series by subtracting the local means. In addition, the MFDMA method can be used to investigate the local fluctuations in timeseries 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 Kepler30 timeseries as obtained by the presearch 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 Kepler30 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 longcadence (data sampling every ~30 min; Jenkins et al. 2010a) and shortcadence (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 presearch data conditioning (PDC) timeseries. The SAP flux is the sum of all calibrated backgroundsubtracted 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 socalled cotrending 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 outoftransit variations on timescales comparable to the stellar rotation, we make use of the longcadence data of Kepler30 publicly available at the MAST archive^{1}. These cover four years and the mean relative accuracy of each data point is ~ 260 ppm.
3 Kepler30 preparation and preprocessing
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 timeseries, although presenting some residual instrumental effects and longterm systematic trends, is not affected by this problem. Therefore, we investigate the variability of Kepler30 by analysing both the SAP and PDC timeseries.
Because of the flux offsets between adjacent quarters and the longterm instrumental trends within each quarter (cf. Jenkins et al. 2010b), the timeseries 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 timestamp 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 crosscorrelation coefficient for zero lag (i.e. when the timeseries are not shifted) of the series before and after removing the outliers anddatapoints with a nonzero quality flag. To this end, we use the MATLAB function crosscorr^{2}. The values of the coefficients found were 0.92 for the PDC data and 0.94 for the SAP one. The gaps in our timeseries 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 thirdorder polynomial was fitted to the data to detrend the timeseries 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 outoftransit 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 outoftransit 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 SanchisOjeda et al. (2012). We also performed a visual inspection and verified that the obtained outoftransit 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 Kepler30d. 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 timeseries 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 timeseries 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^{−4} in relative flux units. The number of data points is N = 60007 and N = 60 018 for the PDCand SAP timeseries, respectively.
The final PDC and SAP timeseries and their Lomb–Scargle periodograms (cf. Lomb 1976; Scargle 1982) with the significance levels are illustrated in Fig. 5. For both timeseries, there are two main peaks at 7.8 and 16.8 days. Nevertheless, for SAP timeseries, 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 timeseries (RTS; i.e., the difference) between theSAP and PDC timeseries 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 timeseries (RTS). Only a few points were lost during the process of subtracting the timeseries. 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 timeevolution 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 Kepler30 PDC timeseries. Top panel: PDC timeseries (black dots) and polynomial thirdorder fits. Bottom panel: final PDC timeseries adjusted by polynomial functions. The signatures of three planets are maintained. 
Fig. 3 Detrending of Kepler30 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 timeseries. Black dots indicate the PDC pipeline of Kepler30 after the first three steps of our procedure. Detrended flux is a combination of noise background and planet transits, indicating the trending of the timeseries potentially due to rotational modulation of the star. 
Fig. 5 Final PDC (first panel) and SAP (third panel) timeseries adjusted according to the steps of Sect. 3. The corresponding Lomb–Scargle periodograms are shown at the bottom. The dasheddotted 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 timeseries has a flat periodogram for periods longer than 20 days confirming that the PDC pipeline tends to remove all the longterm trends. 
Fig. 6 Difference between the SAP and PDC timeseries, 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 Kepler30 timeseries
The small v sin i ~ 1.94 ± 0.25 km s^{−1} of Kepler30(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 illposed problem because it tries to extract a bidimensional map from a onedimensional timeseries. 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 Kepler30 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 (SanchisOjeda 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 closeby transiting planet Kepler30c, 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 Kepler30 virtually equatoron, 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 twodimensional (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 nonuniqueness of lightcurve inversion.
We include the effect of solarlike 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 Kepler30 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 solarlike 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 lightcurve inversions with the real distributions of the sunspot groups and the spots detected by occultation during planetary transits, respectively (SilvaValio & 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 Kepler30. In consideration of the systematic differences introduced by the Kepler pipeline between the SAP and the PDC timeseries, 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 timeseries
5.1 Basic concepts and definitions
The light variations in the timeseries of Kepler30 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 quasiperiodic and nonstationary 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 timeseries 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 selfsimilarity 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 rescaling of the fluctuations by a factor a^{H}, that is, Δf(t) → a^{H}Δf(at), in order for the rescaled fluctuations to follow the same statistical distribution as the initial ones. A timeseries that obeys this scaling property for an arbitrary factor a is called a selfaffine or selfsimilar timeseries and the exponent H is called the Hurst exponent of the timeseries. Several stochastic fluctuations, including those associated with turbulent convection, satisfy this definition.
The Hurst exponent H characterises the kind of selfsimilarity. The simplest example is the timeseries of the meansquare displacement x^{2} (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∕2}x in order to preserve the statistical distribution of the variable x. In other words, the timeseries x(t) exhibits a selfsimilar 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 selfsimilarity (or scale invariance) is one of the fundamental properties that define fractals.
Selfaffine timeseries 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 selfaffine 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 timeseries 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 selfaffine timeseries x(t) with the variable sampled at equidistant times t_{i}, with i = 1, 2, …, N, we define the timeseries 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 (antipersistence). The degree of persistence can be quantified by the autocorrelation function in the case of a stationary timeseries, 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 timeseries, the autocorrelation is defined as: (1)
where σ^{2} is the variance of the timeseries of the increments Δx_{i}.
When the Δx_{i} are uncorrelated, C(n) = 0 for n > 0. A shortrange 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 longterm correlations, the asymptotic behaviour of C(n), that is, its scaling for large n, can be described by a power law: (2)
where the exponent 0 < γ < 1 characterises the longterm correlations. It can be shown that a longterm correlated, i.e. persistent, behaviour of Δx_{i} leads to a selfaffine 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 timeseries.
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 socalled 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 rootmeansquare (RMS) fluctuation function F_{ν}(n), which is defined by (3)
where ϵ_{ν}(i) is the deviation of the ith datapoint from a suitable moving average introduced to remove the longterm variation in the timeseries (see Appendix A for its precise definition). Starting from the RMS fluctuation function, we define the fluctuation function of order q as (4)
where N_{n} is the number of nonoverlapping segments fora given segment of size n in the timeseries, 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 powerlaw given by (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 selfaffine 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 timeseries can be produced, for example, by simple, nonlinear, recurrence laws, such as the socalled logistic law x_{i+1} = λ_{ll}x_{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 quasiperiodic modulations in the timeseries
The multifractal character of a stochastic timeseries 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 longterm correlations and persistences – that depend on the amplitude of the fluctuations – and/or by their statistical distribution that, even in the absence of longterm 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 socalled surrogates of the original timeseries. 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 longterm correlations, producing a new timeseries with h(q) close to 0.5, that is the Hurst exponent of a random timeseries consisting of uncorrelated points. A second kind of surrogate can be obtained by taking the Fourier transform of the original timeseries 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 timeseries whose points retain most of the correlations present in the original timeseries, but whose distribution is nearly Gaussian, thus effectively eliminating the effect of nonGaussianity 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 timeseries 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 phaserandomised series will be similar. On the other hand, when the dominant source of multifractality is a nonGaussian statistical distribution of the datapoints, the converse will be true.
From a more general point of view, timeseries can result from both deterministic and stochastic processes, the former usually being responsible for longterm 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 longterm trends and shortterm 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 timeseries 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 solarlike 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 timeseries were monofractal. In a subsequent work on the same sample, de Freitas et al. (2016) established that their timeseries 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 nonoverlapping 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 is then defined as the arithmetic mean of all the individual S_{ph} indexes determined from each of the kP_{rot} intervals along the timeseries. 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_{2}(n) ∝ F_{ν}(n) ~ n^{H} which is evaluated with the fluctuations obtained after subtracting a moving average from the timeseries. 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 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 anticorrelation between the two indexes with a Pearson coefficient of − 0.33 when considering the abovementioned 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 timeseries 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 Hindex 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 multiscale 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)^{3} and Tang et al. (2015), is a powerful technique that provides valuable information on the fluctuations of a timeseries. 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 longrange correlations, which is important in our case because multifractality in our timeseries 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 timedependent 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 nonstationarity, including the nonstationarity 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 timeseries Δ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 timeseries 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 timeseries, 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 bidimensional phase space. In Fig. 7, we plot the case corresponding to the optimal value τ = 8.7 days, which shows a repetitive, quasistationary 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 Timedelayed PDC timeseries 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 timeseries. 
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 timeseries of Kepler30 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 Kepler30 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^{−7}) and a standard deviation σ_{0} = 3.78 × 10^{−4} 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^{−6}, but a standard deviation of σ_{0} = 3.746 × 10^{−4} 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 lightcurve 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 Kepler30 (solid black dots) and the unregularised best fit obtained with our spot model (red solid line). Lower panels: residuals of the best fit. 
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 MEregularised 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 Kepler30 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 Kepler30, 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 timeseries (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 falsealarm probability (FAP) of 0.29 according to the analytical formula by Zechmeister & Kürster (2009). The corresponding bestfitting 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 shortterm modulation in the total spotted area of Kepler30, reminiscent of the shortterm cycles found in CoRoT2 (Lanza et al. 2009; Zaqarashvili et al. 2021).
Such a shortterm cycle of ~34 days is relatively close to the synodic period P_{syn} = 35.2 days of the planet Kepler30b as computed with a rotation period P_{rot} = 16.0 days and an orbital period P_{orb} = 29.334 days from . 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.1day period (see Fig. 6) is close to the synodic period of Kepler30c. However, an interpretation in terms of tidal or magnetic starplanet 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^{4}. More precisely, both the tidal torque on the star and the energy density of its coronal magnetic field decrease with distance as , 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 Kepler30. 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 Kepler30 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 Loglog plot of the fluctuation functions F_{q}(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 timeseries, 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 timeseries.
For both the PDC and SAP timeseries, 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 periodone 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 timeseries 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 timeseries 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 Kepler30 because the rotational modulation was remarkably reduced by subtracting the two timeseries 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 timeseries, while it was strongly reduced in the PDC timeseries owing to the detrending applied by the Kepler pipeline (cf. Sect. 3). Therefore, the longterm intrinsic variability stands out in the difference timeseries.
The fluctuation functions of the timeseries 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 timeseries. 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 Sunlike stars (cf. Meunier et al. 2015). The other crossover at ~35 days coincides with the possible shortterm 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 largescale fluctuations is amplified with respect to that of smallscale fluctuations in the computation of the moments of order q > 0. Conversely, smallscale fluctuations will have the most important impact on the fluctuation functions F_{q} (n) when q < 0 because the effect of largescale 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. 15 As in Fig. 13 but for final RTS timeseries. 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 Kepler30. 
Fig. 17 Hurst surface h(q, τ) calculated for final PDC timeseries. 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. 
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 qmoment (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 smallscale fluctuations, the effect of which on h is amplified for q < 0, have a strong degree of persistence, while the largescale 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 largescale fluctuations in the Kepler30 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 timeseries, 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 timeseries 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 activity index for Kepler30 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 timeseries 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 timeseries of Kepler30 occurs for timescales shorter thanhalf the rotation period of the star which corresponds to only 10% of the time interval sampled by the index. Therefore, such an index provides a very limited description of the light variability in the case of Kepler30, while the Hurst exponent surfaces allow a detailed description of the dependence of the variability on the timescale. The shortterm spot cycle of ~ 34 days of Kepler30 has a period shorter that the five rotations adopted to compute the index, therefore this index may not be appropriate to sample the activity timescales characteristic of Kepler30.
Fig. 19 As in Fig. 17 but for RTS data. In this case, random walktype 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 timeseries
As mentioned in Sect. 5.2, multifractality can have two sources: longterm correlations or a fattailed probability distribution of the fluctuations. In the case of Kepler timeseries, 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 Kepler30 by calculating h(q, τ) for the shuffled and phaserandomised surrogates of the three timeseries considered above.
Figure 20 shows the average results of 200 realisations of the shuffled and phaserandomised 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 phaserandomised 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 longrange correlation and nonlinearity in turn due to fattailed probability distributions, but longrange 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 Sunlike star Kepler30 accompanied by a threeplanet system, using both its PDC and SAP timeseries by means of an MFDMAbased multifractality analysis approach in a standard and a new multiscale version. In the latter case, the PDC and SAP timeseries, 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 Kepler30 are significantly shorter than the five rotations adopted to compute the 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 Kepler30 timeseries 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 timeseries, 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 timeseries, while it corresponds to the evolutionary timescale of active regions (~ 23.5 days) for the RTS timeseries, 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 shortterm cycle in the total area of the starspots (see below), respectively;
 (iii)
The multifractality of the Kepler30 timeseries is principally due to the longrange correlations with a minor contribution from a broad nonGaussian probability density distribution of the fluctuations. This result is found by comparing the original timeseries with their shuffled and phaserandomised surrogates. Our h(t, τ) plots also suggest that the smallscale fluctuations that dominate the function F_{q} for q < 0 have a remarkable persistence, while the largescale 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 Kepler30 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 shortterm 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 timeseries 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 shortterm spot cycle of ~34 days is close to the synodic period of 35.2 days of the planet Kepler30b when a mean rotation period of 16 days is adopted for Kepler30. From the Lomb–Scargle periodogram of the RTS timeseries, a 23.1day period emerges that is close to the synodic period of Kepler30c 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 Kepler30 as indicated by the spot modelling. For the most distant planet Kepler30d, 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 Kepler30, an active Sunlike 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 highprecision photometric timeseries. 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 multifractal 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 Kepler30 (SanchisOjeda 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 phaserandomised (right side) versions of timeseries: 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 CNPqPQ2 (grant no. 311578/20187). 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 COVID19 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 timeseries x(t) over time t = 1, 2, 3, …, N, assuming the data points are evenly spaced: (A.1)
where N is the total number of data points in the timeseries and the index i is the time index t.
Step 2.
Calculate the moving average function of Eq. (A.1) in a moving window: (A.2)
where n is the window size, and is the smallest integer that is not smaller than argument (x)^{5};
Step 3.
Detrend the series by removing the moving average function, ỹ (i), and obtain the residual sequence, ϵ(i), through: (A.3)
where n ≤ i ≤ N. We divide the residual series ϵ(i) into N_{n} = int[N∕n − 1] nonoverlapping 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 rootmeansquare (RMS) fluctuation function, F_{ν} (n), for a segment of size n: (A.4)
Step 5.
Generate the fluctuation function F_{q}(n) of the qth order: (A.5)
for all q≠0, where the qthorder function is the statistical moment (e.g. for q = 2, we have the variance), while for q = 0, (A.6)
The scaling behaviour of F_{q}(n) follows the relationship (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: (A.8)
Finally, the singularity strength function, α(q), and the multifractal spectrum, f(α), are obtained via a Legendre transform: (A.9)
By definition, for a monofractal signal, h is the same for all values of q and is equal to α giving a singlevalue 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 timeseries with more than 10^{4} 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 α_{0}
The α_{0} 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 lefthand 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 timeseries 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.
Multifractal indicators extracted from the standard MFDMA method for PDC, SAP, and RTS.
Degree of asymmetry (A)
This index is defined as the ratio: (A.11)
where α_{0} 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: rightskewed (A > 1), leftskewed (0 < A < 1), or symmetric (A = 1).
Degree of multifractality (Δα)
This index represents the measure of the interval: (A.12)
where α_{max} and α_{min} are defined above. A low value of Δα indicates that the timeseries 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 timeseries 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 phaserandomised series (blue and green lines, respectively).
For PDC, SAP, and RTS from Kepler30, 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 timeseries, 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 timeseries are denoted with an (O) superscript, while those associated with the shuffled and the phaserandomised series with an (S) or (P) superscripts, respectively.
Fig. A.1 Multifractal analysis for PDC timeseries following steps 5 and 6 presented in Sect. 4. Top panel: original (in red), shuffled (green), and phaserandomised (blue) data are based on the procedure mentioned in Sect. 4. Left top: multifractal fluctuation function F_{q} (n) obtained from MFDMA method for only q = 2, indicated henceforth as a big circle. Right top: qorder Hurst exponent (h(q)) as a functionof qparameter. 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 timeseries, respectively. 
Fig. A.4 Deviations of scaling exponent for (top left) PDC, (top right) SAP, and (bottom) RTS data. 
For all the series and qdomains, the multifractality due to correlation is stronger than that due to the nonGaussianity of the distribution of the fluctuations with fat tails as indicated by ‘1’ and the positive values of . This value does not imply that there are only longterm correlations in fluctuations, but it implies that nonlinearities due to fattailed distributions are very weak. As a consequence, a timeindependent structure with h(q) ≈ 0.5 is shown in the shuffled timeseries (see Figs. A.1, A.2 and A.3). After shuffling, all timeseries 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 phaserandomised 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 nonlinearities 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 α_{0}, H, Δα and A following the behaviour of their calculated values from the surrogate series.
Firstly, the relation between the values of α_{0} in the three timeseries show a clear difference among them. It is interesting to note that SAP data have larger values than the PDC timeseries, 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 phaserandomised 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 timeseries of RTS indicates such data may promote some values of the fluctuations, making the signal structure richer. Nevertheless, there is an important detail. As 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 timeseries 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 rightskewed (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 timeseries, 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 α_{0} 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 timeseries is slightly more positive. Some differences between both timeseries 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 timeseries, the source of multifractality due to the contribution of longrange 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 rescaled range analysis, the socalled 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 lowamplitude periodicities. The combination of these contributions should make the simulated timeseries closer to the observed ones.
Fig. B.1 Loglog plot of fluctuation functions F_{q}(n) vs. n for a sinusoidal timeseries (shown in the subplot) with a period of ~6.3 days. 
Appendix C Spot modelling
To map the photosphere of Kepler30, 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, solarlike 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 (SanchisOjeda 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: (C.1)
where a_{p}, b_{p}, and c_{p} are derived from the limbdarkening coefficients given by SanchisOjeda 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 timeseries, 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 Kepler30 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 FisherSnedecor 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 threespot 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 Kepler30, 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 Kepler30 (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 Kepler30c 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^{−6} 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 Kepler30. 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).
Parameters adopted in the spot modelling of Kepler30.
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 Kepler30. 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 χ^{2} between theobserved and the modelled fluxes, we minimise a linear combination of the χ^{2} and the regularising functional S, as: (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 χ^{2} and residuals that are symmetrically distributed around a zero mean μ_{res}. When we increase λ_{ME}, we increase the χ^{2} 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 Kepler30 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 χ^{2} 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 σ_{0} of the residuals of the unregularised best fit and the number of data points M in each interval Δt_{f} as . 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 timeseries 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
 Agarwal, S., Del Sordo, F., & Wettlaufer, J. S. 2017, AJ, 153, 12 [Google Scholar]
 Ballot, J., García, R. A., & Lambert, P. 2006, MNRAS, 369, 1281 [Google Scholar]
 Basri, G., & Shah, R. 2020, ApJ, 901, 14 [Google Scholar]
 Belete, A. B., Bravo, J. P., Canto Martins, B. L., et al. 2018, MNRAS, 478, 3976 [Google Scholar]
 Bonomo, A. S., & Lanza, A. F. 2012, A&A, 547, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Collier Cameron A. 2018, Handbook of Exoplanets, eds. H. Deeg, & J. A. Belmonte, (Switzerland: Springer Nature), 3, 1791 [Google Scholar]
 de Franciscis, S., PascualGranado, J., Suárez, J. C., et al. 2019, MNRAS, 487, 4457 [Google Scholar]
 de Freitas, D. B., & De Medeiros, J. R. 2009, Europhys. Lett., 88, 19001 [Google Scholar]
 de Freitas, D. B., Leão, I. C., Lopes, C. E. F., et al. 2013a, ApJ, 773, L18 [Google Scholar]
 de Freitas, D. B., França, G. S., Scherrer, T. M., Vilar, C. S., & Silva, R. 2013b, Europhys. Lett., 102, 39001 [Google Scholar]
 de Freitas, D. B., Nepomuceno, M. M. F., de Moraes Junior, P. R. V., et al. 2016, ApJ, 831, 87 [Google Scholar]
 de Freitas, D. B., Nepomuceno, M. M. F., de Moraes Junior, P. R. V., et al. 2017, ApJ, 843, 103 [Google Scholar]
 de Freitas, D. B., Nepomuceno, M. M. F., Alves Rios, L. D., Das Chagas, M. L., & De Medeiros, J. R., 2019a, ApJ, 880, 151 [Google Scholar]
 de Freitas, D. B., Nepomuceno, M. M. F., Cordeiro, J. G., Das Chagas, M. L., & De Medeiros, J. R., 2019b, MNRAS, 488, 3274 [Google Scholar]
 De Medeiros, J. R., Lopes, C. E. F., Leão, I. C., et al. 2013, A&A, 555, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabrycky, D. C., Ford, E. B., Steffen, J. H., et al. 2012, ApJ, 750, 114 [Google Scholar]
 Gieratowski, J., Zebrowski, J., & Baranowski, R. 2012, Phys. Rev. E, 85, 021911 [Google Scholar]
 Gilliland, R. L., Chaplin, W. J., Jenkins, J. M., Ramsey, L. W., & Smith, J. C. 2015, AJ, 150, 133 [Google Scholar]
 Gu, G.F., & Zhou, W.X. 2010, Phys. Rev. E, 82, 011136 [Google Scholar]
 Henry, B., Lovell, N., & Camacho, F. 2001 Nonlinear dynamics time series analysis, ed. A., Citeseer, 5, 10 [Google Scholar]
 Hippke, M., David, T. J., Mulders, G. D., & Heller, R., 2019, AJ, 158, 143 [Google Scholar]
 Hurst, H. E. 1951, Trans. Am. Soc. Civ. Eng., 116, 770 [Google Scholar]
 Ihlen, E. A. F. 2012, Front. Physiol. 3, 141 [Google Scholar]
 Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al. 2010a, ApJ, 713, L87 [NASA ADS] [CrossRef] [Google Scholar]
 Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al. 2010b, ApJ, 713, L120 [NASA ADS] [CrossRef] [Google Scholar]
 Kantelhardt, J. W., 2015, Fractal and multifractal time series, in Encyclopedia of Complexity and Systems Science (New York: Springer Science & Business Media) [Google Scholar]
 Kawaler, S. D., 1988, ApJ, 333, 236 [Google Scholar]
 Kinemuchi, K., Barclay, T., Fanelli, M., et al. 2012, PASP, 124, 963 [Google Scholar]
 Kraft, R. P., 1967, ApJ, 150, 551 [Google Scholar]
 Krzyszczak, J., Baranowski, P., Zubik, M., et al. 2018, Theor. Appl. Climatol., 137, 1811 [Google Scholar]
 Lanza, A. F. 2012, A&A, 544, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F. 2013, A&A, 557, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F. 2016, Lecture Notes in Physics, 914, 43, Berlin Springer Verlag; https://doi.org/10.1007/9783319241517_3 [Google Scholar]
 Lanza, A. F., Rodono, M., Pagano, I., Barge, P., & Llebaria, A. 2003, A&A, 403, 1135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Bonomo, A. S., & Rodonò, M. 2007, A&A, 464, 741 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Pagano, I., Leto, G., et al. 2009, A&A, 493, 193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Das Chagas, M. L. & De Medeiros, J. R. 2014, A&A, 564, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lanza, A. F., Netto, Y., Bonomo, A. S., et al. 2019, A&A, 626, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
 Luo, A. C. J., & Han, R. P. S., 1992, Chaos Solitons Fractals, 2, 335 [Google Scholar]
 McQuillan, A., Aigrain, S., & Mazeh, T. 2013, MNRAS, 432, 1203 [Google Scholar]
 Makowiec, D., & Fuliński, A. 2010, Acta Phys. Pol. B, 41, 1025 [Google Scholar]
 Mathur, S., Salabert, D., Garcia, R. A., & Ceillier, T. 2014, J. Space Weather Space Clim., 4, 15 [Google Scholar]
 Maxted, P. F. L. 2016, A&A, 591, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meunier, N., Lagrange, A.M., Borgniet, S., et al. 2015, A&A, 583, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Panichi, F., Gozdziewski, K., Migaszewski, C., Szuszkiewicz, E., 2018, MNRAS, 478, 2480 [Google Scholar]
 Radick, R. R., Lockwood, G. W., Henry, G. W., et al. 2018, ApJ, 855, 75 [Google Scholar]
 SanchisOjeda, R., Fabrycky, D. C., Winn, J. N., et al. 2012, Nature, 487, 449 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
 Skumanich, A., 1972, ApJ, 171, 565 [NASA ADS] [CrossRef] [Google Scholar]
 SilvaValio, A., & Lanza, A. F. 2011, A&A, 529, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
 Strassmeier, K. G. 2009, A&ARv, 17, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
 Suyal, V., Prasad, A., & Singh, H. P. 2009, Sol. Phys., 260, 441 [Google Scholar]
 Tang, L., Lv, H., Yang, F., & Yu, L. 2015, Chaos, Solitons Fractals, 81, 117 [Google Scholar]
 Tanna, H. J., Pathak, K. N., 2014, Astrophys. Space Sci., 350, 47 [Google Scholar]
 Thompson, S. E., Christiansen, J. L., Jenkins, J. M., et al. 2013, Kepler Data Release 23 Notes (KSCI19063001) [Google Scholar]
 Van Cleve, J. E., & Caldwell, D. A. 2009, Kepler Instrument Handbook, KSCI19033 [Google Scholar]
 Van Cleve, J. E., Jenkins, J. M., Caldwell, D. A., et al. 2010, Kepler Data Release 5 Notes, KSCI19045001 [Google Scholar]
 Walkowicz, L. M., & Basri, G. S., 2013, MNRAS, 436, 1883 [Google Scholar]
 Wang, J., Shang, P., & Cui, X., Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2014, 89, 032916 [Google Scholar]
 Zaqarashvili, T. V., Albekioni, M., Ballester, J. L., et al. 2021, Space Sci. Rev., 217, 15 [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
See https://in.mathworks.com/help/econ/crosscorr.html, for more details.
MATLAB codes for MFDMA analysis can be found in the arXiv version of the Gu & Zhou (2010) paper: https://arxiv.org/pdf/1005.0877v2.pdf
It is possible to adopt different kinds of moving averages instead of that defined by Eq. (A.2). For example, instead of considering data points that precede the given time t, it is possible to consider points that follow it or that are taken half before and half after the time t (see Gu & Zhou 2010).
All Tables
Multifractal indicators extracted from the standard MFDMA method for PDC, SAP, and RTS.
All Figures
Fig. 1 Kepler30 PDC timeseries. Top panel: PDC timeseries (black dots) and polynomial thirdorder fits. Bottom panel: final PDC timeseries adjusted by polynomial functions. The signatures of three planets are maintained. 

In the text 
Fig. 2 As in Fig. 1 but for SAP timeseries. 

In the text 
Fig. 3 Detrending of Kepler30 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 timeseries. Black dots indicate the PDC pipeline of Kepler30 after the first three steps of our procedure. Detrended flux is a combination of noise background and planet transits, indicating the trending of the timeseries 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) timeseries adjusted according to the steps of Sect. 3. The corresponding Lomb–Scargle periodograms are shown at the bottom. The dasheddotted 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 timeseries has a flat periodogram for periods longer than 20 days confirming that the PDC pipeline tends to remove all the longterm trends. 

In the text 
Fig. 6 Difference between the SAP and PDC timeseries, 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 Timedelayed PDC timeseries 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 timeseries. 

In the text 
Fig. 8 Upper panels: PDC light curve of Kepler30 (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 Kepler30. 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 Kepler30 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 Loglog plot of the fluctuation functions F_{q}(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 timeseries. 

In the text 
Fig. 15 As in Fig. 13 but for final RTS timeseries. 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 Kepler30. 

In the text 
Fig. 17 Hurst surface h(q, τ) calculated for final PDC timeseries. 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 walktype 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 phaserandomised (right side) versions of timeseries: 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 timeseries following steps 5 and 6 presented in Sect. 4. Top panel: original (in red), shuffled (green), and phaserandomised (blue) data are based on the procedure mentioned in Sect. 4. Left top: multifractal fluctuation function F_{q} (n) obtained from MFDMA method for only q = 2, indicated henceforth as a big circle. Right top: qorder Hurst exponent (h(q)) as a functionof qparameter. 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 timeseries, respectively. 

In the text 
Fig. A.2 As in Fig. A.1 but for SAP timeseries. 

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 Loglog plot of fluctuation functions F_{q}(n) vs. n for a sinusoidal timeseries (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 Kepler30. 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 Kepler30 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.