© The Authors 2023

1 Introduction

The stellar rotation period, P_rot, is a fundamental quantity in stellar astrophysics because it is closely linked to the star’s activity level and its age. Skumanich (1972) first demonstrated that the average equatorial rotational velocity and the emission luminosity in the cores of the Ca II H and K lines both decrease with stellar age, t, according P_rot ~ t^1/2. In the years following, it was shown that, on average, young stars rotate faster and are more active, whereas old stars rotate more slowly and are less active (e.g., Noyes et al. 1984). The pioneering work of Skumanich (1972) inspired the notion of age-dating the star using its rotation period. This semi-empirical method, now known as gyrochronology (Barnes 2003, 2007), calibrates the relation between the stellar mass, rotation period, and age. Hence, knowing the stellar rotation period is essential for estimating the stellar age – a fundamental quantity of the star that cannot be measured directly.

Most commonly, stellar rotation periods are measured by observing stellar brightness variations over time and searching for repeatable patterns in long-term (photometric) time series caused by star spots rotating in and out of view. Owing to the Kepler space telescope’s almost uninterrupted photometric observations of ~150 000 main-sequence stars for a period of 4 yr, the rotation periods have been measured for several tens of thousands of stars (McQuillan et al. 2013a,b, 2014; Reinhold et al. 2013; Walkowicz & Basri 2013; Nielsen et al. 2013; do Nascimento et al. 2014; García et al. 2014; Reinhold & Gizon 2015; Ceillier et al. 2016, 2017; Santos et al. 2019, 2021).

Among these studies, one of the largest samples of rotation periods was provided by McQuillan et al. (2014, hereafter McQ14), measuring periodic brightness variations in 34 030 Kepler stars, which is still one of the largest collections of rotation periods today and used in numerous studies covering a wide range of topics, from constraining stellar dynamo theories to understanding the evolution of our Galaxy (see, e.g., van Saders et al. 2019 for one of the most recent examples). Despite this huge number, McQ14 could not unambiguously detect periods in an even larger sample of 99 000 stars.

Recently, Santos et al. (2019, 2021, hereafter S21) reanalyzed the full Kepler archive and significantly increased the number of detected rotation periods to 55232 out of 159442 targets. Both studies (McQ14; S21) identified a decrease of the period detection rate with increasing effective temperature. M dwarfs have detection fractions of 70–80% and K dwarfs around 50%, whereas the detection fraction drops to ~30% or less for F and G dwarfs. This observation can be explained by the fact that the variability pattern changes from cooler to hotter stars. The light curves of M dwarfs show very regular periodicity over many rotation periods, whereas G-type stars exhibit more irregular variability, showing hardly any periodicity over long sections of the observations. It has been suggested that the source of this behavior is the fact that the spot lifetimes of these stars are often shorter than the stellar rotation period, which leads to irregularities in the light curves (Giles et al. 2017; Basri et al. 2022).

As a consequence, the rotation periods of many stars around solar spectral type remain undetected in an automated period search because the rotational periodicity is not stable enough to generate a significant peak in the frequency analysis (Reinhold et al. 2021). This implies that the conclusions drawn from many studies of near-solar rotators (e.g., van Saders et al. 2019; Reinhold et al. 2020; Okamoto et al. 2021) might be strongly biased toward the behavior of more active stars, for which rotation periods could be determined. In particular, the relatively small number of stars with known rotation periods and variabilities similar to that of the Sun conveys a false picture of the Sun being unusually quiet compared to other stars with detected rotation period (Reinhold et al. 2020).

The main goal of our study is to make use of recent developments in understanding stellar brightness variability and utilizing new methods to determine rotation periods of a larger sample of Kepler stars than ever before. This extended sample of stars with determined rotation periods should at least partly remove these biases. Possible applications range from the comparison of observed period distributions in the Kepler field to predictions of Galactic evolution models (van Saders et al. 2019), comparing solar and stellar variabilities (Reinhold et al. 2020), as well as in the search for superflares on solar-like stars (Okamoto et al. 2021; Vasilyev et al. 2022).

To achieve this goal, more rotation period measurements of stars with small, solar-like variabilities are needed. Recently, Shapiro et al. (2020) showed that the correct rotation period of stars with irregular variability can reliably be detected by a novel method that considers the gradient of the (global wavelet) power spectrum (GPS), instead of the power spectrum itself.

The GPS method has been successfully applied to measure the solar rotation period (Amazo-Gómez et al. 2020a). Furthermore, it shows good agreement with the previously reported periods of Kepler stars (Amazo-Gómez etal. 2020b). We emphasize that, in contrast to classical period analysis methods, the GPS method does not require a repeatable spot pattern in the time series but is sensitive to the typical dip durations of spots crossings (see Sect. 3). Consequently, it even works in cases when the magnetic features live shorter than the stellar rotation period such that no reoccurring transits of the same magnetic features are required. Hence, the GPS method is ideally suited to measure rotation periods for stars where classical methods have previously failed to detected reliable periods (Reinhold et al. 2022).

2 Data and sample selection

2.1 Kepler data

In this work, we analyze the long-cadence light curves processed by the latest version of the Kepler pipeline (Data Release 25).

The data¹ are released in “quarters” with lengths of ~90 days, with exceptions for quarters Q0, Q1, and Q17, which have shorter observing times of 10–33 days. In the following, we use all available quarters except for Q0, Q1, and Q17 because these are significantly shorter than the other quarters, which becomes important in the period analysis.

The Kepler data are known to suffer from various instrumental effects acting on different time scales, and affecting each observing quarter differently strongly. One of the most severe effects are drifts of stars across the detector, leaving long-term up- and downward trends in the light curves. These long-term signals can mimic the variability of slow rotators, and must be treated with caution. Previous attempts cleaned the data from instrumental signals by searching for shared signals across the detector. These so-called cotrending basis vectors are removed from the data by subtracting a linear combination of them from the time series (Kinemuchi et al. 2012; Stumpe et al. 2012; Smith et al. 2012). It was found that this approach bears the risk of underfitting because instrumental signals were not fully removed. An updated version of the Kepler pipeline separates the instrumental systematics by frequency (Stumpe et al. 2014). This approach, however, at times overcorrects the data and removes true astrophysical signals.

