Issue 
A&A
Volume 638, June 2020



Article Number  A69  
Number of page(s)  10  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202037666  
Published online  12 June 2020 
Shapes of stellar activity cycles
^{1}
Department of Physics, University of Helsinki, PO Box 64, 00014 Helsinki, Finland
email: teemu.willamo@helsinki.fi
^{2}
Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
^{3}
Department of Computer Science, Aalto University, PO Box 15400, 00076 Aalto, Finland
Received:
5
February
2020
Accepted:
17
April
2020
Context. Magnetic activity cycles are an important phenomenon both in the Sun and other stars. The shape of the solar cycle is commonly characterised by a fast rise and a slower decline, but not much attention has been paid to the shape of cycles in other stars.
Aims. Our aim is to study whether the asymmetric shape of the solar cycle is common in other stars as well, and compare the cycle asymmetry to other stellar parameters. We also study the differences in the shape of the solar cycle, depending on the activity indicator that is used. The observations are also compared to simulated activity cycles.
Methods. We used the chromospheric Ca II H&K data from the Mount Wilson Observatory HK Project. In this data set, we identified 47 individual cycles from 18 stars. We used the statistical skewness of a cycle as a measure of its asymmetry, and compared this to other stellar parameters. A similar analysis has been performed for magnetic cycles extracted from direct numerical magnetohydrodynamic simulations of solartype convection zones.
Results. The shape of the solar cycle (fast rise and slower decline) is common in other stars as well, although the Sun seems to have particularly asymmetric cycles. Cycletocycle variations are large, but the average shape of a cycle is still fairly well represented by a sinusoid, although this does not take its asymmetry into account. We find only slight correlations between the cycle asymmetry and other stellar parameters. There are large differences in the shape of the solar cycle, depending on the activity indicator that is used. The simulated cycles differ in the symmetry of global simulations that cover the full longitudinal range and are therefore capable of exciting nonaxisymmetric largescale dynamo modes, and wedge simulations that cover a partial extent in longitude, where only axisymmetric largescale modes are possible. The former preferentially produce positive and the latter negative skewness.
Key words: stars: activity / stars: chromospheres / Sun: activity
© ESO 2020
1. Introduction
The shape of the 11year sunspot cycle is not perfectly symmetric, but characterised by a faster rise from minimum to maximum and a slower decline from maximum to minimum (Waldmeier 1935). Another common feature that deviates from a sinusoid shaped cycle is the typical double peak, known as the Gnevyshev gap (Gnevyshev 1963). Gnevyshev (1967, 1977) suggested that the solar cycle generally consists of two waves of activity. There is an asymmetry in solar activity between the northern and southern hemisphere (e.g. Newton & Milsom 1955; Deng et al. 2016). Norton & Gallagher (2010) studied the solar cycle separately on each hemisphere and concluded that differences in the hemispheres cannot explain the Gnevyshev gap, but a mechanism must be producing it for both hemispheres. One factor that might affect it is the complexity of active regions. Simple active regions, with unipolar or bipolar sunspot groups, on average appear earlier in the solar cycle than more complex active regions (Nikbakhsh et al. 2019). Thus the simple regions dominate the first peak of the maximum, and the complex regions only have a notable effect on the latter peak. Feminella & Storini (1997) found that the activity dip in the Gnevyshev gap is more evident in highenergy phenomena, such as the occurrence of longlasting energetic flares, while the occurrence of flares and other phenomena with lower energies tend to follow the simple 11year cycle. The cycle amplitude and length of the rising phase are also anticorrelated. This is known as the Waldmeier effect (Waldmeier 1935, 1939).
There is no reason for stellar analogues of the solar cycle to be perfectly symmetric either, but they are usually fitted with simple sinusoids, and not much attention has been paid to their shape. Reinhold et al. (2017) showed that cycles derived from the variability of Kepler stars deviate from simple sinusoids, the average shape showing a sharp maximum and flattened minimum. The authors also hypothesised that this effect might have a temperature dependence because it was weak for the coolest stars.
The solar cycle has been modelled with many different mathematical formulations accounting for their asymmetry (Nordemann 1992; Elling & Schwentek 1992; Hathaway et al. 1994; Volobuev 2009; Du 2011). Takalo & Mursula (2018) applied the principal component analysis to the solar cycle and divided it into two components, an average cycle component, which always has the same shape, with varying period and amplitude, and one component that varies from cycle to cycle.
One parameter that has been used to measure asymmetries of solar cycles is the skewness. This is a measure of asymmetry that is commonly used in statistics. Ramaswamy (1977) reported a relation between the ratio of the maximum sunspot number of the following cycle to the current cycle μ and the skewness γ of the current cycle as
Lantos (2006) improved the correlation by separately considering even and odd cycles, and derived the following formulae:
Stellar cycles, however, have not been modelled as extensively. Garg et al. (2019) found the Waldmeier effect in stars from Mount Wilson observatory data. They also studied stellar cycle asymmetries by fitting similar functions as for the solar cycle, and they calculated the skewness. Pipin & Kosovichev (2016) found from numerical meanfield simulations for solartype stars that magnetic cycles of a higher amplitude are more asymmetric, until at some amplitude, the asymmetry becomes saturated.
The methods that are commonly used to study stellar cycles are not capable of taking cycle asymmetries into account because usually the cycles are assumed to have sinusoidal form. The LombScargle periodogram (Lomb 1976; Scargle 1982) is a commonly used method, but it assumes a strict periodicity, which is generally not the case in stellar activity cycles. The duration of the solar cycle, for instance, varies from about 8 years to 14 years. The use of quasiperiodic models allows the cycles not to be strictly periodic (Olspert et al. 2018). Here, we study each cycle individually to account for cycletocycle differences in the duration and shape of the cycles.
2. Data
2.1. Mount Wilson data
We used the publicly available Ca II H&K Sindex measurements from the Mount Wilson (MW) Observatory, a programme started by Wilson (1978). The data set, including almost 2300 stars, was gathered between 1966 and 1995, with additional data for 35 stars extended to 2001. The Sindex, defined as
is a sensitive indicator of chromospheric magnetic activity (e.g. Egeland et al. 2017). Here H and K indicate flux integrated over narrow passbands centred around the Ca II H and K line cores, and V and R are broad continuum bands on the violet and red sides of the Ca lines. α is a calibration factor that is determined for each night from standard lamp and standard star observations.
Baliunas et al. (1995) determined the periodicity of the MW stars with LombScargle periodograms, and divided the stars with cyclic variations into four different categories based on the falsealarm probability (FAP), the probability that a peak as strong as the observed peak would randomly occur in the LombScargle periodogram, assuming purely Gaussian noise. These categories are labelled “excellent”, “good”, “fair”, and “poor”, corresponding to FAP ≤ 10^{−9}, 10^{−9} < FAP ≤ 10^{−5}, 10^{−5} < FAP ≤ 10^{−2}, and 10^{−2} < FAP ≤ 10^{−1}, expressed in percent, respectively. The authors note, however, that because of variations due to the growth and decay of active regions, for instance, which is nonGaussian noise, the FAP should not be taken too literally.
Olspert et al. (2018) compared the cycle periods in Baliunas et al. (1995) and periods derived with quasiperiodic methods. They found that the results were similar for the “excellent” stars, while the resemblance weakens gradually for the “good”, “fair”, and “poor” stars. Some of the differences, however, can be explained by their use of additional data from the extended 2001 data set, and by the higher significance level.
In our sample we included all the stars defined as “excellent” or “good” by Baliunas et al. (1995), with the exception of HD 78366, HD 201092, and HD 156206, which are labelled “good”. HD 78366 is left out because it is not clear where its minima are because there are multiple secondary minima in the data. HD 201092 was excluded because its minimum around JD2444000 = 2500 is very difficult to define; there seems to be a local maximum where the minimum should be according to the 11.7year cycle reported by Baliunas et al. (1995). They found no secondary shorter cycle in HD 201092, although visual inspection indicates that this would be the case. HD 156206, on the other hand, does not have data to cover any cycle completely (from minimum to minimum). Whithout these stars, our sample consists of 18 stars, all with fairly clear cycles. All the stars in our analysis have also been found to be cyclic by Olspert et al. (2018).
Most of our stars are mainsequence stars, but we also include three giants. The MW database also includes Ca II H&K measurements of the Sun. They were made by measuring the Moon because the lunar spectrum for the H&K lines is effectively just reflected sunlight. Because the Mount Wilson data include only one full cycle for the Sun, we extended our data for the Sun by including Sacramento Peak (SP) Ca II K observations, which were scaled to the same level as the MW Sindex as S_{SP} = 2.61K_{SP} − 0.0647, as was done by Olspert et al. (2018). This combined data set includes three full solar cycles.
The series for the Sun was even further extended back to 1907, including data from solar cycles 15 to 24 by Egeland et al. (2017), who also added Ca II K plage index measurements from the Kodaikanal Observatory in India and calibrated them to the MW scale. However, we only used the data from MW and SP observatories, as was done in Olspert et al. (2018).
2.2. Sunspot numbers
To compare the stellar cycles to the solar cycle, we also analysed sunspot data in addition to the solar chromospheric measurements. We compared the MW+SP data to the classical Wolf sunspot number (WSN)^{1} and to the group sunspot number (GSN), which is recalibrated for different observers with the active day fraction method by Willamo et al. (2017). Reaching as far back as 1610, the sunspot series is much longer than any time series of other active stars. We used the data for sunspot cycles 9–23, from 1843.5 to 2008.9, where multiple of these series are available (MW+SP, WSN, and GSN).
3. Methods
3.1. Defining times of minima and maxima
We defined the times for minima and maxima of the stellar activity cycles individually for each cycle. To define the exact time, we fitted a parabola to the data around the minimum or maximum, and the interval of data included varied depending on the specifics of the cycle. When the cycle was very asymmetric around the minimum or maximum, only a short interval could be used when a symmetric function was fitted, whereas with a poorly covered cycle, a longer interval had to be used to obtain enough data for a reliable fit. The times of minima and maxima defined by this method along with the intervals we used are listed in the appendix (Table A.1). One cycle is then defined as the time between two consecutive minima.
For the dates of minima and maxima for the Sun, we used the minimum and maximum value of the 13month mean value of the sunspot number. This is a commonly used definition of solar minima (see e.g. Hathaway 2015). This is the minimum of the sunspot number cycle, and the chromospheric emission need not be at its minimum at the same time – there are indeed differences of even several years in the timing of the solar minima between different activity indicators, such as the sunspot number, sunspot area, and 10.7 cm radio flux (Hathaway 2015). When the same minima times are used for different solar activity indicators, however, the analysis for the MW cycles of the Sun is comparable to that for the sunspot cycles in Sect. 4.4. For other stars we have only MW data, therefore they are not directly comparable in this sense to the solar cycle.
3.2. Skewness
Skewness is a statistical measure of the asymmetry of a probability distribution, which has been used to measure asymmetries of solar cycles (Ramaswamy 1977; Lantos 2006; Du 2011). The skewness γ, or third moment, of a set of data points x_{i} is defined as
where N is the number of data points, is the sample mean, and σ_{x} is the standard deviation of the sample. A positive skewness indicates a distribution leaning to the left, or in the case of a stellar cycle, a cycle with faster rise time and slower declining time. A negative skewness indicates a leaning to the right, or longer rise time and shorter declining time. A symmetric distribution has γ = 0, although zero skewness does not always mean that the distribution is symmetric. For instance, a distribution with a long and thin tail on the one side and short but thick on the other could also have γ = 0.
In order to calculate the skewness of an activity cycle, the cycle has to be transformed into a onedimensional distribution. We did this by dividing the cycle into ten bins of equal length, where the centre of the bin is at t_{bin}. In the cases when the gaps in the data were too long and some bins would have no data points at all, we reduced the number of bins into the largest number that still included data points in each bin. Then we calculated the mean value of the data points in each bin, and built the final distribution, emulating the cycle, by multiplying this mean value by 10 000 in order to derive an integer value n from data with four decimals, and added n occurrences of t_{bin} to the distribution.
In order to compare the skewness of stellar cycles to the solar cycle, their zerolevels must be comparable. The sunspot cycle approaches zero at solar minimum, but the Sindex of active stars does not. To correct for this, we shifted all the bins of a cycle with a constant value, so that the bin with the lowest value reached S_{min} = 0.001 (corresponds to the t_{bin} appear n = 10 times in the distribution). This was performed similarly for each cycle. An example of this type of distributions emulating the cycles of HD 81809 is shown in Fig. 1. We repeated the same analysis with the same shift of the zerolevel also for the sunspot cycle. For each star, we calculated the skewness for each cycle, and used the average of these cycles as a measure of the average cycle asymmetry for this star.
Fig. 1.
Cycles of HD 81809. The crosses are the original calibrated MW data, and the histograms are the distributions built from these. Vertical lines show the times of minima, dividing the data set into three complete cycles. The zerolevels of the histograms are defined individually for each cycle, but here they are plotted on the same level. The correct individual shifts for each cycle are therefore missing in the visualisation for simplicity. The value on the yaxis, S × 10 000 (which has been shifted in the ydirection), equals the number of data points in a bin, n. 
4. Results
4.1. Rise and decline times of cycles
A simple way to estimate the asymmetry of a cycle is to compare the duration of the rising and the declining phases of the cycle. In the Sun the rising phase is typically shorter.
For each star we calculated the ratio of the average duration for the rising phase ⟨t_{r}⟩ and average duration of the declining phase ⟨t_{d}⟩ of a cycle. Figure 2 shows the relation of this ratio to the average skewness of the stellar cycles. As both are a measure of asymmetry, the almost linear relation is expected. The values of ⟨t_{r}⟩/⟨t_{d}⟩ are also listed in Table 1.
Fig. 2.
Ratio of average rise time and average decline time vs. average skewness. Blue dots represent mainsequence stars, and red dots giants. The Sun is shown in yellow. The continuous line shows the best linear fit. 
Our sample of Mount Wilson stars.
For the calculation of the ⟨t_{r}⟩/⟨t_{d}⟩ parameter for the Sun we used sunspot data for solar cycles 1–23 for better statistics. With any other star, the maximum number of cycles is six.
As the main measure of the cycle asymmetry we used the skewness of the cycle (see Sect. 4.2), but the correlation of the skewness and ⟨t_{r}⟩/⟨t_{d}⟩ confirms that both are usable parameters to measure cycle asymmetries. We calculated a Pearson correlation coefficient r = −0.78 between these two parameters, which indicates a fairly strong negative correlation. Assuming linearity, we derived the relation between them as
The times of minima and maxima are listed for each star in the appendix (Table A.1). When available, we used times of maxima of incomplete cycles to obtain better statistics. For instance, for HD 32147 only one complete cycle is available, but three times of maxima.
4.2. Average skewness of MW cycles
The average skewness ⟨γ⟩ for the cycles of each star and its standard deviation σ for those stars with multiple cycles is shown in Table 1. Figure 3 shows the distribution of the skews of all cycles of all stars. Most cycles (34 of a total number of 47 = 72%) have a positive skew. The peak values are between 0.1 and 0.2. The Sun has a considerably high asymmetry, with a mean skew from MW+SP data of 0.394. Taking all 47 cycles into account, we obtain an average skewness of 0.13, with a standard deviation 0.26.
Fig. 3.
Distribution of skews of all cycles of all stars. 
We compared the average skewness for each star to other stellar parameters; cycle period P_{cyc}, rotation period P_{rot}, effective surface temperature T_{eff}, and activity index , in Figs. 4–7. The figures show the mean value and standard deviation of the variation of γ (and P_{cyc}) for stars with multiple cycles (σ and the standard deviation of P_{cyc} in Table 1). Values for P_{rot} and are from Olspert et al. (2018). The T_{eff} values are from Gaia DR2 (Gaia Collaboration 2016, 2018; Andrae et al. 2018), except for the Sun.
Fig. 4.
Average skewness plotted against P_{cyc}. Blue dots represent mainsequence stars, and red dots giants. The Sun is shown in yellow. The error bars represent the cycletocycle variations for stars with multiple cycles. The vertical line represents γ = 0. 
The three giants and the Sun have a considerably high skewness; the Sun has the highest skewness of the stars in our sample. This is mainly due to the third cycle (solar cycle 23), which is the second most positively skewed cycle of any star in our sample; solar cycles 21 and 22 are much more symmetric. Sunspot data also give a much lower skewness for cycle 23 than MW+SP data (see Sect. 4.4). The value ⟨t_{r}⟩/⟨t_{d}⟩ = 0.625 for the Sun, which was calculated from sunspot cycles 123, is also a very asymmetric value, but two stars have an even larger asymmetry in the rise and decline times. The skewness of the solar cycles from MW+SP data might thus be slightly biased as a result of an overrepresentation of very asymmetric cycles. We recall, however, that the rise and decline times were calculated from sunspot data, which might behave differently than chromospheric data.
We calculated the Pearsons correlation coefficients r between ⟨γ⟩ and the other parameters. These are shown in Table 2 along with their pvalues. ⟨γ⟩ and P_{cyc} or P_{rot} show at best a very weak positive correlation. There might be a slightly stronger positive correlation between ⟨γ⟩ and T_{eff}, but the most relevant is the negative correlation between ⟨γ⟩ and (r = −0.67). The less active stars might thus have more asymmetric cycles in general. This is plausible because young, active stars are known to have more irregular cycles than older, less active stars (Baliunas et al. 1995). Irregular cycles might be skewed in either direction and then be averaged close to symmetric cycles with zero skewness, if enough cycles are included. This correlation is unclear, however. More data are required to be analysed before this can be claimed with some certainty. This differs from the simulated results of Pipin & Kosovichev (2016), who found stronger cycles to be more asymmetric in the regime of weak cycles in their meanfield simulations.
Correlation coefficients between γ and other parameters.
There is much variation in the skewness of the cycles for individual stars, which is seen as high values of σ. Table 1 shows that σ is generally high compared to γ. This is probably not only due to the limited amount of cycles because the sole star with 6 cycles (HD 149661) has the second highest value of σ. The large cycletocycle variations are expected because this is the case in the Sun as well (see Sect. 4.4).
We also compared our values for the average skewness to those of Garg et al. (2019), who studied the same data. This is shown in Fig. 8. There are some large differences in the values. This might be due to the definition of the zerolevel or the binning, which are not described in detail in Garg et al. (2019) because they focused more on the Waldmeier effect than on the cycle skews.
Fig. 8.
Our average skewness for each star compared to that of Garg et al. (2019). The orange line is y = x, where these values would be equal. 
One source of uncertainty is the number of bins, which might affect the skewness of the cycle. In most cases, the data are abundant and regular enough to allow us to divide them into ten bins, but in some stars the number of bins is reduced, in the worst case, to five bins. This is inevitable because these data contain large gaps. The number of bins used for each star is listed in Table 1.
We tested the effect of the binning with the Sun, for which n_{bin} = 10, and average skewness ⟨γ⟩ = 0.394. When the number of bins is reduced to n_{bin} = 8, we obtain ⟨γ⟩ = 0.402, and with n_{bin} = 5, ⟨γ⟩ = 0.329. For stars with poorer data quality, this effect can be expected to be even stronger. It appears to be evident, however, that the stars with n_{bin} = 10 are most comparable to each other.
Because of the gaps in the data, long cycles are probably more reliable than short cycles because the seasonal gaps affect a proportionally smaller part of the cycle, and the shape of the cycle can be identified with more certainty.
4.3. Average cycle shape
We also tried to combine all 47 individual cycles of all stars into an average cycle. The cycle amplitudes were normalised with the same binning as was used in the calculation of γ, so that the lowest bin has the value 0 and the highest bin 1. We added all the data points from individual cycles, scaled between 0 and 1, to the combined cycle without any averaging. The scaling of individual cycles was, however, made with the mean values of the bins to avoid that extreme data points set the scale for the cycles. The cycle duration was normalised to a phase between 0 and 1.
For the resulting average cycle, γ was calculated similarly as for individual cycles, except that now we divided the data to 20 bins instead of 10 because the data are much more abundant. We obtain the value for the skewness as γ = 0.078, which is slightly lower than the average skewness for individual cycles.
We fitted a sinusoid of the form
to the averaged cycle with 20 data points (same as in the binning when we calculated γ). The cosine function has its minimum or maximum at 0, which is defined as the cycle minimum, so that it is forced to the same phase as the cycle. The fitted cosine function, with the fitted parameters a = −0.386 and c = 0.489, is shown in Fig. 9. The fit is plausible, even though individual cycles can be very irregular. Quantitatively, we obtain the chisquared statistics between the data points and the fitted sinusoid as χ^{2} = 11.43, and the pvalue p = 0.91. The cosine curve, however, is not able to take the asymmetry into account because the actual cycle rises to its maxium faster than the cosine. The bestfit cosine also has its maximum at a similar level as the actual average cycle, but its minimum is not as deep. This feature is different than the feature described by Reinhold et al. (2017) for Kepler stars, where the maximum was sharper, and the minimum flatter than the sine curve. They used the amount of photometric variability as a proxy of magnetic activity, however. The variability should be highest around activity maximum, but the details of the cycles might still be different than those found from the Sindex, and furthermore, the span of the Kepler data allow the detection of cycle periods only up to around six years.
Fig. 9.
All cycles phased and normalised together. The blue crosses show the resulting cycle, and the orange curve is a cosine fitted to the cycle. 
In Fig. 9 there might also be an indication of a double peak, as is commonly seen in the Sun, with the Gnevyshev gap in between. The feature is rather weak, however, therefore based on our data, we do not claim the existence of double peaks and the Gnevyshev gap in other stars.
4.4. Comparison to sunspot cycles
The stellar data can be compared to sunspot data. We used the same method to calculate the skewness for monthly values of both the classical WSN and the group sunspot number (GSN) series. The GSN ignores individual spots and only counts the number of spot groups, which reduces observational errors and makes the old observations more reliable. We used the minimum values of the 13month average number of sunspots as the times of solar minima: 1843.5, 1855.9, 1867.2, 1878.9, 1890.2, 1902.0, 1913.5, 1923.6, 1933.7, 1944.1, 1954.3, 1964.8, 1976.2, 1986.7, 1996.3, and 2008.9 (see e.g. Hathaway 2015).
Our values for the skewness of solar cycles 9–23 are shown in Table 3. We also compare our values to the skewness for the WSN published by other authors. Our values agree well with those of Lantos (2006), but not so well with those of Du (2011). The skewness of the GSN is very similar to the skewness of the classical WSN.
Skewness of the solar cycles.
The are some notable differences in the MW data and sunspot data. Especially for cycle 23 do the MW data give a very high skewness of γ = 0.614, whereas sunspot data give γ = 0.282. If the MW cycles for the Sun are not comparable to the sunspot cycle, then cycles for other stars cannot be expected to be directly comparable to the sunspot cycle either.
5. Comparison to simulations
To compare our observational results to numerical simulations, we used the direct numerical magnetohydrodynamic (MHD) simulations of convective dynamos in solarlike stars, described in Viviani et al. (2018), Warnecke (2018), and Warnecke & Käpylä (2019). Some of the simulations, presented in Viviani et al. (2018), are global MHD simulations, ranging between 0.7–1.0 R in the radial direction, and only omitting the polar regions, modelling the star between latitudes −75° to +75°, and the full longitudinal range. A few of the runs in Viviani et al. (2018) and all the runs in Warnecke (2018) and Warnecke & Käpylä (2019) are wedges in the azimuthal direction, covering only the longitudes from 0 to π/2. These are labelled with the superscript “W” in Table 4. A few of the global simulations are run in higher resolution. These are marked with the superscript “a”. The higher resolution runs are slightly more realistic because they are more turbulent than their lower resolution counterparts. In addition to comparing them with the observed MW cycles, we also investigated whether the differing geometry of the simulation setup affects the results.
Skewness of the simulated cycles.
In all of these runs, turbulent convection under the influence of rotation generates differential rotation and largescale dynamo action. As a result, dynamically significant dynamo modes at the system scale are generated and maintained by the flow.
The radial magnetic field at 0.98 R is decomposed into spherical harmonics, where m = 0 mode contains the axisymmetric part of the radial magnetic field, m = 1 is the first nonaxisymmetric mode, m = 2 is the second mode, and so on. We studied the evolution of the dominating dynamo mode in each simulation (found in Table 4 in Viviani et al. 2018), which is m = 0 or m = 1 in all runs. In all the wedge runs m = 0 is the dominating mode, containing most of the magnetic energy on large scales. We note here that a substantial amount of magnetic energy in all runs comes from the smallscale nonaxisymmetric field, but for the comparison with observed cycles, only the largescale magnetic field is relevant.
We chose the runs where cycles for the dominating mode could be defined for a closer study; this includes 20 runs in total. We chose only runs where more than one cycle could be identified in order to obtain some estimate for the cycletocycle variability. We point out that the simulations do not always produce strictly cyclic dynamo solutions, which is likely due to the competition of different dynamo modes in the simulated system. Therefore defining the cycle minima was more challenging from the models than from the MW data. Thus, the results may not be as reliable for the simulated data.
We built the distributions emulating the cycles similarly as for the MW data by multiplying the value of each data point so that we obtained an integer number, and added this many occurrences of this time point to the distribution. We then fitted a parabola around the minimum to define its exact location. Then we calculated the skewness of each cycle similarly as with the MW data.
Figure 10 shows a histogram of the distribution of the skews of the simulated cycles for all cycles together and including only the global runs or only the wedge runs. There is a visible difference between the global and wedge runs: while the histogram including all cycles is centred around zero, with a mean skewness 0.00 and standard deviation 0.32, the histogram including only the global runs has a mean skewness of 0.06 and standard deviation 0.31, and the one including only wedge runs has a mean of −0.06 and standard deviation 0.31. It would thus seem that global simulations produce more positively skewed cycles than wedge runs, although in both cases the cycletocycle variation is large, as it is in real stars as well.
Fig. 10.
Distribution of skews of individual cycles for all the runs (left), only the global runs (centre), and only the wedge runs (right). 
In both global and wedge runs, the deviation (0.31 in both cases) is larger than the difference between them (0.06 − ( − 0.06) = 0.12). To investigate if the difference is significant, we additionally calculated the standard error σ_{⟨γ⟩} of the mean of the distribution:
where n is the sample size. We obtain σ_{⟨γ⟩,global} = 0.04 and σ_{⟨γ⟩,wedge} = 0.03. These are smaller than the difference, which indicates that it is significant and not noise caused by a small sample size. However, we note that for the global runs, a significantly large fraction of the cycles (14 of 48) are from the run K1, which has higher average skewness than most of the runs, and might induce a bias to the result. Nevertheless, the difference between the global and wedge runs, although small, is probably real.
The wedge assumption forces the largescale dynamo to be axisymmetric, whereas in global simulations, nonaxisymmetry is also allowed. These simulations therefore not only allow us to study the cycle asymmetry as a function of rotation or cycle period, but also to study the effect of the degree of nonaxisymmetry on it. By comparing skewness and axisymmetry of global simulations to observational data, it might be deduced whether cycles in real stars are dominated by axisymmetric or nonaxisymmetric modes. Although the parameters in the simulations are still far removed from the real stellar conditions, this may provide a diagnostic tool in the future to further classify observational data into axis and nonaxisymmetric modes.
We note, however, that even the global runs have a lower average skewness than the observed MW cycles. When we assume that the observed chromospheric emission is directly proportional to the magnetic field strength, it is thus plausible that some ingredient is still missing in the simulations, which causes the asymmetry in the observed cycles. The simulations are, for example, still in a parameter regime that is too mildly turbulent, and they do not include realistic photospheres or chromospheres. The other alternative is that cycles are more symmetric for more rapidly rotating stars (for which there is a weak correlation in the MW data). In this case, the different parameter regime of the observations and simulations might explain their difference because the rotation was much faster in most of the simulated runs than in the observed stars.
Similarly to the observational data, we also compared the mean skewness of each run to other stellar parameters. Table 4 shows the mean skewness of all these runs, and the rotation rate of the simulated star, normalised to the solar rotation rate . The rotation rate is transformed into the rotation period by , where P_{⊙} = 26.09 d is the rotational period of the Sun. ⟨γ⟩ is plotted against P_{rot} in Fig. 11, and against P_{cyc} in Fig. 12. Global and wedge simulations are separated by colour in the figures, as are the higher resolution global runs. We calculated Pearson correlation coefficients between ⟨γ⟩ and P_{rot}, and ⟨γ⟩ and P_{cyc} for all simulations together and separately for the global and wedge runs. These are shown in Table 5. The strongest correlation is r = −0.57 for P_{rot} for the global simulations, although this is fairly weak. Moreover, the correlation is positive for the wedge runs. For P_{cyc} the correlations are even weaker. We draw no other conclusions from this, except for the lack of strong correlations between cycle asymmetry and other parameters, as was the case with observed cycles as well.
Fig. 11.
Mean skewness of the simulated cycles as a function of P_{rot}. The error bars represent the standard deviation of the cycles in the run. Global runs are shown in red, with the highresolution runs in orange and the wedge runs in blue. The horizontal line represents γ = 0. 
Correlation coefficients between ⟨γ⟩ and other parameters in the simulated cycles.
We note that the cycle period, which we defined from the times of minima, was determined differently by Viviani et al. (2018). The authors counted how many times the mean magnetic energy level is crossed. The cycle period is also different in Warnecke (2018), who determined the period using power spectra.
It must also be noted that the rotation rate, although the most relevant parameter, is not the only parameter that was varied between the simulations. Other input parameters that were changed between the runs are the grid resolution, the fluid, subgridscale, and magnetic Prandtl numbers, the Taylor number, and the Rayleigh number. We did not analyse, however, how these affect the cycle asymmetry because these parameters are not known for real individual stars.
Table 6 summarizes the main features of the simulated cycles, both including all cycles and when the global and wedge runs are separated. Despite the low number of statistics and the different parameter space between the observations and simulations, the comparison is deemed useful. Firstly, this gives a realistic view on the current state of direct numerical simulations. Furthermore, comparing the observed and simulated trends in skewness may serve as an additional tool for deducing what type of dynamo is operating in a star.
Average skewness and its standard deviation of the observed and simulated cycles.
6. Conclusions
We draw the following conclusions from our study. A fast rise and slower decline is common for stellar activity cycles. The Sun has particularily asymmetric cycles. More active stars might have less asymmetric cycles, but the correlation between the skewness and other parmeters is mostly unclear. Individual cycles might have very irregular shapes, but the average cycle shape is fairly well represented with a sinusoid. The average cycle still reaches its maximum before the sinusoid because it is asymmetrical. The chromospheric and sunspot cycles do not have exactly the same shape. This means that MW cycles for other stars can probably not be directly compared to the sunspot cycle.
The numerically simulated cycles, with shorter rotation periods than the observed real stars, have on average more symmetric cycles, with a distribution in the skewness values centred very close on zero. Perhaps the simulations miss something that makes the cycles asymmetric in real stars. This might indicate that the physics that is still not captured by these models, such as the missing photosphere and chromosphere, is crucial for creating the cycle asymmetries. Other explanations for this might be a difference in the cycles between slow and fast rotators, for which there is some support from the weak correlation between the skewness and the rotation period, and the stronger anticorrelation between the skewness and in the MW data. The simulation geometry affects the asymmetry of the simulated cycles, with the wedge simulations having on average more negatively skewed cycles than the global simulations.
Source: WDCSILSO, Royal Observatory of Belgium, Brussels; available at http://www.sidc.be/silso/datafiles
Acknowledgments
The HK_Project_v1995_NSO data derive 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 HarvardSmithsonian 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. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. TW acknowledges the financial support from the Alfred Kordelin Foundation, and thanks Angie Breimann for our discussion during the BCool meeting in Exeter. TH acknowledges the financial support from the Academy of Finland for the project SOLSTICE (decision No. 324161). MJK and NO acknowledge the support of the Academy of Finland ReSoLVE Centre of Excellence (Grant No. 307411). This project has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (project “UniSDyn”, grant agreement n:o 818665). MV was enrolled in the International Max Planck Research School for Solar System Science at the University of Göttingen.
References
 Andrae, R., Fouesneau, M., Creevey, O., et al. 2018, A&A, 616, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Cox, A. N. 2000, Allen’s Astrophysical Quantities (New York: AIP Press) [Google Scholar]
 Deng, L. H., Xiang, Y. Y., Qu, Z. N., & An, J. M. 2016, AJ, 151, 70 [Google Scholar]
 Du, Z. 2011, Sol. Phys., 273, 231 [Google Scholar]
 Egeland, R., Soon, W., Baliunas, S., et al. 2017, ApJ, 835, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Elling, W., & Schwentek, H. 1992, Sol. Phys., 137, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Feminella, F., & Storini, M. 1997, A&A, 322, 311 [NASA ADS] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [Google Scholar]
 Garg, S., Karak, B. B., Egeland, R., Soon, W., & Baliunas, S. 2019, ApJ, 886, 132 [Google Scholar]
 Gnevyshev, M. N. 1963, Sov. Astron., 7, 311 [Google Scholar]
 Gnevyshev, M. N. 1967, Sol. Phys., 1, 107 [Google Scholar]
 Gnevyshev, M. N. 1977, Sol. Phys., 51, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H. 2015, Liv. Rev. Sol. Phys., 12, 4 [Google Scholar]
 Hathaway, D. H., Wilson, R. M., & Reichmann, E. J. 1994, Sol. Phys., 151, 177 [Google Scholar]
 Lantos, P. 2006, Sol. Phys., 236, 199 [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Newton, H. W., & Milsom, A. S. 1955, MNRAS, 115, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Nikbakhsh, S., Tanskanen, E. I., Käpylä, M. J., & Hackman, T. 2019, A&A, 629, A45 [CrossRef] [EDP Sciences] [Google Scholar]
 Nordemann, D. J. R. 1992, Sol. Phys., 141, 199 [Google Scholar]
 Norton, A. A., & Gallagher, J. C. 2010, Sol. Phys., 261, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Olspert, N., Lehtinen, J. J., Käpylä, M. J., Pelt, J., & Grigorievskiy, A. 2018, A&A, 619, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pipin, V. V., & Kosovichev, A. G. 2016, ApJ, 823, 133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramaswamy, G. 1977, Nature, 265, 713 [Google Scholar]
 Reinhold, T., Cameron, R. H., & Gizon, L. 2017, A&A, 603, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Takalo, J., & Mursula, K. 2018, A&A, 620, A100 [Google Scholar]
 Viviani, M., Warnecke, J., Käpylä, M. J., et al. 2018, A&A, 616, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Volobuev, D. M. 2009, Sol. Phys., 258, 319 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Waldmeier, M. 1935, Astronomische Mitteilungen der Eidgenössischen Sternwarte Zurich, 14, 105 [Google Scholar]
 Waldmeier, M. 1939, Astronomische Mitteilungen der Eidgenössischen Sternwarte Zurich, 14, 470 [Google Scholar]
 Warnecke, J. 2018, A&A, 616, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Warnecke, J., & Käpylä, M. J. 2019, A&A, submitted [arXiv:1910.06776] [Google Scholar]
 Willamo, T., Usoskin, I. G., & Kovaltsov, G. A. 2017, A&A, 601, A109 [Google Scholar]
 Wilson, O. C. 1978, ApJ, 226, 379 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Minima and maxima of individual MW cycles
Times of minima and maxima for the MW stars and the time intervals we used to derive these as upper and lower index, and period and skewness for each cycle.
All Tables
Correlation coefficients between ⟨γ⟩ and other parameters in the simulated cycles.
Average skewness and its standard deviation of the observed and simulated cycles.
Times of minima and maxima for the MW stars and the time intervals we used to derive these as upper and lower index, and period and skewness for each cycle.
All Figures
Fig. 1.
Cycles of HD 81809. The crosses are the original calibrated MW data, and the histograms are the distributions built from these. Vertical lines show the times of minima, dividing the data set into three complete cycles. The zerolevels of the histograms are defined individually for each cycle, but here they are plotted on the same level. The correct individual shifts for each cycle are therefore missing in the visualisation for simplicity. The value on the yaxis, S × 10 000 (which has been shifted in the ydirection), equals the number of data points in a bin, n. 

In the text 
Fig. 2.
Ratio of average rise time and average decline time vs. average skewness. Blue dots represent mainsequence stars, and red dots giants. The Sun is shown in yellow. The continuous line shows the best linear fit. 

In the text 
Fig. 3.
Distribution of skews of all cycles of all stars. 

In the text 
Fig. 4.
Average skewness plotted against P_{cyc}. Blue dots represent mainsequence stars, and red dots giants. The Sun is shown in yellow. The error bars represent the cycletocycle variations for stars with multiple cycles. The vertical line represents γ = 0. 

In the text 
Fig. 5.
Same as Fig. 4, but for P_{rot}. 

In the text 
Fig. 6.
Same as Fig. 4, but for T_{eff}. 

In the text 
Fig. 7.
Same as Fig. 4, but for . 

In the text 
Fig. 8.
Our average skewness for each star compared to that of Garg et al. (2019). The orange line is y = x, where these values would be equal. 

In the text 
Fig. 9.
All cycles phased and normalised together. The blue crosses show the resulting cycle, and the orange curve is a cosine fitted to the cycle. 

In the text 
Fig. 10.
Distribution of skews of individual cycles for all the runs (left), only the global runs (centre), and only the wedge runs (right). 

In the text 
Fig. 11.
Mean skewness of the simulated cycles as a function of P_{rot}. The error bars represent the standard deviation of the cycles in the run. Global runs are shown in red, with the highresolution runs in orange and the wedge runs in blue. The horizontal line represents γ = 0. 

In the text 
Fig. 12.
Same as Fig. 11, but for P_{cyc}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.