Multiple and changing cycles of active stars

Z. Kolláth; K. Oláh

doi:10.1051/0004-6361/200811303

Home

All issues

Volume 501 / No 2 (July II 2009)

A&A, 501 2 (2009) 695-702

Full HTML

Free Access

Issue		A&A Volume 501, Number 2, July II 2009


Page(s)		695 - 702
Section		Stellar atmospheres
DOI		https://doi.org/10.1051/0004-6361/200811303
Published online		05 May 2009

I. Methods of analysis and application to the solar cycles

Z. Kolláth - K. Oláh

Konkoly Observatory of the Hungarian Academy of Sciences, PO Box 67, 1525, Hungary

Received 6 November 2008 / Accepted 1 April 2009

Abstract
Context. Long-term observational data have information on the magnetic cycles of active stars and that of the Sun. The changes in the activity of our central star have basic effects on Earth, such as variations in the global climate, so that understanding the nature of these variations is extremely important.
Aims. The observed variations related to magnetic activity cannot be treated as stationary periodic variations, therefore methods like Fourier transform or different versions of periodograms only give partial information on the nature of the light variability. We demonstrate that time-frequency distributions provide useful tools for analysing the observations of active stars.
Methods. We tested and used different methods, such as short-term Fourier transform, wavelet, and generalised time-frequency distributions, for analysing temporal variations in timescales of observational data.
Results. With test data we demonstrate that the observational noise has practically no effect on the determination in the long-term changes of time-series observations of active stars. The rotational signal may modify the determined cycles, therefore it is advisable to remove it from the data. Wavelets are less powerful in recovering complex long-term changes than other distributions that are discussed. By applying our technique to the sunspot data we find a complicated, multi-scale evolution in the solar activity.

Key words: Sun: activity - methods: data analysis

1 Introduction

It is well known that solar cycles cannot be represented by single periodic variations. Besides the 11 year quasi-periodicity (Schwabe 1843) of solar activity, lots of different rhythms have been discovered or at least suggested by different studies (see e.g. Frick et al. 1997; Forgács-Dajka & Borkovits 2006). One remarkable and well-studied periodicity of the solar variation was found by Gleissberg (1966) about 70 years ago. He described this variability as a ``long-periodic fluctuation of the sun-spot numbers, the period being about seven or eight sun-spot cycles''. Frick et al. (1997) derived the length of the Gleissberg cycle as $\approx$ 100 years from wavelet analysis of a 384 years long dataset (Hoyt et al. 1993). The non-stationarity in the long-term variation of active stars has already been found in active stars as well (Oláh et al. 2007).

Time-frequency representation of a non-stationary signal yields information about characteristics of the dataset in the time-frequency plane. Then transient events or frequency/amplitude changes can be easily monitored with these tools. In the literature, a frequently used method is the wavelet-map and some other techniques based on wavelets. However, wavelet is only one of the several time-frequency representations, and other methods may give better results than wavelets. The key problem is the balance between temporal and frequency resolution. The frequency-time uncertainty principle gives strict constraints for the time duration ( $\Delta t$ ) and the bandwidth ( $\Delta \omega = 2\pi \Delta f$ ) of a signal: $\Delta t * \Delta \omega \ge 1/2$ . These measures of broadness are determined by the standard deviation of frequency and time, and the uncertainty principle with these quantities applies to any time-frequency distribution. The linear transforms (such as wavelet or windowed Fourier-transform) have a limit for joint resolution of time and frequency for local structures; for example, in general it is not possible to determine the frequency precisely at high temporal resolution. Bilinear time-frequency transformation, called the kernel method (Cohen 1995) can break the above principle locally, but its cost is paid by artifacts at other frequencies.

While wavelets and other methods based on Fourier transforms are exclusively used to analyse stellar activity, more sophisticated methods have been available for decades: the Wigner distribution (Wigner 1932) and its extensions by Cohen (1995). These methods have also been proven to be useful in analysing of variable star data of chaotic origin or non-stationary behaviour (see Buchler et al. 1995; Kolláth & Buchler 2001). In this paper we present a method that can be used to study changes of stellar activity cycles in time using data with yearly gaps and with short-term signals (rotational modulation).

As example we reanalysed the available sunspot and radio records with modern time-frequency methods to give a complete picture of the beating of the Sun in the period range from a few years to 200 years. In the companion paper (Oláh et al. 2009, Paper II) we apply the method described in the present paper to 20 active stars, and discuss the time variability of their cycles.

2 Time-frequency distributions

A simple but still powerful method is the ``short-term Fourier transform'' (STFT hereafter). In this case the dataset is multiplied by a Gaussian window centred on a given epoch t₀. If this windowed part of the signal is then Fourier-transformed, the frequency spectrum of the data around the time t₀ can be generated. By shifting the centre of the Gaussian window in time, a sequence of Fourier spectra is obtained, which gives the time-frequency representation of the data. If the sampling of t₀ is dense enough, then a two-dimensional map of the energy distribution on the time-frequency plane can be plotted.