Even though the data used here were reduced with the latest pipeline, visual inspection showed that the reduction was far from being perfect. Many quarters still contain instrumental trends, which showed an increased variability compared to the other quarters for a given star. Such instrumental trends are often found every fourth quarter because the Kepler telescope rolls by 90 degrees every quarter, such that a certain target falls on the same CCD every four quarters. To identify such cases in an automated way, we used a common metric that characterizes the light curve variability: the variability range, R_var (Basri et al. 2010, 2011). This measure computes the difference between the 95th and 5th percentiles of the sorted differential intensities.

Here, we computed the variability range from the three-hour binned time series for each quarter individually, which we denote as R_var,q, and we also computed the median of all quarters R_{var, med}. After trying different thresholds, we found that all quarters with variabilities R_{var q} > 3 · R_{var med} should be discarded from the analysis. The full Table B.1, with all removed quarters, is available at the CDS.

2.2 Sample

The GPS method was originally developed and calibrated to measure periods of stars with near-solar effective temperatures, including the solar rotation period (see Sect. 3). However, it was found that the method also yields reliable periods for stars of later spectral type (Amazo-Gómez et al. 2020b). Thus, we used the revised stellar properties catalog of Mathur et al. (2017) and selected main-sequence stars with effective temperatures of T_eff ≤ 6500 K and surface gravities of log g ≥ 4.0. We further discarded all stars matching the Kepler eclipsing binary cata-log² by Kirk et al. (2016), as well as nine stars with residual instrumental systematics (with KIC numbers 6063291, 6126271, 7627042, 7800157, 11393439, 11414728, 11515679, 11805150, and 11808713). This selection leaves 142168 stars in total. These rather loose criteria should ensure that the targets lie on the main sequence (or close to it).

Recently, Berger et al. (2020, hereafter B20) published an updated catalog of fundamental parameters of ~ 186 000 Kepler stars taking into account Gaia DR2 parallaxes. The fundamental parameters of both catalogs clearly show some deviation. In particular, the B20 temperatures are roughly 200 K cooler and the B20 surface gravities are on average 0.2 dex smaller than the values given in Mathur et al. (2017). We find that 113 867 of the selected 142168 stars fulfill the chosen criteria for the parameters given in the catalog of B20. We note that the main goal of this study is to measure the rotation periods for as many stars as possible, while refrain from making assumptions on the accuracy of the fundamental parameters, as this is beyond the scope of this study.

Fig. 1 Example of a variable Kepler star and different period analysis methods applied to the data. Top panel: original (black dots) and 3 h-binned (red line) light curve of the star KIC 7831394. Second panel: global (black) and local (red) auto-correlation functions (ACFs). The best ACF period is found at 26.28 days. Third panel: global wavelet power spectrum (black). The vertical blue line indicates the position of the inflection point. The ratio between the power at the inflection point and the minimum is defined as signal-to-noise ratio (S/N), indicated by the orange arrow between the dotted lines. Bottom panel: gradient of the power spectrum (GPS). The maximum is found at the inflection period PIP = 4.76 days, marked by the blue line and the red asterisk.

3 Methods

3.1 The GPS method

First, we prepared the final light curves used in the following analysis. The flux of each observing quarter Q2 to Q16 is divided by its median and subtracted by unity, then outliers are removed, namely, those that exceed six times the median absolute deviation. The time series are appended and the resulting light curve is binned to three hours, forming the final light curves used in this analysis. The binning reduces the (photon + granulation) noise in the light curve and 3–6 h is a typical granulation timescale that also dominates the variability in this time interval. We note that rotational variability (i.e., our main focus) starts to dominate on timescales of 6 h and longer.

An example light curve of the star KIC 7831394 is shown in the top panel of Fig. 1. The light curve clearly shows variability on rotational time scales (especially after binning over 3 h), making it a promising candidate for our period analysis. In the second panel of Fig. 1, we compute the auto-correlation function (ACF) using the IDL function A_CORRELATE. We further subtracted the ACF minimum and normalized it by dividing out the maximum such that all ACF values lie between zero and one. For stars with periodic variability patterns, the ACF has proven to be a good tool for measuring the rotation period. The ACF of the full time series is shown in black. The highest peak is found at a period of P_rot,ACF = 26.28 days. As goodness measure of the period, we compute the local peak height (LPH) as the difference between the highest peak (red asterisk) and the mean of the two troughs on either side (see, e.g., Reinhold et al. 2021 for details).

Additionally, we computed the ACF for each quarter Q2 to Q16 individually and compute the mean ACF power in 0.1 day bins. This function is referred to as the local ACF, in contrast to the global ACF described above. The local ACF is shown in red on top of the global ACF in Fig. 1. For quite periodic light curves, both functions are very similar. However, for stars with more irregular variability, the local and global ACF often differ. As a consequence, the highest peak of the one does not match the one of the other, and the period is unclear. We return to this point in Sect. 3.2.

Then, we computed the wavelet power spectrum of the full time series (as seen in the third panel in Fig. 1). The highest power is found at a period of ~15.17 days. As we go on to see later, this period is likely an artifact of the high-pass filtering of the Kepler data. The bottom panel of Fig. 1 shows the gradient of the power spectrum (GPS). The gradient is computed using Eq. (3) in Shapiro et al. (2020; for details, we refer to Shapiro et al. 2020). The maximum of the gradient corresponds to the position of the inflection point (IP), namely, the point where the curvature of the high-frequency tail of the power spectrum changes its sign. The period at this inflection point, P_IP, is linked to the stellar rotation period, P_rot, via the simple equation:

(1)

where α is a calibration factor (see Sect. 4). The main idea behind the GPS method is that the high-frequency tail of the power spectrum is much less affected by the evolution of magnetic features than the power spectrum peak associated with the rotation period (see Fig. 3 in Shapiro et al. 2020). In this case, the inflection period is found at P_IP = 4.76 days, indicated by the blue line and the red asterisk. This period is linked to the rotation period by Eq. (1). Using the calibration factor α = 0.217 derived by Reinhold et al. (2022) yields a rotation period P_rot,GPS = 21.94 days.

