Timescales of starspot variability in slow rotators

There is an intriguing proximity between the turnover time τ MLT of the standard mixing length theory of the Sun and the timescale τ lam of solar activity patterns at the space scale of giant laminar convection assumed in deep layers of the Sun. To verify the reliability of this correspondence, we analyzed the light curves of 637 slowly rotating stars, observed by the Kepler mission, with periods from 16 to 30 days. The proximity τ MLT ≈ τ lam is conﬁrmed. The performed study also conﬁrms the manifestation of large scale turbulence in the dynamics of surface activity such as that in the Sun. These results open a new way to measure the key astrophysical parameter τ MLT and to study deep convection that has been undetected with asteroseismology.


Introduction
The mixing of plasma in deep layers of the solar convection zone (hereafter deep-mixing) plays a key role in calibration of the standard mixing length theory (MLT), which is a standard in modern astrophysics. In particular, the convective turnover time τ MLT at the bottom of convection zone plays an important role in stellar modeling and rotation-activity relations. At the same time, this parameter cannot be measured directly. Only the relative variations of the convective turnover time with stellar color index or mass were found for main sequence stars using the rotation-activity relation. The absolute measure of the parameter requires a calibration with the theoretically calculated value of τ MLT using the solar MLT model (e.g., Noyes et al. 1984;Wright et al. 2011).
However, the appropriateness of the MLT approach, based on the abstraction of identical convection cells (Böhm-Vitense 1958), is not obvious in a turbulent cascade of stellar convection (e.g., Li & Yang 2007). Nevertheless, the outcomes of such an approach look reasonable in the context of observational studies and are widely accepted in stellar physics. This paradox suggests a certain physical reality in the background of the MLT model with the supposed identical convective cells. In this paper we focus on an interesting formal proximity between τ MLT and the timescale τ lam of solar surface activity-pattern at the space scale of laminar convection (Arkhypov et al. 2013). 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 et al. 2013). However, the solar proximity τ MLT ≈ τ lam needs extensive verification for other stars with various τ MLT . The main sequence stars with rotation periods P < 16 days in the Kepler photometric survey are known to show manifestations of deep mixing (Arkhypov et al. 2015a(Arkhypov et al. , 2016.
The subject of the present report is the verification of our previous results using another stellar set, consisting of rotators with the solar-type periods 16 < P < 30 days. Such slow rotators seem best indicators of deep mixing effects because of their relatively slight magnetism and, hence, the shortest lifetimes of starspots, which appear as a masking factor in this work.
The presented analysis became possible owing to an improved method explained below in Sect. 2. The obtained results are described in Sect. 3 and concluded in Sect. 4.

Method and stellar set
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 only the key ideas. We analyzed the rotational modulation of the stellar radiation flux F, which reflects the longitudinal distribution of spots (Fig. 1). The data used were already preprocessed to correct for the instrumental and environmental effects (PDCSAP_FLUX from the Kepler mission archive 1 ). As in our previous studies (Arkhypov et al. 2015a(Arkhypov et al. ,b, 2016, we considered the amplitudes A 1 and A 2 of the light-curve rotational Fourier harmonics with periods P and P/2, respectively, where P is the period of stellar rotation. In Arkhypov et al. (2016) it was shown that the squared amplitudes A 2 1 and A 2 2 of these harmonics are statistically proportional to the sunspot number. In view of that, we used the values A 2 1 and A 2 2 as activity indexes and found the timescales of their stochastic 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 chronological series of the corresponding activity indexes (i.e., A 2 1 or A 2 2 , respectively) at the time lag of one rotational period P. This approach is based on a simple 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.
In fact, this approach is an adaptation of our method applied to the analysis of solar synoptic charts (Arkhypov et al. 2013). We used the same preprocessing of light curves as in our previous study in Arkhypov et al. (2016), aside from the improved later procedure for removing flares. In this work we approximated the one-period light curve, which covers a time interval of one rotation period P, with a sixth-order polynomial F p = a p δt 6 + b p δt 5 + c p δt 4 + d p δt 3 + e p δt 2 + f p δt + g p , where δt is the light curve running time and a p , b p , c p , d p , e p , f p , and g p are the fitted coefficients. This polynomial is used for the replacement of the measured flux F = F p if |F − F p | > 2σ p , where σ p = (F − F p ) 2 1/2 is the standard deviation in the considered one-period time interval.
To exclude the artificial (e.g., instrumental) drifts and long period stellar nonrotational variability (i.e., with timescale P), we removed the linear trend in the one-period light curve using the following improved (as compared the method introduced in Arkhypov et al. 2015a) method. In the new algorithm we linearly approximated the first 0.05 P (i.e., the beginning) and last 0.05 P (i.e., the end) parts of the one-period light curve to be F = a b + b b δt and F = a e + b e δt, respectively. Here a b , b b , a e , and b e are fitted coefficients. Using these approximations, one can estimate the corresponding boundary fluxes F b = a b and F e = a e + b e P for δt = 0 and P, respectively, and calculate the drift correction, i.e., F − [(F e − F b )/P](δt − P/2). As result, we obtained the light curve F(δt) with F(0) ≈ F(P), which we further analyzed to study its rotational Fourier harmonics. We note that this null trend rule is justified statistically. Indeed, although in reality the starspot evolution could generate slight accidental differences in F(0) and F(P), the average difference F(0) − F(P) should be negligible. Therefore, the main sense of the trend removing is a prevention of the residual instrumental trends. Figure 2 shows the result of application of the abovementioned procedures to a one-period light curve of the star KIC 5435923, including the flux interpolation in short gaps. Such prepared one-rotation light curves are further analyzed to study their rotational Fourier harmonics only.
For the study we selected the light curves of 637 slow rotators. This set includes the stars with rotation periods 16 < P < 30 days (according to measurements in Nielsen et al. 2013 andMcQuillan et al. 2014) and effective temperatures 3337 ≤ T eff ≤ 7117 K from the Mikulski Archive for Space Telescopes (MAST) 2 ). All the selected stars belong to the main sequence (surface gravity log(g[cm s −2 ]) > 4 in the catalogs used). These stars have high quality light curves without interferences (i.e., with no detectable short period pulsations or double periodicity from companions). The details of the star selection procedure are described in Arkhypov et al. (2015bArkhypov et al. ( , 2016. Particularly, the value P, measured mainly in McQuillan et al. (2014), was visually reverified for the light curve of each selected star. In fact, this check does not reveal any misestimated periods by a factor 2. This result confirms our previous conclusion (Arkhypov et al. 2015b) on the reliability of the P catalogs (i.e., Nielsen et al. 2013 andMcQuillan et al. 2014).

Results
According to Arkhypov et al. (2015aArkhypov et al. ( , 2016, a kind of a gradientlike function is a useful diagnostic tool for starspot variability, that is, For example, the Kolmogorov theory of turbulence (e. g., Lang 1974) predicts a universal relation between the characteristic size of turbulent eddies (L) and the timescale of their variability (τ L ), i.e., τ L ∝ L 2/3 . Taking into account that L ∝ m −1 , we obtain, for the timescale of rotational harmonics, τ m ∝ m −2/3 , where m is the harmonic number. Substitution of this proportionality in Equation (1) gives β 12 = −2/3 = 0.6(6). Figure 3 shows β 12 estimate histograms for individual stars from the considered set of slow rotating objects. One can see the proximity of the histogram maxima to the predicted values -2/3. In particular, Fig. 3b is constructed with the bin width 0.09 localizing the histogram maximum at β = −0.68 ± 0.045, i.e., with 7% precision. The significance of this main peak is verified considering the binomial distribution of the near-maximum estimates in the limited interval −1.8 < β 12 < 0.4 (Fig. 3c). We calculated the probability of the main histogram peak, or a higher one, supposing a random distribution of β 12 estimates Fig. 2. Example of one-period light curve for the star KIC 5435923 before (a) and after (b) the flare removing procedure. In both curves, the trend is also removed and the flux is normalized F/F o , where F o is the average flux for the considered time interval P. The interpolated gap is seen at 8.0 < δt < 8.8 days. where n = 26 is the maximal number of β 12 estimates in one bin of the histogram,N = 367 is the total number of the considered β 12 estimates, and q = 1/25 is the probability for one estimate to be in a certain bin from the set of 25 bins. Since the obtained probability W ≥n is one chance of 1/0.0075 ≈ 133, the considered histogram maximum at β 12 ≈ −0.67 is apparently not a statistical fluctuation, but is an indication of the turbulence manifestation. This is consistent with the analogous results for the stars with 4 < P < 16 days (Arkhypov et al. 2016) and an outcome of the corresponding study for the Sun in Arkhypov et al. (2013) for the higher order harmonics 1 m 20. The turbulent convection generates irregular subphotospheric flows at the local height scale (m >> 1). 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, Fig. 3 provides an argument for the turbulence manifestation (β 12 ≈ −0.67) at m = 1 and 2. Hence, the global (m = 1 and 2) character of the revealed turbulence contradicts to the small-scale photospheric and subphotospheric turbulent vortices. However, the deep large-scale convection is predicted in numerical modeling (e.g., Miesch et al. 2008) and found experimentally as a regular drift of supergranules (Hathaway et al. 2013). It is natural to connect the global turbulent effect in the starspot pattern with the giant convection in deep layers of a star.
This probable connection could be additionally verified by a comparison between the found timescales τ 1 , τ 2 and the theoretical timescale of deep mixing τ MLT of the MLT. Figure 4 shows that these values are different. But there is certain proximity between τ MLT and the reductions of τ 1 and τ 2 to the spatial scale D lam ≈ 2 3/2 H cz (Dibaj & Kaplan 1976) of laminar convection based on the confirmed values for the deep turbulence time scaling, τ m ∝ m −2/3 (see Fig. 3), that is, where m lam = 2πR star /D lam ; R star is the stellar radius from MAST and H cz is the convective zone depth (van Saders & Pinsonneault 2012). This proximity is shown in Fig. 5.

Discussion and conclusions
Although the turbulent convection generates irregular subphotospheric flows smaller than supergranules at m 1, nevertheless the manifestation of turbulence (β 12 ≈ −2/3) at global scales (m = 1 and 2) is confirmed also in stars with the solar-type rotational periods 16 < P < 30 days (Sect. 3). This means that the emersion of magnetic tubes in such stars including solarlike objects could be modulated by global turbulent vortices, which have vanishingly small direct manifestation in the solar photosphere. At the same time such global motions are predicted in deep layers of the convective zone (e.g., Miesch et al. 2008). Alternative interpretations of the starspot pattern dynamics in terms of magnetic dynamo or starspot decaying predict other values of the parameter β 12 , which are not supported by the histograms in Fig. 3. In particular: Fig. 4. Timescales of variability of A 2 1 and A 2 2 , averaged in temperature bins (a): log(τ 1 ) and (b): log(τ 2 ) indicated by X symbols with error bars for the stars with 16 < P < 30 days vs. effective temperature T eff . For comparison, the estimates of the τ MLT turnover time by Noyes et al. (1984), are depicted by a dashed line and Wright et al. (2011) as error bar crosses. (c) In highly magnetized subphotospheric plasma, such as porous media, or in a network flow, the dependence of an average squared displacement of magnetic element x 2 on time t can differ from the linear relation. In a general case it may be represented as x 2 ∝ t α , where α < 1 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 timescale τ m . In this case τ m ∝ m −2/α and, according Eq.
(1), β 12 = −2/α. For example, the subdiffusion of magnetic elements with α = 0.6 ± 0.2 was found in the solar photoshere 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). (d) 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 timescale of the considered feature blurring (i.e., A 2 m damping) is τ m ∝ m −1 , which corresponds to β 12 = −1 in Eq. (1). However, the corresponding peak at n = 19 in Fig. 3c is rather insignificant because its W ≥n = 0.22 according to Eq. (2). The aforementioned alternative predictions (i.e., β 12 ≈ 0, −2, < −2, and −1, respectively) significantly (no less than than a half-width of the histogram bin) deviate from the histogram maxima at the prediction β 12 ≈ −2/3 of the Kolmogorov's turbulence theory. This fact is an argument for the deep-convective modulation of starspot pattern.
Such modulation is observable if its timescale is longer than the typical lifetime of starspots. If it were otherwise, the starspot pattern dynamics are controlled by the long lifetime of spots, which smooths the shorter living effects from deep mixing. In this case, the maximum of β 12 -histogram should correspond to the spot decay process. Detailed studies (Arkhypov et al. 2016(Arkhypov et al. , 2018 have confirmed that stars with saturated magnetism have the histogram maximum at β 12 ≈ −2, corresponding to the diffusive decay of long-living spots. However, the slow rotators considered in this study have nonsaturated magnetism (see Fig. 2a in Arkhypov et al. 2018), sufficiently short-living spots, and the most frequent β 12 = −2/3, corresponding to the giant turbulence of deep-mixing.
Another possible manifestation of deep mixing appears to be the confirmed proximity of our timescale τ lam for laminar giant convective cells to the turnover time τ MLT (Fig. 5). In the successful MLT this turnover time is attributed to the standard (i.e., laminar) convective cells at the bottom of the convection zone (e.g., Gilliland 1985). The decaying giant (∼45 o ) cells of laminar convection in turbulence were modeled in the solar convection zone (e.g., Brun & Toomre 2002) and indirectly detected in solar activity patterns (Arkhypov et al. 2012(Arkhypov et al. , 2013 and in the regular drift of supergranules (Hathaway et al. 2013). However, the directly observed solar photospheric granulation and even subphotospheric supergranules have much smaller diameters ( 0.4 o ) and shorter lifetimes ( 1.8 day; Rieutord & Rincon 2010) than the scales τ MLT and τ lam1,2 (>10 days) in Fig. 5. Hence, the obtained scales τ lam1,2 ≈ τ MLT support the hypothesis regarding the starspot manifestation of plasma mixing in deeper layers of convective zones.
Since we confirmed the deep-mixing manifestations in the solar-like stars, such phenomena are additionally justified in the Sun as well. In particular, the found proximity β 12 ≈ −2/3 suggests that some deep global flows generate the turbulent cascade at maximal scales. This raises a question on the source of such flows.
The standard MLT describes mixing in a range of situations, from planetary interiors to stars and the interstellar medium. However, its applicability to turbulence in astrophysical problems is often questionable. Nevertheless, MLT is acceptable for turbulence in stars. The confirmed proximity τ MLT ≈ τ lam1,2 suggests a way to resolve this paradox in terms of laminar contribution in the deep mixing.
The confirmed proximity τ MLT ≈ τ lam1,2 opens a new way to measure the key astrophysical parameter τ MLT , which has so far only been calculated theoretically or estimated semiempirically using the MLT paradigm. Correspondingly, our results seem promising for related fields such as stellar activity and evolution.