We demonstrate this method on the yearly sunspot numbers - a dataset that is well known and analysed by several authors. Figure 1 displays the STFT of the solar data. The Fourier spectra of the individually windowed subsets of the data are also plotted with thick lines on the time-frequency distribution map. By changing the width of the Gaussian window, the balance between time and frequency resolution can be modified. All the locations on the map, of course, satisfy the uncertainty principle $\Delta t * \Delta \omega \ge 1/2$ , and the ratio of time and frequency resolution is constant on the map.

$\begin{figure} \par\includegraphics[width=9cm,clip]{1303f1} \end{figure}$	Figure 1: Short-term Fourier transform ( bottom) of the sunspot number time series ( top). Fourier-transforms of the subsets of the dataset are also plotted. The colour scale of relative variations in the distributions is presented in Fig. 7.
Open with DEXTER

While in this demonstration we used a time-centred Gaussian window to plot the curves of Fourier transforms, in the implementation of STFT, we use windowing in Fourier space:

$\begin{displaymath}T(t,f) = 2 \int_0^{+\infty} F(\nu) \exp \left( -{1\over 2} ... ...r \sigma^2(f)} \right) \exp ({\rm i}2\pi t \nu) {\rm d}\nu , \end{displaymath}$

(1)

where $F(\nu)$ denotes the Fourier transform of the signal. The inverse Fourier transform is only calculated for the positive frequencies. This method has several advantages. At low frequencies it reduces the distortion due to the components at negative frequencies; moreover, when $\sigma$ is large, the transform reduces to the analytic function (see Buchler & Kolláth 2001). Depending on the functional form of $\sigma(f)$ , Eq. (1) reduces to STFT or a wavelet. With $\sigma=f_0/\alpha$ , where f₀ is a constant frequency, the equivalent temporal window function of STFT is given by

$\begin{displaymath}h(\tau) = \exp \left( {-\tau^2 f_0^2 \over 2 \alpha^2} \right)\cdot \end{displaymath}$

(2)

Usually it is a good choice to set f₀ to the maximum or the central frequency of the studied period domain. Then $\alpha \sim 1$ provides a naturally well balanced temporal and frequency resolution. It can be seen from Eq. (2) that the half width of the time-centred windowing function is $\alpha/f_0$ . Then, if $\sigma$ depends on frequency, the width of the window also varies with frequency. Especially if $\sigma(f)= f/\alpha$ then T(t,f) is equal to wavelet, with the Morlet kernel, since the width of the temporal windowing function decreases with 1/f.

$\begin{figure} \par\includegraphics[width=8.0cm,clip]{1303f2a}\vspace{4mm} \in... ...303f2b}\vspace{4mm} \includegraphics[width=8.0cm,clip]{1303f2c} \end{figure}$	Figure 2: Wavelet ( top and middle) and short-term Fourier transform (STFT) distribution ( bottom) of the yearly sunspot numbers.
Open with DEXTER

In the standard Morlet wavelet $\alpha$ is fixed to 1. The introduction of a free parameter ( $\alpha$ ) that controls the balance between frequency and temporal resolution has several advantages as already described by Kolláth & Szeidl (1993) and Soon et al. (1999).

As can be seen from the above definitions, a wavelet with Morlet-kernel is very similar to STFT; however, for wavelets the balance between temporal and frequency resolution depends on the frequency. This feature is demonstrated in the top panels of Fig. 2. The parameter (f₀) is defined in such a way, that for a given $\alpha$ , the wavelet and STFT have the same frequency/time resolution at the frequency f₀. The top and middle panels of Fig. 2 show wavelets where the lower frequency part (top) and the high-frequency part (middle) was matched to the STFT, respectively. Resolution in time increases with frequency, since the same number of cycles is covered on a shorter timebase, at higher frequencies. The disadvantage of this representation is that harmonic structures, which display the same characteristics in STFT, are distorted. Figure 2 thus demonstrates that, while the time-frequency distribution in the lower panel clearly show the synchronised variation of the $\approx$ 11 year cycle with its double frequency harmonic, the wavelet destroys this parallel structure.

We used a different scaling of the amplitudes at different frequencies. In time-frequency distributions the lower amplitude parts are easily washed out by the contributions from power of other components. Far from our topics but similar in mathematical tools, it is well known in speech analysis that the obscured parts of the frequency distribution can be recovered by applying a pre-emphasis filter. A multilevel pre-emphasis Fourier filter is applied to visualise the whole studied period range equally. For the solar data relative to the highest amplitude $\approx$ 11 year periodicity (0.07-0.14 c/y), the signal in different frequency ranges were amplified by the following factors ( $R_{\rm e}$ ): $R_{\rm e}=1.75$ for 0.0 < f < 0.05, and $R_{\rm e}=3.5$ above f=0.16 (f is given in c/y). A smooth change in $R_{\rm e}$ is used in the transition regions. With this pre-emphasis filtering there is a factor of 3.5 difference in the amplitudes to get similar grey levels around 11 and 5.5 year periods.