This result is in good agreement with the rotation period of 22.11 days derived by S21. However, the period, P_rot,GPS, shows some discrepancy to the ACF period, and is very different from the one at the highest peak of the power spectrum. As we go on to see later, the ACF method often fails to detect the correct rotation period, the more irregular the variability gets.

3.2 Goodness of the GPS periods

For each light curve, the GPS method returns an inflection point period. However, it is not always clear if this period can be associated with the rotation of active features over the stellar surface. Thus, we define different goodness metrics for the derived GPS and ACF periods. Eventually, these metrics are combined to a point system that assigns a certain number of points to each star to assess the period reliability, namely: the higher the number of points, the more periodic the signal.

Most commonly, classical period analysis methods define the highest peak of the power spectrum (or the first ACF peak) as the strongest periodicity in the data. It was found by visual inspection of many different Kepler stars that the highest GPS peak nicely scales with the periodicity in the light curve. If the highest peak lies in the range 0.5–10 days, we save the inflection point period, P_IP, and the associated peak height, h_IP. The lower limit of 0.5 days is basically determined by the 3 h binning of the data (the Nyquist frequency would be 1/6 h), whereas the upper limit of 10 days should help deflect issues related to the data reduction (see Fig. 8 and subsequent discussion). According to Eq. (1), the considered range of inflection periods enables us to detect rotation periods in the range of ≈2.3–46.7 days.

To compute the power spectra, we used the IDL function WV_CWT, which returns the continuous wavelet transform, and we set the keyword dscale=1/32, which affects the peak height values. In this normalization, the peak height distribution h_IP ranges from 1 to 1.15, with a median height of 1.05. Strictly periodic stars have large peak heights h_IP > 1.06 to which we assign 1 point. Less periodic but still variable stars have peak heights 1.04 < h_IP < 1.06, to which we assign 0.5 points. Light curves that are completely dominated by noise have even smaller peak heights and thus get 0 points.

Even though it was shown that the highest power spectrum peak itself is not necessarily a good measure of the rotation period, we can still use the power spectrum to define another goodness measure: we call the ratio between the power at the inflection point and the minimum power of the spectrum the signal-to-noise ratio S/N (see third panel in Fig. 1). Similarly to h_IP, this quantity also scales with the periodicity. Based on a visual inspection of many light curves and the S/N distribution of all stars, we assigned 1 point if S/N > 50, 0.5 points if 10 < S /N < 50, and 0 points otherwise.

For the ACF, we already have a definition of the LPH as a goodness metric (see, e.g., Reinhold et al. 2021). These authors found that strong periodicity is usually found for LPH > 0.2 (1 point). Less periodic time series still reach values 0.1 < LPH < 0.2 (0.5 points) and purely noisy stars exhibit small LPH < 0.1 (0 points). We note that the LPH used here always refers to the global ACF. Additionally, we compared the global and the local ACF period. If these two periods agree within 10%, we add another 0.5 points.

By visual inspection, it was found that many light curves are dominated by noise and hardly show any variability and even less periodicity. Nevertheless, all the metrics above will return some values since they respond to any signal (even to pure noise) in the time series. However, these very quiet stars can be identified by comparing the variabilities of the unbinned and binned time series. In Fig. 2, we show the relative fluxes of the unbinned (black) and 6-hour binned (red) data as a histogram. The left panel of Fig. 2 shows that both flux histograms nicely overlap (also see top panel in Fig. 1), which means that the light curve of the star KIC 7831394 is dominated by (rotational) variability. On the contrary, the binned and unbinned flux distributions of the star KIC 11802969 (right panel) look very different. The 6-hour binning reduced the noise in the light curve and the remaining variability is small, which means that the photometric variability of this star is completely dominated by noise.

Instead of looking at flux distributions, we can simply compute the variability ranges R_var of the unbinned and 6-h binned time series and compute their ratio. If the R_var,6h/R_var ratio is close to unity, the light curve variability is dominated by rotation. By visual inspection, we assign 1 point if R_var,6h/R_var > 0.6, 0.5 points if 0.4 < R_var,6h/R_var < 0.6, and 0 points otherwise. We note that also the 3-hour binned value of R_var,3h could have been used instead, but the 6h-binning effectively reduced the noise even more.

Combining all metrics defined above yields points in the range 0 (pure noise) to 4.5 (highest periodicity) for each star. For the sake of clarity, the point metric is summarized in Table 1. We will return to this point system in Sect. 4.

4 Results

For the 142 168 stars in our sample, we could measure an inflection point period within 0.5–10 days for 141151 stars. We note that a fraction of this sample was assigned zero points in the end (compare Fig. 5). We now compare our results to a sample of stars with previously determined periods. This comparison will show whether the calibration factor α = 0.217 that was previously determined for models of solar-like stars (Reinhold et al. 2022) still holds for real data. In Fig. 3, we show the measured inflection periods P_IP against the rotation periods derived by McQ14 (upper panel) and S21 (lower panel) for the stars in common. The solid black line shows the relation given in Eq. (1) with the calibration factor α = 0.217. Both panels clearly show the linear dependence between both periods for the vast majority of stars. However, a second branch with (slightly more than) twice the inflection period, and consequently twice the rotation period, is visible in both panels. The origin of this “double period branch” is discussed in Sect. 4.4. We further note that the agreement of our period measurements (P_IP/α) becomes weaker for P_rot,Santos > 30 days. This is caused by the fact that the PDC-MAP pipeline does not preserve variability on these timescales, so the periods retrieved by classical periods analysis tools (as used in McQ14 and S21) are less reliable (also see Sect. 4.1).

Figure 3 reveals that there is indeed a linear dependence between the inflection period, P_IP, and the rotation period, P_rot. However, this relation is accompanied by large scatter. As mentioned in Shapiro et al. (2020), and later shown in detail by Reinhold et al. (2022), the calibration factor a has an intrinsic uncertainty of ~25%. In Fig. 4, we show the ratio of the inflection period to the rotation period derived by McQ14 (left panel) and S21 (right panel). The distributions have a Gaussian shape centered at 〈α〉 = 0.212 (left) and 〈α〉 = 0.213 (right) with standard deviations σ_α = 0.023 (left) and σ_α = 0.029 (right). These values are in good agreement with the value α = 0.217 derived by Reinhold et al. (2022). The small bump at twice these α values is associated with the upper branch in Fig. 3.

