© ESO 2018

1. Introduction

In the mid 1960s, O. C. Wilson and collaborators started flux measurements at the centers of the Ca II H+K lines “[...] for the purpose of initiating a search for stellar analogues of the solar cycle.” (Wilson 1968). Ten years later, Wilson (1978) presented long-term Ca II H+K flux measurements of 91 main sequence stars, providing the first evidence for cyclic activity in some stars. Over the years an empirical activity index describing the emission in the H and K line cores was established (Vaughan et al. 1978), now referred to as the Mount Wilson S-index. This index can be converted into the chromospheric H and K line surface flux (Middelkoop 1982; Rutten 1984). Subtracting the photospheric flux from this quantity, one arrives at the currently most commonly used activity indicator $\log R_{\text{HK}}^{\prime }$, where $R_{\text{HK}}^{\prime }$ denotes the chromospheric flux in the Ca II H+K lines normalized by the bolometric flux (Linsky et al. 1979; Noyes et al. 1984).

Most stellar activity cycles known these days rely on periodicities measured in long-term observations of the S-index (Baliunas et al. 1995; Saar & Brandenburg 1999; Böhm-Vitense 2007). Constraining the basic physical parameters that determine the nature of the activity cycle, especially its period, is a matter of current research. Various authors (Saar & Brandenburg 1999; Böhm-Vitense 2007; Lehtinen et al. 2016) identified a strong dependence of the cycle period on the rotation period, and, as a consequence thereof, the existence of certain sequences (active and inactive sequences) in the P_rot − P_cyc plane. Böhm-Vitense (2007) suggested that the two distinct sequences might indicate the existence of two dynamos operating at different depths in the star. In contrast, studies by Reinhold et al. (2017) and Saikia et al. (2018) strongly question the existence of the two sequences, especially for the active sequence. Moreover, Metcalfe et al. (2016) and Metcalfe & van Saders (2017) proposed a secular increase of the cycle period once the star reaches a critical Rossby number Ro ∼ 2.

The periodic changes of the Ca II H+K emission originate from a change in the magnetic field generated by the stellar dynamo. Concurrently, the emergence of active regions in the photosphere (such as dark spots and bright faculae) is expected during activity maximum, as is observed in the Sun. Over the course of the solar cycle, variability is seen both in the total solar irradiance (TSI) data (e.g., Shapiro et al. 2016), as well as in the S-index (e.g., Egeland et al. 2017). In the Sun the photometric (TSI) and chromospheric time series (disk-integrated Ca II K-line) are in phase (Preminger et al. 2011). At activity maximum the Sun is brightest due to the presence of bright faculae regions, slightly overcompensating the contributions of dark spots.

Active stars are known to show photometric variability caused by dark spots and bright faculae, as well as enhanced chromospheric emission in the Ca II H+K line cores. Whether or not all active stars exhibit an underlying activity cycle, that is, whether the surface activity changes periodically or randomly over time, remains unclear. We aim to understand the physical mechanism that drives long-term photometric and chromospheric brightness variations. To reach this goal, we search for cycle periods in contemporaneous photometric and chromospheric time series. Assuming a common origin of the periodicity, we fit the same cycle period to both time series, which allows us to measure photometric and chromospheric amplitudes, and the phase difference between the variations in both time series. We show that the phase difference reveals the dominant type of stellar activity, that is, dark spots or bright faculae. Moreover, the phase difference shows a strong dependence on the stellar activity level, the rotation period, and the photometric amplitude. These correlations may help to understand stellar dynamos as a function of activity level, rotation period, and age. In the following we study a sample of 30 stars, for which chromospheric activity cycles have been reported for 18 stars (17 in Baliunas et al. 1995, 14 in Saikia et al. 2018, and 6 in Olspert et al. 2018), and photometric activity cycles for 6 stars (Messina & Guinan 2002). Taking advantage of the simultaneous analysis of both time series, we detect activity cycles in 27 stars in total. In 8 stars thereof, activity cycles have never been measured before.

2. Data

In the current study, data from three different sources have been used. Long-term V band and Strömgren b and y photometric time series have successfully been requested from Messina & Guinan (2002) and Lockwood et al. (2007), respectively. For all stars, chromospheric emission data from the Mount Wilson survey are publicly available¹. In total our sample contains 30 stars. The Julian date range, number of data points, mean time stamp, brightness, and standard deviation of each observing season for both the photometric and chromospheric time series are available in electronic form at the CDS. The basic stellar parameters of our sample are listed in Table 1. Effective temperatures, surface gravities and metallicities are taken from the references given. The spectral types have been adopted from Messina & Guinan (2002) and Lockwood et al. (2007). The different data sets are briefly described in the following.

2.1. Photometric data

Messina & Guinan (2002, 2003) studied the long-term magnetic activity of five young solar analogues (HD 1835, HD 20630, HD 72905, HD 129333, HD 206860), and a young K dwarf (HD 82443). These authors compiled photometric observations available from the literature (with different band passes) and new observations taken in Johnson V band into one time series covering up to ∼16 years. The authors state that the typical standard deviation of the V band magnitudes is equal to 0.007 mag, which we adopt as photometric uncertainty of the measurements. For all stars in their sample, Messina & Guinan (2002, 2003) detected a photometric activity cycle. For details, we refer the reader to the above publications.

Lockwood et al. (2007) presented up to 20 years of differential Strömgren b and y photometry of 32 stars with contemporaneous observations of the Mount Wilson S-index (see also Lockwood et al. 1997; Radick et al. 1998, 2018) with the purpose of studying long-term photometric and chromospheric variability. We received data for 26 out of the 32 stars in their sample. Photometric uncertainties σ_by have been adopted from Table 2 in Lockwood et al. (2007). In total, our sample consists of 30 stars, as two stars (HD 1835 and HD 129333) are present in both photometric samples.

2.2. Chromospheric emission data