The optimal balance of temporal and frequency resolution depends on the signal itself. A question arises as to whether it is possible to construct a distribution, where the resolution in both coordinates is optimal, but some other disadvantage is allowed. The answer is more general group of methods for signal processing introduced by Cohen (see e.g. Cohen 1995). The generalised time-frequency distribution is given by

C(t,f)	=	$\displaystyle {1\over 2\pi} \int\int\int \exp(-{\rm i}\xi t - 2\pi~ {\rm i}\tau~ f - {\rm i}\xi\theta)$
		$\displaystyle \times \Phi(\xi,\tau) s^*(t-\tau/2) s(t+\tau/2) {\rm d}\theta {\rm d}\tau {\rm d}\xi,$	(3)

where s(t) is the analysed time series and $\Phi(\xi,\tau)$ the kernel of the distribution that determines the specific properties of the distribution. The simplest time-frequency distribution is the Wigner-Ville transformation with $\Phi(\xi,\tau)=1$ . This distribution only works well for special datasets, otherwise the resulting map is heavily contaminated by cross terms of the different components. A bi-Gaussian kernel provides a distribution that, in the limit of the kernel width, is equivalent to STFT. This so-called pseudo-Wigner distribution (PWD hereafter) is widely used in sound processing. For comparison purposes, we also use the Choi-Williams distribution (CWD) (see Choi & Williams 1989), which applies an exponential kernel.

The top panel of Fig. 6 shows the PWD of the yearly sunspot numbers. It can be seen that the resolution in both coordinates is improved compared to the result with STFT displayed in Fig. 2 (bottom panel); however, there are artificial features visible on the PWD distribution that have no real meaning in the observed signal. The reason for these cross terms is the non-linearity of the method; e.g., if two sinusoidal waves with different frequencies coexist in the signal, then false power arises between the two ridges representing the cycles.

3 Preprocessing observational data

Stellar observations in reality have lots of deficiencies, like observational noise and imperfect sampling. Then time-frequency methods cannot be applied directly to the observational data, and the analyses should be preceded by a sequence of preprocessing steps.

In the case of active stars, the situation is even more complicated, because different timescales of variability coexist, which interfere with the sampling in different ways. The rotational periods fall into the range of days to weeks, very close to the order of the sampling time. In contrast, long-term changes have a characteristic time above a year, so averaging and resampling on this level is needed to obtain information on these components. Thus, the preprocessing of the data requires extreme care.

To avoid false averages due to the undersampling of the rotational component, it is advisable to whiten all observations with the rotational period. Because the rotational component is modulated by e.g. the differential rotation, prewhitening cannot be performed simultaneously for the whole dataset.

The main problem is that the light-curves are unevenly sampled, and they contain seasonal gaps. Although methods exist to calculate e.g. wavelets of this kind of data directly (like the adaptive wavelet), according to our previous experiments we prefer to interpolate the data when it is possible. When the sampling is denser than the dominant periodicities in the data, then a simple moving averaging can provide a continuous signal. Next, the averaged data can be smoothed with a cubic smoothing spline (Reinsch 1967). This operation provides a continuous function S(t) based on the discrete time series ( $\{t_i,x_i\}$ i=1,n) by minimising the integral of S''(t)², with the following constraint:

$\begin{displaymath}{1\over n} \sum_{i=1}^n (x_i-S(t_i))^2 \le \sigma_{\rm s}, \end{displaymath}$

(4)

where $\sigma_{\rm s}$ , the smoothing factor, is a free parameter, which controls the smoothness of the spline. In theory its value should be chosen in the confidence interval corresponding to the left side of Eq. (4), thus $\sigma_{\rm s}$ should be of the same order of amplitude as the standard deviation of the observational noise. In practice, we usually manually check the resulting spline-smoothed curves to fine tune of the averaging and spline smoothing. As a final step, S(t) is resampled to generate an equally spaced time series.

4 Testing the methods

4.1 Artificial test data sets

We created artificial datasets to test the methods for a general signal representing the cycles of active stars. Signal ``I'' consists of a periodic signal with a constant frequency of 0.002 c/d, a component with linearly increasing frequency (chirp) with a modulation between f=0.0007 and f=0.0011 c/d, and a wave packet at f=0.0003 c/d with a Gaussian amplitude modulation: $\exp(-(t-2500)^2/2000^2)$ ). The amplitudes of the constant f=0.002 c/d signal, the chirp, and the wave packet are 0.02, 0.04, and 0.04, repectively.