We now turn to the point system defined at the end of Sect. 3. In Fig. 5, we show the distribution of points allocated to each star. The sample of McQ14 is shown in blue and the stars in common with S21 are shown in red. It is obvious that the number of stars with previously determined periods steeply increases with the number of points and that the vast majority of those stars has the highest possible number of points. That means that these stars exhibit very periodic light curves, where the periodicity is picked up easily with standard tools.

As with almost every frequency analysis tool, setting thresholds for the detected peaks or the period significance is quite subjective. Visual inspections of many different quiet and active stars have led us to count all stars with ≥3 points as period detections. This threshold requires that at least one of the metrics has 1 point assigned. It is therefore not very conservative but a reasonable choice. We further note that below this threshold only very few rotation periods have been found in the surveys of McQ14 and S21, which makes it unlikely that many stars with measurable period have been missed. Although we cannot completely rule out the possibility that there exist stars with measurable periods below 3 points, lowering the threshold would make the periods less reliable. In total, the threshold is satisfied for 67 515 stars, which we focus on in the following.

In Fig. 6, we show the distribution of 67 515 GPS rotation periods of all stars with ≥3.0 points. Here, we used α = 0.213, derived in Fig. 4, to retrieve the rotation period. The periods determined by McQ14 and S21 are shown as blue and red curves, respectively. Compared to the previous surveys, the GPS method retrieves periods of a larger number of stars.

In particular, 20 397 new periods were detected that have not been reported before. Most of the newly determined periods are longer than ~28 days. Furthermore, the median variability of these 20 397 stars with newly detected periods equals R_var,3h = 0.085%, which is very close to the solar value R_var,Sun = 0.07% (compare Reinhold et al. 2020), and so much smaller than the average variability of all 67 515 stars (R_var,3h = 0.17%). We emphasize that the detection of these stars with near-solar rotation periods and variabilities is a clear benefit of the GPS method.

The reasons why these periods have been missed in previous surveys are manifold: McQ14 used quite conservative thresholds to detect periods, which removes many less periodic slow rotators. Moreover, these authors only considered Q3–Q14 data. Since that time, the Kepler pipeline has also changed several times. S21 analyzed light curves reduced with their own pipeline as well as those reduced with the latest version (DR25). These authors also combined different period analysis tools and used a machine-learning approach to finally detect periodicity. The GPS method could detect even more (and longer) periods because it does not require a repeatable spot pattern (as the other techniques do). The very short periods, however, cannot be accessed because we set a lower limit of 0.5 days to the inflection period, which translates into a lower rotation period limit of ~2.3 days.

We also checked the cases where a period was reported by McQ14 and/or S21 but not detected in this survey. There exist 11094 periods in these surveys that do not have a period reported here. The majority of these stars were initially not considered in this study because they either had effective temperatures greater than 6500 K or log g < 4.0 dex. Only 3313 out of the 11094 stars were considered in this study but had a point number smaller than 3.0. Of those, 1064 stars had a reported period outside the considered range of inflection periods between 0.5 and 10 days, and thus could not have been detected. However, the point distribution of the remaining 2249 stars continuously increases toward our lower limit of 3.0, with a mean number of 2.0 points. This result indicates that the imposed point limit might be lowered, which would eventually lead to even more period detections.

The better performance of the GPS method can also be seen in Fig. 7. Here, we show the detection rate as a function of stellar variability. In general, the detection rate increases with variability (without any variability, nothing can be detected). Interestingly, the steepest rise happens shortly after the solar mean variability (gray dashed line). This observation shows that the spot signals – compared to the photon noise in the light curves – start to dominate at this variability. We further see that the black and the red curves show very similar qualitative behavior: both curves steeply increase at small, solar-like variabilities R_var > 0.1%, and level off at a detection rate of ≈93% for variabilities R_var > 1%. However, it is obvious that the GPS method detects much more periods for stars with smaller variability.

The detection rate of the McQ14 sample (blue curve) shows a slightly different behavior for variabilities R_var > 0.3%, where it breaks through the black and red curves, and eventually reaches almost 100% for the most variable stars. As mentioned above, a different pipeline has been used in McQ14 that better preserved stellar variability. However, this cannot be the reason for the higher detection rate because all values used in Fig. 7 have been computed from the latest pipeline used in this study. We attribute the even higher detection rate of the variable stars to an extensive visual light curves inspection of McQ14, in contrast to the purely automated approach applied by S21 and in this study.

Fig. 2 Distribution of the 6 h binned (red) and unbinned (black) flux values of the stars KIC 7831394 (left) and KIC 11802969 (right).

Fig. 3 Inflection point period, PIP, vs. rotation period derived by McQ14 for 30 582 stars (upper panel) and S21 for 48 264 stars (lower panel). The solid black line shows the relation given in Eq. (1) with α = 0.217. Both panels exhibit an upper branch where the GPS method detects (slightly more than) the double of the inflection period.

Fig. 4 Distribution of α = PIP/Prot for the same stars, as shown in Fig. 3 for the rotation periods Prot of McQ14 (left panel) and S21 (right panel). The red curve shows a Gaussian fit with mean 〈α〉 = 0.212 and standard deviation σα = 0.023 (left panel) and 〈α〉 = 0.213 and σα = 0.029 (right panel).

Fig. 5 Distribution of the points metric for all 141 151 stars with measured inflection period (black) and the ones in common with McQ14 (blue) and S21 (red).

Fig. 6 Rotation period distribution of the 67 515 stars with points ≥3.0. The rotation periods derived by McQ14 and S21 are shown in blue and red, respectively.

Fig. 7 Detection rate as a function of the variability range, Rvar. The black curve shows the detection rate of the GPS method and the detection rates of the methods employed by McQ14 and S21 are shown in blue and red, respectively.

4.1 GPS vs. ACF periods