The Mount Wilson survey has taken chromospheric emission data from 1966 to 1995 for the majority of the stars in their sample, and for 35 stars until the year 2001. Thereof, 24 stars belong to our sample and have been observed for more than 33 years. The remaining six stars of our sample have shorter observing time spans between 9 and 28 years. Recently, this unprecedented data set became publicly available¹. The Mount Wilson S-index is the main quantity of interest describing the amount of chromospheric emission in the centers of the Ca II H+K lines. A definition of this quantity can be found in Vaughan et al. (1978). From the year 1977 onwards, an observational weight, WT, accounting for the photon noise is supplied. This weight can be converted into a relative S-index uncertainty σ_S/S = WT^−1/2 (Duncan et al. 1991). These uncertainties were used in our computations whenever available.

3. Detection of activity cycle period and phase

In Fig. 1 we show an example of contemporaneous photometric (top panel) and chromospheric (bottom panel) observations. Both time series show short-term (seasonal) variability, partly caused by active regions on the stellar surface rotating in and out of view, and long-term brightness changes over several seasons, possibly caused by an underlying activity cycle. However, detecting a cycle period (and in a next step, the phase difference between the two time series) is challenging because the scatter in the individual observing seasons hampers the detection of long-term brightness changes. In the following we present a model capable of detecting stellar rotation signals. This model is applied to each observing season and subtracted from the data in order to enhance the visibility of the signal from the activity cycle. Afterwards we search for a cycle period present in both time series in order to compute the phase difference between the two time series.

Fig. 1. Photometric (top panel) and chromospheric (bottom panel) time series of the star HD 1835. The red box indicates the observing season shown in Fig. 2.

3.1. The Gaussian process model

To account for the variation due to stellar rotation in each observing season, we apply a quasi-periodic Gaussian process² (GP) model that accounts for periodic and long-term trends in the data. The GP was implemented using the python package George (Ambikasaran et al. 2015). A quasi-periodic GP model has successfully been applied to Kepler data to infer stellar rotation periods (Angus et al. 2018). The underlying kernel function of this GP model reads

$$\begin{array}{cc}\hfill k(x_i,x_j)& =h^2\underset{\mathrm{SE}}{\underbrace{exp(-\frac{{(x_i-x_j)}^2}{2l^2})}}·\underset{\mathrm{ESS}}{\underbrace{exp(-\mathrm{\Gamma }{sin}^2\left[\frac{\pi }{P}(x_i-x_j)\right])}}\hfill \end{array}$$(1)

$$\begin{array}{cc}& +{\sigma }^2{\delta }_{\mathit{ij}}.\hfill \end{array}$$(2)

In general, a kernel function k(x_i, x_j) describes the covariance between arbitrary data points x_i and x_j. The quasi-periodic kernel is composed of two well-established kernels, namely the squared exponential (SE) and the exponential sine squared (ESS) kernel. The SE kernel contains the length scale parameter l, which accounts for the covariance between observations taken at times x_i and x_j. The ESS kernel contains two parameters: the period P to account for periodicities in the data, and the parameter Γ to control the harmonic behavior of the periodicity. The parameter h is sometimes called the “amplitude” of the covariance because it is linked to the variance in the time series. The second term, σ² δ_ij, accounts for the white noise in the time series, with δ_ij being the Kronecker delta. For all stars in our sample, the rotation periods have been measured (Baliunas et al. 1996), and we use them as initial period guesses P in the GP prediction. Finally the full set of so-called hyper-parameters {σ, h, l, Γ, P} is estimated by fitting the above model to the data assuming Gaussian-distributed uncertainties.

We applied uniform priors to all parameters, with the period being limited to ±30% of the initial value. It is worth noting that these limits do not exclude possible harmonics of this period needed to fit a double-dip shape in the time series arising from spots on opposite sides of the star. The parameter Γ in the GP accounts for the harmonic behavior of the fit, given an initial period. Additionally, we checked the initial values reported by Baliunas et al. (1996) and Messina & Guinan (2002), and can confirm all rotation periods within ±20%.

The parameter estimation is performed using the affine-invariant ensemble sampler emcee (Foreman-Mackey et al. 2013; Goodman & Weare 2010). In total we used 32 walkers and 2000 iterations when running emcee. The walkers were initialized around the solution found by a previous maximum likelihood estimation (MLE) step and in almost all cases this provided a starting point very close to the final parameters found after the MCMC. For each run we calculated the integrated auto-correlation time and thinned the chains accordingly. In addition to computing the effective sample size, we also visually inspected the chains and the marginal posterior distributions to ensure that they looked reasonable. This allows us to estimate the marginalized posterior probability density of the period given the other hyper-parameters. Additionally, the mean flux value of each observing season is predicted by the MCMC sampler.

In Fig. 2 we present data from the observing season indicated by the red box in Fig. 1 and the optimized GP model including 1σ uncertainties. The GP model predicts a rotation period of P_rot = 7.80 d, close to the literature value of P_rot = 7.756 d (Messina & Guinan 2002). The length scale parameter l accounts for long-term trends over the course of the observing season, likely caused by the active regions’ lifetimes.

Fig. 2. Selected observing season of Mount Wilson S-index time series of the star HD 1835 (indicated by the red box in Fig. 1). The GP prediction is shown as solid line with the 1σ uncertainties of the prediction in blue. The GP model predicts a rotation period of Prot = 7.80 d.

The main purpose of using the GP model is to remove the rotational short-term activity rather than measuring the rotation period (a Lomb-Scargle periodogram or a Fourier transform would be sufficient for the latter). Whenever it was not possible to detect a periodic signal (either due to the absence of active regions or their short lifetime compared to the stellar rotation period), the GP model predicts a flat line, which is a fair model in quiet epochs. The GP model presented in Eq. (1) is applied to each observing season in each data set.

3.2. The reduced time series