In addition to the components of signal ``I'', signal ``II'' has a higher frequency variation, demonstrating the rotational period of active stars, and a uniformly distributed noise with amplitude comparable to the rotational component. To test the preprocessing and the time-frequency methods together, the test signal ``II'' was finalised by resampling the data with the time series of the observations of EI Eri (see Paper II) The data were averaged with a time base of 90 days, then the averages were spline smoothed with $\sigma_{\rm s} = 0.005$ and resampled with $\Delta t = 10$ days.

4.2 Observational noise and time-frequency distributions

It is well known that the Fourier transform has a very good noise rejection property in the case of periodic or quasi-periodic signals. The same is true for time-frequency distributions. If the observational noise and the signal are not correlated, then the method provides a good result even in the case of a large amount of noise. A nice example of this nature of time-frequency distributions is the analysis of the semi-regular star V CVn (Buchler et al. 2004), in which two datasets are available for the same period of time. One of the records consists of visual magnitude estimates by amateur astronomers, while the other data are professional measurements by photoelectric photometry. It is clear that even the statistical properties of the noise sources are different in the two cases. Figure 2 in Buchler et al. (2004) clearly demonstrates that the time-frequency plots are almost fully independent of the choice of the dataset, so the effect of observational noise is negligible for signals with good SNR and sampling.

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{1303f3a}\hspace*{5mm} \in... ...{5mm} \includegraphics[width=7.5cm,clip]{1303f3f} \vspace*{0.5mm} \end{figure}$

Figure 3:

Time-frequency distribution of test data. Panel A): PWD of signal ``I''; panel B): PWD of signal ``I'' with white noise; panel C): PWD of signal ``II'' with averaging and spline interpolation; panel D): PWD of signal ``II'' after Fourier filtering of the rotational period in the independent observational seasons and then averaging and spline smoothing; panel E): STFT of the same signal as in panel D); panel F): CWD of the same signal as in panel D).

Open with DEXTER

We demonstrate this noise rejection property of time-frequency distributions on our first test time series. The top panels of Fig. 3 display the PWD of signal ``I'' (A) and the same signal modified with noise, with large relative temporal variance, added to it (B). The amplitude of the noise is seen in the top panel of the figure. While the noise has practically no effect on the lower frequencies, the f=0.002 c/d component is slightly distorted. The results with signal ``I'' are presented only for PWD, but we note that the other time-frequency methods (STFT, wavelet, CWD) have the same properties.

We have the sunspot observations for the sensitivity to noise, too. White or Gaussian noise with the same amplitude as the maximum sunspot number does not significantly distort the time-frequency plots.

4.3 Effect of uneven sampling and noise

Because of the averaging and spline interpolation, the effect of observational noise is mixed with the effects of sampling, so we performed our tests with ill-sampled data. The resulting time-frequency distributions obtained from the gapped signal ``II'' are displayed in the lower panels of Fig. 3, where, on panel C the rotational signal is included, while it is filtered out from the data in panels D, E, and F. Averaging and smoothing spline interpolation introduces a low-pass filtering of the signal, which can be compensated for by the pre-emphasis filtering. In these cases we selected $R_{\rm e}=4.0$ for f> 0.0015. All the components above about 2 years of characteristic time are recovered well with these settings and processing, especially in the case where the rotation period is filtered out. Traces of the F=0.002 c/d (P=1.37 yrs) component can be found, but it is strongly distorted. The reason for this distortion is that these periodicities cannot be exactly recovered because of the sampling. Except for this deficiency, the main time scales of the ill-sampled data can be recovered by the methods described above. This result gives a warning that, because of the gapped nature of stellar data, activity cycles around 1-2 years cannot be easily recovered, if at all.

The above conclusions are valid for all the time-frequency distributions we tested. We display the results for the filtered signal ``II'' for three different transformations. The main components are recovered with all distributions. The cross terms are clearly visible in PWD and CWD (especially in CWD). Closely spaced frequencies, together with their cross terms, tend to form loops (or ``bubbles''), like the ones visible in all PWD figures between $t \approx 2000$ and 4000 days at low frequencies. These are well-defined signatures of structures with two different frequencies, even short-lived ones, and their appearance remains mostly invariant against changes in the preprocessing (e.g. smoothing). STFT with our test data clearly displays the same double frequency structure, however, according to our experience, PWD is superior in some cases to STFT in this property.