The rotation periods in McQ14, and to a large extent also those in S21, have been determined by the auto-correlation function (ACF). Thus, we also computed ACF periods for each star, with an upper period limit of 70 days. In this section, we compare the ACF and the GPS periods with each other to test the performance of both methods, and to show their limitations.

Figure 8 shows the ACF against the GPS periods for all stars where both periods could be measured and the local peak height (LPH) of the ACF peak at least fulfills the mild criterion LPH > 0.1. For periods less than 20 days, there is good overlap between the two methods, as indicated by the 1:1 line (red). The secondary branch at twice the ACF period is also visible.

One particularly striking feature is certainly the pile-up of ACF periods around 15 days, which is most visible at the top histogram. We attribute these periods to the high-pass filtering of the data in the latest data reduction and emphasize that most of these periods are of instrumental origin because such an accumulation of periods is not seen for the GPS periods. It is important to note that there are actually cases where also the GPS returns a period around 15 days and this periodicity is clearly seen in the light curves. We conclude that the ACF periods in the range of 10–20 days cannot be trusted without independent confirmation by another method. The pile-up of periods at around 15 days was also noted and dismissed by Basri et al. (2022), but without the benefit of this independent analysis.

In Fig. 9, we show the same as in Fig. 8 for different LPH thresholds. Additionally, the data are color-coded with the point metric defined above. It is apparent that the majority of stars with ≥3.0 points are located either along the 1:1 or the 2:1 line for all LPH thresholds. This tendency becomes even more evident with increasing the LPH threshold from 0.1 to 0.4 (upper left to lower right panel), which empties most of the other plot regions. Similar to Fig. 4, we compute α = P_IP/P_rot,ACF for the different LPH thresholds, and fit the distributions with a Gaussian. Very similar mean values have been derived, increasing from α = 0.204 ± 0.038 (LPH > 0.1) to α = 0.214 ± 0.023 (LPH > 0.4). We further note that also the ACF peak around 15 days becomes less pronounced as LPH increases.

4.2 Dependence of α on stellar parameters

We have seen that the mean a values are very similar for the different samples considered so far. In this section, we show how a depends on different stellar parameters and we consider what can be learned from that about the stars. In the following, we assume that the derived ACF periods are good measures of the stellar rotation period, at least for stars with a high number of points. In Fig. 10, we show α = P_IP/P_rat,ACF as a function of the ACF period P_rat,ACF. We see that the stars with the largest number of points accumulate around α = 0.213 (indicated by the black horizontal line) and a bit more than twice that value (i.e., the double period branch). The plot aptly shows that a does not show any dependence on rotation up to periods of ≈35 days. Beyond that period, there are much fewer stars with a high point score and, most importantly, the ACF periods become less reliable. Thus, we conclude that the same α value can be used to derive rotation periods of fast and slow rotators.

The following three figures show the dependence of α on the stellar fundamental parameters: effective temperature T_eff (Fig. 11), surface gravity log g (Fig. 12), and metallicity [Fe/H] (Fig. 13). Figure 11 reveals that α shows very little dependence on effective temperature from 4000–6000 K. For cooler stars below 4000 K, α seems to increase. The opposite effect is found for stars hotter than 6000 K, where α decreases. To emphasize this effect, we restrict the main a branch to those stars between 0.1 < α < 0.3 and ≥3.0 points and we overplot the mean alpha value in the 100 K wide temperature bins as violet star symbols.

For the hot stars, the relative decrease of α can be explained by the fact that the inflection period is sensitive to the spot lifetimes. Giles et al. (2017) showed that spots have shorter lifetimes on G- and F-type stars compared to later-type stars. This observation was recently confirmed by Basri et al. (2022), who used a similar approach as Giles et al. (2017) to assess the spot lifetimes. Reinhold et al. (2022) showed that the periods measured by the GPS method are shorter than the rotation period when the spot lifetimes are shorter than two complete rotations. As a consequence, the α values are smaller than the average value for these hot stars.

For stars cooler than 4000 K, however, we attribute the increase of α to another effect. Reinhold et al. (2022) further showed that the inflection point period is sensitive to the duration of a spot crossing. This dip duration depends on the spot latitude and the stellar inclination (neglecting the spot evolution on this comparatively short time scale). Spots at higher latitudes generate more sinusoidal dips in the light curves, and so have a longer dip duration than, for instance, equatorial spots. The same is true for lower latitude spots on a highly inclined star. We cannot break this degeneracy but we can argue that inclination, as a geometrical effect, is independent on effective temperature. Thus, we argue that inclination is partly responsible for the spread of α along the mean value but can be ruled out as explanation for the increased α here.

Instead, we propose that these very cool stars exhibit spots at higher latitudes than warmer stars. This idea is tested by a simple spot model in Sect. 4.4. Additionally, we tested if the center-to-limb variation (CLV) for cooler stars changes such that α might show an increase (see the appendix). However, this is not the case and can be ruled out as explanation here.

Figure 12 shows the dependence of α on log g. We can see that there is also very little dependence on surface gravity, except for the high and low gravity ends. These dependencies, however, are the same as those in the previous plot because of the dependence between surface gravity and effective temperature on the main sequence. For instance, the stars at log g > 4.8 with the highest α values have almost exclusively temperatures below 4000 K.

In Fig. 13, we show the dependence of α on metallicity. Along the main branch (black line), no dependence on metallicity is visible. This result is consistent with the latest tests of the GPS method on simulated data: Reinhold et al. (2022) found that α does not show any dependence on metallicity between −0.4 < [Fe|H] ≤ 0.4 dex for simulated time series of solar-like stars.

Fig. 8 ACF vs. GPS periods for more than 90 000 stars with LPH > 0.1, with the associated histograms to their sides. The GPS periods have been calculated using the same α = 0.213 as in Fig. 6. The red solid line shows the 1:1 identity.

Fig. 9 ACF vs. GPS periods for different LPH thresholds indicated at the top of each panel. The data are color-coded with the points system. The black line shows the 1:1 identity.

Fig. 10 ACF periods vs. α = PIP/Prot,ACE. The data are color-coded with the points metric. The black horizontal line indicates the mean value α = 0.213.