To remove the rotation signature from each observing season, we subtract the GP model from the data and add the mean flux value of the season as predicted by the MCMC sampler. This preserves the seasonally averaged mean flux values and pronounces the signature of the activity cycle. The reduced time series of the star HD 1835 is shown in Fig. 3. Comparing the reduced time series to the original data in Fig. 1 shows that the scatter in each season is reduced significantly. Alternatively to using all data points, the time series can be simplified by compressing each season to the predicted mean value, and using the standard deviation of the residuals (i.e., data minus GP model) as uncertainties. Since these uncertainties reflect the residual activity of each season, and not the uncertainty of the predicted mean value, we decided to use the full time series in the following. We note that we also searched for cyclic behavior of the rotation periods inferred from the GP analysis as done by Messina & Guinan (2003). These authors interpreted (cyclic) changes of the rotation period in terms of a possible butterfly diagram. However, our analysis of the “rotation period time series” does not show a significant cycle period.

Fig. 3. Reduced photometric (top panel) and chromospheric (bottom panel) time series of Fig. 1. Sine fits to the data are shown in red with their cycle periods given at the top of each panel. In the lower panel a sine wave with a long period of 28.1 yr has been subtracted (see Table B.1).

The reduced time series are used to search for cycle periods in the following way. Using Lomb-Scargle periodograms (Zechmeister & Kürster 2009) we first search for periodicity in the photometric and chromospheric time series individually. In a second step we search for a single cycle period fitting both data sets at the same time (see Sect. 3.3). However, finding the correct cycle period requires visual inspection of each time series. The data often show long trends on the order of the observing time span, appearing as the highest peak in the periodogram. Whenever periods longer than 75% of the observing time span are detected, we subtract a long sine wave from the data and compute the periodogram again. Since it cannot be decided whether such a trend is caused by a very long activity cycle or if it is an instrumental artifact, we report these long periods in Table B.1 as potential long activity cycles.

Another problem arises for the Mount Wilson data before the year 1980. At these early epochs the data sampling is much sparser due to a different observing strategy. It seems that including these early observations generates a long trend over the full observing time span, although no systematic offset or a common trend could be identified for all stars. Therefore, for 18 stars in our sample, early Mount Wilson data have been excluded from the further analysis. All Mount Wilson observing seasons used in the following can be found at the CDS.

3.3. Detection of a single cycle period and phase difference

After the subtraction of potential long activity cycles we search for periods present in both the photometric and chromospheric time series. Computing the Lomb-Scargle periodogram and selecting the highest peak usually yields a good estimate of the true cycle period. In Fig. 3 the Lomb-Scargle periodogram returns a cycle period of P_cyc = 6.4 yr for the photometric time series, and a cycle period of P_cyc = 7.1 yr for the chromospheric time series. To compute the phase difference between the two time series one has to fit the same period in both data sets simultaneously. Thus, a sine fit was applied to fit both time series at the same time using a single period, with individual amplitudes, phases and offsets for the two data sets. As an initial guess of the true cycle period, we use the mean of the cycle periods inferred from the individual time series. The result of the fit to both data sets of HD 1835 simultaneously is shown in Fig. 4. To combine the two different types of data into a dimensionless quantity the time series have been normalized by their sine amplitudes. The simultaneous sine fit returns a single cycle period of P_cyc = 7.07 ± 0.02 yr. From this fit we can compute the phase difference between the two time series to check for possible phase shifts. These results are presented in Sect. 4.2. The sine fit was performed using the Levenberg-Marquardt least-squares minimization method. We also tried fitting the sine wave using an MCMC approach to retrieve more realistic uncertainties. This approach did not converge unless we restricted the priors of the parameters to a very small region around the best values returned by the least-squares fit. Given the nature of the model, where interpretability in terms of phase differences is preferred over being optimal in terms of model structure (i.e., more activity-cycle-based), the application of MCMC was not ideal and a simpler least-squares approach was favored. Thus, we report the formal uncertainties of the fit parameters, which should be considered as lower limits. We also tried to model the cycle behavior with the same GP model that was used for the rotation, and to use cross correlation to determine the relative phases; however, this was found to be unreliable. One reason for not trying to add a phase difference into the kernel of the GP is that understanding the behavior of the resulting covariance matrix is not straight forward.

Fig. 4. Photometric (blue) and chromospheric (red) time series of the star HD 1835. The simultaneous sine fit to the data (same colors, respectively) yields a cycle period of Pcyc = 7.07 ± 0.02 yr. Both data sets were normalized by their respective sine amplitudes.

4. Results

4.1. Cycle periods

Cycle periods could be detected for 27 stars in total. These periods are given in Table B.1. Following the cycle classification as introduced by Baliunas et al. (1995), we provide a quality flag for each cycle period that has been assigned by visual inspection of both the photometric and the chromospheric time series individually. A higher number indicates stronger periodicity of the data, with integer values between 1 (weak periodicity), and 3 (excellent periodicity). Two example light curves are shown in Appendix A. We identified a subset of eight stars exhibiting an excellent chromospheric activity cycle. We return to these stars in Sects. 4.2 and 5. To quantify the goodness of the periodicity we computed the reduced chi-square values of the fit to each data set individually (${\chi }_{\text{phot}}^2$ and ${\chi }_{\text{chrom}}^2$), and also to both data sets simultaneously (${\chi }_{\text{sim}}^2$). The values are reported in the last three columns of Table B.1. However, we found that these values do not correlate with the periodicity of the data. This is owing to the fact that both the period and the amplitude of the activity cycle can vary from cycle to cycle. For instance, the varying cycle amplitude of the star HD 201091 in the Mount Wilson time series creates quite a high value for ${\chi }_{\text{chrom}}^2$, although the periodicity is evident in this star (flag_chrom = 3). Moreover, the shape of the activity cycle does not need to be sinusoidal (as assumed by the sine fit), even when the period and amplitude are constant over time.

