Time-scales of stellar rotational variability and starspot diagnostics

The difference in stability of starspot distribution on the global and hemispherical scales is studied in the rotational spot variability of 1998 main-sequence stars observed by Kepler mission. It is found that the largest patterns are much more stable than smaller ones for cool, slow rotators, whereas the difference is less pronounced for hotter stars and/or faster rotators. This distinction is interpreted in terms of two mechanisms: (1) the diffusive decay of long-living spots in activity complexes of stars with saturated magnetic dynamos, and (2) the spot emergence, which is modulated by gigantic turbulent ﬂows in convection zones of stars with a weaker magnetism. This opens a way for investigation of stellar deep convection, which is yet inaccessible for asteroseismology. Moreover, a subdiffusion in stellar photospheres was revealed from observations for the ﬁrst time. A diagnostic diagram was proposed that allows differentiation and selection of stars for more detailed studies of these phenomena.


I N T RO D U C T I O N
Spots in stellar photospheres are extensively studied as an indicator of physical processes in stars and as a proxy of sunspot phenomenology. Nevertheless, particular features in evolution of starspot patterns and their physical drivers in stars of various types are insufficiently explored. For example, the measurements of spot lifetime lead to unexpectedly high estimates of magnetic diffusivity in some stars, supposing an anomalous magnetic diffusion related to the famous giant convection cells (Bradshaw & Hartigan 2014). Besides magnetic cycles, the shorter time-scales of starspot variability traditionally were interpreted in terms of diffusive decay and spot disruption by the shearing of stellar differential rotation (e.g. Bradshaw & Hartigan 2014 and therein). However, the distinction between these effects rarely was a subject of practical study (Hall & Henry 1994). Mainly both the mechanisms were considered in stellar modelling without distinction and clear diagnostics (e.g. Isik, Schüssler & Solanki 2007 and therein). Additionally, the spot emergence frequency was recognized as the third factor affecting the lifetime of starspot pattern.
The importance of the emergence modulation for the evolution of spot pattern was confirmed in our solar and stellar studies (Arkhypov, Antonov & Khodachenko 2013;Arkhypov et al. 2015aArkhypov et al. ,b 2016. In particular, a diagnostic method has been pro-E-mail: oleksiy.arkhypov@oeaw.ac.at (OVA); maxim.khodachenko@ oeaw.ac.at (MLK) posed there to identify the manifestations of the aforementioned mechanisms of starspot variability. In this Letter, we describe several important express-results of application of this method to an extended set of main-sequence stars, which opens for research community the new ways to probe the stellar physics.

M E T H O D A N D DATA S E T
Detailed descriptions of the used spectral-autocorrelation method can be found in our previous papers (Arkhypov et al. 2015a(Arkhypov et al. ,b, 2016. Therefore, we briefly describe here only its key ideas. We analyse the rotational modulation of the stellar radiation flux F (PDCSAP_FLUX from the Kepler mission archive 1 ), which reflects the longitudinal distribution of spots.
Like in our previous studies (Arkhypov et al. 2015a(Arkhypov et al. ,b, 2016, we consider here the squared amplitudes A 2 1 and A 2 2 of the lightcurve rotational harmonics with periods P and P/2, respectively, where P is the period of stellar rotation. Following this approach, we found the time-scales of their variability τ 1 = −P/ln [r 1 (P)] and τ 2 = −P/ln [r 2 (P)]. Here, r 1 (P) and r 2 (P) are the autocorrelation functions of the chronological series of the A 2 1 or A 2 2 , respectively, at the time lag of one rotational period P. This method is based on the simplest approximation of the logarithm of autocorrelation function at the shortest lag: ln (r m ) ≈ − t/τ m (for details, see in  Arkhypov et al. 2016), where m = 1 or 2 is the harmonic number, and t = P is the lag. The squared amplitudes of harmonics are used in the analysis because of their statistical proportionality to the solar spot number (Arkhypov et al. 2016). Analogously, Messina et al. (2003) used as an activity index the maximum amplitude (A max ) of rotational variations of stellar light curve, which is a proxy of our amplitude A 1 of the fundamental harmonic. They found that their index of various main-sequence stars (0.1 < P < 20 d; M4V to F5V) is related to normalized X-ray luminosity: L x /L ∝ A b max , where b ≈ 2, and L is the stellar bolometric luminosity. The ratio L x /L is widely used as a standard activity index (e.g. Wright et al. 2011), related to spots (Wagner 1988;Ramesh & Rohini 2008). Since A 2 1 and A 2 2 are naturally the major contributors to A 2 max , the approximate relation L x /L ∝ A 2 max supports applicability of A 2 1 and A 2 2 as stellar activity indexes for the non-solar-type stars too. Moreover, we found that our activity index A 2 1 is related to the Rossby number similarly to the X-ray ratio R x ≡ L x /L (see table 3 in Arkhypov et al. 2016), supporting the general statistical proportionality L x /L ∝ A 2 1 . Note that the measurements of time-scales τ 1 and τ 2 are insensitive to this proportionality. For example, the variation cycle of generalized index A γ m with harmonic number m = 1, 2, and the constant γ has the same duration (i.e. the same time-scale) at any γ .
According to Arkhypov et al. (2015aArkhypov et al. ( , 2016, the effects of starspot variability can be distinguished using the 'gradient' ratio In fact, this ratio measures the extent to which the spot distribution evolves faster on smaller scale compared to the largest scales. For example, the Kolmogorov's theory of turbulence (e.g. Lang 1974) predicts a universal relationship between the characteristic size of turbulent eddies (L) and the time-scale of their variability (τ L ): τ L ∝ L 2/3 or τ m ∝ m −2/3 (taking into account that L ∝ m −1 ). Substituting this proportionality into equation (1), one can obtain β 12 = −2/3 −0.67.
The horizontal diffusion of magnetic elements in photosphere effectively decreases our activity indexes A 2 m (m = 1, 2), when the average displacement of magnetic elements is about the longitudinal period of harmonic 2π/m. For the normal diffusion, this happens during the characteristic time τ m ≈ (2π/m) 2 η −1 (where η is the diffusion coefficient), i.e. τ m ∝ m −2 . In this case, equation (1) gives β 12 = −2.
In highly magnetized sub-photospheric plasma, like porous media, or in a network flow, the dependence of a squared displacement of magnetic element x 2 on time t, can differ from the linear one. In general case, it may be represented as x 2 ∝ t α , where α is a constant. Therefore, the squared displacement x 2 ≈ (2π/m) 2 = ητ α m corresponds to the noticeable decrease of the activity index A 2 m during the harmonic time-scale τ m . In this case, τ m ∝ m −2/α and, according equation (1), β 12 = −2/α. For example, the subdiffusion of magnetic elements with α = 0.6 ± 0.2 was found in the solar photosphere at the spatial scale of supergranulation (∼10 4 km), which is related to the material network flow that traps these elements in the conjunction points (Iida 2016).
The differential rotation of a star stretches an activity region, or a complex of active regions, over the longitudinal harmonic scale m = 2π/m during the time τ m = m / , where is a typical variation of angular velocity over the latitudinal extension of the activity region or complex. Hence, the time-scale of the considered feature blurring (i.e. A 2 m damping) is τ m ∝ m −1 that corresponds to β 12 = −1 in equation (1).
A stellar activity cycle modulates the total spot number with the identical period or time-scale for all rotational harmonics of a light curve. This means β 12 = 0.
The aforementioned predictions are tested below by measuring the parameter β 12 for the main-sequence stars. For this purpose, we combined our previously analysed data set (Arkhypov et al. 2016) that contains the light curves of 1361 main-sequence stars observed by the Kepler space observatory, with the light curves of additionally selected 637 slow rotators. In summary, the extended sample includes 1998 stars with 0.5 < P < 30 d (according to measurements in Nielsen et al. 2013;McQuillan, Mazeh & Aigrain 2014) and the effective temperatures 3227 ≤ T eff ≤ 7171 K (according to Huber et al. (2014) in the Mikulski Archive for Space Telescopes 2 ). All selected stars belong to the main sequence (surface gravity log (g[cm s −2 ]) > 4 in the used catalogues). The availability of a high-quality light curve without any interferences (i.e. no detectable short period pulsations or double periodicity from companions) was a special criterion for the compiling of analysed sample of stars. Further details on the star selection and light curve preparing (i.e. removing of gaps, flares, artefacts, trend) and processing are described in Arkhypov et al. (2015bArkhypov et al. ( , 2016.

D I AG N O S T I C D I AG R A M
After processing of the light curves, we have measured the parameter β 12 for the selected stars. Fig. 1(a) shows the P − T eff distribution of smoothed values β 12 as a result of β 12 averaging over individual stars with log (P c ) − 0.2 < log (P) < log (P c ) + 0.2 and log (T c ) − 0.05 < log (T eff ) < log (T c ) + 0.05 in a sliding window with the central values of stellar rotation period P c and effective temperature T c . In average, 117 (up to 400) individual estimates of β 12 appeared within this sliding window. The average standard error of β 12 is σ = 13 per cent. The errors σ = [ (σ − β 12 ) 2 ] 1/2 in the individual windows of averaging does not exceed 15 per cent in the majority (83 per cent) of cases for β 12 . The total distribution of σ / β 12 is shown in Fig. 1(c). One can see that the found pattern of β 12 is mainly reliable, excluding the corners of the diagram with a depleted stellar population there (Fig. 1d) and increased relative errors (Fig. 1c). Fig. 1(a) reveals an interesting pattern, consisting of the two regions with an increased (dark colour) and decreased (brighter tint) value of β 12 . Fig. 1(b) depicts these regions in terms of predictions for different mechanisms of starspot variability. The used intervals for schematics of β 12 correspond to ±σ . Since the difference between β 12 ≈ −2/3 in the 'turbulent' region (black colour in Fig. 1b) and the 'sub-diffusion' region with β 12 < −2.3 (grey) is more than 12.5 σ , these diagram details are significant. Fig. 2(a) demonstrates that the discussed regions in the diagnostic diagram are separated approximately along the lines, corresponding to the Rossby number Ro = P/τ MLT = 0.13, which was found as a border between the saturated and unsaturated stellar magnetism (Wright et al. 2011). To construct these border lines, we used two versions of the turnover time τ MLT in the mixing length theory (MLT): (a) the blue line is based on the classical definition τ MLT (B − V) by Noyes et al. (1984), where the standard colour index B − V is related to T eff (Flower 1996); (b) the red line corresponds to the commonly used τ MLT (m * ) by Wright et al. (2011), where m * is the stellar mass in solar unit related to T eff (Habets & Heintze 1981). Fig. 2(b) illustrates a similar pattern as in Fig. 2(a), but for an analogously smoothed with the sliding-window time-scale log (τ 1 ) . Its standard error σ τ = [log (τ 1 ) − log (τ 1 ) ] 2 1/2 is shown in Fig. 2(c). While the error, averaged over the whole plot, is σ τ = 0.03, the value log (τ 1 ) in sliding windows varies from 1.28 to 2.88, corresponding to 19 ≤ τ 1 ≤ 765 d in the most reliable region of the diagram at log (P) > 0. Since this variation of log (τ 1 ) is about 53 σ τ , the found pattern log (τ 1 ) in Fig. 2(b) is significant. The time-scale maximum clearly corresponds to the region with minimal β 12 in Fig. 2(a), or the sub-diffusion region in Fig. 1(b). This τ 1 ↔ β 12 correspondence confirms the reality of stellar sub-diffusion, which prolongs the lifetime of spots due to their suppressed decay.
However, the regions of diffusion (β 12 ≈ −2) as well as the differential rotation (β 12 ≈ −1) in Fig. 1(b) could be artefacts of smoothing in the transition area between the diagram features with values β 12 < −2.3 and β 12 ≈ −2/3. To test this possibility, we consider the distributions of β 12 values, measured for individual stars. Fig. 3(a) shows such distribution for the stars with certainly saturated magnetism at Ro = P/τ MLT < 0.13, where τ MLT is taken in the conservative version by Noyes et al. (1984). One can see the histogram maximum at the predicted diffusion value β 12 = −2, confirming the reality of starspot decay by the normal diffusion process. Fig. 3(b) demonstrates the β 12 -histograms for the stars with formally unsaturated magnetism at Ro > 0.13. The most probable values β 12 ≈ −0.67 correspond to the prediction of the turbulent manifestation (dashed lines). The significance of the main peak is demonstrated using various bins in Figs 3(b) and (c). Considering the binomial distribution of the near-maximum estimates in the limited interval −1.75 < β 12 < 0.4 (Fig. 1c), we calculated the probability of the main histogram peak, or a higher one, supposing a random distribution of β 12 estimates where n = 108 is the maximal number of β 12 estimates in one bin of the histogram; N = 1062 is the total number of the considered β 12 estimates; q = 1/13 is the probability for one estimate to be in a certain bin from the set of 13 bins. Since the obtained probability W ≥n is one chance of 344, the considered histogram maximum at β 12 ≈ −0.67 is rather not a statistical fluctuation, but an indication of the turbulence manifestation. It is worth to note that the absence in Figs 2(b) and (c) of any hint on the peak at β 12 = −1 predicted for the differential rotation. Therefore, the differential rotation effect remains unconfirmed. Apparently, the difference of selected stars in typical light-curve amplitude and measured rotational periods lead to the dispersion of the autocorrelation accuracy, hence, to the scattering of the estimated β 12 in Fig. 3a,b. For example, the more noisy low-amplitude light curves of the hot and slow rotators could give somewhat decreased time-scales τ 1 and τ 2 caused by decrease of the related autocorrelation coefficients r 1 (P) and r 2 (P), respectively. Since generally A 2 < A 1 , this effect could decrease τ 2 more than τ 1 , decreasing therefore the gradient ratio β 12 . However, one can see in Fig. 1(a) Wright et al. 2011, 'blue': τ MLT by Noyes et al. 1984) on the β 12 pattern; (b) the value log (τ 1 ) averaged in the same sliding windows as in (a); (c) the relative error σ τ / log (τ 1 ) of estimates log (τ 1 ) in (a). observed difference in β 12 for hot-slow and cold-fast rotators is not a result of the selection effect. The quasi-symmetric core of the distribution of β 12 estimates in Fig. 3(b) and the correspondence of the histogram peaks to the predicted value β 12 = −2/3 demonstrate insignificance of the selection effect for our conclusions. Since the decrease of β 12 due to the light-curve noise is statistically insignificant in the considered stars with unsaturated magnetism, this effect must be weaker and negligible in the high-amplitude light curves of saturated stars. Hence, the decreased β 12 for cold rotators can be interpreted as real manifestation of the sub-diffusion decay of starspot patterns in stars with saturated magnetism.

D I S C U S S I O N A N D C O N C L U S I O N S
In general, we found that the global spot distribution evolves slower than the hemispherical one. Particularly, the largest scale is much more stable than the smaller one for cold-slow rotators, whereas the difference is less pronounced for hotter stars and/or faster rotators. The corresponding gradient ratio β 12 < 0 appears an applicable indicator for a dominating mechanism of starspot pattern evolution. Thus, the constructed diagnostic diagram in Fig. 1(b) opens the way to differentiate stars with respect to the processes which control the stochastic variability of starspot pattern. For example, it is first found that the sub-diffusion (−3.0 < β 12 < −2.3 in Fig. 1(a), hence, 0.7 < α < 0.9) is a dominating process in the decay of activity complexes in stars with saturated magnetism and certain values T eff and P. Note that the solar sub-diffusion is a result of deceleration of the normal diffusion process in the converging sub-photospheric plasma flows between supergranules at the scale ∼10 4 km (Iida 2016). Analogously, the sub-diffusion at much larger scales up to ∼10 6 km (m = 1 and 2) requires regular mega-flows MNRASL 473, L84-L88 (2018) Downloaded from https://academic.oup.com/mnrasl/article-abstract/473/1/L84/4554392 by guest on 28 July 2018 in the sub-photospheric plasma. Such flows, converging in active regions, are found in the Sun (Hindman, Haber & Toomre 2009).
The turbulent convection generates non-regular sub-photospheric flows at the local height scale (m 1) in the standard MLT approach. For example, in the solar photosphere the turbulent cascade has been described at scales smaller than that of supergranules, i.e. m > 200 (Abramenko et al. 2001;Stenflo 2012). However, Figs 1(b), 3(b) and (c) argue for the turbulence manifestation (β 12 ≈ −0.67) at m = 1 and 2 in solar-type stars (4500 T eff 6900 K and 6 P 30 d). This incommensurability of the photospheric and found global turbulence means that the last is the manifestation of plasma mixing in deep layers of the stellar convection zones. The connection between starspot pattern and the deep convection was predicted in numerical models of magnetic tube emergence which take into account the modulation effect from the convective flows (e.g. Weber, Fan & Miesch 2013). However, the deep-mixing effect is masked when its time-scale becomes shorter than the typical spot lifetime in the stars with saturated magnetism.
Figs 1(a) and (b) do not reveal any signs of dominating manifestation of magnetic cycles. Note, the period (>10 3 d) of a typical activity cycle, which standardly described in terms of MFT (Brandenburg & Subramanian 2005), must be much longer than the time-scales of the aforementioned processes. Our autocorrelation method for the estimation of τ m is focused on the shortest time-scale, i.e. on the activity complex decay or turbulence manifestation.
Therefore, the described diagnostic approach appears as a useful instrument for future studies of starspot phenomenology, magnetic diffusion and flows at stellar photospheres as well as in deep layers of convection zones.