Fig. 11 Dependence of α on Teff. The colors and the horizontal line are the same as in Fig. 10. The violet star symbols show the mean α value in the 100 K wide temperature bins.

Fig. 12 Dependence of α on log g. The colors and the horizontal line are the same as in Fig. 10.

Fig. 13 Dependence of α on [Fe|H]. The colors and the horizontal line are the same as in Fig. 10.

4.3 Dependence of α on activity

In the previous section, we demonstrated that α shows very little dependence on rotation and the fundamental stellar parameters. Here, we test the dependence of α on two well-known measures of stellar activity. In Fig. 14, we show α as a function of the S -index, which is defined as the ratio of the flux in the Ca IIH and K lines, normalized to the flux in the R and V bands (see Vaughan et al. 1978 for details). This activity indicator was chosen because it is a well-established measure of stellar chromospheric activity (see, e.g., Noyes et al. 1984). Furthermore, it is independent of the Kepler data, in contrast to other photometric activity indices (such as the index S_ph used by Mathur et al. 2014 or the measure MDV used by Basri et al. 2013). The S-indices are taken from the catalog of Zhang et al. (2022) using the calibration to the Mount Wilson scale (Eq. (6) in Zhang et al. 2022). We find 18 797 matches between the LAMOST catalog and the stars in our sample. Figure 14 shows that a does not strongly depend on S-indices greater than 0.3. However, the spread of a becomes much stronger toward smaller S-indices. Since the S-index itself is not corrected for its dependence on effective temperature, this spread cannot solely be attributed to smaller activity.

Taking a look at the LAMOST S-index distribution shows that the vast majority of stars exhibits rather small S -indices between 0.1 and 0.3. This result is likely a selection effect because all these stars have effective temperatures greater than ≈4800 K. Moreover, the S-index distribution is quite narrow around the mean value of 0.2 (even for different effective temperature bins). We attribute this to the low resolution of the LAMOST instrument, which makes it more difficult to assess the true stellar activity level and, in particular, the true dependence of α on activity.

A well-known quantity closely related to activity is the photometric variability. Here, we use the quantity R_var,3h as a measure of the rotational variability, and show α as a function of R_var,3h in Fig. 15. This figure clearly shows that a is almost constant down to variabilities of R_var,3h = 0.2%, which nearly equals the solar maximum variability (Reinhold et al. 2020), as well as demonstrating that these stars have a high number of points. Down to smaller variabilities around R_var,3h = 0.1%, the a values start to show large spread. However, the horizontal branches extend down to very low variabilities with a moderate point number greater than 2 (green dots). This result clearly shows that the GPS method is able to detect the correct rotation period even for stars with very small variabilities. At the same time, also the ACF method yields the correct rotation period for these greenish stars (recalling that we defined α = P_IP/P_rot,ACF). For the cloud of blue dots below R_var,3h = 0.1%, likely the ACF method returns a wrong period because we did not impose any LPH threshold here.

Fig. 14 Dependence of α on the LAMOST S-index. The colors and the horizontal line are the same as in Fig. 10.

4.4 The double period branch

To better understand the origin of the double period branch, we employed a simple spot model. In this toy model, circular spots of a certain size and contrast can be placed on a sphere. The sphere can then be rotated, either as rigid body or differentially, and viewed from different inclination angles i. Since the GPS method is sensitive to the spot profile, we chose a simplified model with a single spot of fixed radius 5° sitting at random latitudes and longitudes (both uniformly distributed) and inclination angles (uniform in cos i). We arbitrarily chose a rotation period of 10 days and simulated five complete revolutions of rigid rotation, namely, a time series spanning 50 days.

We computed 500 models and applied the GPS method to them. The result is shown in Fig. 16. Here, we plot the inclination of the model star against the spot latitude because both quantities affect the spot profile in the light curve, and we color-coded each point with the derived a value. There is a sharp separation between blue points in the lower right half, where the correct value of α ≈ 0.21 is derived, and the yellow-greenish dots in the upper left part of the diagram, where roughly twice the correct a value is measured. We note that there are more dots in the lower right half, which means that we can measure the correct period for a large number of possible spot and inclination angle alignments, and that the transition between the two regimes is rather discrete (also compare Fig. A.1).

This result again confirms that the GPS method is sensitive to the spot profile in the light curve (Reinhold et al. 2022). It is known that both higher latitude spots and/or strongly inclined stars render the spot profiles more sinusoidal, which is equivalent to an increase of the dip duration in the light curve. Since the inflection point is proportional to the dip duration, the inflection periods become larger, reaching roughly twice the correct α value when the dip duration equals one full rotation period.

Thus far in the study, we have mainly promoted the novel GPS method, especially when stellar variability is low and more irregular. However, also this method has its shortcomings, for instance, it sometimes detects twice the correct rotation period. A visual inspection of those light curves where GPS detects twice the period derived by McQ14 or S21 shows that this mostly happens in cases when the variability is very periodic and so the LPH values are large. Consequently, in such cases, the ACF method yields the correct period with high confidence.

To better quantify this observation, we once again consider the rotation period sample of S21. For this purpose, we define the number of stars in the main branch, N_main, as those stars where the GPS and the S21 period differ by less than 10%. Similarly, the number of stars in the double period branch, N_double, are defined such that the GPS and twice the S21 period values differ by less than 10%. In the left panel of Fig. 17, we show the “double fraction” N_double/(N_main + N_double) as a function of the LPH. We see that this fraction is rather flat for LPH > 0.2 (the increase in the last bin can safely be attributed to the much smaller number of stars in both branches). This so-called “double floor” with a median of ≈6% accounts for the cases with very symmetric light curves. Assuming the employed toy model correctly distinguishes between the main and double period branch, this percentage must be considered as the fraction of stars with rotation axes sufficiently inclined towards Earth that their light curves look rather sinusoidal even for spots at low latitude.