We compare our cycle periods to measurements from the literature in Fig. 5. Panel a shows the measurements by Baliunas et al. (1995) that contain all 30 stars from our sample. These authors found periodicity in 17 out of the 30 stars. We also detect periodicity in these 17 stars, although no periodicity could be detected in the photometric time series of the star HD 201092. In panel a–c this star is shown as an open circle to indicate that the periodicity has only been measured in the chromospheric time series. Taking advantage of up to ten years of additional Mount Wilson data, we detect periodicity in 11 stars where Baliunas et al. (1995) did not report a cycle. Comparing the period measurements to one another generally shows good agreement within ±20% (dotted lines), although some outliers are apparent. In HD 120136, Baliunas et al. (1995) detected a “poor” cycle of ∼12 yr, which can be disproved by considering the additional data. In HD 149661, these authors detected a “good” cycle of ∼17.4 yr. Again we cannot confirm this long period, which might be due to the sparse data of the early observations. For the star HD 190007, Baliunas et al. (1995) detected a “fair” ∼13.7 yr cycle. This period might be a modulation of the potentially very long cycle of ∼38 yr that has been subtracted in our analysis. The two stars HD 1835 and HD 115404 lie close to the upper 20% limit. The periods determined by Baliunas et al. (1995) are likely affected by trends in the early observations.

Fig. 5. Comparison of our cycle period measurements to results from (a) Baliunas et al. (1995), (b) Saikia et al. (2018), (c) Olspert et al. (2018), and (d) Messina & Guinan (2002). The open circle indicates the star HD 201092 for which only a chromospheric cycle period was measured. Crosses indicate secondary periods connected by vertical dotted lines to the primary periods. The dashed line shows the 1:1 identity, and the dotted lines to either side mark the ±20% range. We set the cycle period to zero for stars in common for which no cycle period was reported.

At first glance, the measurements by Saikia et al. (2018) in panel (b) of Fig. 5, corresponding to the 14 stars that we have in common, seem to show greater deviations from our results. These overall greater deviations might arise from the fact that these latter authors also considered the much sparser data of the early seasons in their analysis. For HD 1835, HD 82443, and HD 115383, long periods have been subtracted in our analysis, which explains the large deviation. HD 20630 shows a long trend in the early data, which is no longer visible in the later seasons. Our cycle periods for HD 1835 and HD 20630 have been identified as secondary cycles in Saikia et al. (2018). For HD 149661, Saikia et al. (2018) claimed a primary cycle of 15.3 yr and a secondary cycle of 7.7 yr, respectively. Our analysis does not confirm any of these periods, not even when considering all available “raw” data. Furthermore, we cannot reproduce the period claimed for HD 18256, although the deviation from our results is smaller in this case. For HD 190007, Saikia et al. (2018) claimed a cycle period of 5.4 yr. This period is visible in the computed periodograms with a smaller peak height. Therefore, we expect this period to be a modulation of the longer one. Although most deviations can be explained by a different data analysis, some period detections remain mysterious.

A comparison of our results to those by Olspert et al. (2018) in panel c of Fig. 5 shows remarkably good agreement. In total, 17 stars are contained in both samples. All six stars for which both authors detected periodicity almost lie on the 1:1 line. As mentioned above, for the star HD 201092 we do not detect a photometric cycle period. All cycle period detections in panels a–c are based solely on the Mount Wilson data set. In contrast to Baliunas et al. (1995) and Saikia et al. (2018), Olspert et al. (2018) use a similar approach as described here based on GP models. Additionally, these authors account for long-timescale trends in the time series, which explains the overall good agreement to our results.

Finally, we compare our measurements to the cycle periods detected by Messina & Guinan (2002) in panel d of Fig. 5. In contrast to the measurements in panels a–c, these periods are solely based on photometric data. We find good agreement to our results for five of the six stars within the ±20% range. For HD 72905, two different photometric data sets have been combined in the year 1993 (compare Table 5 and Fig. 6 in Messina & Guinan 2002). Comparing the photometric and the chromospheric time series of this star raises the question of whether this combination may have introduced an offset between the different photometric data sets because no such offset is seen in the chromospheric time series. We calculated the median brightness of the data before and after the year 1993, and added the difference to the data taken after the year 1993. This removes the long period found by Messina & Guinan (2002), and yields a cycle period of ∼4.2 yr, which is almost twice the period reported by these authors.

4.2. Activity cycle phase difference

The phase difference Δϕ between the fits to the photometric and chromospheric time series reveals the dominant source of activity on the stellar surface, that is, dark spots or bright faculae regions. A phase difference of Δϕ ≈ 0 indicates that the star is brightest at activity maximum. This is the case in the Sun. The opposite case of minimum brightness during activity maximum, that is, Δϕ ≈ 0.5, demonstrates that the star is dominated by dark spots. In Fig. 6 we show the phase difference Δϕ against the mean stellar magnetic activity level $\log R_{\text{HK}}^{\prime }$.

Fig. 6. Activity cycle phase difference Δϕ against the magnetic activity index $\log R_{\text{HK}}^{\prime }$. Dots show dwarf stars (luminosity class V) and diamonds show (sub-) giants (luminosity class III, IV, IV-V). The colors show different spectral types: F stars (blue), G stars (orange), and K stars (red). The vertical dashed line indicates the approximate location of the Vaughan–Preston gap, with its boundaries indicated by the dotted lines, adopted from Fig. 2 of Noyes et al. (1984). The horizontal dashed lines at Δϕ = 0.5 and Δϕ = 0 mark the region of spot and faculae domination, respectively. The yellow star indicates the position of the Sun, adopting the mean value logR′HK,  Sun = −4.94 from Egeland et al. (2017).

The maximum phase difference between two time series equals |Δϕ| = 0.5. In Fig. 6 we plot Δϕ between −0.25 and 0.75 to allow for a scatter around Δϕ = 0 and Δϕ = 0.5. We find two distinct groups of stars: one scattering around Δϕ = 0.5 and the other scattering around Δϕ = 0. The two groups are clearly separated by their mean activity level $\log R_{\text{HK}}^{\prime }$. The transition from Δϕ = 0.5 to Δϕ = 0 occurs around $\log R_{\text{HK}}^{\prime }\hspace{0.17em}\approx -4.75$, which is the location of the Vaughan–Preston (VP) gap (Vaughan & Preston 1980; Noyes et al. 1984). This clear separation tells us that active stars (above the VP gap) are spot-dominated whereas inactive stars (below the VP gap) are faculae-dominated.

