LONG-TERM TREND OF SUNSPOT NUMBERS

Some of the most important indicators of solar activity reach all time lows during the solar minimum between solar cycles (SCs) 23 and 24, and the level of solar activity during SC 24 is relatively lower. Thus, previous studies on the long-term trends of solar activity have concentrated on whether the level of solar activity would be reduced significantly during the consequent SCs, somewhat resembling the Maunder Minimum.

Penn & Livingston (2006, 2010) proposed that the sunspot magnetic field strength shows a decreasing trend while it waxes and wanes with SC from 1998 to 2010. Their results imply that there would be very few or virtually no sunspots during SC 25. Then papers by Pevtsov et al. (2011, 2014) suggested that there is not a statistically significant long-term trend in sunspot magnetic field strength from 1920 to 1958, from 1957 to 2011, or from 1920 to 2012. However, from 1874 to 2012, the magnetic field proxy shows a weaker long-term trend: the value gradually increases and reaches a broad maximum in 1920–1960, then decreases gradually (Pevtsov et al. 2014).

Gleissberg (1939) noted a periodicity of about 100 years, the so-called Centennial Gleissberg Cycle, in the maximum yearly sunspot numbers of SCs from 1750 to 1928. Similar findings were proposed by Garcia & Mouradian (1998), Javaraiah et al. (2005), and Hathaway (2010). Feynman & Ruzmaikin (2011) speculated that the low level of solar activity during SC 23/24 minimum favors the minimum of the Centennial Gleissberg Cycle. Reports by Solanki & Krivova (2011), Nielsen & Kjeldsen (2011), Zolotova & Ponyavin (2014), and Javaraiah (2015) confirm their conclusion.

The Hilbert–Huang transform (HHT), designed specifically for analyzing nonlinear and nonstationary data, can extract the intrinsically cyclical component from a sequence and has unprecedented prowess in revealing hidden physical meanings in data (Kataoka et al. 2009; Chen et al. 2010; Li et al. 2012; Deng et al. 2014; Gao & Xu 2016). Thus, in this paper, we investigated the long-term trend of yearly mean total sunspot numbers in the time interval of 1700–2015, which come from the World Data Center—the sunspot Index and Long-term Solar Observations (WDC-SILSO), using the HHT method.

In the next section, we describe briefly the data set and the methodology employed in this study. Then, Section 3 gives our results. Finally, conclusions and discussions are revealed in Section 4.