A different behavior is seen for small LPH values where the double fraction increases. In this regime (LPH < 0.2), the light curve variability is less periodic, and the GPS method is superior to standard methods (Reinhold et al. 2022). As a consequence, here the “double period” measured by GPS is likely to be the correct rotation period. We also note that in this regime much fewer stars are contained in both branches because it is more difficult to detect periods in light curves with shallow periodicity (using standard methods such as S21). This fact is accounted for by the larger error bars assuming for each branch.

In complement to the LPH dependence, we show the double fraction as a function of the GPS rotation period in the right panel of Fig. 17. Also, different regimes can be observed here: as mentioned previously, the GPS method is not sensitive to very short periods, and so the increase of the double periods for P_rot,GPS< 5 days can be ignored. Over the broad period range of between 5–30 days, the double fraction shows a shallow decrease, which may be consistent with a decrease of spot latitude with rotation period but this is rather speculative. For periods greater than 30 days, a steep increase of the double period is observed (although the number of stars also strongly decreases beyond 40 days). This period range is exactly the regime where GPS is superior to classical methods and likely returns the correct period, similarly to the LPH < 0.2 range in the left panel of this figure.

Fig. 15 Dependence of a on the variability range Rvar. The colors and the horizontal line are the same as in Fig. 10.

Fig. 16 GPS method applied to model light curves. The dots show the inclination and the spot latitude of the 1-spot models. The color-coding shows the derived a value.

Fig. 17 Fraction of stars in the double period branch as a function of the LPH (left panel) and GPS rotation period (right panel) for the stars in common with S21.

4.5 Selection of final rotation period

Figure 17 nicely demonstrates the region of validity of either method, and so helps to define a “final” rotation period. As the final rotation period, P_rot,fin, we used the ACF period if 0 < P_rot,ACE≤ 10 days and LPH ≥ 0.1. This decision is based on the fact that the ACF finds very accurate periods for fast rotators with rather sinusoidal shape and moderate to high peak heights.

In the period range of 10 < P_rot,ACE≤ 20 days, we see that the Kepler pipeline induces an artificial pile-up (see Fig. 8). We then compared the ACF and GPS periods in this range and used the ACF period if they differed by less than 10% and the ACF fulfilled the minimal requirement of LPH ≥ 0.1; otherwise, the GPS period was used.

For P_rot,ACE > 20 days, the periodicity becomes weaker, the LPH values smaller and so, the GPS method superior to the ACF. Thus, we used the GPS periods here, as well as so for the few cases where no ACF period could be measured. The chosen parameters for the definition of the final period are summarized in Table 2.

We note that these period thresholds are highly subjective (as with most thresholds), but the result depends on the expertise of the authors with various kinds of light curves and frequency analysis methods. Based on the period and LPH thresholds, the number of final rotation periods determined by the GPS and the ACF methods varies. Using our thresholds, we detected 67 163 final rotation periods with 17 246 ACF and 49 917 GPS periods. We note that there are 352 stars (67 515 – 67 163) with ≥3.0 points that have a measured GPS period but no final period assigned. These stars mostly do not satisfy the very mild LPH ≥ 0.1 thresholds. A parameter table for all stars with measured GPS period and ≥3.0 points assigned is given in Table C.1.

5 Summary and conclusions

In this study, we applied the novel GPS method to the light curves of tens of thousands main-sequence stars observed by the Kepler telescope. Although this huge data set has previously been combed for rotation periods (e.g., see Reinhold et al. 2013; McQuillan et al. 2014; Santos et al. 2021), we show that the GPS method was able to measure 20 397 periods that have not been detected before. One reason for that is that the GPS method is superior to standard frequency analysis methods (such as the ACF) for detecting rotation periods of stars with small and irregular variability.

Another important result seen for these 20 397 new periods was that their average rotation period was found to be ~28 days, while their variability is R_var,3h = 0.085%, so both the rotation period and the variability amplitude are very close to the solar values. The detection of these “solar-like” rotation periods is certainly an advantage of the GPS method.

Another benefit of the GPS method is that the high-frequency tail of the power spectrum, namely, the period regime in which GPS searches for the inflection point is much less sensitive to instrumental residuals. In stark contrast, the ACF method reveals a pile-up of periods between 10–20 days (for low to mid LPH, i.e., 0.1 < LPH < 0.2). We conclude that one cannot trust ACF periods in this period and LPH range without independent confirmation by the GPS method. Another option would be to apply a customized data reduction that does not undergo a high-pass filter (likely responsible for the period pile-up at ~15 days) and then to apply a standard method such as ACF.

Furthermore, our work suggests a temperature dependence of the GPS calibration factor α for T_eff < 4000 K. Employing a simplified spot model, we argue that this increase of the inflection point periods is caused by spots located at higher latitudes for these cool (likely early M dwarf) stars. We emphasize that this information could not be extracted by other methods before. Thus, we conclude that GPS rotation periods should be combined with other spectroscopic measurements to better constrain potential spot locations on the stellar surface.

In total, we were able to measure 67 163 final rotation periods by combining the ACF and the GPS methods. Compared to previous surveys, we find that 86.2% of the McQ14 and 77.4% of the S21 periods agree within 20% with our final rotation period, P_rot,fin. For periods lower than 20 days, the difference mostly originates from the double branch, whereas for periods greater than 20 days the intrinsic uncertainty of α dominates. We note that the difference in this period regime is smaller for the McQ14 sample and attribute this result to their more conservative threshold of LPH > 0.3. In the period range from 10–20 days, the GPS method is superior to classical ones (especially ACF) because it is not affected by the data reduction that would otherwise lead to the unphysical pile-up of ACF periods around 15 days (see top panel in Fig. 8). We note that the ACF method usually returns more accurate periods for very periodic light curves; but the shallower the periodicity gets, the more reliable the GPS period becomes, with respect to the period regime of P_rot,ACF > 20 days.

In summary, this study deals with the largest set of rotation periods known to date. This work clearly demonstrates the power of novel methods, such as GPS, for detecting new rotation periods even in large data sets that have been trawled through for periods before. Other promising methods might be Gaussian processes, possibly equipped with a suitable kernel function that does not a priori require a periodicity parameter. It is solely such methods that may reveal the true rotation periods of less active stars and will thereby help to improve on previous solar-stellar comparison studies.

Appendix A Center-to-limb variation