As stars age they spin down because angular momentum is lost due to stellar winds. This spin-down is accompanied by a decrease of the magnetic activity (Skumanich 1972). Hence, we expect all stars to transition from spot-dominated to faculae-dominated surface activity, depending on their activity level and age. Taking advantage of the close connection between activity and rotation (Fig. 8 in Noyes et al. 1984), we show the phase difference versus the Rossby number (instead of the activity level on the abscissa) in Fig. 7.

Fig. 7. Activity cycle phase difference Δϕ against the Rossby number Ro = Prot/τ. Colors, symbols, and horizontal dashed lines have the same meaning as in Fig. 6. The vertical dashed line indicates the transition around Ro = 1. The yellow star indicates the position of the Sun, RoSun = 2.16, adopting the values B − VSun = 0.653 (Ramírez et al. 2012), and the synodic rotation period Prot,  Sun = 26.3 d.

The Rossby number Ro = P_rot/τ equals the rotation period divided by the convective turnover time τ. The latter has been computed following Noyes et al. (1984). As in Fig. 6, the separation into a spot-dominated and a faculae-dominated group is evident with the active stars at Ro ≲ 1, and the inactive stars at Ro ≳ 1. Below Ro = 1 rotation is the driving source for surface activity; above Ro = 1 stars rotate slower, and convection plays a more important role. The advantage of the Rossby number over the rotation periods is that the convective turnover time accounts for the mass dependence of the rotation period.

The presence of two distinct groups in Fig. 8 provides further evidence for our interpretation that stars around Δϕ = 0 are faculae-dominated and stars around Δϕ = 0.5 are spot-dominated. Stars of the spot-dominated group have photometric amplitudes that are on average six times larger than those of the faculae-dominated group. This result is consistent with the work of Shapiro et al. (2014), who showed that with the increase of magnetic activity the spot contribution to the photometric variability increases faster than the faculae contribution. Consequently, less active stars with smaller photometric variability are faculae-dominated, while more active stars with larger photometric variability are spot-dominated. The total photometric amplitude of 0.37 mmag for the solar cycle has been computed by Shapiro et al. (2016). Such a small value was explained in Witzke et al. (2018) by the incident combination of solar fundamental parameters and the spectral location of the Strömgren b and y pass bands. Furthermore, from Table B.1, it can be seen that 8 of the 27 stars with measured cycle periods exhibit a chromospheric activity cycle classified as excellent (flag_chrom = 3). Concurrently, their photometric amplitudes are shallow. All these stars belong to the group around Δϕ = 0, that is, they are faculae-dominated.

Fig. 8. Activity cycle phase difference Δϕ against the photometric sine amplitude. Colors and symbols have the same meaning as in Fig. 6. The photometric amplitude of the solar cycle (yellow star) was adopted from Shapiro et al. (2016). These authors calculated a total amplitude of 0.37 mmag, which results in a sine amplitude of ampphot,  Sun = 0.185 mmag.

5. Discussion

Taking advantage of contemporaneous photometric and chromospheric time series, we have determined activity cycle periods and the corresponding phase differences for 27 stars in our sample. For two stars in our sample, no activity cycle could be detected. HD 124570 and HD 143761 are both highly inactive and show almost flat variability both in the photometric and the chromospheric time series. HD 201092 shows an activity cycle of ∼10.8 yr in the chromospheric time series, whereas no periodicity could be detected in the photometric time series.

In Figs. 6–8 we showed that faculae-dominated stars are inactive slow rotators with small photometric amplitudes, compared to their active spot-dominated counterparts with shorter rotation periods and larger amplitudes. The transition from spot- to faculae-dominated activity occurs at an activity level of $\log R_{\text{HK}}^{\prime }\hspace{0.17em}\approx -4.75$ (the VP gap), and around a Rossby number Ro ≈ 1. This result is consistent with the work of Montet et al. (2017). These authors searched for photometric activity cycles by analyzing Kepler full-frame images (FFIs). They correlated the long-term brightness variations of the FFIs with photometric short-term variability measured by the magnetic activity proxy S_ph from the Kepler light curves (see Mathur et al. 2014). These authors found that most stars with rotation periods shorter than 15–25 days show short-term photometric variability amplitudes that are anti-correlated with the seasonal brightness in the FFIs. This supports the interpretation that the activity cycles of fast rotating stars are spot-dominated. For stars with longer rotation periods, a correlation between short-term variability amplitudes and average brightness suggests that their activity cycles are faculae-dominated. In their sample of Sun-like stars with 0.6 ≲ B − V ≲ 0.7, the rotation period range between 15 and 25 days marks the transition from spot- to faculae-dominated activity. Using activity-rotation relations from Mamajek & Hillenbrand (2008), this period range coincides with the VP gap, where we find the transition from spot- to faculae-dominated activity. Using the same relations, this period range can be transferred to Rossby numbers: Ro(P_rot = 15 d) = 1.25 and Ro(P_rot = 25 d) = 2.09 using the mean color B − V_mean = 0.65.

To further support the finding of a transition from spot- to faculae-dominated activity, we attempted a similar analysis to Montet et al. (2017) by correlating the following photometric time series: the first one consists of the seasonal mean brightness values, and the second one contains the variability range R_var of each observing season (see Basri et al. 2010, 2011). This quantity measures the rotation-induced variability in each season similar to the quantity S_ph used by Montet et al. (2017). We correlated these two time series to see whether the long- and short-term variability is in or out of phase. Our results are, however, inconclusive given the relative imprecision of the ground-based data.

