© T. Asikainen & J Mantere, Published by EDP Sciences 2023

1 Introduction

Prediction of the sunspot number has been an everlasting interest in the space science community since the discovery of the sunspot cycle by Schwabe (1844). The sunspot number is an indirect indicator of many different solar phenomena, e.g., total and spectral solar radiation (e.g. Krivova et al., 2011; Frohlich, 2012), coronal mass ejections (e.g. Richardson & Cane, 2012), solar flares and magnetic active regions (e.g. Toriumi et al., 2017). Its cyclic variation can even be used as a pacemaker to time different aspects of solar activity, solar wind, and resulting geomagnetic variations (Chapman et al., 2021; Leamon et al., 2022). Therefore, there is considerable practical interest in predicting the evolution of future sunspot cycle(s). This is especially true in today’s technological society where space hazards pose a significant threat, e.g., to satellites, communications, and electric grids on ground (e.g. Lanzerotti, 2001). Another interest in predicting sunspots arises from the relatively recently recognized influences of variable solar radiation and solar wind activity on Earth’s climate system (Gray et al., 2010; Ward et al., 2021; Salminen et al., 2020).

Over the period of about the last 100 years, a vast array of different methods ranging from statistical methods to intensive physical simulations have been developed for predicting sunspots. As an unbiased introduction to all of them would be a futile effort here, the interested reader is referred to several excellent reviews on the subject by Hathaway (2009), Pesnell (2012), and Petrovay (2020). According to the classic solar dynamo theory, the poloidal solar magnetic field in the solar minimum gets stretched to the toroidal magnetic field of the next cycle, which then produces sunspots and magnetic active regions that ultimately make up the poloidal field in the next solar minimum (Parker, 1955; Babcock, 1961; Leighton, 1969; Charbonneau, 2020). Physically motivated solar cycle predictions are based on numerical modeling of the solar dynamo process and the transport of magnetic flux on the solar surface (Charbonneau, 2020; Nandy, 2021; Karak, 2023). However, some of the most successful, yet much simpler prediction methods are based on precursors that serve as indicators for the strength of the coming solar cycle. The so-called polar field precursor methods have found a good correlation between the magnetic field observed at the solar polar region up to a few years before the sunspot minimum and the amplitude of the next sunspot cycle (Schatten et al., 1978; Petrovay, 2020; Kumar et al., 2021, 2022). As the polar field reflects the strength of the poloidal phase of the solar magnetic field the precursor methods are deeply rooted in the core idea of the dynamo theory.

Because the polar field has been systematically measured only since the 1970s some long-running proxy measures for the polar field have also been used. The most successful polar field proxies are based on geomagnetic activity measures close to the solar minimum (e.g., Ohl, 1966; Du et al., 2012). It has been shown that certain geomagnetic activity indices, e.g., the aa index correlate quite well with the open solar magnetic flux carried by the solar wind (e.g. Lockwood et al., 1999, 2014). Close to the solar minimum, the geomagnetic activity is dominantly driven by high-speed solar wind streams (e.g. Richardson & Cane, 2012), which emanate from coronal holes that eventually form the solar polar field in the declining phase of the cycle (Krieger et al., 1973; Bame et al., 1976). During these times the geomagnetic activity is most directly connected to the polar field and, thereby acts as a good proxy for the amplitude of the next sunspot cycle.

In accordance with the solar dynamo theory the sunspot cycle prediction methods often consider only the predictability of the cycle based on the previous cycle. Most of these methods do not typically take into account the 22-year Hale cycle of solar magnetism, which is well-known in the solar cycle (e.g. Gnevyshev & Ohl, 1948; Takalo & Mursula, 2018; Leamon et al., 2022) and geomagnetic phenomena (e.g. Chernosky, 1966; Takalo, 2021; Chapman et al., 2021). However, some recent studies have considered the even/odd cycle parity in the context of sunspot cycle prediction (e.g. Du, 2020; Kakad & Kakad, 2021; Penza et al., 2021; Du, 2022b; Nagovitsyn & Ivanov, 2023). Here we study the prediction of sunspot cycles accounting for the 22-year Hale cycle by considering the differences in the even- and odd-numbered sunspot cycles. In Section 2, we first present our data and then in Section 3 proceed to show that the time evolution and shape of all sunspot cycles are extremely well-described by a simple parameterized mathematical expression. In Section 4, we discuss the mutual dependencies of the parameters describing the sunspot cycles and in Section 5, we show that the parameters can be predicted using precursors found partly from past sunspot number values and partly from geomagnetic activity, which is often used as a precursor proxy for solar polar magnetic field in the sunspot minimum. Most importantly, though, we find that even and odd sunspot cycles obey different statistics therefore implying Hale-cyclicity in their predictability. Separately these statistical relations are stronger and more significant than those based on combining even and odd cycles together. Using these statistics we construct in Section 6 a new method to separately predict even and odd sunspot cycles and apply it to the coming solar cycle 25. We also find an interesting connection between consecutive odd-even sunspot cycle pairs, which allows us to make rough early predictions for cycle 26 as well. In Section 7, we discuss the results and give our conclusions.

2 Data

In this work we use monthly values of the version 2 Sunspot Number (SSN) obtained from the SILSO World Data Center (1749–2023). At the time of writing this paper the SSN v2 series covers time from 1749 to July 2023 therefore covering full sunspot cycles 1–24 and about 3 years from the start of sunspot cycle 25. The monthly values of SSN were first smoothed with a 13-month window, where the first and last month are given a weight of 0.5 and all other months a weight of 1. Using this smoothed SSN curve we identified the times of sunspot minima and maxima.

In addition to the sunspot data we use in this work the geomagnetic aa index. Recently Lockwood et al. (2018a,b) presented a homogenized version of the aa index, which inter-calibrates the observations of the different stations used to compose the aa index. They also correct the index for the secular drift in the location of the auroral oval (due to secular changes in Earth’s magnetic field) in relation to the observing stations. The homogenized aa index was obtained from the supplementary data of Lockwood et al. (2018b), which offers the 3-hourly values for the years 1868–2017.