To give a quantitative test on the combined effect of sampling and noise in a statistical sense, we performed Monte-Carlo simulation of the effect of observational noise on signal ``II''. We calculated 10⁴ realisations of the added Gaussian noise sequence with different amplitude (standard deviation). Then for each realisation we calculated the same procedure as for the original test data. The time-frequency distribution of the difference of the noisy and the noiseless smoothed data was obtained and the statistics of the amplitude of the largest feature (peak or ridge) was estimated. Figure 4 shows the false-alarm probability FAP(A), i.e., the probability that a structure with amplitude larger than A appears in the transform due to observational noise. This test was performed with different standard deviations ( $\sigma_n$ ) of the added Gaussian noise and for all the time-frequency distributions we use. In the figure from left to right the curves belong to $\sigma_n= 0.0025$ , 0.005, 0.01, 0.02, 0.04, and 0.08 mag, representing a wide range of observational noise. The different distributions behave similarly in this test, with only slight differences. The results with STFT (black) and CWD (grey) are displayed in the figure, the curves obtained with PWD are located between the two plotted ones. It is clear that the location of the curves are scaled linearly with the amplitude of the added noise. Only very high noise, above the displayed ones, distorts this trend, as the spline interpolation becomes unstable. The amplitude of the observed features of the distributions of the original signal is displayed by a black rectangle in Fig. 4. It is clearly demonstrated that even a lot of noise is negligible compared to the real signal in the time-frequency distribution. We have to note that this result depends on the sampling and the signal content of the observations, so it should be repeated for datasets with different characteristics. The nature of noise also influences the above results. We performed our tests with a broad band (white) noise. For correlated, e.g. colour noise, with the same power, false structures appear with higher probability within the frequency range of the noise. Since the power density increases linearly with the inverse of the bandwidth of the noise (with the same total power), a similar shift is expected in the FAP(a) curves. However, as our experience suggests, the effect of a realistic noise turns out to be negligible in most applications.

The most probable time-frequency distribution due to observational noise can be calculated during the Monte-Carlo simulation. For pure Gaussian noise, the expected value of the time-frequency distribution $\langle T(t,f)\rangle$ is the same for all time-frequency pairs. The averaging and spline smoothing interpolation, however, acts as a low-frequency filter, so $\langle T(t,f)\rangle$ depends on f. Because of the uneven sampling, this filtering also depends on the time. Therefore, $\langle T(t,f)\rangle$ of the noise components can be used as a combined time-frequency transfer function of the preprocessing of the data. Figure 5 displays this transfer function for the added noise of $\sigma_n=0.04$ for STFT and PWD. With PWD the most probable range of spurious components is narrower in time than that with STFT. This test confirms the above findings that the amplitude of the signal is decreased in the high frequency part of the distribution (around 1 year). For times when the sampling is very poor, in addition, an amplification is possible at low frequencies. This amplification is located on the time-frequency plane at periods that are comparable to the length of the longest gaps. However, according to our experience, this effect can be drastically reduced with a thorough interactive control of the parameters. In the Monte-Carlo test, of course, it is not possible to do the same fine tuning, and the structures are visible as result in Fig. 5.

$\begin{figure} \par\includegraphics[width=9cm,clip]{1303f4} \end{figure}$	Figure 4: False-alarm probability that a structure appears in STFT (black) and CWD (grey) due to observational noise. The amplitude of the added Gaussian noise ( from left to right) are 0.0025, 0.005, 0.01, 0.02, 0.04, and 0.08. The black rectangle represents the amplitudes of the visible real structures in Fig. 3.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=9cm,clip]{1303f5a}\vspace*{3mm} \includegraphics[width=9cm,clip]{1303f5b} \end{figure}$	Figure 5: Time-frequency transfer function of the processing. Top: STFT; bottom: PWD.
Open with DEXTER

5 Application to solar cycles

Long-term solar records have been analysed by several authors to date. A review of solar cycle evolution is presented by Usoskin & Mursula (2003), with an impressive list of references on the subject. The long-term behaviour of the sunspot group numbers have been analysed using wavelet technique by Frick et al. (1997) who plotted the changes of the Schwabe cycle (length and strength) and studied the grand minima. However, there are still features remaining that can be revealed by the use of time-frequency methods. Here we present some of these characteristics from sunspot numbers and solar radio data.

5.1 Variations related to the Schwabe cycle

$\begin{figure} \par\includegraphics[width=7.95cm,clip]{1303f6a}\vspace{3mm} \i... ...1303f6b}\vspace{3mm} \includegraphics[width=7.95cm,clip]{1303f6c} \end{figure}$	Figure 6: PWD of the yearly sunspot number time series ( top) and the same of the datasets generated from the January ( middle) and July ( bottom) monthly averages.
Open with DEXTER

It is well known that the sunspot number is not the best physical indicator of solar activity; however, it is the only direct measure we have for a relatively long time span. Moreover, the sunspot number has an intrinsic fluctuation because of the arbitrariness of the number of visible spots and groups. If we assume that the averages of e.g. the January data are of independent measures those of the July data (i.e., the ``noise'' for the averages of a given month are mostly independent from that of the month half a year before and after), then it is possible to construct two independent datasets: one for the January and one for the July data. The comparison of the results for the two sets gives an idea of the errors in the time-frequency distributions. Figure 6 (middle and bottom panels) exhibits the time-frequency distributions for periods longer than 2.5-yr for both data subsets. The PWD of the two subsets of the data is hardly distinguishable for periods longer than $\approx$ 4 years.