Moreover, we find that the majority of faculae-dominated stars in our sample exhibit a chromospheric activity cycle with a strong sinusoidal shape (flag_chrom = 3). All eight of these flag_chrom = 3 stars lie below the VP gap. Activity-age relations from Mamajek & Hillenbrand (2008) reveal ages greater than 2.7 Gyr for these eight stars, suggesting that the transition from complex to sinusoidal or smooth cycles is complete at this age. This confirms earlier results by Oláh et al. (2016) who analyzed activity cycles of Mount Wilson stars, finding that slowly rotating stars show smooth activity cycles, whereas more rapidly rotating stars exhibit more complex cycles, with a clear transition from complex to smooth cycles at ages between ∼2 and 3 Gyr.

For six stars in our sample, the magnetic field topology has been measured from Zeeman Doppler imaging (See et al. 2016). Their results suggest that spot-dominated stars show toroidal fields, whereas faculae-dominated stars exhibit mostly poloidal fields. A toroidal field has indeed been detected in HD 20630 and HD 131156A, for which we detect spot-dominance. HD 10476 and HD 201091 are faculae-dominated stars with smooth chromospheric activity cycle (flag_chrom = 3) showing strong poloidal fields and almost no toroidal fields. The spot-dominated star HD 206860 does not fall in either group, showing 50% of each field component, and HD 120136 shows a slightly more toroidal than poloidal field, although our results suggest faculae-dominance for this star. Unfortunately, the number of stars with observed magnetic field topology that overlap with our sample is too small to either confirm or reject the scenario suggested by See et al. (2016).

5.1. Can the VP gap explain the dearth of intermediate rotation periods?

McQuillan et al. (2013) first discovered a lack of rotation periods around ∼25 days in the rotation period distribution of the Kepler M dwarfs. It has been shown that this shortage extends to hotter stars covering a rotation period range of ∼15–25 days (Reinhold et al. 2013; McQuillan et al. 2014; Reinhold & Gizon 2015). McQuillan et al. (2013) explained the detected bimodality of the rotation periods by an overlap of two stellar populations of different age in the Kepler field, with the younger stars rotating faster, on average, than the older stars. Davenport (2017) also discussed the bimodality seen in the Kepler rotation period distribution. By using absolute G-band magnitudes from the Gaia mission, this author demonstrated that the McQuillan et al. (2014) sample might be strongly contaminated by subgiants. Because of their deeper convective envelopes, their rotation period distribution is significantly different from dwarf stars, which might explain the previous non-detection of the bimodality in these stars. Removing those evolved stars from the sample, Davenport (2017) showed that the bimodality persists for stars hotter than 5000 K. This author further finds that the fast and slow rotators in his sample also exhibit a different distribution of the total proper motion. Hence, the existence of two kinematically separate groups would favor the explanation of two epochs of stars formation in the Kepler field. This explanation is further supported by the work of Davenport & Covey (2018), showing that the bimodality also correlates with Galactic height, which is assumed to be related to stellar age.

Over the last few years, astrophysicists have debated whether or not the lack of intermediate active stars at the VP gap can explain the lack of intermediate rotation periods. To test this hypothesis, we revisit the rotation period distribution of the Kepler sample measured by Reinhold & Gizon (2015) in Fig. 9. The rotation periods are plotted against B − V color, and the data points are color-coded with the variability range R_var (we return to this quantity in Sect. 5.2). The dearth of rotation periods becomes visible at ∼15 days at B − V ≈ 0.9 and increases to ∼25 days at B − V ≈ 1.5. We find that the Kepler rotation periods are well represented by gyrochronology relations from Barnes (2010): the blue and red lines show the gyrochronology isochrones (gyrochrones) for ages of 300 Myr and 2550 Myr, respectively. The color-dependent lack of rotation periods follows a gyrochrone with an age of ∼800 Myr (dotted black line).

Fig. 9. Kepler rotation periods Prot against B − V color taken from Reinhold & Gizon (2015). The data points are color-coded with the variability range Rvar. Black and yellow symbols show spot- and faculae-dominated stars, respectively, for the current sample. As before dots show dwarf stars (luminosity class V) and diamonds show (sub-) giants (luminosity class III, IV, IV-V). The solid blue, dotted black, and solid red lines show 300 Myr, 800 Myr, and 2550 Myr isochrones, respectively, using gyrochronology relations from Barnes (2010).

The age of 2550 Myr has been selected because it refers to the activity level of $\log R_{\text{HK}}^{\prime }\hspace{0.17em}\approx -4.75$ at the VP gap. In Fig. 9 we over-plot the spot- and faculae-dominated stars of the current sample as black and yellow symbols, respectively. For colors B − V >  0.5, where gyrochronology relations are well calibrated (Barnes 2007), we find all spot-dominated stars below the 2550 Myr isochrone, and all the faculae-dominated stars above this age. In the color range B − V <  0.5 where gyrochronology is not well calibrated we use the same activity-age relations as before to infer the ages of the two F stars HD 18256 and HD 120136. This analysis yields ages of 2.65 Gyr and 2.45 Gyr, respectively, consistent with the age of 2.55 Gyr associated with the activity level at the VP gap. Hence, the expected age of stars at the VP gap clearly separates the spot- and faculae-dominated stars in our sample. This result demonstrates that the observed lack of stars with ages of ∼800 Myr cannot be explained in terms of the VP gap.

5.2. Alternate explanation for the dearth region

The Kepler stars in Fig. 9 are color-coded by their variability range R_var, which describes the amplitude of the rotational variability. Interestingly, the dearth region around ∼800 Myr is dominated by stars with small variability amplitudes (mostly blue points). Moreover, there seems to be a variability gradient towards the dearth region from either side, that is, stars younger and older than ∼800 Myr show larger variability. Based on the result that faculae are the dominant ingredient of the surface activity for the activity cycle after 2.55 Gyr, we expect that the faculae surface coverage may increase from an earlier age. To explain both the dearth of stars and the decreased variability, we propose the idea of a period determination bias, and suggest the following scenario. Fast rotating stars are able to generate strong magnetic fields, which in turn create large and long-lived spots. As the stars slow down, the magnetic field topology changes such that the magnetic field is no longer able to maintain large spots, and instead generates smaller active regions surrounded by bright faculae. This process diminishes the variability amplitudes and continues until faculae become so prominent on the surface that bright and dark active regions partially cancel each other. The stars in the dearth region might have smaller and more homogeneously distributed spots, which simultaneously would decrease the variability, and render period determination more difficult. Hence, many stars with small variability amplitudes will not be detected due to certain amplitude constraints.