We first extended this series forward in time from January 2018 to December 2021 by calculating the monthly means of the homogenized aa index and by calibrating the raw aa index (obtained from ISGI: http://isgi.unistra.fr/) against the corresponding homogenized values. The calibration was found by fitting a regression line to the logarithmic monthly averaged homogenized aa (aa_H) and raw aa values using the data between 1980 and 2017 when the aa index is based on the latest pair of observatories (Hartland in England and Canberra in Australia). The best-fit line was

$$\log \left(a{a}_{\mathrm{H}}\right)=1.049\left(\pm 0.004\right)\times \log \left(\mathrm{aa}\right)-0.257(\pm 0.011).$$(1)

The uncertainties of the fit parameters have been indicated in parentheses. Note that this scaling uses logarithmic values, since in the logarithmic scale the residuals of the fit are closely homoscedastic (constant variance) while in the linear scale they display large heteroscedasticity thereby compromising the basic assumptions of the least-squares fit. The 3-hourly raw aa values since January 2018 were then scaled with equation (1) and the resulting data was appended to the homogenized aa index time series.

We also extended the aa index backward in time using the daily magnetic declination based Ak(D) indices recorded at the Helsinki geomagnetic observatory since 1844 (Nevanlinna, 2004) and obtained from https://space.fmi.fi/MAGN/K-index/. These values have previously been used successfully to extend the aa index time series backward from 1868 to 1844 (e.g. Lockwood et al., 2014).

Unlike Lockwood et al. (2014), who used annual averages, we calibrated here the daily Ak(D) values against the simultaneous homogenized aa index values by Lockwood et al. (2018b) for the overlapping time period from 1.1.1868 to 31.12.1879, where the Ak(D) data is rather continuous. Starting from 1880 the Ak(D) data series has large data gaps. We found that the daily Ak(D) values can be scaled to the corresponding daily homogenized aa index values by the following equation

$$\log \left({\mathrm{aa}}_{\mathrm{H}}\right)=1.199\left(\pm 0.012\right)\times \log \left(\mathrm{Ak}\left(D\right)+5 \mathrm{nT}\right)-0.97(\pm 0.04).$$(2)

Also here the logarithmic scale ensures homoscedasticity of the fit residuals. The correlation between the scaled Ak(D) values and the homogenized aa index values is 0.84 indicating a rather reliable scaling.

Figure 1 shows the annual averages of homogenized aa index (blue), scaled raw aa index (red) (shown since 1980), and scaled Helsinki Ak(D) values (shown for years 1844–1879). The times when the different curves overlap indicate a very good correspondence between the different datasets. The extended aa index time series is formed by the scaled Ak(D) data from 1844 to 1867, by the homogenized aa index from 1868 to 2017, and by the scaled raw aa index since 2018. For the purposes of this study we will use monthly averages of the extended aa index composite.

Figure 1 Annual averages of homogenized aa index (blue), scaled raw aa index (red) (shown since year 1980) and scaled Helsinki Ak(D) values (shown for years 1844–1879).

3 Parameterization of the sunspot cycle

Sunspot cycles are often fitted with a parameterized curve (e.g. Stewart & Panofsky, 1938; Hathaway et al., 1994; Volobuev, 2009). Stewart & Panofsky (1938) showed that the sunspot cycles could be described roughly by curves of the form $c{\left(t-t_0\right)}^a{e}^{-b(t-t_0)}$, where a, b, c, and t₀ are free parameters. Hathaway et al. (1994) used a slightly modified version

$$f\left(t\right)=\frac{a{\left(t-t_0\right)}^3}{\exp \left(\left(t-t_0\right)/b^2\right)-\mathrm{c}},$$(3)and fitted this model curve to sunspot cycles 1–21. They also showed that many of the parameters correlated rather well with each other thereby offering a possibility to reduce the effective number of parameters in the curve down to one free parameter. Many studies have since used similar parameterizations. However, to be useful in predicting the sunspot cycles the parameters of the future cycle should be predicted by some means. Several studies have found relatively good correlations between different precursors and the amplitude of the sunspot cycle or the maximum SSN during the cycle. The parameterizations described above (e.g., Eq. (3)) may not be optimal in light of these correlations, because for those the amplitude (maximum) of the sunspot cycle depends on a combination of several parameters.

Therefore, we formulated a new parameterization for the sunspot curve, where the parameters are perhaps better interpreted in terms of well-known solar cycle properties (amplitude, rise time, asymmetry etc.). We use here an asymmetric Gaussian curve of form

$$f\left(t\right)=A\exp \left(-\frac{{\left(t-B\right)}^2}{{\left(C\times g\left(t,B\right)\right)}^2}\right),$$(4)where A is the sunspot cycle maximum, B is the time of the sunspot maximum measured from the sunspot minimum beginning the cycle (i.e., B is the cycle rise time), C is the time scale of the rising phase (cf. standard deviation of a Gaussian), and function g(t, B) is defined as

$$g\left(t,B\right)=\{\begin{array}{c}1, \mathrm{if} t\le B,\\ D,\mathrm{if} t<B.\end{array}$$(5)

Therefore, the parameter D appearing in g(t, B) describes the asymmetry of the time scales in the declining and rising phases. The more positive the D is the longer the declining phase is compared to the rising phase.

We then fitted equation (4) to the 13-month smoothed SSN of each sunspot cycle separately. Each cycle was defined from the time of the sunspot minimum to the next minimum. The 4-parameter fit was done with the non-linear Levenberg-Marquardt optimization implemented in Matlab software. Although the consecutive sunspot cycles are known to be overlapping (so that the cycle already starts before the minimum SSN time), the fitting of our parametric model is largely dependent on the whole cycle and results are not strongly sensitive to the data around sunspot minima. This was tested either by leaving out up to 1 year of data from the beginning of the cycles or by extending the previous fitted cycle and subtracting it from the next cycle. In either case, the fitted parameter values remained practically the same.

Figure 2 shows the time series of the 13-month smoothed SSN in black, the 4-parameter model fits for cycles 1–24 in red and the Hathaway et al. (1994) fit of equation (3) in blue. One can see that the fitted curves describe all the individual cycles reasonably well, although some of the detailed structure in the SSN cycles cannot be described by a smooth asymmetric Gaussian. Such structures are for example the very sharp-peaked cycles like cycles 1–4 and 8–11 and the often seen double peaks (see e.g., Karak et al., 2018), which are quite prominent, e.g., in cycles 22–24. However, the rising and declining phases of each cycle are well-captured by the fit and the cycle amplitudes are quite close to the real cycle amplitudes. The average R²-value of the 4-parameter fit over all cycles is 0.973, indicating that over 97% of the variability in the SSN cycles is captured by the model curves. The average root-mean-squared error for all the cycles is 8.7 and there is no statistically significant difference in this between even and odd numbered sunspot cycles. We also note that for most cycles the asymmetric Gaussian function used here is very close to the function (Eq. (3) used by Hathaway et al. (1994). However, in some cycles (3, 4, 8, 10) there is a clear difference with the asymmetric Gaussian providing a somewhat better fit.

Figure 2 Time series of 13-month smoothed sunspot number (black), parameterized fits for each sunspot cycle (red curves), and the fit by Hathaway et al. (1994) for comparison (blue curves). The sunspot cycle numbers are indicated on the plot as well. The top panel shows cycles 1–8, middle panel cycles 9–16, bottom panel cycles 17–24.

4 Relationships between cycle parameters

Let us next consider the relationships between the fitted values of the four parameters over all cycles. Figure 3 shows as a scatter plot the relationship between parameters C and D^−1/2. In Figure 3, the odd-numbered cycles have been indicated with blue dots and even-numbered cycles with red squares.

Figure 3 Relationship between the C and D−1/2 parameters of the 4-parameter fits. Blue dots (red squares) indicate even (odd) numbered cycles. The cycle numbers have been further indicated beside all points. The yellow line depicts the linear fit to all the cycles having a correlation coefficient of 0.94 (p-value is 10−11).

One can see that all the cycles depict quite a strong correlation between C and D^−1/2 and that there is no significant difference between the odd and even cycles. The correlation coefficient between C and D^−1/2 over all the cycles is 0.94 and it is highly significant (p-value is 10⁻¹¹). This relationship indicates that cycles with a steep rising phase (i.e., small C) have a relatively more gradual (i.e., large D and small D^−1/2) declining phase and vice versa. The correlation in Figure 3 is therefore a manifestation of the well-known property of sunspot cycles that small amplitude cycles tend to be more symmetric about the sunspot maximum time than large amplitude cycles (Waldmeier, 1968).

For our purposes the tight relationship between C and D allows us to eliminate the parameter D from the sunspot model curve and thereby reduce the number of free parameters to 3. The parameter D is therefore replaced in equations (4) and (5) by the expression

$$D={\left(0.25\left(\pm 0.05\right)+0.226\left(\pm 0.02\right)\times C\right)}^{-2},$$(6)which corresponds to the linear fit (yellow line) in Figure 3.

After replacing D with equation (6) we repeated the fitting procedure for each sunspot cycle, but now using the remaining 3-parameter model. The correlation between C and D^−1/2 is so high that the 3-parameter model is roughly equally good as the 4-parameter model in describing each sunspot cycle. The average R²-value for the 3-parameter model is 0.959 and therefore not much smaller than 0.973 for the 4-parameter model. Note also that the elimination of D from the model does not significantly alter the values of the remaining parameters A, B, and C. The correlations of the corresponding parameters of the 4-parameter and 3-parameter fits exceed 0.97.

Figure 4 shows the relationship between the parameter B (time of SSN maximum counted from SSN minimum in years) and cycle amplitude A. One can see that there is a general anti-correlation between A and B. This anti-correlation between the cycle rise time and amplitude has been long known as the Waldmeier effect (Waldmeier, 1935). Here one can see that this anti-correlation is somewhat stronger (cc = −0.88, p-value = 10⁻⁴) in even cycles than in odd cycles (cc = −0.69, p-value = 0.014). The correlation between A and B found here is clearly higher than, e.g., the correlation between SSN maximum and cycle rise time found by Karak & Choudhuri (2011) (cc = −0.5). It therefore appears that the parameters of the fitted SSN model curve indicate the Waldmeier effect more robustly than the exact SSN data. This is likely because of the sensitivity of the Waldmeier effect to the timing and height of the sunspot maximum, which can be difficult to determine if the cycle has multiple peaks Karak & Choudhuri (2011).

Figure 4 Relationship between the A and B parameters of the 3-parameter fits. Blue dots (red squares) indicate even (odd) numbered cycles. The cycle numbers have been further indicated beside all points. The blue and red lines depict the linear fit to odd and even cycles respectively.

The linear fits to even and odd cycles depicted by the red and blue lines in Figure 4 are given by equations

$$B_{\mathrm{even}}=5.9\left(\pm 0.4\right)-0.013\left(\pm 0.003\right)\times A_{\mathrm{even}},$$(7)

$$B_{\mathrm{odd}}=7.2\left(\pm 1.0\right)-0.016\left(\pm 0.006\right) \times A_{\mathrm{odd}}.$$(8)

The correlation coefficients and the slopes/intercepts of the best-fit regression lines are different in the even and odd cycles only with a weak statistical significance (the p-values for the differences exceed 0.12). However, the mean squared errors of the linear fits to the even and odd cycles are highly significantly different (0.14 for even cycles and 1.01 for odd cycles and the p-value for the difference is 0.003). This result indicates that the Waldmeier effect is more strongly valid for even cycles than for odd cycles.

Figure 5 shows the relationship between the final parameter C (time scale of rising phase) and cycle amplitude A. Here one can also see a general anti-correlation, which is yet another manifestation of the Waldmeier effect. Actually, there is also a quite good correlation (cc = 0.91, p-value = 8 × 10⁻¹⁰) between the B and C parameters, which explains the two manifestations of the Waldmeier effect. The linear fits to even and odd cycles depicted by the red and blue lines in Figure 4 are given by equations

$$C_{\mathrm{even}}=3.5\left(\pm 0.4\right)-8\left(\pm 2\right){\times 10}^3\times A_{\mathrm{even}},$$(9)

$$C_{\mathrm{odd}}=4.7\left(\pm 0.7\right)-1.3\left(\pm 0.4\right){\times 10}^{-2}\times A_{\mathrm{odd}}.$$(10)

Figure 5 Relationship between the A and C parameters of the 3-parameter fits. Blue dots (red squares) indicate even (odd) numbered cycles. The cycle numbers have been further indicated beside all points. The blue and red lines depict the linear fits to odd and even cycles respectively.

The conclusions about the differences between even and odd cycles are the same as for B vs. A relation in Figure 4. That is, the linear relationships or the correlations are only weakly statistically significantly different, but the difference in the mean squared errors is highly significant (p-value is 0.038).

Overall the results in Figures 4 and 5 indicate that the cycle amplitude A could further be used to reduce the number of parameters in the sunspot cycle fit. Moreover, there are indications that the accuracies of these fits are significantly different in even and odd cycles.

5 Precursors for cycle parameters

Based on the above relations we could further reduce the SSN model curve to a 1-parameter model. However, each simplification of the model makes it less flexible and decreases the model accuracy. Therefore, it is useful to first consider whether the cycle parameters A, B, and C can be directly predicted by some suitable precursor. Cameron & Schüssler (2007) showed that the solar activity level 3 years before the sunspot minimum starting the cycle is a relatively good predictor for the amplitude of the cycle. This result for sunspot numbers was shown by Petrovay (2020) (his Figure 6), who found a correlation of 0.8 between the maximum sunspot number of the cycle and the sunspot number taken 3 years before the sunspot minimum that starts the cycle. A correlation of 0.8 implies that only about 62% variability in cycle amplitudes could be explained by the past sunspot number, thereby offering a rather mediocre accuracy in predicting the coming cycle amplitudes as Petrovay (2020) mentions. However, still motivated by this result we calculated in Figure 6 the correlation coefficient between the cycle amplitude A and the lagged sunspot number as a function of the time lag in months before the SSN minimum. Unlike in past studies, we did the calculation here separately for even and odd cycles and using the 2nd power of SSN (SSN²), which we found to produce slightly larger correlations than SSN (the difference, however, is rather small and not statistically significant).

Figure 6 Correlation coefficient between cycle amplitude A and lagged 13-month smoothed sunspot number (SSN2) as a function of the time lag in months before the SSN minimum. Blue (red) curve indicates the odd (even) cycles and the correspondingly colored regions indicate the 95% confidence limits for the correlation coefficients. The green dot indicates the maximum correlation at 41 months before the SSN minimum for the even cycles.

One can see that for the even cycles the correlation is systematically higher than for the odd cycles. The difference between even and odd cycles is quite large, but because of the rather low number of data points (only 12 even and odd cycles) the 95% confidence limits remain rather large even when the correlation coefficient is high. However, there is a location around the lag of 41 months (about 3.4 years) before sunspot minimum where the correlations do differ from each other statistically significantly. This lag of optimal correlation is quite close to the 3 years found by Cameron & Schüssler (2007). Note, however, that the good correlation is not specific to exactly the optimum of 41 months but is seen over a broad range of lags from 34 to 44 months.

We chose the lag of 41 months for a closer inspection in Figure 7, which displays all the three parameters A, B, and C^1/2 as a function of the sunspot number taken 41 months before the sunspot minimum (SSN(41)). One can see that the SSN(41)² correlates quite well not only with cycle amplitude A but with all three parameters in the even cycles. In the odd cycles, there is evident correlation too, but it is clearly much lower than for the even cycles due to the larger scatter. In particular the SSN(41)² vs. cycle amplitude A resembles the plot in Figure 6 of Petrovay (2020), but now reveals a large difference between even and odd cycles. Recently Du (2020) also found that the sunspot number 39 months before the sunspot minimum is a precursor for the maximum SSN of the next cycle and this relation is stronger for even cycles. A similar finding was done by Nagovitsyn & Ivanov (2023).

Figure 7 Dependence of A, B, and C parameters on 13-month smoothed sunspot number 41 months before the sunspot minimum (SSN(41)2). Blue dots (red squares) indicate odd (even) numbered cycles. The cycle numbers have been further indicated beside all points. The red lines depict the linear fits to even numbered cycles.

We shall later discuss the reasons for the better correlation in even cycles but for now we concentrate on the fact that SSN(41) can be used as a quite accurate predictor of the three cycle parameters for the even cycles. For cycle amplitude, the correlation coefficient is 0.984 (p-value = 8 × 10⁻⁹) indicating that about 96.9% of the variation in cycle amplitude may be predicted with SSN(41). For cycle rise time B the correlation is slightly lower −0.86 (p-value = 0.0003) and for C^1/2 (time scale of rising phase) the correlation is −0.89 (p-value = 10⁻⁴).

The linear fits to the even cycles indicated by the yellow lines in Figure 7 are given by equations

$$A_{\mathrm{even}}=88\left(\pm 6\right)+0.0103\left(\pm 0.0006\right)\times \mathrm{SSN}{\left(41\right)}^2,$$(11)

$$B_{\mathrm{even}}=4.8\left(\pm 0.3\right)-1.3\left(\pm 0.3\right)\times {10}^{-4}\times \mathrm{SSN}{\left(41\right)}^2,$$(12)

$$C_{\mathrm{even}}={\left(1.70\left(\pm 0.05\right)-3.1\left(\pm 0.5\right)\times {10}^{-5}\times \mathrm{SSN}{\left(41\right)}^2\right)}^2.$$(13)

We also studied the correlation between cycle amplitude and geomagnetic activity, which is known to be a good precursor for the sunspot cycle amplitude. We tested various different quantities calculated from the extended homogenized geomagnetic aa index. We found that the best predictor for the cycle amplitude for odd sunspot cycles is given by the maximum average September–October aa³ value within the 3-year window extending backward from the sunspot minimum. Figure 8 shows the relationship between the cycle amplitude and this geomagnetic activity measure separately for the odd cycles (blue) and even cycles (red). The yellow line indicates the best linear fit of

$$A=130\left(\pm 7\right)+0.0036\left(\pm 0.0004\right)\times {\max }_{\left[\mathrm{0,3}\mathrm{y}\right]}\left({\mathrm{aa}}_{\mathrm{Sep}-\mathrm{Oct}}^3\right).$$(14)

Figure 8 Relationship between cycle amplitude A and maximum September–October aa3 from a 3-year window before SSN minimum. The blue points indicate odd cycles and red squares the even cycles. The blue line indicates the linear fit to the odd cycles. The green circle indicates the estimated amplitude for cycle 25.

The linear correlation coefficient between the odd cycle amplitude and maximum ${\mathrm{aa}}_{\mathrm{Sep}-\mathrm{Oct}}^3$ is 0.981 (p < 9 × 10⁻⁵). One can also see from Figure 8 that the corresponding relation for the even cycles is much worse, but still statistically significant (correlation 0.76, p = 0.027). It should be noted though, that trying different quantities calculated from the aa index one can obtain higher correlations for the even cycles as well, but it seems that none of such correlations exceeds the extremely good correlation between even cycle amplitude and SSN(41) found above.

We also evaluated the correlation between the cycle amplitude and the geomagnetic precursor defined above as a function of time lag counted from the sunspot minimum preceding the cycle. That is, instead of taking the maximum ${\mathrm{aa}}_{\mathrm{Sep}-\mathrm{Oct}}^3$ in a 3-year window ending at the cycle minimum the window was shifted backward in time by varying lags. The resulting correlations are shown in Figure 9 separately for odd (blue) and even (red) cycles. One can see that the highest correlations for odd cycles are indeed found up to 16 months before the cycle minimum. There is a drop in correlation between 16 and 55 months before minimum but the correlation then rises again and remains rather higher (about 0.8) between 55 and 100 months before the SSN minimum. Note however, that this does not imply a longer lead time for prediction, since the timing of the SSN minimum is still needed to determine the precursor. Figure 9 also shows that the correlation for even cycles is clearly much lower than for the odd cycles at practically all lags. However, because of the somewhat limited number of data points the statistical 95% confidence limits for the correlations are rather large. Even though there are indications for systematically higher correlation for the odd cycles, there is no firm evidence that the correlations are statistically significantly different for the even and odd cycles.

Figure 9 Correlation coefficient between cycle amplitude A and lagged maximum monthly September–October aa3 within a 3 year window. Blue (red) curve indicates the odd (even cycles) and the correspondingly colored regions indicate the 95% confidence limits for the correlation coefficients.

The specific recipe for calculating the geomagnetic activity-based precursor for odd cycles was found by trial and error, manually testing a few tens of different (rather random) possibilities. This raises the question of whether the found precursor is statistically significantly better than some other, perhaps more common, choice like a mere average aa index in the year of sunspot minimum. We tested the possibility of randomly obtaining a geomagnetic activity-based precursor as good as the one found above. This was done by generating 10⁴ random geomagnetic precursor time series by varying the calendar months from which aa index is taken (2 random calendar months), the length of the time window (randomly selected between 1 and 5 years), and the exponent assigned for the aa index (randomly selected between 1 and 3). In addition, we randomly varied whether we take the maximum of the yearly values within the time window (as we did for the chosen precursor) or the average. For each such randomly generated precursor we calculated the correlation between the precursor and the following cycle amplitude as a function of lag up to 11 years and found the maximum correlation. This procedure simulates the act of randomly choosing a recipe for determining the precursor and finding the maximal correlation over all these lags. Finally, we randomly grouped the remaining correlations into groups of 10 values and determined the maximum correlation in each group. This simulated the act of testing 10 random precursors and choosing the one that gives the maximal correlation at some lag. Now, comparing the correlation coefficient (0.981) for the precursor used above in Figures 8 and 9 to the maximal correlations of these randomly generated precursors in shows that there is only a probability of less than 4.5% to randomly obtain a correlation higher than 0.981. This indicates that the recipe for the geomagnetic precursor we have chosen for the odd cycles is indeed statistically significant by a 95% significance level. We also compared the geomagnetic precursor found above to a more commonly used geomagnetic precursor, which is the annual average of the aa index at the solar minimum, for which the correlation with the following cycle amplitude is 0.834. The correlation for our precursor was 0.981 and the difference of it from 0.834 is statistically highly significant (p-value for the difference is about 0.0007).

The length of the sunspot cycle is also an interesting and important quantity. However, it is difficult to estimate from the predicted sunspot cycle curve because typically the SSN does not drop to zero at the end of the cycle. We studied the association of the cycle length and the four-cycle parameters but found no strong relationships, which would allow one to estimate the length of the cycle based on the parameters of the same cycle. However, we found that the length of the cycle seems to significantly correlate with the ratio of D and C parameters (i.e., D/C) of the previous cycle. This relationship is shown in Figure 10. The correlation coefficient between D/C of the previous cycle and the length of the current cycle is 0.66 (p = 5.7 × 10⁻⁴) and is statistically very significant. There are three cycles (6, 13, and 23) which are clear outliers. Excluding these cycles yields an even higher and more significant correlation of 0.88 (p = 2.6 × 10⁻⁷). We note that these correlations are both significantly higher than, e.g., the correlation between the cycle length and the amplitude of the same cycle (about −0.5) (Wolf, 1861; Petrovay, 2020).

Figure 10 Relationship between the cycle length and the D/C ratio of the previous cycle. The numbers of each data point refer to the current cycle. The blue dots (red squares) correspond to odd (even) cycles and the green dot to the estimated length of the cycle 25. The yellow line indicates the linear fit excluding the three cycles 6, 13, and 23.

The yellow line in Figure 10 shows the linear fit excluding cycles 6, 13, and 23. The equation for this fit is

$$L_i=9.2\left(\pm 0.4\right)+1.5\left(\pm 0.3\right)\frac{D_{i-1}}{C_{i-1}},$$(15)where L_i is the length of cycle i in years and D_i−1 and C_i−1 refer to the parameters of the previous cycle i − 1.

6 Prediction of sunspot cycles

6.1 Cross-validated predictions for past cycles and prediction for cycle 25

Based on the discussion above we can predict the sunspot cycle once the time of the sunspot minimum starting the new cycle is known. We use the 3-parameter description (A, B, and C) for the sunspot cycle of equations (4) and (5). Parameter D has been eliminated using equation (6). For the even sunspot cycles we directly estimate the A, B, and C parameters with SSN(41)² as shown in Figure 7.

For the odd cycles we first estimate the cycle amplitude A from the geomagnetic precursor as shown in Figure 8 and then estimate B and C using the linear relationships for the odd cycles given depicted in Figures 4 and 5, respectively.

Using this approach we predicted past even sunspot cycles starting from cycle 2 and each odd sunspot cycle starting from cycle 11. We predicted each cycle separately with the so-called leave-one-out cross-validation method. This means that when predicting the ith cycle all the fits between different parameters and precursor values discussed above were determined by using data from all other cycles except the ith cycle. Therefore, the fitted relationships between the parameters and precursors (Eqs. (8), (10)–(14)) change from one predicted cycle to the next. This variability in the model parameters together with the residual variability incorporates the total prediction uncertainty of the model. It is important to note that the numerical values in equations (11), (12), and (13) for A, B, and C of even cycles and in equations (14), (8), and (10) for A, B, and C of odd cycles as well as their standard errors correspond to the fitted values when all available sunspot cycles are used in the fit. These values are therefore appropriate for the prediction of sunspot cycles from cycle 25 onward.

An important step in applying the cross-validation method is to estimate the prediction error of the model. Therefore, when proceeding through all the past cycles 1–24 (excluding odd cycles 1, 3, 5, 7, and 9, for which no geomagnetic precursor could be determined) we obtained the residuals of the 13-month smoothed SSN^3/4 values and the corresponding predicted values each time neglecting the cycle to be predicted as

$$r={\mathrm{SSN}}^{3/4}-{\mathrm{SSN}}_{\mathrm{pred}}^{3/4}.$$(16)

The exponent of 3/4 in the above equation was used because the residuals calculated this way were more homoscedastic, i.e., the residual variance was more uniform over different values of SSN compared to regular residuals in linear scale (exponent of 1 in Eq. (16)). We then used the collection of residuals from equation (16) to generate a spread of possible sunspot number predictions for the ith predicted cycle. This was done by bootstrapping 10⁵ values for the residuals (i.e., randomly resampling with replacement from the collection of residuals) for each predicted monthly sunspot number value in the cycle. The residuals were added to the ${\mathrm{SSN}}_{\mathrm{pred}}^{3/4}$ values and the result was then converted back to the linear scale by the transformation

$${\mathrm{SSN}}_k={({\mathrm{SSN}}_{\mathrm{pred}}^{3/4}+r_k)}^{4/3},$$(17)where k = 1, 2, …, 10⁵, r_k indicates the kth bootstrapped residual and SSN_k indicates the kth predicted SSN value for a particular month. These 10⁵ values form an ensemble of SSN predictions for each monthly SSN value in the predicted cycle.

Figure 11 shows the predicted past cycles (green curves) and their uncertainty ranges (blue shading) together with the 13-month smoothed sunspot number (black curve) and the optimal 4-parameter model curves (red curves) as in Figure 2. The figure also shows the predicted curve for cycle 25 (magenta curve) and as a comparison the cycle 25 prediction offered by the Space Weather Prediction Center (https://www.swpc.noaa.gov/products/solar-cycle-progression). One can see that the predicted past cycles are rather close to both the optimal 4-parameter model curves and the 13-month smoothed sunspot number. Particularly noteworthy is the fact that the amplitude of all of the predicted past cycles quite accurately matches the real amplitude of the sunspot cycles. This is true for the very small cycles 6, 12, 14, and 24 and also for the largest cycles, e.g., cycle 19. Furthermore, by definition about 95% of the monthly values of the 13-month smoothed sunspot number are within the 2-standard deviation uncertainty range of the predicted curve. This extremely good performance of the past predictions gives strong reasons to believe that the prediction of the sunspot cycle 25 is reliable as well.

Figure 11 Predicted past sunspot cycles as well as the prediction for cycle 25. The black curve shows the original 13-month smoothed SSN. The red curves show the optimal 4-parameter fits to the SSN cycles. The green curves show the sunspot cycles predicted using the precursor values. The dark and light blue shadings show the 1- and 2-standard deviation ranges based on the ensemble of SSN predictions. The magenta curve shows the prediction for cycle 25 and the dashed purple curve (and yellow uncertainty range) shows the cycle 25 prediction according to the international sunspot prediction panel obtained from the Space Weather Prediction Center.

Figure 12 shows a closeup of cycle 25 prediction in the same format as in Figure 11, but with the addition of the monthly unsmoothed values of SSN also included as the thin black curve. Our predicted curve peaks in May 2024 and the peak value of the cycle is 171 ± 23 (1 standard deviation uncertainty). Using equation (15) and the D and C parameters of cycle 24 the predicted length of cycle 25 is about 9.7 ± 1.9 years (2 standard deviation confidence interval). This indicates that the end of cycle 25 is attained likely in September 2029 with the 2-standard deviation uncertainty range extending from October 2027 to July 2031.

Figure 12 Closeup of the sunspot prediction for cycle 25. The thick (thin) black curve shows the original 13-month smoothed (unsmoothed monthly) SSN. The red curves show the optimal 4-parameter fits to the SSN cycles. The green curves show the sunspot cycles predicted using the precursor values. The dark and light blue shadings show the 1- and 2-standard deviation ranges based on the ensemble of SSN predictions. The magenta curve shows the prediction for cycle 25 and the dashed purple curve (and yellow uncertainty range) shows the cycle 25 prediction according to the international sunspot prediction panel obtained from the Space Weather Prediction Center

At the time of writing this paper, cycle 25 has already begun and the recorded SSN is already significantly higher than the prediction offered by the SWPC. On the other hand, the recorded SSN follows quite closely our predicted curve, which indicates that cycle 25 will be considerably stronger than cycle 24 and roughly the same size as cycle 23. There have also been a range of other predictions for cycle 25. For example, Upton & Hathaway (2018) used advective flux transport modeling to predict that cycle 25 amplitude will be about 95% of cycle 24 amplitude which would make cycle 25 the weakest cycle in 100 years. Bhowmik & Nandy (2018), on the other hand, used solar dynamo modeling to predict that cycle 25 would be slightly stronger (about 14%) than cycle 24. Pesnell & Schatten (2018) used the SODA index (based on a combination of solar polar magnetic field and solar F10.7 index) to predict an amplitude of 135±25 for cycle 25, i.e., slightly larger than cycle 24. Kumar et al. (2021) used various polar field precursors extracted, e.g., from solar magnetograms, to predict a cycle amplitude of 126 ± 3 for cycle 25. Kumar et al. (2022) used the correlation between the rise rate of the polar field and the amplitude of the next cycle to predict an amplitude of 137 ± 23 for cycle 25. Sarp et al. (2018) used a non-linear time series modeling approach to predict a clearly stronger cycle 25 with an amplitude of 154 ± 12 peaking in early 2023. Li et al. (2018) used statistical modeling to reach a similar prediction with a predicted amplitude of 168.5 ± 16.3 and peak of the cycle in October 2024. Both Sarp et al. (2018) and Li et al. (2018) predictions are fairly close, but slightly lower than our prediction. A quite different prediction was given by McIntosh et al. (2020), who used the timing of termination of toroidal bands of solar activity to predict a rather strong cycle 25 with an amplitude of 233 ± 21. However, this prediction was recently revised to an amplitude of 184 ± 17 (McIntosh et al., 2023), which is in agreement with our prediction when considering the ranges of uncertainty. Du (2022b) used the correlation between the cycle rising rate and cycle amplitude to predict cycle 25 based on the 2 years of data from the beginning of the cycle. They predicted the cycle amplitude to be 135.5 ± 33.2, which is somewhat lower than our prediction, although not in complete disagreement given the uncertainty range. Penza et al. (2021) found a correlation between the parameters describing the sunspot cycle shape of even and subsequent odd cycles. Based on this they predicted cycle 25 to be similar or slightly larger than cycle 24. While the studies discussed above are only a subset of all predictions made for cycle 25 it seems that our prediction is clearly above the average in predicted cycle amplitude and also clearly above the SWPC prediction issued by the Solar Cycle Prediction Panel.

6.2 Attempt at predicting the cycle 26

Above we found that for even numbered cycles the 13-month smoothed sunspot number evaluated 41 months before the start of the cycle provides an extremely good estimate for the amplitude and other parameters of the cycle. This leads to an interesting question: How accurately could one use the predicted cycle 25 SSN curve to provide a prediction of cycle 26? Evidently the uncertainty of such a prediction would be rather large but we shall here attempt to make one also for cycle 26. The first step is to evaluate how well the predicted SSN 41 months before the sunspot minima actually correspond to the true values, which were used as a precursor for even cycles.

Figure 13 shows the relationship between the real and predicted 13-month smoothed SSN evaluated 41 months before the minima that end each sunspot cycle. Both odd and even cycles seem to adhere to the same overall linear relationship, which has a high correlation of 0.832 (p = 10⁻⁵). The linear relationship is given by the equation

$${\mathrm{SSN}\left(41\right)}_{\mathrm{real}}=20\left(\pm 10\right)+0.72\left(\pm 0.12\right)\times {\mathrm{SSN}\left(41\right)}_{\mathrm{pred}}.$$(18)

Figure 13 Relationship between the real and predicted 13-month smoothed SSN evaluated 41 months before the minima that end each sunspot cycle. The blue dots (red squares) indicate odd (even) cycles and the yellow line is a linear fit to all cycles.

Ideally there would be a one-to-one relationship between the real and predicted values, but the modeled SSN curves seem to systematically slightly underestimate (overestimate) small (large) SSN(41) values. We can use this fit and its 95% prediction error limits (see Fig. 13) to predict the SSN 41 months before the sunspot minimum that starts cycle 26. The timing of this minimum is determined by the timing of the minimum that starts cycle 25 and the cycle 25 length evaluated above from equation (15). Once the SSN(41) value to be used as a precursor for cycle 26 is known we can use equation (11) to estimate the cycle amplitude.

Because of the uncertainties associated with the length of cycle 25 and the scaling of predicted SSN(41) according to equation (18) we evaluated the spread of possible cycle 26 amplitudes using a Monte Carlo simulation having 10⁴ rounds. In each round, we randomly generated a value for the length of cycle 25 within the range of its uncertainty. This was then used to calculate the timing of 41 months before the end of cycle 25 and the SSN value at that time using the predicted SSN curve for cycle 25 (Fig. 12). This value was then used to calculate a prediction for the cycle amplitude using equation (11). The histogram of the Monte Carlo simulation results is shown in Figure 14. As expected the results indicate a quite a large uncertainty range covering practically all past sunspot cycles. However, despite the large range of uncertainty, some interesting and non-trivial aspects are seen. The median of the results implies that cycle 26 would be even slightly stronger than cycle 25. In fact, based on these results the probability that cycle 26 will be weaker than cycle 25 is only about 19%. The results also imply an even clearer difference to cycle 24, which was the weakest cycle of the last 100 years. The probability that cycle 26 would be weaker than cycle 24 is only 0.8%.

Figure 14 Histogram of predicted cycle 26 amplitudes based on a Monte Carlo simulation with 104 repetitions. For reference the dashed vertical lines show the cycle amplitudes for cycle 19, 24, and predicted 25. The red dashed line indicates the median of the simulated values.

7 Discussion and conclusion

We have studied here the prediction of sunspot cycles paying particular emphasis on the differences between odd and even numbered cycles. Our approach was to first parameterize the sunspot cycle as an asymmetric 4-parameter Gaussian function, which was found to represent each sunspot cycle extremely well. Several of the parameters were found to correlate, which allowed us to simplify the representation of the sunspot cycles. Interestingly it was found that the correlations between the cycle parameters were significantly stronger in even cycles compared to the odd ones. For example, the well-known Waldmeier effect that cycle amplitude and its rise time are inversely correlated was more strictly valid for even cycles. A similar finding for the Waldmeier effect was also reported, e.g., by Du (2022a) and Dikpati et al. (2008). This result implies that for some reason the even-numbered cycles possess a smaller dimensionality than odd-numbered cycles. Such a systematic difference between even and odd cycles may imply a connection to the more fundamental 22-year magnetic Hale cycle of the Sun. Du (2022a) also found significant differences in the Waldmeier effect between the two hemispheres with the southern hemisphere displaying the normal Waldmeier effect, which is stronger in even cycles, while the northern hemisphere displayed an inverse-Waldmeier effect, which is stronger in odd cycles. This implies that there may also be a connection between the Waldmeier effect and the hemispheric asymmetry in the sunspot number.

Our approach to the prediction of sunspot cycles is based on statistical precursors for the cycle parameters. Following Cameron and Schüssler (2007) we found that the 13-month smoothed sunspot number evaluated 41 months before the sunspot minimum that begins a cycle, SSN(41), is on average a fair predictor for the next cycle when considering all sunspot cycles. However, we found that this relation is much stronger for even cycles than for odd cycles, for which the SSN(41) could not be very accurately used as a predictor. The advantage of SSN(41) precursor is that we can predict for even cycles all the cycle parameters (which correlate well with the amplitude) quite accurately. Cameron & Schüssler (2007) explained the tendency of past SSN to correlate with the next cycle amplitude as a result of overlapping sunspot cycles. While the past cycle is still declining the new cycle begins. Because of the Waldmeier effect the stronger the new cycle will be the faster it will rise to the maximum and the earlier the intersection time of the old and the new cycle is attained. The earlier the intersection time, the higher the SSN 41 months earlier in the declining phase of the previous cycle. The fact that this connection is found to be more strictly valid for even cycles arises because of the fact that the Waldmeier effect is also tighter in even cycles. The question of why the Waldmeier effect is different in even and odd cycles is still open and should be investigated in future studies.

For odd cycles we found that the maximum September–October geomagnetic aa index within 3 years preceding the sunspot minimum that begins the cycle is an extremely good precursor for the cycle amplitude. It has been long recognized that geomagnetic activity close to solar minimum reflects the strength of the open solar flux, which is at these times connected to the poloidal magnetic flux extending from the Sun’s polar coronal holes. It is curious, however, that the best precursor for odd cycles was found to be connected to geomagnetic activity close to the fall equinox. Some other past studies have used geomagnetic activity averaged in other ways over different periods of time and found good correlations to the next cycle amplitude. While our September–October aa precursor is probably not statistically significantly better than some of these other measures there might be a physical reason for the preference of September–October season.

It is known that geomagnetic activity has strong seasonal variation, which is largely due to the Russell-McPherron (RMP) effect (Russell & McPherron, 1973; Lockwood et al., 2020). The RMP effect describes the fact that the interplanetary magnetic field (IMF), which is oriented along the solar equatorial plane, projects onto the Z-axis of the GSM coordinate system close to the fall and spring equinoxes and thereby may enhance energy input from the solar wind into the magnetosphere/ionosphere system. The RMP effect is also dependent on the polarity of the IMF so that during the fall (spring) equinox the IMF pointing away from (toward) the Sun enhances solar wind energy input and therefore leads to larger geomagnetic activity as well. Often both polarities of the IMF are seen within a solar rotation but especially close to solar minima the heliospheric current sheet is often flatter, which allows Earth to be more exposed to the magnetic polarity connected to Sun’s northern (southern) pole in fall (spring) due to the Earth’s changing heliographic latitude over the course of the year (Rosenberg & Coleman Jr., 1969). According to the 22-year Hale cycle the northern solar pole has a positive magnetic polarity close to sunspot minima preceding odd cycles, which leads to a dominance of the away sector close to fall equinoxes (Hiltula & Mursula, 2007; Vokhmyanin & Ponyavin, 2012). Furthermore, it has been shown that the heliospheric current sheet is systematically tilted southward in the declining phase (Hiltula & Mursula, 2006), which further enhances the dominance of the away sector in fall and decreases the dominance of the toward sector in spring prior to odd cycles. Therefore, it is expected that the geomagnetic activity portrayed by the aa index is most sensitively proportional to the strength of the IMF (i.e., to open solar magnetic flux connected to solar polar field) at these times. A detailed confirmation of this interpretation is warranted, but out of the scope of this study.

In addition to cycle amplitude and other parameters we found a curious statistical relation between the length of the sunspot cycle and the ratio of D (cycle asymmetry) and C (time scale of rising phase) of the preceding cycle. The fact that such properties of the preceding cycle somehow are connected to the length of the next cycle again highlights the fundamental 22-year Hale cyclicity of solar magnetism. While there is no physical explanation for this relationship at the moment it can be statistically used to estimate the cycle length perhaps a bit more accurately than previous metrics (Petrovay, 2020).

Using the found precursors we used cross-validation to test their prediction accuracy by predicting the past solar cycles. For all past cycles the predictions were very close to the real sunspot cycles thereby giving strong confidence that the prediction of future cycles would be equally successful. We proceeded to predict the odd cycle 25 and found that its amplitude will be 171 ± 23, thus about 1.6 times stronger than cycle 24. There are already clear indications that our prediction closely follows the progression of cycle 25 which has already started at the time of writing this paper. It is also noteworthy that the prediction issued by the Solar Cycle 25 Prediction Panel at the SWPC suggests cycle 25 to be similar to cycle 24, which is already clearly below the current sunspot levels.

Using the predicted cycle 25 and the fact that SSN(41) could be used as a predictor for the even cycles we provided a rough prediction also for cycle 26. As expected the uncertainty range of the prediction was rather large, but based on the results it seems rather likely that cycle 26 will be stronger than both cycles 24 and 25. Therefore, we find no evidence for an imminent drop of solar activity to a grand solar minimum as suggested by several past studies (e.g. Abreu et al., 2008; Owens et al., 2011; Lockwood et al., 2011). Overall these results display the capability to predict even and odd cycles using different precursors with rather high accuracy. The results also clearly indicate a connection between odd and even pairs of sunspot cycles and highlight the 22-year Hale cyclicity. Accordingly, the Hale cyclicity should be considered more carefully also by more physically motivated dynamo and flux transport prediction models of solar activity.

Acknowledgments

We acknowledge the financial support by the Academy of Finland to the PROSPECT (project no. 321440). The sunspot number data was obtained from World Data Center SILSO, Royal Observatory of Belgium, Brussels (https://www.sidc.be/SILSO/home). The editor thanks three anonymous reviewers for their assistance in evaluating this paper.

References

All Figures

	Figure 1 Annual averages of homogenized aa index (blue), scaled raw aa index (red) (shown since year 1980) and scaled Helsinki Ak(D) values (shown for years 1844–1879).
In the text

	Figure 2 Time series of 13-month smoothed sunspot number (black), parameterized fits for each sunspot cycle (red curves), and the fit by Hathaway et al. (1994) for comparison (blue curves). The sunspot cycle numbers are indicated on the plot as well. The top panel shows cycles 1–8, middle panel cycles 9–16, bottom panel cycles 17–24.
In the text

	Figure 3 Relationship between the C and D−1/2 parameters of the 4-parameter fits. Blue dots (red squares) indicate even (odd) numbered cycles. The cycle numbers have been further indicated beside all points. The yellow line depicts the linear fit to all the cycles having a correlation coefficient of 0.94 (p-value is 10−11).
In the text

	Figure 4 Relationship between the A and B parameters of the 3-parameter fits. Blue dots (red squares) indicate even (odd) numbered cycles. The cycle numbers have been further indicated beside all points. The blue and red lines depict the linear fit to odd and even cycles respectively.
In the text

	Figure 5 Relationship between the A and C parameters of the 3-parameter fits. Blue dots (red squares) indicate even (odd) numbered cycles. The cycle numbers have been further indicated beside all points. The blue and red lines depict the linear fits to odd and even cycles respectively.
In the text

Figure 6 Correlation coefficient between cycle amplitude A and lagged 13-month smoothed sunspot number (SSN2) as a function of the time lag in months before the SSN minimum. Blue (red) curve indicates the odd (even) cycles and the correspondingly colored regions indicate the 95% confidence limits for the correlation coefficients. The green dot indicates the maximum correlation at 41 months before the SSN minimum for the even cycles.

In the text

	Figure 7 Dependence of A, B, and C parameters on 13-month smoothed sunspot number 41 months before the sunspot minimum (SSN(41)2). Blue dots (red squares) indicate odd (even) numbered cycles. The cycle numbers have been further indicated beside all points. The red lines depict the linear fits to even numbered cycles.
In the text

	Figure 8 Relationship between cycle amplitude A and maximum September–October aa3 from a 3-year window before SSN minimum. The blue points indicate odd cycles and red squares the even cycles. The blue line indicates the linear fit to the odd cycles. The green circle indicates the estimated amplitude for cycle 25.
In the text

	Figure 9 Correlation coefficient between cycle amplitude A and lagged maximum monthly September–October aa3 within a 3 year window. Blue (red) curve indicates the odd (even cycles) and the correspondingly colored regions indicate the 95% confidence limits for the correlation coefficients.
In the text

	Figure 10 Relationship between the cycle length and the D/C ratio of the previous cycle. The numbers of each data point refer to the current cycle. The blue dots (red squares) correspond to odd (even) cycles and the green dot to the estimated length of the cycle 25. The yellow line indicates the linear fit excluding the three cycles 6, 13, and 23.
In the text

Figure 11 Predicted past sunspot cycles as well as the prediction for cycle 25. The black curve shows the original 13-month smoothed SSN. The red curves show the optimal 4-parameter fits to the SSN cycles. The green curves show the sunspot cycles predicted using the precursor values. The dark and light blue shadings show the 1- and 2-standard deviation ranges based on the ensemble of SSN predictions. The magenta curve shows the prediction for cycle 25 and the dashed purple curve (and yellow uncertainty range) shows the cycle 25 prediction according to the international sunspot prediction panel obtained from the Space Weather Prediction Center.

In the text

Figure 12 Closeup of the sunspot prediction for cycle 25. The thick (thin) black curve shows the original 13-month smoothed (unsmoothed monthly) SSN. The red curves show the optimal 4-parameter fits to the SSN cycles. The green curves show the sunspot cycles predicted using the precursor values. The dark and light blue shadings show the 1- and 2-standard deviation ranges based on the ensemble of SSN predictions. The magenta curve shows the prediction for cycle 25 and the dashed purple curve (and yellow uncertainty range) shows the cycle 25 prediction according to the international sunspot prediction panel obtained from the Space Weather Prediction Center

In the text

	Figure 13 Relationship between the real and predicted 13-month smoothed SSN evaluated 41 months before the minima that end each sunspot cycle. The blue dots (red squares) indicate odd (even) cycles and the yellow line is a linear fit to all cycles.
In the text

	Figure 14 Histogram of predicted cycle 26 amplitudes based on a Monte Carlo simulation with 104 repetitions. For reference the dashed vertical lines show the cycle amplitudes for cycle 19, 24, and predicted 25. The red dashed line indicates the median of the simulated values.
In the text