In Fig. 6 (top) we present the PWD of the whole yearly averaged sunspot data. The STFT (Fig. 2, bottom) shows basically the same structure as the PWD, but with less visibility and resolution. Also, because the STFT does not contain any false cross terms, it validates the higher resolution results of the PWD.

Both the Schwabe and the Gleissberg cycles can be easily identified in Fig. 6. In the following we denote the instantaneous frequencies of these cycles by $f_{\rm S}$ and $f_{\rm G}$ , respectively. The modulation of the Schwabe cycle is well known. The period of the main variation in solar activity fluctuates between 9 and 13 years, which is clear in our results (Figs. 2 and 6). However, the previous investigations neglected to check the harmonics of frequency $f_{\rm S}$ . Since the solar cycle is not a sinusoidal variation, one should find any power at $2f_{\rm S}$ and possibly at $3f_{\rm S}$ . Since the harmonic frequencies are synchronised to $f_{\rm S}$ , investigating these harmonics leads to an independent estimate of the variation in the cycle length. The inspection of the history of the $2f_{\rm S}$ component confirms that around 1780, the cycle had a period of about 9 years, then it increased to 11 years after the minimal activity.

Till the middle of the 20th century, the period decreased to 10 years, then it made a swing to a longer period. It is noticeable that the amplitude of the $2f_{\rm S}$ component is very small when the Schwabe period changes. It looks as if there were some locked phases, with almost constant period, when the power at $2f_{\rm S}$ is high. Another important feature related to the nature of the $2f_{\rm S}$ component is that its amplitude is increased for the solar cycles with high amplitude more than expected from the ratio of the maximum levels. Since the even harmonics indicate the asymmetry of the form of a solar cycle (symmetric waveforms only have odd harmonics), another indication of the Waldmeier effect and related phenomena (see e.g. Cameron & Schüssler 2008).

The ``bubble'' like structures in time-frequency plots indicate that two frequencies coexist close to each other for a finite time interval. The power is clearly split into two frequency parts after 1950, and it is most visible around 1970. The frequency splitting is the same at all three places ( $\delta f=0.025{-}0.03~c/y$ ), which indicates a nonlinear connection between different periodicities. The splitting is also visible in the other time-frequency distributions (STFT, CWD), but most prominent in PWD as already mentioned in Sect. 4.3. As an additional test we calculated the time-frequency distribution from the 10.7 cm radio flux of the Sun, which is available from 1947, and is a good proxy for its activity. The solar radio dataset was recorded at the Dominion Radio Astrophysical Observatory (DRAO) at Penticton, Canada (Tapping & Charrois 1994), three times a day, and are given in solar flux units or Janskys. The result is plotted in Fig. 7, and it is clear that the same structure is found in this diagram like in Fig. 6 from sunspot numbers; however, we cannot resolve the Gleissberg cycle from the radio data because of the short time base. The possible frequencies appearing in the splitting of the power are indicated by horizontal dotted lines in Figs. 6 and 7. This test strongly confirms the physical origin of the structure that we have uncovered. To unfold all the regularly spaced frequencies ( $f_{\rm G}$ , $f_{\rm G}+\delta f$ , $k*f_{\rm S}$ and $k*f_{\rm S}+\delta f$ , k=1, 2), at least one more frequency is needed in addition to the well-known ones ( $f_{\rm G}$ and $f_{\rm S}$ , i.e., the Gleissberg and the Schwabe frequencies, respectively).

$\begin{figure} \par\includegraphics[width=5.2cm,clip]{1303f7a}\hspace*{4.8mm} \includegraphics[width=3.2cm,clip]{1303f7b} \end{figure}$	Figure 7: PWD of the monthly solar radio flux data.
Open with DEXTER

5.2 Long-term variations in the solar cycle

$\begin{figure} \par\includegraphics[width=12cm,clip]{1303f8} \end{figure}$	Figure 8: PWD of the long term envelope of the sunpot data.
Open with DEXTER

The temporal evolution of the Gleissberg cycle can also be seen in the time-frequency distribution of the solar data. The Gleissberg cycle has been found to be as variable as the Schwabe cycle. It has two higher amplitude occurrences: first around 1800 (during the Dalton minimum), and then around 1950. Very interesting behaviour is the continuous decrease in the frequency (increase of period). While near 1750 the cycle length was about 50-yr, it lengthened to approximately 130-yr by 1950. This period variation explains why different works give different periods for the Gleissberg cycle. This systematic change in the frequency does not support the conclusion of Lawrence et al. (1995) that the power at long periods is a fingerprint of a period-doubling cascade. The regular pattern of frequency splitting after 1950 gives the possibility of estimating the value of the Gleissberg period for the last decades, and $\delta f = 0.03~c/y$ can be precisely determined from the splitting at $f_{\rm S}$ and $2f_{\rm S}$ . Similarly, $f_{\rm G}+\delta f = 0.031{-}0.032~c/y$ is well-defined. From these two ranges ( $\delta f=0.025{-}0.03$ and $f_{\rm G}+\delta f=0.031{-}0.032$ , so $f_{\rm G} < 0.007~c/y$ ), we suggest that the period of Gleissberg cycle has already reached a level above 140 years. This result agrees with the rate of frequency decrease from 1800 to 1950.