To better illustrate the idea of spot and faculae cancelation, we applied the Spectral And Total Irradiance REconstruction model (SATIRE, see Fligge et al. 2000; Krivova et al. 2003), which is the state-of-the-art model of solar brightness variations, for synthesizing light curves brought about by the mixture of facular and spot features. Such a model light curve is shown in Fig. 10. We have considered a star with a rotational period of 30 days at various inclination angles. The light curve was synthesized for 1400 days, which is comparable to the total duration of the Kepler observations. The mean surface coverage of spots was chosen to be 0.1% (i.e., the star is slightly more active than the Sun). We used a log-normal distribution of spot sizes from Baumann & Solanki (2005), and assumed a linear decay of spot sizes with time (following Martinez Pillet et al. 1993). We assumed that each spot is surrounded by a facular region, which is a reasonable assumption for modeling rotational stellar variability (see, e.g., Shapiro et al. 2017). The faculae-to-spot coverage ratio at the day of emergence was taken to be 3.4, and the lifetimes of faculae were considered to be three times longer than those of spots, which roughly corresponds to the solar case.

Fig. 10. Top panel: model light curve with a rotation period of 30 days viewed at an inclination angle of 90°. Lower left: Lomb-Scargle periodogram of the time series. Lower right: auto-correlation function of the time series. The vertical dashed line indicates the model rotation period.

The lower left panel in Fig. 10 shows the Lomb-Scargle periodogram of the time series. The power was normalized to unity (i.e., perfect sinusoidal periodicity); none of the peaks are significant, following the threshold of a peak height of 0.1 as used by Reinhold & Gizon (2015). The lower right panel shows the auto-correlation function (ACF) of the time series, which was used by McQuillan et al. (2014) to search for periodicity in the Kepler sample. The peaks indicate a periodicity with multiples of 10 days, although with low power. The correct rotation period cannot be recovered either by the periodogram or by the ACF method. Thus, we conclude that 1) this star would have been considered as non-periodic using the periodogram approach, or 2) would have been considered as periodic using the ACF method, however measuring an incorrect rotation period. This model light curve has a variability range of R_var = 0.56%; if a significant period were detected, it would appear as a blue dot in Fig. 9. We note that a similar spot and faculae cancelation happens in the solar case. Shapiro et al. (2017) showed that the 27-day peaks in the spot and facular component of the solar brightness variability (calculated for the period 1996–2015) cancel each other, so that the power spectrum of solar brightness variations is almost flat around the rotational period.

Figure 9 further shows that the variability increases after an age of ∼800 Myr, suggesting that spots again become the dominant source of activity. This scenario is supported by the presence of the four spot-dominated stars between the ∼800 Myr and the 2.55 Gyr isochrone. It is worth mentioning that Lockwood et al. (2007) also detected spot dominance for these four stars. Once the stars slow down further and pass the VP gap, faculae become the dominant source of surface activity on activity cycle time scales. This is supported by the large abundance of stars with small variability amplitudes above the transitionary age of 2.55 Gyr. Although the transition from spot to faculae domination has been inferred from variability on activity cycle time scales, this observation suggests that at least some rotation periods have been inferred from long-lived faculae rather than short-lived spots, as seen in the quiet Sun. Based on observations of solar activity this idea has already been proposed by Pettersen et al. (1992). During solar activity maximum bright and dark features cancel each other, hampering a reliable period detection, whereas rotation can be inferred from faculae during activity minimum. To decide on the dominant surface feature, spot-crossings during transits (e.g., Morris et al. 2018) or Doppler images would be needed.

We emphasize that classifying a star as spot- or faculae-dominated refers to the type of feature that causes it to appear overall fainter or brighter during the peaks of its activity cycle. On rotational timescales, spots are usually the dominant source of variability seen in the light curve because they are more localized than faculae (Shapiro et al. 2016). Moreover, faculae-dominated stars older than 2.55 Gyr are still able to maintain spots, as is observed in the Sun. Applying the simultaneous sine fit to the individual observing seasons, and searching for the same rotation period instead of the same cycle period in both data sets, we indeed find spot dominance for most of the stars where a rotation period could be detected. That implies that a star can be spot-dominated on rotational timescales, whereas it can be faculae-dominated on activity cycle timescales.

6. Conclusion

In this paper we have determined activity cycle periods for 27 out of the 30 stars in our sample. By computing the phase difference between the photometric and the chromospheric time series we have shown that stellar activity transitions from spot to faculae domination at an activity level of $\log R_{\text{HK}}^{\prime }\hspace{0.17em}\approx -4.75$. This value is equal to the activity level at the Vaughan–Preston gap. Using activity-rotation-age relations this activity level is associated with a Rossby number Ro ≈ 1 and an age of ∼2.55 Gyr. Furthermore, the transitional age of ∼2.55 Gyr is consistent with measurements of a change of the activity cycle shape and the global magnetic field topology from the literature (Oláh et al. 2016; See et al. 2016). These results provide observational constraints on the underlying dynamo. The connection between the dominant source of surface activity and the magnetic field topology sounds promising to either infer the type of surface activity from magnetic field measurements or vice versa.

We further conclude that the lack of intermediate active stars at the Vaughan–Preston gap cannot explain the lack of intermediate rotation periods discovered by the Kepler telescope. Nevertheless, the results presented in this work strongly suggest that the importance of faculae contribution to the big picture of stellar activity has so far been underestimated. We propose a scenario to explain the lack of intermediate rotation periods based on the cancelation of the spot and faculae contribution at rotation time scales. The proposed idea has been tested by simulating light curves with a solar faculae-to-spot ratio. We found that neither the periodogram nor the auto-correlation function were able to detect the correct rotation period. We plan to refine the light-curve model and extend it to stars of various spectral types in a forthcoming publication.