Fig. A.1 GPS method applied to model light curves for three different combinations of photosphere and spot temperatures. The α distributions shift toward the smaller values when the temperatures decrease.

We tested another possible explanation for the increase of α toward cooler stars, namely, the center-to-limb variation (CLV). This effect describes the decline of the surface intensity from the center of the star (brightest) toward the limb (darkest). This is not only true for the quiet photosphere but also affects (in our model) cooler active region. For instance, a spot appears darker at the disc center because its contrast to the quiet photosphere is larger at the center than at the limb. In turn, the CLV might have an effect on a because the inflection point period is sensitive to the profile of individual spot crossing events (Reinhold et al. 2022).

The CLV has been extensively studied for the Sun (e.g., Claret 2000). Now, model grids of stars of various effective temperatures, surface gravities, and metallicities are available (Kostogryz et al. 2022) that were computed using the state-of-the-art radiative transfer code MPS-ATLAS (Witzke et al. 2021).

Using the simple spot model employed in Sect. 4.4, we computed different sets of light curves (500 realizations each) for three different combinations of quiet photosphere and spot temperatures. Here, we fixed the surface gravity at log g = 4.4 and metallicity at M/H = 0.0 dex. In each of the three sets, the same inclination and spot position (longitude and latitude) have been used to focus on the temperature dependence.

Figure A.1 shows the result of this small study. We chose the three combinations of photosphere and spot temperatures as T_phot = 5800 K & T_spot = 4800 K (black), T_phot = 4800 K & T_spot = 3800 K (orange) and T_phot = 4000 K & T_spot = 3600 K (red). These combinations were chosen arbitrarily, but do appear reasonable since the difference between the spot and photosphere temperatures becomes smaller with decreasing effective temperature (see Fig. 7 in Berdyugina 2005).

As expected, one can see that the majority of the distributions scatters around α = 0.2 and a smaller fraction around roughly twice that value. However, the models with smaller temperatures are shifted to smaller values of α, in contrast to our initial assumption. Although the origin of this effect is not entirely clear, we conclude that CLV can be ruled out as an explanation for the increase of α toward cooler stars. The exact reason for the observed shift would require a more thorough analysis, which is beyond the scope of this study.

Appendix B Bad quarters

Appendix C Rotation period table

References

All Tables

All Figures

Fig. 1 Example of a variable Kepler star and different period analysis methods applied to the data. Top panel: original (black dots) and 3 h-binned (red line) light curve of the star KIC 7831394. Second panel: global (black) and local (red) auto-correlation functions (ACFs). The best ACF period is found at 26.28 days. Third panel: global wavelet power spectrum (black). The vertical blue line indicates the position of the inflection point. The ratio between the power at the inflection point and the minimum is defined as signal-to-noise ratio (S/N), indicated by the orange arrow between the dotted lines. Bottom panel: gradient of the power spectrum (GPS). The maximum is found at the inflection period PIP = 4.76 days, marked by the blue line and the red asterisk.

In the text

	Fig. 2 Distribution of the 6 h binned (red) and unbinned (black) flux values of the stars KIC 7831394 (left) and KIC 11802969 (right).
In the text

	Fig. 3 Inflection point period, PIP, vs. rotation period derived by McQ14 for 30 582 stars (upper panel) and S21 for 48 264 stars (lower panel). The solid black line shows the relation given in Eq. (1) with α = 0.217. Both panels exhibit an upper branch where the GPS method detects (slightly more than) the double of the inflection period.
In the text

	Fig. 4 Distribution of α = PIP/Prot for the same stars, as shown in Fig. 3 for the rotation periods Prot of McQ14 (left panel) and S21 (right panel). The red curve shows a Gaussian fit with mean 〈α〉 = 0.212 and standard deviation σα = 0.023 (left panel) and 〈α〉 = 0.213 and σα = 0.029 (right panel).
In the text

	Fig. 5 Distribution of the points metric for all 141 151 stars with measured inflection period (black) and the ones in common with McQ14 (blue) and S21 (red).
In the text

	Fig. 6 Rotation period distribution of the 67 515 stars with points ≥3.0. The rotation periods derived by McQ14 and S21 are shown in blue and red, respectively.
In the text

	Fig. 7 Detection rate as a function of the variability range, Rvar. The black curve shows the detection rate of the GPS method and the detection rates of the methods employed by McQ14 and S21 are shown in blue and red, respectively.
In the text

	Fig. 8 ACF vs. GPS periods for more than 90 000 stars with LPH > 0.1, with the associated histograms to their sides. The GPS periods have been calculated using the same α = 0.213 as in Fig. 6. The red solid line shows the 1:1 identity.
In the text

	Fig. 9 ACF vs. GPS periods for different LPH thresholds indicated at the top of each panel. The data are color-coded with the points system. The black line shows the 1:1 identity.
In the text

	Fig. 10 ACF periods vs. α = PIP/Prot,ACE. The data are color-coded with the points metric. The black horizontal line indicates the mean value α = 0.213.
In the text

	Fig. 11 Dependence of α on Teff. The colors and the horizontal line are the same as in Fig. 10. The violet star symbols show the mean α value in the 100 K wide temperature bins.
In the text

	Fig. 12 Dependence of α on log g. The colors and the horizontal line are the same as in Fig. 10.
In the text

	Fig. 13 Dependence of α on [Fe|H]. The colors and the horizontal line are the same as in Fig. 10.
In the text

	Fig. 14 Dependence of α on the LAMOST S-index. The colors and the horizontal line are the same as in Fig. 10.
In the text

	Fig. 15 Dependence of a on the variability range Rvar. The colors and the horizontal line are the same as in Fig. 10.
In the text

	Fig. 16 GPS method applied to model light curves. The dots show the inclination and the spot latitude of the 1-spot models. The color-coding shows the derived a value.
In the text

	Fig. 17 Fraction of stars in the double period branch as a function of the LPH (left panel) and GPS rotation period (right panel) for the stars in common with S21.
In the text

	Fig. A.1 GPS method applied to model light curves for three different combinations of photosphere and spot temperatures. The α distributions shift toward the smaller values when the temperatures decrease.
In the text