We extended the Schove (1955) series that lists activity maxima, minima, and amplitudes, based mostly on aurora records before 1610 and using optical records afterwards, by the envelope of the recent sunspot data. The dates and values of the maxima were spline-smoothed and interpolated to get an equally sampled envelope curve. The last part of the time-frequency distribution of the envelope displays the same structure as the low-frequency part of the PWD of the recent data (Fig. 6). Figure 8 exhibits the smoothed Schove series together with the variation in the Gleissberg cycle for 500 years.

The Gleissberg cycle had a long period during the Maunder minimum. After 1700 the period jumped to a lower value (50 years) and started a slow increase. If we accept that the Gleissberg cycle gives the primary indication of the long term behaviour of sunspot numbers, we can interpret the 500-year-long history of solar activity as follows. The long period of the Gleissberg cycle is a good indicator of an extended Maunder minimum, while the subsequent, sudden increase in the Gleissberg period (after 1700) coincides with its termination. An increase in the amplitude of the Gleissberg cycle around 1800, when its period was relatively short, resulted in a shorter (Dalton) minimum. During the last century, the Gleissberg period became long, but now indicates a long lasting active stage of the Sun. Can we derive any prediction about the possible future of solar cycles from the extension of this history? The answer depends on the nature of the newly found splitting of the frequencies. If it is a sign of the restart of the Gleissberg cycle at a shorter period, then we can expect faster changes in the long-term activity level, perhaps a sudden termination of the high activity period. If the 30 year cycle (i.e., $f_{\rm G}+\delta f \approx 0.032$ , $\sim$ 30-yr) is just a temporary one, then a very slow decrease in activity level is expected. To rule out or confirm any of the two possibilities we have to wait for further observations or construct a physical model that is able to produce all the observed behaviour and unfold the nature and connection of all the different periodicities.

Unfortunately, at present we are very far from the construction of a proper physical model of the Sun, if it is at all possible. Bushby & Tobias (2007) discuss the possibilities of predicting the solar cycle through mean-field models, in which the cycle modulations originated from stochastic or deterministic processes. They found that good prediction even of the nearest cycles is impossible in both cases. However, time-series analysis techniques, which find a mapping function for the data, may give better results at least for a short-term prediction.

6 Conclusions

It is demonstrated that time-frequency distributions are useful tools for analysing nonstationary time series, because the short-term Fourier transform, the Wigner distribution, and its extensions are more useful than wavelets. Long-term changes (cycles)of active stars can be correctly identified with the ill-sampled observational data after careful preprocessing. We found an extremely complicated multi-scale evolution in the solar activity. All the observed features in the time-frequency history of the Sun should provide strong constraints on modelling the solar magnetism.

Whether it is possible to predict the long-scale future of solar activity from the similarities of the structures that have emerged in the last decades and those in the early history of the Sun is uncertain. Its verification would be an important future project both in theory and data analysis, because of the probable connection between the long-term changes in our climate and solar activity (see Velasco & Mendoza 2007, and references therein, for more).

Acknowledgements

We are grateful to A. F. Lanza and W. Soon for useful discussions during the preparation of the manuscript. Our thanks go to the referee, Dr. F. Baudin, whose suggestions helped to improve this paper, as well as the companion paper. K.O. acknowledges support from the Hungarian Research Grants OTKA T-048961 and T-068626.

References