As we know, sunspot counts are available in the two main series: the International Sunspot Number (defined as ${R}_{z}=k(10g+f)$, where k is a correction factor for each observer, g is the number of sunspot groups, and f is the number of individual sunspots) created by R. Wolf in 1849, and the Group sunspot Number (defined as ${R}_{g}=\tfrac{12.08}{N}{\sum }_{i=1}^{N}({k}_{i}{g}_{i})$, where g_i represents the number of sunspot groups observed by the ith observer, k_i represents the correction factor for ith observer, and N represents the number of observers) initiated by Hoyt & Schatten (1998a, 1998b). However, the two sunspot series do not match by various aspects. Thus, through diagnosing and correcting flaws and biases affecting the two sunspot series, on 2015 July 1, a new version of the sunspot number is announced by WDC-SILSO, Royal Observatory of Belgium, Brussels (Clette et al. 2014). It is a milestone event in the history of the sunspot number. The yearly mean total sunspot numbers for the past 316 years (1700–2015) used in our study can be freely downloaded from the web site of SILSO (http://sidc.oma.be/silso/datafiles).

2.1. Hilbert–Huang Transform

HHT is proposed by Huang et al. (1996, 1998, 1999) and Wu & Huang (2009). It consists of two data analysis tools: the well-known Hilbert spectral analysis (HSA) and the recently developed ensemble empirical mode decomposition (EEMD), which is the key part of HHT. EEMD is a recursive "sifting" algorithm that locally extracts the robust and statistically significant cyclic components—defined as intrinsic mode functions (IMFs)—presented in a signal. IMFs satisfy the following conditions. (1) The number of extrema is equal to that of zero crossings or they differ at most by one. (2) The mean of the upper and lower envelopes (by connecting the maxima and minima respectively) is zero (Huang & Wu 2008). The effective algorithm of EEMD is summarized below (Huang & Wu 2008; Wu & Huang 2009).

• 1.  
  Add a white noise series to the time series $x(t)$ with equal spacing $\bigtriangleup t$.
• 2.  
  Find all the local extrema of $x(t)$ with added white noise.
• 3.  
  Connect the maxima (minima) by a cubic spline line as the upper (lower) envelope.
• 4.  
  Calculate the local mean of the two envelopes m₁.
• 5.  
  ${h}_{1}=x(t)-{m}_{1}$. If h₁ satisfies the definition of an IMF, we obtain IMF1: C1. Otherwise, iterate on h₁ through steps 2–4 until we obtain IMF1.
• 6.  
  ${R}_{1}=x(t)-C1$ is treated as the new time series in the next iteration steps 2–5 to obtain IMF2, IMF3, and so on. When no more IMF can be extracted from the residue, R_n, the decomposition process stops. R_n can be either the adaptive trend or a constant.
• 7.  
  Repeat steps 1–6 with different white noise series.
• 8.  
  Average the ensemble of corresponding IMFs and residue as the final result.

That is to say, EEMD can decompose the signal into a finite number of monocomponent IMFs, which are the intrinsic cycles of the original signal, and a residue, which can be either the adaptive trend or a constant.

Traditional frequency-space analysis methods, such as the Fourier transform, assume that the underlying processes are stationary in time. They have a prior basis, which does not necessarily represent the variety of underlying physical processes. Thus, when they are applied to the analysis of real-world processes, the results produced may contain artificial information. In contrast, EEMD is based on the local characteristics of the signal and has a posteriori adaptive basis. Compared to the Fourier transform, the advantages of EEMD is that it is designed specifically for analyzing nonlinear and nonstationary data.

The Fourier transform decomposes a signal into the frequencies. For EEMD, $x(t)={\sum }_{i=1}^{n}{C}_{i}+{R}_{n}$, where C_i is the ith IMF, which is extracted from the data $x(t)$,  and R_n is the residue of the data $x(t)$ after n IMFs are extracted. EEMD is more appropriate for removing short-term fluctuations and reveals long-term trends. This is another advantage of EEMD. Thus, in this study, applying the EEMD method, we investigate the long-term trend of yearly mean total sunspot numbers.

HHT applies the HSA method to obtain instantaneous frequency of IMFs (Huang et al. 1998, 1999; Wu et al. 2011):

$$\begin{eqnarray}&&{\omega }_{i}(t)=\displaystyle \frac{d\left(\arctan \tfrac{{y}_{i}}{{C}_{i}}\right)}{{dt}},\end{eqnarray} \tag{ 1 }$$

where C_i is the ith IMF and y_i is the Hilbert transform (Gabor 1946; Van der Pol 1946) of the ith IMF. Then, we can calculate the mean periods of the IMF.

To ensure that an IMF contains a true signal, Wu & Huang (2004, 2005) proposed a method to test the statistical significance of IMFs.

• 1.  
  Calculate the energy of IMFs

  $$\begin{eqnarray}&&{{NE}}_{n}=\displaystyle \sum _{j=1}^{N}| {C}_{n}(j){| }^{2},\end{eqnarray} \tag{ 2 }$$

  where ${C}_{n}(j)$ represents the nth IMF and N represents the number of data points.
• 2.  
  The first IMF is used to obtain the relative energy for other IMFs.
• 3.  
  Select the confidence-limit level, e.g., 95%, and calculate the spread function of 95th percentile of white noise.
• 4.  
  If the energy level of any IMF lies above the spread line for the 95th percentile, the IMF is statistically significantly at 95% confidence levels.

In this study, the code of HHT prepared by Z. Wu (http://rcada.ncu.edu.tw/research1.htm) is utilized to investigate the long-term trend of yearly mean total sunspot numbers.

First, we plot the yearly mean total sunspot numbers as a function of time for the interval 1700–2015, as shown in Figure 1. Then, applying the EEMD method, the yearly mean total sunspot numbers are decomposed into seven IMFs and an adaptive trend. Figure 2 shows the seven IMFs of yearly mean total sunspot numbers. We also calculate the mean periods of seven IMFs and their errors, which represent the uncertainty in the mean ($\sigma /{n}^{1/2}$, where σ and n refer to the standard deviation and the number of data points, respectively), which are collected in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to ensure that an IMF contains a true signal, we test the statistical significance of the seven IMFs (Figure 3). Figure 3 shows that IMFs 1–6 are statistically significant at 95% confidence levels, which indicate that periodicities of 3.71 ± 0.46 (Quasi-Triennial Oscillations; QTOs), 10.53 ± 0.11 (Schwabe Cycle), 21.06 ± 0.60 (Hale Cycle), 51.41 ± 1.16, 106.13 ± 8.24 (Centennial Gleissberg Cycle), and 161.34 ± 8.24 years are statistically significant in the yearly mean total sunspot numbers.

Zoom In Zoom Out Reset image size

The adaptive trend of the yearly mean total sunspot numbers is shown in Figure 4. From the adaptive trend, we can find that the value gradually increases during the time period 1700–1975 (from 58 to 91), then decreases gradually from 1975 to 2015 (from 91 to 90).

Zoom In Zoom Out Reset image size

EEMD is a recursive "sifting" algorithm. IMF7 is below the 95% confidence level line, suggesting that a periodicity of 325.65 ± 12.37 years is not statistically meaningful. Thus, in order to further determine the trend of yearly mean total sunspot numbers in the time interval of 1700–2015, we then remove IMFs 1–6 (they are statistically significant) and reveal the overall trend (the sum of the IMF7 and the adaptive trend) of the yearly mean total sunspot numbers, also shown in Figure 4. From Figure 4, we can find that the overall trend is similar to the adaptive trend and, since the 1970s, the speed of decrease in the overall trend is faster than that in the adaptive trend.

The maximum and minimum yearly mean total sunspot numbers of SC in the time interval of 1700–2015 are shown in Figure 5. We can find that the maximum and minimum yearly mean total sunspot numbers wax and wane, consistent with the IMF5 of the yearly mean total sunspot numbers (Figure 5), the so-called Centennial Gleissberg Cycle. That is to say, the Centennial Gleissberg Cycle is extracted from the yearly mean total sunspot numbers, which is the reason why there is a periodicity of about 100 years in the maximum and minimum yearly mean total sunspot numbers, supplementing the results obtained by Gleissberg (1939), Garcia & Mouradian (1998), Javaraiah et al. (2005), and Hathaway (2010).

Zoom In Zoom Out Reset image size

The IMF5 (Centennial Gleissberg Cycle) reach minima during Dalton Minimum and Gleissberg Minimum. Recently, IMF5 (Centennial Gleissberg Cycle) reached minimum once again. From the IMF5 of the yearly mean total sunspot numbers, the maximum and minimum yearly mean total sunspot numbers of SC, we can infer that (1) the new grand minimum is in progress and (2) the Dalton Minimum, the Gleissberg Minimum, and the low level of solar activity during SC 24 (the part of the new grand minimum) all can be understood as minima of the Centennial Gleissberg Cycle, which confirms the conclusions of Feynman & Ruzmaikin (2011), Solanki & Krivova (2011), Nielsen & Kjeldsen (2011), Zolotova & Ponyavin (2014), and Javaraiah (2015).

We then remove IMFs 1–4 (they are statistically significant and their mean periods are shorter than 100 years) and reveal the 100-year and longer timescale trend (the sum of IMFs 5–7 and the adaptive trend) of the yearly mean total sunspot numbers. The 100-year and longer timescale trend of the yearly mean sunspot numbers is shown in Figure 5. The maximum and minimum yearly mean total sunspot numbers of SC from 1700 to 2015 are also shown in Figure 5.

SCs 5–7 are regarded as the Dalton Minimum and SCs 12–16 are regarded as the Gleissberg Minimum by Zolotova & Ponyavin (2014), who proposed a similarity between SC 23 and SCs 4 and 11 on the eve of the Dalton Minimum and the Gleissberg Minimum. Thus we consider SC 24 as the part of the new grand minimum.

Figures 4 and 5 show that the adaptive (overall) trend of the yearly mean total sunspot numbers is consistent with the variations of the minimum values of the 100-year and longer timescale trends of the yearly mean total sunspot numbers. In addition, the variations of the maximum and minimum yearly mean total sunspot numbers of SC during the Dalton Minimum, the Gleissberg minimum, and SC 24 (part of the new grand minimum) are consistent with the adaptive (overall) trend of the yearly mean total sunspot numbers (Figure 4). Based on the adaptive (overall) trend, the 100-year and longer timescale trend of yearly mean total sunspot numbers, we can infer that the level of solar activity during the new grand minimum may be close to that during the Gleissberg Minimum and slightly higher than that during the Dalton Minimum.

That is to say, the level of solar activity is relatively low during the next several SCs. However, from the adaptive (overall) trend, the 100-year and longer timescale trend of the yearly mean total sunspot numbers, we can infer that the level of solar activity during the new grand minimum is significantly higher than that during the Maunder Minimum.

In this paper, we investigate the long-term trend of yearly mean total sunspot numbers from 1700 to 2015, which is a new version of sunspot number announced by WDC-SILSO. First, using the HHT method, the yearly mean total sunspot numbers are decomposed into seven IMFs and an adaptive trend. Then, we investigate the adaptive (overall) trend, the Centennial Gleissberg Cycle (IMF5), and the 100-year and longer timescale trend of the yearly mean total sunspot numbers. The main results of our study can be summarized as follows. (1) From the adaptive trend, we can find that the value gradually increases during the time period 1700–1975, then decreases gradually from 1975 to 2015. (2) The Centennial Gleissberg Cycle is extracted from the yearly mean total sunspot numbers and confirms that the Dalton Minimum, the Gleissberg Minimum, and the low level of solar activity during SC 24 can all be understood as minima of the Centennial Gleissberg Cycle. That is to say, a new grand minimum is in progress. (3) Based on the adaptive (overall) trend, the 100-year and longer timescale trend of yearly mean total sunspot numbers, we can infer that the level of solar activity during the new grand minimum may be close to that during the Gleissberg Minimum, slightly higher than that during the Dalton Minimum, and significantly higher than that during the Maunder Minimum.

There is the completeness of the sunspot record (sunspots have been observed every day) from SC 10 to present, the so-called modern era sunspot SCs which are considered to be the most reliably known SCs (Wilson 1988; Li et al. 2015). Thus, we also investigate the long-term trend of yearly mean total sunspot numbers from SC 10 to present (from 1856 to 2015) and can also obtain similar results.

From linear least-squares fitting of magnetic field strength series, Penn & Livingston (2006, 2010) found that a decreasing trend in the magnetic field strength from 1998 to 2010. Based on second-degree polynomial approximation, Pevtsov et al. (2014) showed a long-term trend: the value of magnetic field proxy increases gradually from 1874 to 1920 and reaches a broad maximum in 1920–1960, then decreases gradually from 1960 to 2012. In this study, using EEMD, the yearly mean total sunspot numbers are decomposed into seven IMFs and an adaptive trend. From the adaptive trend, we can find that the value gradually increases from 1700 to 1957, then decreases gradually from 1957 to 2015. In order to further determine the trend of yearly mean total sunspot numbers in the time interval of 1700–2015, IMFs 1–6, which are statistically significant at 95% confidence levels are removed from the yearly mean total sunspot numbers and the overall trend is obtained. The overall trend is similar to the adaptive trend. The adaptive (overall) trend is consistent with the trends in the magnetic field strength from 1998 to 2010 and from 1874 to 2012 (Penn & Livingston 2006, 2010; Pevtsov et al. 2014).

The mean periods of statistically significant IMFs are 3.71 ± 0.46, 10.53 ± 0.11, 21.06 ± 0.60, 51.41 ± 1.16, 106.13 ± 8.24, and 161.34 ± 8.24 years. The periods of 10.53 ± 0.11 years and 21.06 ± 0.60 years correspond to the most prominent periodicities: the 11-year sunspot number cycle (Schwabe 1844), the so-called Schwabe Cycle, and the 22-year Hale magnetic cycle (Hale 1924), respectively. The period of 3.71 ± 0.46 years corresponds to the so-called QTO. QTO behavior has been reported in many solar and geomagnetic activity indices (Gonzalez & Gonzalez 1987; Kane 1997, 2005a, 2005b; Gao et al. 2012; Qu et al. 2015).

Numerous authors noted the Centennial Gleissberg Cycle in the sunspot cycle amplitudes (Gleissberg 1939; Garcia & Mouradian 1998; Javaraiah et al. 2005; Hathaway 2010). The low level of solar activity during SC 23/24 minimum and SC 24 is speculated to favor the minimum of the Centennial Gleissberg Cycle (Feynman & Ruzmaikin 2011; Nielsen & Kjeldsen 2011; Solanki & Krivova 2011; Zolotova & Ponyavin 2014; Javaraiah 2015). Using EEMD, the Centennial Gleissberg Cycle is extracted from the yearly mean total sunspot numbers. Furthermore, from the Centennial Gleissberg Cycle, we can infer that the new grand minimum is in progress and the Dalton Minimum, the Gleissberg Minimum, and the low level of solar activity during SC 24 (part of the new grand minimum) all can be understood as minima of the Centennial Gleissberg Cycle.

Penn & Livingston (2006, 2010) proposed that, if this trend—a decrease in the sunspot magnetic field strength during SC 23 and part of the rise phase of SC 24 (1998–2010)—continues, then there would be very few or virtually no sunspots during SC 25, indicating the onset of a Maunder-Minimum-type Grand Minimum. In this study, the adaptive (overall) trend in the time interval of 1998–2010 is consistent with the trend in the magnetic field strength (Penn & Livingston 2006, 2010; Pevtsov et al. 2014). However, from the adaptive (overall) trend, the 100-year and longer timescale trend, which have been extracted from yearly mean total sunspot numbers for the past 316 years (1700–2015), we can infer that the level of solar activity during the new grand minimum may be close to that during the Gleissberg Minimum, slightly higher than that during the Dalton Minimum, and significantly higher than that during the Maunder Minimum. Our results do not support the suggestion that a new grand minimum, somewhat resembling the Maunder Minimum, is in progress.

P.X.G. thanks the referee and the editor for their careful reading and constructive comments, which improved the original version of the manuscript. The sunspot numbers sequence in this study is announced by WDC-SILSO, Royal Observatory of Belgium, Brussels. This work is supported by the National Natural Science Foundation of China (Grant no. 11673061 and 11633008), the Applied Basic Research Foundation of Yunnan Province, China (grant no. 2014FB190), the Foundation of Key Laboratory of Solar Activity of National Astronomical Observatories of Chinese Academy of Sciences (KLSA201507), and the Chinese Academy of Sciences.