Acknowledgments

The research leading to the presented results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (2007–2013)/ERC grant agreement no 338251 (StellarAges). A.I.S. acknowledges the support by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement no. 715947). We thank S. Messina and G. Lockwood and collaborators for kindly sharing their data with us. We also like to thank T. Metcalfe, H. Schunker, and R. H. Cameron for valuable discussion. The HK_Project_v1995_NSO data derived from the Mount Wilson Observatory HK Project, which was supported by both public and private funds through the Carnegie Observatories, the Mount Wilson Institute, and the Harvard-Smithsonian Center for Astrophysics starting in 1966 and continuing for over 36 years. These data are the result of the dedicated work of O. Wilson, A. Vaughan, G. Preston, D. Duncan, S. Baliunas, and many others.

References

Appendix A: Example light curves

Fig. A.1. Top panel: chromospheric time series of the star HD 81809 flagged as excellent periodicity (flagchrom = 3). The red curve shows the best sine fit to the data. Bottom panel: Lomb-Scargle periodogram of the time series.

Fig. A.2. Top panel: chromospheric time series of the star HD 13421 flagged as weak periodicity (flagchrom = 1). The red curve shows the best sine fit to the data. Bottom panel: Lomb-Scargle periodogram of the time series.

Appendix B: Additional table

All Tables

All Figures

	Fig. 1. Photometric (top panel) and chromospheric (bottom panel) time series of the star HD 1835. The red box indicates the observing season shown in Fig. 2.
In the text

	Fig. 2. Selected observing season of Mount Wilson S-index time series of the star HD 1835 (indicated by the red box in Fig. 1). The GP prediction is shown as solid line with the 1σ uncertainties of the prediction in blue. The GP model predicts a rotation period of Prot = 7.80 d.
In the text

	Fig. 3. Reduced photometric (top panel) and chromospheric (bottom panel) time series of Fig. 1. Sine fits to the data are shown in red with their cycle periods given at the top of each panel. In the lower panel a sine wave with a long period of 28.1 yr has been subtracted (see Table B.1).
In the text

	Fig. 4. Photometric (blue) and chromospheric (red) time series of the star HD 1835. The simultaneous sine fit to the data (same colors, respectively) yields a cycle period of Pcyc = 7.07 ± 0.02 yr. Both data sets were normalized by their respective sine amplitudes.
In the text

Fig. 5. Comparison of our cycle period measurements to results from (a) Baliunas et al. (1995), (b) Saikia et al. (2018), (c) Olspert et al. (2018), and (d) Messina & Guinan (2002). The open circle indicates the star HD 201092 for which only a chromospheric cycle period was measured. Crosses indicate secondary periods connected by vertical dotted lines to the primary periods. The dashed line shows the 1:1 identity, and the dotted lines to either side mark the ±20% range. We set the cycle period to zero for stars in common for which no cycle period was reported.

In the text

Fig. 6. Activity cycle phase difference Δϕ against the magnetic activity index $\log R_{\text{HK}}^{\prime }$. Dots show dwarf stars (luminosity class V) and diamonds show (sub-) giants (luminosity class III, IV, IV-V). The colors show different spectral types: F stars (blue), G stars (orange), and K stars (red). The vertical dashed line indicates the approximate location of the Vaughan–Preston gap, with its boundaries indicated by the dotted lines, adopted from Fig. 2 of Noyes et al. (1984). The horizontal dashed lines at Δϕ = 0.5 and Δϕ = 0 mark the region of spot and faculae domination, respectively. The yellow star indicates the position of the Sun, adopting the mean value logR′HK,  Sun = −4.94 from Egeland et al. (2017).

In the text

Fig. 7. Activity cycle phase difference Δϕ against the Rossby number Ro = Prot/τ. Colors, symbols, and horizontal dashed lines have the same meaning as in Fig. 6. The vertical dashed line indicates the transition around Ro = 1. The yellow star indicates the position of the Sun, RoSun = 2.16, adopting the values B − VSun = 0.653 (Ramírez et al. 2012), and the synodic rotation period Prot,  Sun = 26.3 d.

In the text

	Fig. 8. Activity cycle phase difference Δϕ against the photometric sine amplitude. Colors and symbols have the same meaning as in Fig. 6. The photometric amplitude of the solar cycle (yellow star) was adopted from Shapiro et al. (2016). These authors calculated a total amplitude of 0.37 mmag, which results in a sine amplitude of ampphot,  Sun = 0.185 mmag.
In the text

Fig. 9. Kepler rotation periods Prot against B − V color taken from Reinhold & Gizon (2015). The data points are color-coded with the variability range Rvar. Black and yellow symbols show spot- and faculae-dominated stars, respectively, for the current sample. As before dots show dwarf stars (luminosity class V) and diamonds show (sub-) giants (luminosity class III, IV, IV-V). The solid blue, dotted black, and solid red lines show 300 Myr, 800 Myr, and 2550 Myr isochrones, respectively, using gyrochronology relations from Barnes (2010).

In the text

	Fig. 10. Top panel: model light curve with a rotation period of 30 days viewed at an inclination angle of 90°. Lower left: Lomb-Scargle periodogram of the time series. Lower right: auto-correlation function of the time series. The vertical dashed line indicates the model rotation period.
In the text

	Fig. A.1. Top panel: chromospheric time series of the star HD 81809 flagged as excellent periodicity (flagchrom = 3). The red curve shows the best sine fit to the data. Bottom panel: Lomb-Scargle periodogram of the time series.
In the text

	Fig. A.2. Top panel: chromospheric time series of the star HD 13421 flagged as weak periodicity (flagchrom = 1). The red curve shows the best sine fit to the data. Bottom panel: Lomb-Scargle periodogram of the time series.
In the text