Buchler, J. R., & Kolláth, Z. 2001, in Nonlinear Studies of Stellar Pulsation, ed. M. Takeuti, & D. Sasselov, Atrophysics and Space Science Library, 47, 185 (In the text)
Buchler, J. R., Kolláth, Z., Serre, T., & Mattei, J. 1995, ApJ, 462, 489 [NASA ADS] [CrossRef] (In the text)
Buchler, J. R., Kolláth, Z., & Cadmus, R. R. 2004, ApJ, 613, 532 [NASA ADS] [CrossRef] (In the text)
Bushby, P. J., & Tobias, S. M. 2007, ApJ, 661, 1289 [NASA ADS] [CrossRef] (In the text)
Cameron, R., & Schüssler, M. 2008, ApJ, 685, 1291 [NASA ADS] [CrossRef] (In the text)
Choi, H., & Williams, W. J. 1989, IEEE. Trans. Acoustics, Speech, Signal Processing, 37, 862 (In the text)
Cohen, L. 1995, Time-Frequency Analysis (Englewood Cliffs: Prentice-Hall) (In the text)
Forgács-Dajka, E., & Borkovits, T. 2006, MNRAS, 374, 282 [NASA ADS] [CrossRef] (In the text)
Frick, P., Galyagin, D., Hoyt, D. V., et al. 1997, A&A, 328, 670 [NASA ADS] (In the text)
Gleissberg, W. 1939, The Observatory, 62, 158 [NASA ADS] (In the text)
Hoyt, D. V., Schatten, K. H., & Nesme-Ribes, E. 1994, Geophys. Res. Lett., 21, 2067 [NASA ADS] [CrossRef] (In the text)
Kolláth, Z., & Buchler, J. R. 2001, in Nonlinear Studies of Stellar Pulsation, ed. M. Takeuti, & D. Sasselov, Astrophysics and Space Science Library, 47, 29 (In the text)
Kolláth, Z., & Szeidl, B. 1993, A&A, 277, 62 [NASA ADS] (In the text)
Lawrence, J. K., Cadavid, A. C., & Ruzmaikin, A. A. 1995, ApJ, 455, 366 [NASA ADS] [CrossRef] (In the text)
Oláh, K., Strassmeier, K. G., Granzer, T., Soon, W., & Baliunas, S. L. 2007, AN, 328, 1072 (In the text)
Oláh, K. Kolláth, Z., Granzer, T., et al. 2009, A&A, 501, 703 [NASA ADS] [CrossRef] [EDP Sciences] (Paper II)
Reinsch, C. H. 1967, Numer. Math., 10, 177 [CrossRef] (In the text)
Schove, D. J. 1955, J. Geophys Research, 60, 127 [NASA ADS] [CrossRef] (In the text)
Schwabe, H. 1843, reprinted in Early Solar Physics, ed. A. J. Meadows (London: Pergamon), 95 (In the text)
Soon, W., Frick, P., & Baliunas, S. 1999, ApJ, 510, L35 [CrossRef] (In the text)
Tapping, K. F., & Charrois, D. P. 1994, Sol. Phys., 150, 305 [NASA ADS] [CrossRef] (In the text)
Usoskin, I. G., & Mursula, K. 2003, Sol. Phys., 218, 319 [NASA ADS] [CrossRef] (In the text)
Velasco, V. M., & Mendoza, B. 2007, Adv. Space Res., 42, 866 [NASA ADS] [CrossRef] (In the text)
Wigner, J. 1932, Phys. Rev., 40, 749 [NASA ADS] [CrossRef] (In the text)

All Figures

$\begin{figure} \par\includegraphics[width=9cm,clip]{1303f1} \end{figure}$	Figure 1: Short-term Fourier transform ( bottom) of the sunspot number time series ( top). Fourier-transforms of the subsets of the dataset are also plotted. The colour scale of relative variations in the distributions is presented in Fig. 7.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.0cm,clip]{1303f2a}\vspace{4mm} \in... ...303f2b}\vspace{4mm} \includegraphics[width=8.0cm,clip]{1303f2c} \end{figure}$	Figure 2: Wavelet ( top and middle) and short-term Fourier transform (STFT) distribution ( bottom) of the yearly sunspot numbers.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{1303f3a}\hspace{5mm} \in... ...{5mm} \includegraphics[width=7.5cm,clip]{1303f3f} \vspace{0.5mm} \end{figure}$	Figure 3: Time-frequency distribution of test data. Panel A): PWD of signal ``I''; panel B): PWD of signal ``I'' with white noise; panel C): PWD of signal ``II'' with averaging and spline interpolation; panel D): PWD of signal ``II'' after Fourier filtering of the rotational period in the independent observational seasons and then averaging and spline smoothing; panel E): STFT of the same signal as in panel D); panel F): CWD of the same signal as in panel D).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=9cm,clip]{1303f4} \end{figure}$	Figure 4: False-alarm probability that a structure appears in STFT (black) and CWD (grey) due to observational noise. The amplitude of the added Gaussian noise ( from left to right) are 0.0025, 0.005, 0.01, 0.02, 0.04, and 0.08. The black rectangle represents the amplitudes of the visible real structures in Fig. 3.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=9cm,clip]{1303f5a}\vspace*{3mm} \includegraphics[width=9cm,clip]{1303f5b} \end{figure}$	Figure 5: Time-frequency transfer function of the processing. Top: STFT; bottom: PWD.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=7.95cm,clip]{1303f6a}\vspace{3mm} \i... ...1303f6b}\vspace{3mm} \includegraphics[width=7.95cm,clip]{1303f6c} \end{figure}$	Figure 6: PWD of the yearly sunspot number time series ( top) and the same of the datasets generated from the January ( middle) and July ( bottom) monthly averages.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=5.2cm,clip]{1303f7a}\hspace*{4.8mm} \includegraphics[width=3.2cm,clip]{1303f7b} \end{figure}$	Figure 7: PWD of the monthly solar radio flux data.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=12cm,clip]{1303f8} \end{figure}$	Figure 8: PWD of the long term envelope of the sunpot data.
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.