Solar cycle prediction

TL;DR

This paper reviews various methods for predicting solar cycle amplitude and timing, highlighting the superior performance of precursor methods like polar field predictions and the emerging role of model-based forecasts.

Contribution

It provides a comprehensive review of solar cycle prediction techniques, comparing their performance and discussing recent advances and future prospects.

Findings

Precursor methods outperform extrapolation methods in accuracy.

Polar field precursor has been consistently reliable.

Model-based forecasts are gaining confidence despite limited history.

Abstract

A review of solar cycle prediction methods and their performance is given, including early forecasts for cycle 25. The review focuses on those aspects of the solar cycle prediction problem that have a bearing on dynamo theory. The scope of the review is further restricted to the issue of predicting the amplitude (and optionally the epoch) of an upcoming solar maximum no later than right after the start of the given cycle. In their overall performance during the course of the last few solar cycles, precursor methods have clearly been superior to extrapolation methods. One method that has yielded predictions consistently in the right range during the past few solar cycles is the polar field precursor. Nevertheless, some extrapolation methods may still be worth further study. Model based forecasts are quickly coming into their own, and, despite not having a long proven record, their predictions are received with increasing confidence by the community.

Tables

Table 1. Table 1: Predictors based on magnetic field measurements and their forecasts. Forecast limits have been calculated adding/subtracting the scatter (as calculated nominally from the 4 data points) to/from the value given by the fit.

Cycle	Amplitude	Previous WSO	WSO field	Previous	Dipole coef. $D$
	(SSN v2)	field maximum	at minimum	$D$ maximum	at minimum
21	232.9	$\left[1.07\right]$	$1.05$	$4.19$	$4.10$
22	212.5	$1.31$	$1.28$	$4.23$	$3.98$
23	180.3	$1.13$	$1.00$	$3.96$	$3.01$
24	116.4	$0.63$	$0.52$	$1.95$	$1.33$
25		$>0.72$		$>1.93$
Fit coef.		$177.4$	$188.8$	$51.3$	$57.5$
Scatter		$25.9$	$24.8$	$16.8$	$21.9$
Cycle 25
forecast		$>102$		$>82$

Table 2. Table 2: A selection of early forecasts for Cycle 25

Category	Minimum	Maximum	Peak amplitude	Reference
Internal precursors	2019.9	2023.8	175 [154–202]	Li et al (2015)
External precursor
polar precursor			117 $\pm$ 15	Table 1 here
polar precursor			136 $\pm$ 48	Pesnell and Schatten (2018)
helicity			117	Hawkes and Berger (2018)
SoDA		2025.2 $\pm$ 1.5	120 $\pm$ 39	based on Pesnell and Schatten (2018)
rush-to-the-poles	2019.4	2024.8	130	Petrovay et al (2018)
Model-based: SFT
SFT			124 $\pm$ 31	Jiang et al (2018)
AFT	2020.9		110	Upton and Hathaway (2018)
Model-based: dynamo
$2\times 2$D	2020.5 $\pm$ 0.12	2027.2 $\pm$ 1.0	89${}_{-14}{}^{}{}_{}{}^{+29}$	Labonville et al (2019)
Truncated	2019–20	2024 $\pm$ 1	90 $\pm$ 15	Kitiashvili (2016)
Spectral
wavelet decomposition tree		2023.4	132	Rigozo et al (2011)
Attractor analysis
simplex projection analysis		2024.0 $\pm$ 0.6	103 $\pm$ 25	Singh and Bhargawa (2017)
simplex proj./time-delay		2023.2 $\pm$ 1.1	154 $\pm$ 12	Sarp et al (2018)
Neural networks
neuro-fuzzy		2022	90.7 $\pm$ 8	Attia et al (2013)
spatiotemporal		2022–23	57 $\pm$ 17	Covas et al (2019)
Cycle 24 [comparison]	2008.9	2014.3	116

Equations

$$R_Z=k(10g+f),$$

$$R=\frac{1}{24}\left(R_{\mathrm{m},-6}+2\sum _{i=-5}^{i=5}{R}_{\mathrm{m},i}+R_{\mathrm{m},6}\right),$$

$$R_{\pm }=100{\left(R_B/83\right)}^{3/2}$$

$$⟨R_{\max }{⟩}_{\mathrm{G}}^{(n)}=\frac{1}{8}\left(R_{\max }^{(n-2)}+2R_{\max }^{(n-1)}+2R_{\max }^{(n)}+2R_{\max }^{(n+1)}+R_{\max }^{(n+2)}\right).$$

$$\log R_{\max }^{(n)}=C_1+C_{2}f$$

$$R_{\max }^{(n)}=C_1-C_2\left(t_{\max }^{(n)}-t_{\min }^{(n)}\right),$$

$$R_{\max }^{(n)}\propto \left[1/\left(t_{\max }^{(n)}-t_{\min }^{(n)}\right)\right].$$

$$\log R_{\max }^{(n)}=C_1-C_2\left(t_{\max }^{(n)}-t_{\min }^{(n)}\right),$$

$$\log R_{\max }^{(n)}\propto \left[1/\left(t_{\max }^{(n)}-t_{\min }^{(n)}\right)\right].$$

$$R_{\max }=114.3+6.1R_{\min }.$$

$$R_{\max }=79+1.52R(t_{\min }-3\mathrm{yr}).$$

$$D(t)=\frac{3}{2}{\int }_0^{\pi }\overline{B}(\theta ,t)\cos \theta \sin \theta \mathrm{d}\theta .$$

$$D_{NS}(t)$$

$$\frac{\partial B}{\partial t}$$

$$\delta D_{\mathrm{BMR}}=\frac{3}{4\pi R_{\odot }^2}Fd\sin \alpha \sin \theta ,$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{U}\times \mathbf{B}+\alpha \mathbf{B})-\nabla \times ({\eta }_T\times \nabla \mathbf{B}).$$

$$\mathbf{U}={\mathbf{U}}_c+r\sin \theta {\Omega }_0{\mathbf{e}}_{\varphi },$$

$$\mathbf{\Omega }=r\sin \theta \nabla {\Omega }_0,\Omega =-\mathrm{sgn}\frac{\mathrm{d}{\Omega }_0}{\mathrm{d}r}\cdot |\mathbf{\Omega }|,$$

$$\frac{\partial A}{\partial t}=\alpha B-({\mathbf{U}}_c\cdot \nabla )A-(\nabla \cdot {\mathbf{U}}_c)A+{\eta }_T{\nabla }^{2}A,$$

$$\frac{\partial B}{\partial t}=\Omega \frac{\partial A}{\partial x}-({\mathbf{U}}_c\cdot \nabla )B-(\nabla \cdot {\mathbf{U}}_c)B+{\eta }_T{\nabla }^{2}B,$$

$$\dot{A}=\alpha B-A/\tau ,$$

$$\dot{B}=(\Omega /L)A-B/\tau ,$$

$$\ddot{B}=\frac{D-1}{{\tau }^2}B-\frac{2}{\tau }\dot{B},$$

$$\ddot{B}=\frac{D_0-1}{{\tau }^2}B-\frac{2}{\tau }\dot{B}-\frac{D_0-1}{{\tau }^2}Bf(B).$$

$$\ddot{B}=-{\omega }^{2}B-\mu (\xi B^2-1)\dot{B}-\gamma B^3$$

$$\ddot{\xi }=-\xi +\mu (1-{\xi }^2)\dot{\xi },$$

$$R_n=R_0+\sum _{i=1}^p{c}_{n-i}{R}_{n-i}+{ϵ}_n,$$

$$R_n=R_0+\sum _{i=1}^p{c}_{n-i}{R}_{n-i}+{ϵ}_n+\sum _{i=1}^q{d}_{n-i}{ϵ}_{n-i}.$$

Full text

Prediction methods form three main groups. Precursor methods rely on the value of some measure of solar activity or magnetism at a specified time to predict the amplitude of the following solar maximum. The choice of a good precursor often implies considerable physical insight: indeed, it has become increasingly clear that the transition from purely empirical precursors to model-based methods is continuous. Model-based approaches can be further divided into two groups: predictions based on surface flux transport models and on consistent dynamo models.

The implicit assumption of precursor methods is that each numbered solar cycle is a consistent unit in itself, while solar activity seems to consist of a series of much less tightly intercorrelated individual cycles. Extrapolation methods, in contrast, are based on the premise that the physical process giving rise to the sunspot number record is statistically homogeneous, i.e., the mathematical regularities underlying its variations are the same at any point of time, and therefore it lends itself to analysis and forecasting by time series methods.

In their overall performance during the course of the last few solar cycles, precursor methods have clearly been superior to extrapolation methods. One method that has yielded predictions consistently in the right range during the past few solar cycles is the polar field precursor. Nevertheless, some extrapolation methods may still be worth further study. Model based forecasts are quickly coming into their own, and, despite not having a long proven record, their predictions are received with increasing confidence by the community.

††journal: Living Reviews in Solar Physics

1. 1 Introduction

2. 1.0.1 What’s new in this edition?

3. 1.1 The sunspot number (SSN)

4. 1.1.1 Version 1.0

5. 1.1.2 Revision

6. 1.1.3 Alternating series and nonlinear transforms

7. 1.2 Other indicators of solar activity

8. 1.2.1 Group sunspot number (GSN)

9. 1.2.2 Other sunspot data

10. 1.2.3 Other activity indices

11. 1.3 The solar cycle and its variation

12. 1.3.1 Secular activity variations

13. 1.3.2 Does the Sun have a long term memory?

14. 1.3.3 Waldmeier effect and amplitude–frequency correlation

15. 1.4 Approaches to solar cycle prediction

16. 2 Precursor methods

17. 2.1 Cycle parameters as precursors and the Waldmeier effect

18. 2.2 External empirical precursors

19. 2.3 Polar precursor

20. 2.3.1 Polar magnetic field data

21. 2.3.2 Proxy reconstructions of the polar magnetic field

22. 2.3.3 Extending the range of the polar precursor: early forecasts for Cycle 25

23. 2.4 Geomagnetic and interplanetary precursors

24. 2.5 The quest for a precursor of the polar precursor

25. 2.5.1 Photospheric flow variations

26. 2.5.2 Active region tilts

27. 2.5.3 Alternatives to tilt quenching: inflows and latitude quenching

28. 3 Model-based predictions

29. 3.1 Surface flux transport models

30. 3.1.1 Parameter optimization

31. 3.1.2 Nonlinear feedback effects

32. 3.1.3 Extending the polar precursor

33. 3.1.4 The importance of fluctuations: rogue active regions

34. 3.1.5 Explaining the end of the Modern Maximum

35. 3.1.6 Forecasting Cycle 25

36. 3.2 The solar dynamo: a brief summary of current models

37. 3.3 Is model-based cycle prediction feasible?

38. 3.4 Explicit models

39. 3.4.1 Axisymmetric models

40. 3.4.2 Nonaxisymmetric models

41. 3.5 Truncated models

42. 3.6 The Sun as an oscillator

43. 4 Extrapolation methods

44. 4.1 Linear regression

45. 4.2 Spectral methods

46. 4.2.1 The 11-year cycle and its harmonics

47. 4.2.2 The even–odd (a.k.a. Gnevyshev–Ohl) rule

48. 4.2.3 The Gleissberg cycle

49. 4.2.4 Supersecular cycles

50. 4.3 Nonlinear methods

51. 4.3.1 Attractor analysis and phase space reconstruction: the pros …

52. 4.3.2 …the cons …

53. 4.3.3 …and the upshot

54. 4.3.4 Neural networks

55. 5 Summary evaluation

56. 6 Forecasts for Cycle 25 [NEW!]

1 Introduction

Solar cycle prediction is an extremely extensive topic, covering a very wide variety of proposed prediction methods and prediction attempts on many different timescales, ranging from short term (month–year) forecasts of the runoff of the ongoing solar cycle to predictions of long term changes in solar activity on centennial or even millennial scales. As early as 1963, Vitinsky published a whole monograph on the subject, later updated and extended (Vitinsky 1963, 1973). More recent overviews of the field or aspects of it include Hathaway (2009), Kane (2001), Pesnell (2008), and the first edition of this review (Petrovay, 2010). In order to narrow down the scope of the review, we constrain our field of interest in two important respects.

Firstly, instead of attempting to give a general review of all prediction methods suggested or citing all the papers with forecasts, here we will focus on those aspects of the solar cycle prediction problem that have a bearing on dynamo theory. We will thus discuss in more detail empirical methods that, independently of their success rate, have the potential of shedding some light on the physical mechanism underlying the solar cycle, as well as the prediction attempts based on solar dynamo models.

Secondly, we will here only be concerned with the issue of predicting the amplitude (and optionally the epoch) of an upcoming solar maximum no later than right after the start of the given cycle. This emphasis is also motivated by the present surge of interest in precisely this topic, prompted by the unusually long and deep recent solar minimum and by sharply conflicting forecasts for the maximum of the incipient solar Cycle 24.

As we will see, significant doubts arise both from the theoretical and observational side as to what extent such a prediction is possible at all (especially before the time of the minimum has become known). Nevertheless, no matter how shaky their theoretical and empirical backgrounds may be, forecasts must be attempted. Making verifiable or falsifiable predictions is obviously the core of the scientific method in general; but there is also a more imperative urge in the case of solar cycle prediction. Being the prime determinant of space weather and space climate, solar activity clearly has enormous technical, scientific, and financial impact on activities ranging from space exploration to civil aviation and everyday communication. Political and economic decision makers expect the solar community to provide them with forecasts on which feasibility and profitability calculations can be based. Acknowledging this need, around the time of solar minimum the Space Weather Prediction Center of the US National Weather Service does present annually or semiannually updated “official” predictions of the upcoming sunspot maximum, emitted by a Solar Cycle Prediction Panel of experts. The unusual lack of consensus in the early meetings of this panel during the previous minimum (SWPC 2009), as well as the concurrent, more frequently updated but wildly varying predictions of a NASA MSFC team (MSFC 2009) made evident the deficiencies of prediction techniques available at the time. In view of this, Cycle 24 provided us with some crucial new insights into the physical mechanisms underlying cyclic solar activity. In preparation to convene the Solar Cycle 25 Prediction Panel a new call for predictions was issued in January 2019 (Biesecker and Upton 2019).

While a number of indicators of solar activity exist, by far the most commonly employed is still the smoothed relative sunspot number $R$; the “Holy Grail” of sunspot cycle prediction attempts is to get $R$ right for the next maximum. We, therefore, start by briefly introducing the sunspot number and inspecting its known record. Then, in Sects. 2, 4, and 3 we discuss the most widely employed methods of cycle predictions. Section 5 presents a summary evaluation of the past performance of different forecasting methods, while Sect. 6 finally collects some early forecasts for Cycle 25 derived by various approaches.

1.0.1 What’s new in this edition?

For readers familiar with the 1st edition of this review, who would prefer to go through the “new stuff” only, here I briefly list the new or thoroughly rewritten sections.

The revision of the sunspot number series that took place in 2015 is one topic that had to be discussed in detail. The new Sect. 1.1.2 is devoted to this subject but other parts of Sects. 1.1 and 1.2 have also been subjected to a major revision.

For reasons explained in Sect. 1.4, the overall structure of the review has been given a major overhaul: the section on extrapolation methods has now been placed after the section on model-based approaches.

Researches into the origins of the sudden change in the behaviour of our Sun from Cycle 23 to 24 led to important, although still poorly understood realizations, now discussed in a dedicated subsection (Sect. 2.2). And the stellar rise in popularity of the polar precursor led to such an amount of exciting new research that Sect. 2.3, discussing this method, had to be completely rewritten and expanded; furthermore, it gave rise to another section on “The quest for a precursor of the polar precursor” (Sect. 2.5), containing mostly new or thoroughly updated material.

The section on model-based predictions now includes a presentation of the approach based on surface flux transport models (Sect. 3.1). In the field of dynamo-based cycle prediction the major novelty in this decade was the development of nonaxisymmetric models capable to account for the emergence of individual active regions: a new subsection is now devoted to this topic (Sect. 3.4.2).

Finally, the updated summary evaluation and the overview of early predictions for Cycle 25 obviously cover mostly new results (Sect. 5 and 6).

1.1 The sunspot number (SSN)

Despite its somewhat arbitrary construction, the series of relative sunspot numbers constitutes the longest homogeneous global indicator of solar activity determined by direct solar observations and by methods that were, until recently, perceived to be carefully controlled. For this reason, their use is still predominant in studies of solar activity variation.

As aptly noted by Clette et al (2014), until recently the sunspot number series was “assumed to be carved in stone, i.e., it was considered largely as a homogeneous, well-understood and thus immutable data set. This feeling was probably reinforced by the stately process through which it was produced by a single expert center at the Zürich Observatory during 131 years.”

This perception has now been shattered by the major revision of the official sunspot number series that took place in 2015, and opened the way to further periodic revisions.

In what follows, the original series, the revision, and transformed versions of the series will be discussed in turn.

1.1.1 Version 1.0

As defined originally by Wolf (1850), the relative sunspot number is

[TABLE]

where $g$ is the number of sunspot groups (including solitary spots), $f$ is the total number of all spots visible on the solar disc, while $k$ is a correction factor depending on a variety of circumstances, such as instrument parameters, observatory location, and details of the counting method. Wolf, who decided to count each spot only once and did not count the smallest spots, used $k=1$. He also introduced a hierarchical system for the determination of $R_{Z}$ where data from a list of auxiliary observers were used for days when the primary observer (Zürich) could not provide a value.

The counting system employed in Zürich was changed by Wolf’s successors to count even the smallest spots, attributing a higher weight (i.e., $f>1$) to spots with a penumbra, depending on their size and umbral structure. (The exact details and timing of these changes are incompletely documented and controversial, see discussion in the next subsection.) As the changes in the counting and the regular use of a larger telescope naturally resulted in higher values, the Zürich correction factor was set to $k=0.6$ for subsequent determinations of $R_{Z}$ to ensure continuity with Wolf’s work. (Waldmeier 1961; see also Izenman 1983, Kopecký et al 1980, Hoyt and Schatten 1998, Friedli 2016 for further discussions on the determination of $R_{Z}$).

In addition to introducing the relative sunspot number, Wolf (1861) also used earlier observational records available to him to reconstruct its monthly mean values since 1749. In this way, he reconstructed 11-year sunspot cycles back to that date: hence, the now universally used numbering of solar cycles starts with the first complete cycle in the monthly $R_{Z}$ series. In a later work he also determined annual mean values for each calendar year going back to 1700. This reconstruction and calibration work took place in several steps, so the $R_{Z}$ record was very much a project in the making until the end of the 19th century (see Clette et al 2014). It was only from the early 20th century that the series came to be regarded as “carved in stone”.

In 1981, the observatory responsible for the official determination of the sunspot number changed from Zürich to the Royal Observatory of Belgium in Brussels. The website of the SIDC111http://sidc.oma.be (originally Sunspot Index Data Center, recently renamed Solar Influences Data Analysis Center) is now the most authoritative source of archive sunspot number data. The department of SIDC formally responsible for the sunspot number series is WDC-SILSO (World Data Centre for Sunspot Index and Long-term Solar Observations). It has become customary to denote the original Zürich series with $R_{Z}$ (“the Zürich sunspot number”), and its continuation by the SIDC from 1981 to 2015 with $R_{i}$ (International Sunspot Number, ISN). The new, revised series is conventionally denoted by $S_{N}$.

It must be kept in mind that since the middle of the 20th century, the sunspot number is also regularly determined by other institutions: the most widely used such variants are informally known as the American sunspot number (collected by AAVSO and available from the National Geophysical Data Center222http://www.ngdc.noaa.gov/ngdc.html) and the Kislovodsk Sunspot Number (available from the web page of the Kislovodsk Mountain Astronomical Station of Pulkovo Observatory333http://en.solarstation.ru/).

Given that $R_{Z}$ is subject to large fluctuations on a time scale of days to months, it has become customary to use annual mean values for the study of longer term activity changes. To get rid of the arbitrariness of calendar years, the standard practice444Alternative proposals have been put forward by Muñoz-Jaramillo et al (2012) and Podladchikova et al (2017); who pointed out that Eq. (2) does a poor job of filtering out high-frequency variations, and use better or even optimized weight factors for the $R_{\mathrm{m},i}$’s instead. is to use 13-month boxcar averages of the monthly averaged sunspot numbers, wherein the first and last months are given half the weight of other months:

[TABLE]

$R_{\mathrm{m},i}$ being the mean monthly value of the daily sunspot number values for ith calendar month counted from the present month. It is this running mean $R$ that we will simply call “the sunspot number” throughout this review and what forms the basis of most discussions of solar cycle variations and their predictions.

In what follows, $R_{\mathrm{max}}^{(n)}$ and $R_{\mathrm{min}}^{(n)}$ will refer to the maximum and minimum value of $R$ in cycle $n$ (the minimum being the one that starts the cycle). Similarly, $t_{\mathrm{max}}^{(n)}$ and $t_{\mathrm{min}}^{(n)}$ will denote the epochs when $R$ takes these extrema.

1.1.2 Revision

The process that led to the 2015 revision was started by Leif Svalgaard (Stanford) who pointed out a number of inhomogeneities in the series, rooted in changes in the base data and processing techniques. Starting from 2011, at Svalgaard’s initiative, a series of four workshops on the sunspot number were held by the community involved. The 2015 revision is the result of this process. The motivation for and the detailed process of the revision was described by Clette et al (2014) and Clette and Lefèvre (2016) and discussed in a number of papers in a topical issue of Solar Physics (vol. 291, issue 9–10; Clette et al 2016a).

The revision included dropping the $k=0.6$ scaling factor traditionally applied to the Zürich data, so all values increased by a factor of $\textstyle\frac{5}{3}$. In addition to this trivial rescaling and some other minor changes, three major corrections were implemented.

(a) The Locarno drift after 1981. The determination of the International Sunspot Number $R_{i}$ by the SIDC did not follow Wolf’s hierarchical system, taking into account observations from all network stations and only dropping outliers. Nevertheless, in order to ensure continuity, the Locarno solar observatory (Zürich’s successor) still had a special role as a pilot station, all other observers being calibrated to Locarno’s scale. A slow time-varying drift in the Locarno data came to light during the revision process and has been corrected in the new series. This change is apparently uncontroversial and was made with the full consensus of all actors involved (Clette et al 2016b).

(b) The “Waldmeier jump” from 1947. Plotting the original Zürich sunspot numbers against other sunspot-related indices such as sunspot areas or group numbers, or even against non-weighted sunspot numbers determined by non-Zürich observers, a jump was discovered which was suggested to originate in the introduction of a new weighting method (in use in Locarno until the present day) by the new Zürich director, Max Waldmeier and his largely new staff. Under current solar conditions the weighting results in a 15–20% inflation of the sunspot numbers (Svalgaard et al 2017). This assumed inflation of the series has now been corrected.

(c) The Schwabe–Wolf transition (1849–64) The upward correction of 14% applied in this period relies primarily on a comparison of the original sunspot number series with group sunspot numbers (the result being apparently insensitive to which group number series is used). The presumed cause of the discrepancy is that in this period the sunspot number was determined by Wolf using small portable telescopes, while Schwabe also continued his observations. For days not covered by his own observations Wolf used Schwabe’s data without marking these out. It was only in 1861 that, upon cross-correlating their data Wolf determined a correction factor $k=1.25$ for Schwabe, which he also applied retrospectively to the pre-1849 observations of Schwabe (and, by inference, of earlier observers calibrated to Schwabe). However, the correction factor was apparently not considered for his own observations (mixed with Schwabe’s) in the period before the early 1860s.

The revised series, introduced from 1 July 2015, is now considered version 2.0 of the sunspot number series. Further corrections, with proper version tracking, are expected as early data may contain other inconsistencies, and the corrections applied in v2.0 were somewhat crude. In particular, recomputation of the whole series from observational data, wherever available, is planned. The process has now been placed under the ægis of the IAU, with a dedicated Working Group “Coordination of Synoptic Observations of the Sun” focusing on the validation and accreditation process of new SSN versions.

A further possible bias in the series that remains to be corrected may concern the counting of sunspot groups. While in earlier parts of the series physical closeness of spots was considered a sufficient criterion, since the mid-20th century evolutionary information is also taken into account, sometimes resulting in the division into several groups of what would have been considered as a single group by early observers. Svalgaard et al (2017) estimate that this effect may have inflated Waldmeier’s sunspot numbers by 4–5% relative to earlier counts, while the effect on the late 18th century sunspot reconstruction of the SSN by Wolf based on the drawings of Staudacher may be even larger, reaching 25% (Svalgaard 2017). On the other hand, Izenman (1983) notes that Waldmeier’s authoritative 1965 edition of the $R_{Z}$ series does contain slight corrections also to the data published previously, in 1925 by Wolfer. (One might also consider this as a surreptitious earlier minor “revision” of the SSN.) This shows that Waldmeier himself was very much concerned with the long-term homogeneity of the data, already taking some measures to homogenize the data processing.

Amidst all the revision fervour, some caveats may still be in order. On the one hand, for the current generation of solar physicists it is abundantly clear that our Sun is capable of rather sudden unexpected changes and that a varying ratio between different activity indices (even just those related to sunspots) can be a real, physical feature (Georgieva et al 2017; see also Sect. 2.2 below). Suggestions for revisions based purely on variations or jumps in the ratio of the sunspot number to other indices are therefore to be treated with strong caution as it may potentially bias understanding of the physical processes involved. Second, having superior instruments does not entitle us to think we are smarter than our predecessors were, or that we know their data better than they themselves did. Information lost and unknown considerations may well explain practices that, in retrospect, seem incorrect to us. Even in the case of the already implemented corrections (b) and (c) some nagging questions do remain and need further exploration (Clette and Lefèvre 2016).

The significant disagreements between determinations of the SSN by various observatories, observers and methods are even more pronounced in the case of historical data, especially prior to the mid-19th century. In particular, the controversial suggestion that a whole solar cycle may have been missed in the official sunspot number series at the end of the 18th century is taken by some as glaring evidence for the unreliability of early observations. Note, however, that independently of whether the claim for a missing cycle is well founded or not, there is clear evidence that this controversy is mostly due to the very atypical behaviour of the Sun itself in the given period of time, rather than to the low quality and coverage of contemporary observations. These issues will be discussed further in Sect. 4.2.2.

1.1.3 Alternating series and nonlinear transforms

Instead of the “raw” sunspot number series $R(t)$ many researchers prefer to base their studies on some transformed index $R^{\prime}$. The motivation behind this is twofold.

(a) The strongly peaked and asymmetrical sunspot cycle profiles strongly deviate from a sinusoidal profile; also the statistical distribution of sunspot numbers is strongly at odds with a Gaussian distribution. This can constitute a problem as many common methods of data analysis rely on the assumption of an approximately normal distribution of errors or nearly sinusoidal profiles of spectral components. So transformations of $R$ (and, optionally, t) that reduce these deviations can obviously be helpful during the analysis. In this vein, e.g., Max Waldmeier often based his studies of the solar cycle on the use of logarithmic sunspot numbers $R^{\prime}=\log R$; many other researchers use $R^{\prime}=R^{\alpha}$ with $0.5\leq\alpha<1$, the most common value being $\alpha=0.5$.

(b) As the sunspot number is a rather arbitrary construct, there may be an underlying more physical parameter related to it in some nonlinear fashion, such as the toroidal magnetic field strength $B$, or the magnetic energy, proportional to $B^{2}$. It should be emphasized that, contrary to some claims, our current understanding of the solar dynamo does not make it possible to guess what the underlying parameter is, with any reasonable degree of certainty. In particular, the often used assumption that it is the magnetic energy, lacks any sound foundation. If anything, on the basis of our current best understanding of flux emergence we might expect that the amount of toroidal flux emerging from the tachocline should be $\int|B-B_{0}|\,dA$ where $B_{0}$ is some minimal threshold field strength for Parker instability and the surface integral goes across a latitudinal cross section of the tachocline (cf. Ruzmaikin 1997; indeed, a generalized linear $R$–$B$ link involving a threshold field strength has now also been used in the dynamo models of Pipin and Sokoloff 2011 and Pipin et al 2012).

As, however, the lifetime of any given sunspot group is finite and proportional to its size (Petrovay and van Driel-Gesztelyi 1997; Henwood et al 2010), instantaneous values of $R$ or the total sunspot area should also depend on details of the probability distribution function of $B$ in the tachocline. This just serves to illustrate the difficulty of identifying a single physical governing parameter behind $R$.

One transformation that may still be well motivated from the physical point of view is to attribute an alternating sign to even and odd Schwabe cycles: this results in the the alternating sunspot number series $R_{\pm}$. The idea is based on Hale’s well known polarity rules, implying that the period of the solar cycle is actually 22 years rather than 11 years, the polarity of magnetic fields changing sign from one 11-year Schwabe cycle to the next. In this representation, first suggested by Bracewell (1953), usually odd cycles are attributed a negative sign. This leads to slight jumps at the minima of the Schwabe cycle, as a consequence of the fact that for a 1 – 2 year period around the minimum, spots belonging to both cycles are present, so the value of $R$ never reaches zero; in certain applications, further twists are introduced into the transformation to avoid this phenomenon.

After first introducing the alternating series, in a later work Bracewell (1988) demonstrated that introducing an underlying “physical” variable $R_{B}$ such that

[TABLE]

(i.e., $\alpha=2/3$ in the power law mentioned in item (a) above) significantly simplifies the cycle profile. Indeed, upon introducing a “rectified” phase variable555The more precise condition defining $\phi$ is that $\phi=\pm\pi/2$ at each maximum and $\phi$ is quadratically related to the time since the last minimum. $\phi$ in each cycle to compensate for the asymmetry of the cycle profile, $R_{B}$ is a nearly sinusoidal function of $\phi$. The empirically found 3/2 law is interpreted as the relation between the time-integrated area of a typical sunspot group vs. its peak area (or peak $R_{Z}$ value), i.e., the steeper than linear growth of $R$ with the underlying physical parameter $R_{B}$ would be due to the larger sunspot groups being observed longer, and therefore giving a disproportionately larger contribution to the annual mean sunspot numbers. If this interpretation is correct, as suggested by Bracewell’s analysis, then $R_{B}$ should be considered proportional to the total toroidal magnetic flux emerging into the photosphere in a given interval. (But the possibility must be kept in mind that the same toroidal flux bundle may emerge repeatedly or at different heliographic longitudes, giving rise to several active regions.)

1.2 Other indicators of solar activity

1.2.1 Group sunspot number (GSN)

Reconstructions of $R$ prior to the early 19th century are increasingly uncertain. In order to tackle problems related to sporadic and often unreliable observations, Hoyt and Schatten (1998) introduced the Group Sunspot Number (GSN) as an alternative indicator of solar activity. In contrast to the SSN, the GSN only relies on counts of sunspot groups as a more robust indicator, disregarding the number of spots in each group. Furthermore, while $R_{Z}$ was determined for any given day from a single observer’s measurements (a hierarchy of auxiliary observers was defined for the case if data from the primary observer were unavailable), the GSN uses a weighted average of all observations available for a given day.

A further advantage is that, in addition to the published series, all the raw data upon which it is based are made public.

The GSN series published by Hoyt and Schatten (1998) remained unchanged until the 2010s when it was taken under revision in concert with the SSN revision discussed above. As in the case of the GSN there is no generally accepted responsible organization, i.e., no “official” series, revisions were undertaken by several teams, leading to conflicting results.666Note that Hoyt and Schatten (1998) included a coefficient of 12.08 in the definition of their GSN in order to bring it the same scale as the SSN; this coefficient has been omitted in the more recent reconstructions, sometimes leading to some confusion regarding the scale applied.

The common denominator of all efforts to reconstruct the GSN is (or should be) the common set of observations upon which the construction of the series is based. This observational data set has been greatly extended in the past two decades thanks to the discovery and/or publication of many previously inaccessible historical sources. An update of the database providing a good basis for subsequent efforts to construct a GSN series was compiled by Vaquero et al (2016). This archive is now available from the SILSO site and from the Historical Archive of Sunspot Observations (HASO)777http://haso.unex.es/. (In principle, a regular upgrade is planned, with version numbers, v1.0 referring to Hoyt & Schatten, but the project is currently stuck at version 1.12 dated May 2016.)

The original method of Hoyt and Schatten (1998) was subject to a random drift of the mean group number over long timescales. While consistent use of the Greenwich Photoheliograph Results, available from 1874, helped to avoid such a drift in the 20th century, a drift already appears from the late 19th century back, owing to the still evolving techniques of photography used at Greenwich. As a result, group numbers before this period may be systematically lower than what the SSN would suggest. This was the main issue motivating a revision of the GSN series.

New GSN series were compiled by Svalgaard and Schatten (2016) and by Usoskin et al (2016) using two alternative methods: the backbone method and a method based on active day fraction (ADF) statistics. The backbone method resulted in significantly elevated GSN values before about 1900, while the ADF method resulted in a series closer to the original Hoyt & Schatten values. Both of these methods have been subject to criticisms (Cliver 2016; Willamo et al 2017, 2018). Finally, Chatzistergos et al (2017) came up with a variety of the backbone method with an improved methodology for the fitting of successive backbones, resulting in an intermediate series. At the time of writing, this “ultimate backbone” GSN series, available at CDS888http://vizier.cfa.harvard.edu/viz-bin/VizieR?-source=J/A+A/602/A69

(or via the paper’s ADS link)., seems to be the most recommendable version for further analysis.

The GSN series has been reproduced for the whole period since 1611 and it is generally agreed that for the period 1611 – 1818 it is a more reliable reconstruction of solar activity than the relative sunspot number. Yet there have been relatively few attempts to date to use this data series for solar cycle prediction.

1.2.2 Other sunspot data

The classic source of sunspot area and position data is the Greenwich Photoheliographic Results (GPR) catalogue999https://www.ngdc.noaa.gov/stp/solar/solardataservices.html, covering the period 1874–1976. The official continuation of the GPR is the Debrecen Photoheliographic Data (DPD) catalogue101010http://fenyi.solarobs.csfk.mta.hu/en/databases/Summary/, commissioned by the IAU and containing data from 1973 (Baranyi et al 2016; Győri et al 2017). The Debrecen database also includes a revised/enriched version of the GPR in the same format as the DPD. Another GPR extension, with USAF/NOAA data, covering the period 1874–2019 is available from a website maintained by NASA MSFC staff.111111http://solarcyclescience.com/activeregions.html Sunspot data from many other observatories are also available at the NGDC site.

Recent years have seen a surge in the digitization and processing of sunspot drawings made before the photographic era. A major role in this work has been played by a team in Potsdam led by Rainer Arlt (Arlt 2008; Arlt et al 2013; Diercke et al 2015). As a result sunspot positions (butterfly diagrams) have now been reconstructed for the period 1826–1880 from drawings by Schwabe and Spörer (Leussu et al 2016; Leussu et al 2017); for the period 1749–1796 from drawings by Staudacher (Arlt 2009); for the period 1670-1711 from scattered information (Vaquero et al 2015b; Neuhäuser et al 2018); and for the period 1611–1631 from drawings by Scheiner and Galileo (Arlt et al 2016; Vokhmyanin and Zolotova 2018).

Instead of the sunspot number, the total area of all spots observed on the solar disk might seem to be a less arbitrary measure of solar activity. However, these data have been available since 1874 only, covering a much shorter period of time than the sunspot number data. In addition, the determination of sunspot areas, especially farther from disk center, is not as trivial as it may seem, resulting in significant random and systematic errors in the area determinations. Area measurements performed in two different observatories often show discrepancies reaching $\sim$ 30% for smaller spots (cf. the figure and discussion in Appendix A of Petrovay et al 1999). Despite these difficulties, attempts at reconstructing sunpot areas have also been made (Carrasco et al 2016; Senthamizh Pavai et al 2015a), and Muraközy et al (2016) recently proposed a new activity index based on a calibration of the emerged magnetic flux to sunspot areas.

1.2.3 Other activity indices

A number of other direct indicators of solar activity have become available from the 20th century

(see Ermolli et al 2014 for a recent review). These include, e.g., various plage indices or the 10.7 cm solar radio flux – the latter is considered a particularly good and simple to measure indicator of global activity (see Fig. 2). As, however, these data sets only cover a few solar cycles, their impact on solar cycle prediction has been minimal.

A promising exception from this is the nearly three centuries long record of the solar EUV flux, recently reconstructed from the diurnal variation of the geomagnetic field by Svalgaard (2016).

Of more importance are proxy indicators such as geomagnetic indices (the most widely used of which is the aa index), the occurrence frequency of aurorae or the abundances of cosmogenic radionuclides such as 14C and 10Be in natural archives. For solar cycle prediction uses such data sets need to have a sufficiently high temporal resolution to reflect individual 11-year cycles. For the geomagnetic indices such data have been available since 1868, while an annual 10Be series covering 600 years has been published by Berggren et al (2009). Attempts have been made to reconstruct the epochs and even amplitudes of solar maxima during the past two millennia from oriental naked eye sunspot records and from auroral observations (Stephenson and Wolfendale 1988; Nagovitsyn 1997), but these reconstructions are currently subject to too many uncertainties to serve as a basis for predictions. Isotopic data with lower temporal resolution are now available for up to 50 000 years in the past; while such data do not show individual Schwabe cycles, they are still useful for the study of long term variations in cycle amplitude. Inferring solar activity parameters from such proxy data is generally not straightforward. (See review by Usoskin 2017.)

1.3 The solar cycle and its variation

The series of $R$ values determined as described in Sect. 1.1 is plotted in Fig. 3. It is evident that the sunspot cycle is rather irregular. The mean length of a cycle (defined as lasting from minimum to minimum) is 11.0 years (median 10.9 years), with a standard deviation of 1.16 years; actual cycle lengths scatter in the range 9.0–13.6 years. Note that cycle lengths measured between successive maxima show a wider scatter, in the range 7.3 and 16.9 years. This is partly due to the fact that many cycles show a double maximum, the two sub-peaks being separated by 1–2 years.

The mean cycle amplitude in terms of $R$ is 179 (median 183)121212Here and in the rest of this paper $R$ will normally refer to Version 2 values, unless explicitly noted., with a standard deviation of 57. It is this wide variation that makes the prediction of the next cycle maximum such an interesting and vexing issue.

1.3.1 Secular activity variations

Inspecting Fig. 3 one can discern an obvious long term variation. For the study of such long term variations, the series of cycle parameters is often smoothed on time scales significantly longer than a solar cycle: this procedure is known as secular smoothing. One popular method is the so-called Gleissberg filter or 12221 filter (Gleissberg 1967). For instance, the Gleissberg filtered amplitude of cycle $n$ is given by

[TABLE]

The Gleissberg filtered series of solar cycles is plotted in Fig. 4. The most obvious feature of the variation is a cyclic modulation of the cycle amplitudes on a timescale of $\sim$9–10 solar cycles. This so-called Gleissberg cycle will be discussed further in Sect. 4.2.3. The first minimum of this cycle plotted in Fig. 4, known as the “Dalton Minimum”, is formed by the unusually weak cycles 5, 6, and 7. The second secular minimum consists of a rather long series of moderately weak cycles 12 – 16, occasionally referred to as the [last] “Gleissberg Minimum” – but note that most of these cycles are less than $1\,\sigma$ below the long-term average value given at the start of Section 1.3.

Finally the last secular maximum of the cycle comprises the series of strong cycles 17–23 in the second half of the 20th century: the “Modern Maximum”. In addition to this cyclic modulation there is a tendency for an overall secular increase of solar activity in the figure: the Modern Maximum is clearly stronger than previous maxima. However, the strength of this secular increase in the activity level as well as the amount by which the Modern Maximum exceeds previous maxima of the Gleissberg cycle clearly depends on the reconstruction of the measure of activity chosen. The revision of the sunspot numbers has greatly reduced the amount of secular increase compared to Version 1.0, in agreement with the GSN reconstruction by Svalgaard and Schatten (2016). On the other hand the most recent GSN reconstruction (Chatzistergos et al 2017) shows a marked long-term increasing trend again. The cosmogenic record rather unequivocally indicates that the persistently high level of solar activity characterizing the second half of the 20th century had no precedent for thousands of years in the history of solar activity (cf. Table 3 in Usoskin 2017). The currently hotly debated problem of the strength of the Modern Maximum has important implications, e.g., for the understanding of the role of solar forcing in global warming (Lean and Rind 2008; Chylek et al 2014; Nagy et al 2017; Owens et al 2017). In this context it is important to stress that a secular increase in solar activity from the late 19th century (beginning of terrestrial global temperature record) to the mid-20th century is unquestionably present in all solar activity reconstructions. The decreasing trend displayed in the last few decades in the smoothed cycle amplitude series is also potentially important in this respect.

While the Dalton and Gleissberg minima are just local minima in the ever changing Gleissberg filtered SSN series, the conspicuous lack of sunspots in the period 1640 – 1705, known as the Maunder Minimum (Fig. 1) quite obviously represents a qualitatively different state of solar activity. Such extended periods of low activity are usually referred to as grand minima. Ever since the rediscovery of the Maunder Minimum in the late 19th century (Maunder 1894; Eddy 1976) its reality and significance has time to time been brought into question. Recent studies have shown that the 11/22-year solar cycle persisted during the Maunder Minimum, but at a greatly suppressed level (Usoskin et al 2015; Vaquero et al 2015a; Asvestari et al 2017).

A number of possibilities have been proposed to explain the phenomenon of grand minima, including chaotic behaviour of the nonlinear solar dynamo (Weiss et al 1984), stochastic fluctuations in dynamo parameters (Moss et al 2008; Usoskin et al 2009b); a bimodal dynamo with stochastically induced alternation between two stationary states (Petrovay 2007) or stochastic effects like fluctuations in AR tilt (Karak and Miesch 2018) or “rogue” sunspots (Petrovay and Nagy 2018a).

The analysis of long-term proxy data, extending over several millennia further showed that there exist systematic long-term statistical trends and periods such as the so called secular and supersecular cycles (see Sect. 4.2).

1.3.2 Does the Sun have a long term memory?

The existence of long lasting grand minima and maxima suggests that the sunspot number record may have a long-term memory extending over several consecutive cycles. Indeed, elementary combinatorical calculations show that the occurrence of phenomena like the Dalton minimum (3 of the 4 lowest maxima occurring in a row) in a random series of 24 recorded solar maxima has a rather low probability (5%). This conclusion is supported by the analysis of long-term proxy data, extending over several millennia, which showed that the occurrence of grand minima and (perhaps) grand maxima is more common than what would follow from Gaussian statistics (Usoskin et al 2007; Wu et al 2018).

It could be objected that for sustained grand minima or maxima a memory extending only from one cycle to the next would suffice. This intercycle memory explanation of persistent secular activity minima and maxima, however, would imply a good correlation between the amplitudes of subsequent cycles, which is not the case (cf. Sect. 2.1 below). With the known poor cycle-to-cycle correlation, strong deviations from the long-term mean would be expected to be damped on time scales short compared to, e.g., the length of the Maunder minimum. This suggests that the persistent states of low or high activity may be due to truly long term memory effects extending over several cycles.

In an analysis of the GSN series for the period 1799–2011 Love and Rigler (2012) found that the sequence of cycle maxima (and also of time-integrated activity in each cycle), including the Modern Maximum, would not be an unlikely result of the accumulation of multiple random-walk steps in a lognormal random walk of cycle amplitudes where $\ln R$ performs a Gaussian random walk with mean stepsize 0.39 (or 0.28 for the integrated activity). This analysis, however, does not extend to the Maunder Minimum; and in any case, such a random walk should ultimately take the values of $R$ up to arbitrarily high values in sufficiently long time, whereas the cosmogenic record clearly shows that the level of activity is bounded from above.

Further evidence for a long-term memory in solar activity comes from the persistence analysis of activity indicators. The parameter determined in such studies is the Hurst exponent $0<H<1$. Essentially, $H$ is the steepness of the growth of the total range $\cal R$ of measured values plotted against the number $n$ of data in a time series, on a logarithmic plot: ${\cal R}\propto n^{H}$. For a Markovian random process with no memory $H=0.5$. Processes with $H>0.5$ are persistent (they tend to stay in a stronger-than-average or weaker-than-average state longer), while those with $H<0.5$ are anti-persistent (their fluctuations will change sign more often).

Hurst exponents for solar activity indices have been derived using rescaled range analysis by many authors (Mandelbrot and Wallis 1969; Ruzmaikin et al 1994; Komm 1995; Oliver and Ballester 1996; Kilcik et al 2009; Rypdal and Rypdal 2012). All studies coherently yield a value $H=0.85\mbox{\,--\,}0.88$ for time scales exceeding a year or so, and somewhat lower values ($H\sim 0.75$) on shorter time scales. Some doubts regarding the significance of this result for a finite series have been raised by Oliver and Ballester (1998); however, Qian and Rasheed (2004) have shown using Monte Carlo experiments that for time series of a length comparable to the sunspot record, $H$ values exceeding 0.7 are statistically significant.

A complementary method, essentially equivalent to rescaled range analysis is detrended fluctuation analysis. Its application to solar data (Ogurtsov 2004) has yielded results in accordance with the $H$ values quoted above.

The overwhelming evidence for the persistent character of solar activity and for the intermittent appearance of secular cyclicities, however, is not much help when it comes to cycle-to-cycle prediction. It is certainly reassuring to know that forecasting is not a completely idle enterprise (which would be the case for a purely Markovian process), and the long-term persistence and trends may make our predictions statistically somewhat different from just the long-term average. There are, however, large decadal scale fluctuations superposed on the long term trends, so the associated errors will still be so large as to make the forecast of little use for individual cycles.

1.3.3 Waldmeier effect and amplitude–frequency correlation

“Greater activity on the Sun goes with shorter periods, and less with longer periods. I believe this law to be one of the most important relations among the Solar actions yet discovered.”

(Wolf, 1861)

It is apparent from Fig. 3 that the profile of sunspot cycles is asymmetrical, the rise being steeper than the decay. Solar activity maxima occur 3 to 4 years after the minimum, while it takes another 7 – 8 years to reach the next minimum. It can also be noticed that the degree of this asymmetry correlates with the amplitude of the cycle: to be more specific, the length of the rise phase anticorrelates with the maximal value of $R$ (Fig. 5), while the length of the decay phase shows weak or no such correlation.131313Note that Osipova and Nagovitsyn (2017) recently constructed two separate group number series for small and large sunspot groups and they found that the Waldmeier effect applies better to large spots.

Historically, the relation was first formulated by Waldmeier (1935) as an inverse correlation between the rise time and the cycle amplitude; however, as shown by Tritakis (1982), the total rise time is a weak (inverse logarithmic) function of the rise rate, so this representation makes the correlation appear less robust. (Indeed, when formulated with the rise time it is not even present in some activity indicators, such as sunspot areas – cf. Dikpati et al 2008; Ogurtsov and Lindholm 2011.) As pointed out by Cameron and Schüssler (2008), the weak link between rise time and slope is due to the fact that subsequent cycle are known to overlap by 1–2 years, so in steeper rising cycles the minimum will occur earlier, thus partially compensating for the shortening due to a higher rise rate. The effect is indeed more clearly seen when the rate of the rise is used instead of the rise time (Lantos 2000; Cameron and Schüssler 2008) or if the rise time is redefined as the time spent from 20 to 80% of the maximal amplitude (Karak and Choudhuri 2011).

Nevertheless, when coupled with the nearly nonexistent correlation between the decay time and the cycle amplitude, even the weaker link between the rise time and the maximum amplitude is sufficient to forge a weak inverse correlation between the total cycle length and the cycle amplitude (Fig. 5). This inverse relationship was first noticed by Wolf (1861).

A stronger inverse correlation was found between the cycle amplitude and the length of the previous cycle by Hathaway et al (1994). This correlation is also readily explained as a consequence of the Waldmeier effect, as demonstrated in a simple model by Cameron and Schüssler (2007); the same probably holds for the correlations reported by Hazra et al (2015).

Note that in a more detailed study Solanki et al (2002) found that the correlation coefficient of this relationship has steadily decreased during the course of the historical sunspot number record, while the correlation between cycle amplitude and the length of the third preceding cycle has steadily increased. The physical significance (if any) of this latter result is unclear.

In what follows, the relationships found by Wolf (1861), Hathaway et al (1994), and Solanki et al (2002), discussed above, will be referred to as “$R_{\mathrm{max}}\mbox{\,--\,}t_{\mathrm{cycle},n}$ correlations” with $n=0,-1$ or $-3$, respectively.

Modern time series analysis methods offer several ways to define an instantaneous frequency $f$ in a quasiperiodic series. One simple approach was discussed in the context of Bracewell’s transform, Eq. (3), above. Mininni et al (2000) discussed several more sophisticated methods to do this, concluding that Gábor’s analytic signal approach yields the best performance. This technique was first applied to the sunspot record by Paluš and Novotná (1999), who found a significant long term correlation between the smoothed instantaneous frequency and amplitude of the signal. On time scales shorter than the cycle length, however, the frequency–amplitude correlation has not been convincingly proven, and the fact that the correlation coefficient is close to the one reported in the right hand panel of Fig. 5 indicates that all the fashionable gadgetry of nonlinear dynamics could achieve was to recover the effect already known to Wolf. It is clear from this that the “frequency–amplitude correlation” is but a secondary consequence of the Waldmeier effect.

Indeed, an anticorrelation between cycle length and amplitude is characteristic of a class of stochastically forced nonlinear oscillators and it may also be reproduced by introducing a stochastic forcing in dynamo models (Stix 1972; Ossendrijver et al 1996; Charbonneau and Dikpati 2000). In some such models the characteristic asymmetric profile of the cycle is also well reproduced (Mininni et al 2000; Mininni et al 2002). The predicted amplitude–frequency relation has the form

[TABLE]

where $f\sim(t_{\mathrm{min}}^{(n+1)}-t_{\mathrm{min}}^{(n)})^{-1}$ is the frequency.

Nonlinear dynamo models including some form of $\alpha$-quenching also have the potential to reproduce the effects described by Wolf and Waldmeier without recourse to stochastic driving. In a dynamo with a Kleeorin–Ruzmaikin type feedback on $\alpha$, Kitiashvili and Kosovichev (2009) are able to qualitatively reproduce the Waldmeier effect. Assuming that the sunspot number is related to the toroidal field strength according to the Bracewell transform, Eq. (3), they find a strong link between rise time and amplitude, while the correlations with fall time and cycle length are much weaker, just as the observations suggest. They also find that the form of the growth time–amplitude relationship differs in the regular (multiperiodic) and chaotic regimes. In the regular regime the plotted relationship suggests

[TABLE]

while in the chaotic case

[TABLE]

The linear relationship (6) was also reproduced in some stochastically forced nonlinear dynamo models (Pipin and Sokoloff 2011; Pipin and Kosovichev 2011; Pipin et al 2012).

Note that based on the actual sunspot number series Waldmeier originally proposed

[TABLE]

while according to Dmitrieva et al (2000) the relation takes the form

[TABLE]

At first glance, these logarithmic empirical relationships seem to be more compatible with the relation (5) predicted by the stochastic models. These, on the other hand, do not actually reproduce the Waldmeier effect, just a general asymmetric profile and an amplitude–frequency correlation. At the same time, inspection of the the left hand panel in Fig. 5 shows that the data is actually not incompatible with a linear or inverse rise time–amplitude relation, especially if the anomalous cycle 19 is ignored as an outlier. (Indeed, a logarithmic representation is found not to improve the correlation coefficient – its only advantage is that cycle 19 ceases to be an outlier.) All this indicates that nonlinear dynamo models may have the potential to provide a satisfactory quantitative explanation of the Waldmeier effect, but more extensive comparisons will need to be done, using various models and various representations of the relation. In one such exploratory study for instance Nagy and Petrovay (2013) found that solar-like parameter correlations can be obtained in a stochastically forced van der Pol oscillator but only if the perturbations are applied to the nonlinearity parameter rather than to the damping. In another study Karak and Choudhuri (2011) found that in a stochastically forced flux transport dynamo perturbing the poloidal field amplitude is not sufficient to induce solar-like parameter correlations, and perturbations to the meridional flow speed are also needed.

1.4 Approaches to solar cycle prediction

As the SSN series is a time series it is only natural that time series analysis methods have been widely applied in order to predict its future variations, including the amplitude of an upcoming cycle. As a group, however, time series methods have not been particularly successful in attaining this goal. In addition, time series analysis is a purely mathematical tool offering little physical insight into the processes driving cycle-to-cycle variations. In view of this, time series methods (or extrapolation methods) have been relegated to a later section of this review, after dealing with the currently much more lively field of the physically more insightful and more successful alternative approaches: precursor schemes and model-based forecasts.

In the 1st edition of this review the model-based approach, then still very new, was discussed well separated from precursor methods, in a section following the discussion of the time series approach. In the time elapsed since the 1st edition, however, a major surge of activity in surface flux transport (SFT) modelling, new developments in dynamo models and in empirical precursors have made it harder to draw a clear line between the precursor and model based approaches. Indeed, there seem to be at least “five shades of grey” arching between archetypical examples of these two categories:

•

(a) Internal empirical precursors (relying only on the SSN series)

•

(b) External empirical precursors (relying on other activity indicators)

•

(c) Physical[ly motivated] precursors

•

(d) Forecasts based on SFT models

•

(e) Forecasts based on dynamo models

Accordingly, the present edition has been reorganized so the section discussing model-based predictions immediately follows the section on precursors: the two topics have been separated, somewhat arbitrarily, between the classes (c) and (d) above, i.e., the term “model-based” is reserved for methods employing detailed quantitative models rather than empirical or semiempirical correlations based on qualitative physical ideas.

2 Precursor methods

“Jeder Fleckenzyklus muß als ein abgeschlossenes Ganzes, als ein Phänomen für sich, aufgefaßt werden, und es reiht sich einfach Zyklus an Zyklus.”

(Gleissberg, 1952)

In the most general sense, precursor methods rely on the value of some measure of solar activity or magnetism at a specified time to predict the amplitude of the following solar maximum. The precursor may be any proxy of solar activity or other indicator of solar and interplanetary magnetism. Specifically, the precursor may also be the value of the sunspot number at a given time.

In principle, precursors might also herald the activity level at other phases of the sunspot cycle, in particular the minimum. Yet the fact that practically all the good precursors found need to be evaluated at around the time of the minimum and refer to the next maximum is not simply due to the obvious greater interest in predicting maxima than predicting minima. Correlations between minimum parameters and previous values of solar indices have been looked for, but the results were overwhelmingly negative (e.g., Tlatov 2009). This indicates that the sunspot number series is not homogeneous and Rudolf Wolf’s instinctive choice to start new cycles with the minimum rather than the maximum in his numbering system is not arbitrary – for which even more obvious evidence is provided by the butterfly diagram. Each numbered solar cycle is a consistent unit in itself, while solar activity seems to consist of a series of much less tightly intercorrelated individual cycles, as suggested by Wolfgang Gleissberg in the motto of this section.

In Sect. 1.3.2 we have seen that there may be some evidence for a long-term memory underlying solar activity. In addition to the evidence reviewed there, systematic long-term statistical trends and periods of solar activity, such as the secular and supersecular cycles (to be discussed in Sect. 4.2), also attest to a secular mechanism underlying solar activity variations and ensuring some degree of long-term coherence in activity indicators. However, as we noted, this long-term memory is of limited importance for cycle prediction due to the large, apparently haphazard decadal variations superimposed on it. What the precursor methods promise is just to find a system in those haphazard decadal variations – which clearly implies a different type of memory. As we already mentioned in Sect. 1.3.2, there is obvious evidence for an intracycle memory operating within a single cycle, so that forecasting of activity in an ongoing cycle is currently a much more successful enterprise than cycle-to-cycle forecasting. As we will see, this intracycle memory is one candidate mechanism upon which precursor techniques may be founded, via the Waldmeier effect.

The controversial issue is whether, in addition to the intracycle memory, there is also an intercycle memory at work, i.e., whether behind the apparent stochasticity of the cycle-to-cycle variations there is some predictable pattern, whether some imprint of these variations is somehow inherited from one cycle to the next, or individual cycles are essentially independent. The latter is known as the “outburst hypothesis”: consecutive cycles would then represent a series of “outbursts” of activity with stochastically fluctuating amplitudes (Halm 1901; Waldmeier 1935; Vitinsky 1973; see also de Meyer 1981 who calls this “impulse model”). Note that cycle-to-cycle predictions in the strict temporal sense may be possible even in the outburst case, as solar cycles are known to overlap. Active regions belonging to the old and new cycles may coexist for up to three years or so around sunspot minima; and high latitude ephemeral active regions oriented according to the next cycle appear as early as 2 – 3 years after the maximum (Tlatov et al 2010 – the so-called extended solar cycle).

In any case, it is undeniable that for cycle-to-cycle predictions, which are our main concern here, the precursor approach seems to have been the relatively most successful, so its inherent basic assumption must contain an element of truth – whether its predictive skill is due to a “real” cycle-to-cycle memory (intercycle memory) or just to the overlap effect (intracycle memory).

The two precursor types that have received most attention are polar field precursors and geomagnetic precursors. A link between these two categories is forged by a third group, characterizing the interplanetary magnetic field strength or “open flux”. In terms of the classification outlined in Sect. 1.4 above, all these belong to the category (c) of physically motivated precursors. But before considering these approaches, we start by discussing categories (a) and (b): the empirical precursors based on the chance discovery of correlations between certain solar parameters and cycle amplitudes. These parameters involved may also be external to the SSN series (b); but first of all we will focus on the most obvious precursor type: internal empirical precursors (a) — the level of solar activity at some epoch before the next maximum.

2.1 Cycle parameters as precursors and the Waldmeier effect

The simplest weather forecast method is saying that “tomorrow the weather will be just like today” (works in about 2/3 of the cases). Similarly, a simple approach of sunspot cycle prediction is correlating the amplitudes of consecutive cycles. There is indeed a marginal correlation, but the correlation coefficient is quite low (0.35). The existence of the correlation is related to secular variations in solar activity, while its weakness is due to the significant cycle-to-cycle variations.

A significantly better correlation exists between the minimum activity level and the amplitude of the next maximum (Fig. 6). The relation is linear (Brown 1976), with a correlation coefficient of 0.68 (if the anomalous cycle 19 is ignored – Brajša et al 2009; see also Pishkalo 2008). The best fit is

[TABLE]

Cameron and Schüssler (2007) pointed out that the activity level three years before the minimum is an even better predictor of the next maximum. The exact value of the time shift producing optimal results may somewhat vary depending on the base period considered: in the Version 2 SSN series for the whole series of cycles 1–24 the highest correlation results for a time shift of 39 months while considering only cycles 8–24, for which data are sometimes considered more reliable, the best correlation is found at 32 months. Here in the right hand panel of Fig. 6 we simply use the round value of of 3 years (36 months) as originally proposed. The linear regression is

[TABLE]

This method, to be referred to as as “minimax3” for brevity, can only provide an upper estimate for the expected amplitude of Cycle 25 as the minimum has not been reached at the time of writing. As the minimum will take place no earlier than May 2019, the value of the predictor is not higher than 44.8, which would result in a cycle amplitude of 147. This already gives an indication that the upcoming cycle will be weaker than the climatological average.

As the epoch of the minimum of $R$ cannot be known with certainty until about a year after the minimum, the practical use of these methods is rather limited: a prediction will only become available 2 – 3 years before the maximum, and even then with the rather low reliability reflected in the correlation coefficients quoted above. In addition, as convincingly demonstrated by Cameron and Schüssler (2007) in a Monte Carlo simulation, these methods do not constitute real cycle-to-cycle prediction in the physical sense: instead, they are due to a combination of the overlap of solar cycles with the Waldmeier effect. As stronger cycles are characterized by a steeper rise phase, the minimum before such cycles will take place earlier, when the activity from the previous cycle has not yet reached very low levels.

The same overlap readily explains the $R_{\mathrm{max}}\mbox{\,--\,}t_{\mathrm{cycle},n}$ correlations discussed in Sect. 1.3.3. These relationships may also be used for solar cycle prediction purposes (e.g., Kane 2008) but they lack robustness. The forecast is not only sensitive to the value of $n$ used but also to the data set (relative or group sunspot numbers) (Vaquero and Trigo 2008). Similar correlations between the properties of subsequent cycles were used by Li et al (2015) to give a prediction for Cycle 25.

2.2 External empirical precursors

Cycle 24 peaked at an amplitude 35 % lower than Cycle 23. This was the 4th largest intercycle drop in solar activity in the monthly SSN record. The last such occurrence was rather different, following the single anomalously strong cycle 19, while in the case of the drop after Cycle 23 a series of 7 strong cycles was ended. This is illustrated in Fig. 4 where the shaded area marks a $\pm 2\sigma$ band around the mean amplitude of cycles 17–23 comprising the Modern Maximum. Cycle 24 clearly deviates downwards from the Modern Maximum cycles by over $2\sigma$ on the low side.

The first two occurrences of similar intercycle drops in activity (following cycles 4 and 11) are closer analogues of the recent events. These heralded the two previous Gleissberg minima centuries ago. The circumstance that a new Gleissberg minimum was already overdue suggests that this is indeed what we are witnessing. This is in line with current indications that solar activity will remain at Cycle 24 levels also in Cycle 25 (cf. Sect. 5). Yet no work predating 2005 predicted this drop which mostly caught solar physicists by surprise.

This unexpected drop in the level of solar activity from Cycle 23 to 24 has spurred efforts to find previously overlooked earlier signs of the coming change in solar data. A number of such “portents” were indeed identified, as first reviewed and correlated by Balogh et al (2014).

A “seismic portent” was identified by Basu et al (2012). High-frequency solar oscillations, sampling the top of the solar convective zone, have long been known to display frequency variations correlated with the solar cycle. The analysis of Basu et al (2012) showed that for [relatively] lower frequencies the amplitude of the frequency variation was strongly suppressed in Cycle 23, compared with Cycle 22 or with the variation in higher frequency modes (Fig. 7). This suggests that the (presumably magnetically modulated) variations in the sound speed were limited to the upper 3 Mm of the convective zone in Cycle 23, whereas in the previous cycle they extended to deeper layers. Revisiting the issue, Howe et al (2017) confirmed a change in the frequency response to activity during Solar Cycle 23, with a lower correlation of the low-frequency shifts with activity in the last two cycles compared to Cycle 22.

A similar disproportionately strong suppression of Cycle 23 relative to Cycle 22 is seen in the occurrence rate of flares, especially of class X and M (Fig. 8), and also in the variation of the Hα flare index. The suppression is rarely commented on yet clearly seen in the plots of Ataç and Özgüç (2001), Ataç and Özgüç (2006), Hudson et al (2014), or Gao and Zhong (2016).141414Hudson et al (2014) also called attention to the unusal temporal distribution of active region flare productivity during Cycle 23: in the first half of the cycle the number of flares produced by a flaring AR remained at an all-time low, then from the 2003 Halloween events it suddenly rose. It may be worth noting that the relationship between the SSN and the soft X-ray background was also found to differ for different solar cycles by Winter et al (2016).

This “eruptivity portent” is also manifested in the variation of the coronal index (green coronal line emissivity), as seen in the plots of Ataç and Özgüç (2006). A curious disagreement is seen regarding the suppression of the number of C class flares in Cycle 23: while the data of Hudson et al (2014) suggest a significant suppression of the number of these flares, only slightly less than for M and X class flares, Gao and Zhong (2016) found a much less strong suppression for C flares compared to M and X type flares. This is puzzling as both works are based on the same NOAA data, the only apparent difference being the exclusion of flares close to the limb by Hudson et al (2014).

The stronger suppression of larger flares might be interpereted as a relative lack of large active regions harbouring sufficient magnetic energy to produce such flares. Based on the expectation that the magnetic “roots” of larger ARs reach deeper, this would also agree with the seismic portent. Indeed, Howe et al (2017) explicitly speculate that the observed suppression of the low frequency modulation in Cycle 23 is “perhaps because a greater proportion of activity is composed of weaker or more ephemeral regions”.

Apparently in line with the above reasoning, de Toma et al (2013b) reported a strong suppression of the number of very large ($>700$ msh) sunspots and sunspot groups in Cycle 23. But the situation is not so simple as Kilcik et al (2011), Lefèvre and Clette (2011) and Kilcik et al (2014) found an apparently opposite trend: a strong suppression in the number of very small ($<17$ msh) spots or of sunspot groups of Zürich type A and B (pores/pore pairs) while the number of larger, more complex spots/groups is largely unaffected, or even slightly enhanced (Fig. 9). As the contribution to plage areas, radio flux, TSI or disk-integrated magnetic flux density is dominated by these large ARs, no significant suppression of Cycle 23 is detected in these proxies either (Göker et al 2017).

These perplexing findings may also be linked to the apparent decrease of the sunspot magnetic field strengths throughout Cycle 23 (Livingston et al 2012; Nagovitsyn et al 2016). There is, however, as yet no consensus regarding the reality of this trend (de Toma et al 2013a; Watson et al 2014) or regarding to what extent they are cycle related or due to secular trends (Norton et al 2013; Rezaei et al 2012, 2015; Nagovitsyn et al 2017).

Studies pointing to possible interrelationships between the various portents discussed above include Kilcik et al (2018) where a stronger decrease in sunspot count in flaring AR was reported compared to non-flaring regions. While local subsurface flow properties in AR, in particular vorticity, have also been found to correlate with flare productivity (Mason et al 2006; Komm et al 2011, 2015), the apparently only study of the relationship between local disturbances of seismic properties (such as sound speed) in AR and flare index led to inconclusive results (Lin 2014).

2.3 Polar precursor

The polar precursor method, as first suggested by Schatten et al (1978), is based on the correlation between the amplitude of a sunspot maximum with a measure of the amplitude of the magnetic field near the Sun’s poles at the preceding cycle minimum. Its physical background is the plausible causal relationship between the toroidal flux and the poloidal flux that serves as a seed for the generation of toroidal fields by the winding up of field lines in a differentially rotating convective zone.

It is now widely agreed that, beside internal empirical precursor methods based on the Waldmeier rule, the polar precursor method is currently the most reliable way to forecast an upcoming solar cycle. As the first revision of this review concluded, the polar precursor method “has consistently proven its skill in all cycles.” It is now also widely agreed that the polar precursor stands behind the apparent predicting skill of several other forecasting methods, including geomagnetic precursors.

2.3.1 Polar magnetic field data

Observational data on magnetic fields near the Sun’s poles were reviewed by Petrie (2015). Solar magnetograms have been available on a regular basis from Mt. Wilson Observatory since 1974, from Wilcox Solar Observatory (WSO) since 1976, and from Kitt Peak since 1976 (with a major change in the instrument from KPVT to SOLIS in 2003). The most widely used set of direct measurements of the magnetic field in the polar areas of the Sun is from the WSO series (Svalgaard et al, 1978; Hoeksema, 1995). While these magnetograms have the lowest resolution of the three sets, from the point of view of the characterization of the polar fields this is not necessarily a disadvantage, as integrating over a larger aperture suppresses random fluctuations and improves the S/N ratio. As a result of the low resolution, the WSO polar field value is a weighted average of the line-of-sight field in a polar cap extending down to $\sim 55^{\circ}$ latitude on average (with significant annual variations due to the $7^{\circ}$ tilt of the solar axis).

The classic reference on the processing and analysis of WSO polar field data is Svalgaard et al (1978). The inference from their analysis was that assuming a form $B=B_{0}f(\theta)$ with $f(\theta)=\cos^{n}(\theta)$ for the actual mean magnetic field profile ($\theta$ is colatitude) inside the polar cap around minimum, $n=8\pm 1$ while $B_{0}$ was around 10 G for Cycle 21 and the next two cycles, being reduced to about half that value in Cycle 24. While one later study (Petrie and Patrikeeva 2009) points to a possibility that the value of $n$ may be even higher, up to 10, the “canonical” value $n=8$ seems quite satisfactory in most cases (e.g., Fig. 2 in Whitbread et al 2017).

Figure 10 shows the variation of the smoothed amplitude of the WSO polar field, averaged over the two poles. (The presence of undamped residual fluctuations on short time scales illustrates the unsatisfactory nature of the 13-month smoothing, applied here for consistency with the rest of this review. A regularly updated plot of the WSO polar field with a more optimal smoothing (low-pass filter) is available from the WSO web site151515http://wso.stanford.edu/gifs/Polar.gif.)

Also shown is an alternative measure of the amplitude of the poloidal field component, the axial dipole coefficient, i.e., the amplitude of the coefficient of the $Y_{1}^{0}$ term in a spherical harmonic expansion of the distribution of the radial magnetic field strength over the solar disk:

[TABLE]

where $\overline{B}$ denotes the azimuthally averaged radial magnetic field.

This formula assumes the use of the Schmidt quasi-normalization in the definition of the spherical harmonics, widely used in solar physics and geomagnetism (see, e.g., Winch et al 2005). For direct comparison of the amplitudes of harmonics of different degree, a full normalization is sometimes preferred (e.g., in DeRosa et al 2012): this results in a normalized dipole coefficient $\hat{D}=(4\pi/3)^{1/2}D$. While (12) or even (13) are often loosely referred to as the “solar dipole moment”, it should also be kept in mind that the magnetic [dipole] moment, as normally defined in physics, is related to $D$ as $(2\pi R_{\odot}^{3}/\mu_{0})D$ where $R_{\odot}$ is the solar radius and $\mu_{0}$ is the vacuum permeability.

The two curves in Fig. 10 are quite similar even in many of their details: the polar field amplitude follows the variations of the dipole coefficient with a phase lag of about a year. This is hardly surprising as the hemispherically averaged polar field amplitude $|B_{N}-B_{S}|/2$ is clearly proportional to the contribution to $D$ coming from the polar caps:

[TABLE]

As the polar field is formed by the poleward transport of magnetic fields at lower latitudes, it is only to be expected that the variation of the polar cap contribution will follow that of the overall dipole with some time delay.

Indeed, based on the good agreement of the two curves in Fig. 10, $(B_{N}-B_{S})/2$ may simply be used as a simple measure of the amplitude of the dipole (on an arbitrary scale); on similar grounds, $(B_{N}+B_{S})/2$ may be considered a measure of the quadrupole component (e.g., Svalgaard et al 2005; Muñoz-Jaramillo et al 2013a).

We note that the dipole moment in our figure is an average of all available values from different magnetogram data sets for a given date; however, there is a quite good overall agreement among the values from different data sets (see Fig. 9 in Jiang et al 2018).161616With the exception of NSO data (not shown in that plot), which in Cycle 21 significantly deviate from the others.

The behaviour of the curves in Fig. 10 further shows that the times of dipole reversal are usually rather sharply defined. Based on the 4 reversals seen in the plot, the overall dipole is found to reverse $3.44\pm 0.18$ years after the minimum, while the polar contribution to the dipole reverses after $4.33\pm 0.36$ years. (The formal errors given are $1\sigma$.) The low scatter in these values suggests that the cycle phase of dipole reversal may be a sensitive test of SFT and dynamo models.

In contrast to reversal times, maxima of the dipole amplitude are much less well defined (occurring $7.27\pm 1.38$ and $8.33\pm 1.08$ years after minimum for the two curves). The curves display broad, slightly slanting plateaus covering 3 to 5 years (Iijima et al 2017); the dipole amplitude at the time of solar minimum is still typically $84\pm 12\%$ (global dipole) and $90\pm 6\%$ (polar fields) of its maximal value, reached years earlier. This kind of slanting profile is actually good news for cycle prediction as it opens the way to guess the dipole amplitude at the time of minimum, used as a predictor, years ahead. For example, using the rather flat and low preceding maximum in polar field strength, Svalgaard et al (2005) were able to predict a relatively weak Cycle 24 (predicted Version 1 peak SSN value $75\pm 8$ vs. $67$ observed) as early as 4 years before the sunspot minimum took place in December 2008!

The potential use of the dipole amplitude as a precursor is borne out by the comparison with the sunspot number curve in Fig. 10. After our arbitrary rescaling of the SSN, its maxima in each cycle are roughly at level with the preceding plateaus of the solid curves. This certainly seems to indicate that the suggested physical link between the precursor and the cycle amplitude is real.

The flatness of the maxima of the polar field imply that the precursor normally cannot be strongly affected by the exact time when it is evaluated. The cycle overlap effect combined with the Waldmeier relation, affecting the timing of the minimum, is therefore unlikely to explain the predictive skill of the polar precursor: we are here dealing with a real physical precursor (as also argued by Charbonneau and Barlet 2011).

The “polar precursor” may thus be interpreted in four different ways. The precursor may be the value of the global dipole moment or of the contribution of polar fields to this dipole moment only (i.e., the WSO field); and either of these may be evaluated at cycle minimum, or a few years earlier when they reach their respective maxima. Table 1 lists these precursor values for individual solar cycles, compared with the actual cycle amplitude. A homogeneous linear fit with one free parameter to the precursor–cycle amplitude relation yields the fit coefficient values given in the lower part of the table. The nominal random scatter is also indicated. Precursor values for Cycle 25 have been evaluated in late 2018; as it is not yet clear whether the maximum of the dipole moment has passed or when the minimum will take place, the forecasts based on these values are to be interpreted as lower/upper estimates, respectively. Taking into account the given values of the scatter, a combination of these results implies that Cycle 25 should peak in the range 102–133 with a nominal probability of $\sim 70\,$%, i.e. the next cycle should have an amplitude similar to Cycle 24.

2.3.2 Proxy reconstructions of the polar magnetic field

Despite the plausibility of a physical link between polar magnetic fields and cycle amplitude, the shortness of the available direct measurement series represents a difficulty when it comes to finding a convincing statistical link between these quantities. A way to circumvent this difficulty is offered by the availability of proxy data for polar magnetism spanning much longer time scales. Indeed, Schatten’s original suggestion of a polar field precursor (Schatten et al 1978) on a generic physical basis was supported by such proxy studies. It is remarkable that despite the very limited available experience, proxy-based forecasts using the polar field method proved to be consistently in the right range for Cycles 21, 22, and 23 (Schatten and Sofia, 1987; Schatten et al, 1996).

The types of proxies of solar polar magnetism were reviewed by Petrie et al (2014). In the present subsection we focus on polar faculae, which are currently considered the best photospheric proxy data for the reconstruction of polar magnetic field/flux. The interplanetary and geomagnetic precursors discussed in the following sections, however, may also be interpreted as indirect proxies of the solar polar magnetic field.

High resolution observations by the Hinode space observatory confirm that, like all other magnetic fields in the solar photosphere, polar fields are highly intermittent, nearly all of the flux being concentrated in isolated strong magnetic field concentrations (Tsuneta et al, 2008). These magnetic elements are observed as bright facular points in white light and in some spectral lines. The number density of these polar faculae is then related to the intensity of the polar magnetic field, while their total number above a certain latitude is related to the total magnetic flux (Sheeley 1964). This conclusion was indeed confirmed by Li et al (2002) and, more recently, by Tlatov (2009).

Perhaps the most carefully analyzed polar facular data set is the series of observations in the Mt. Wilson Observatory. These data were validated against MDI observations and then calibrated to WSO and MDI magnetic measurements by Muñoz-Jaramillo et al (2012). This resulted in a time series of the solar polar magnetic flux (poleward of $70^{\circ}$ latitude) for each hemisphere, covering the period 1906–2014. Owing to the varying tilt of the solar axis, the data are available with an annual cadence and with a time shift of 0.5 year between the hemispheres. This reconstructed polar flux was subsequently correlated with cycle amplitudes, confirming the usefulness of the polar precursor method (Muñoz-Jaramillo et al 2013a; 2013b). Hemispherically averaged data show a highly significant correlation, albeit with a large scatter ($r=0.69$ at 96%). This imperfect correlation is reflected in Fig. 11 in the wildly varying ratio between the maxima of the black curve and the subsequent maxima of the green dotted curve.

The poor correlation can be improved by separating the hemispheres, but some outlier points appear, apparently obeying an alternative relationship. Arguing that outliers correspond to cycles where the hemispheric asymmetry of polar fields exceeds a threshold, Muñoz-Jaramillo et al (2013a) finally arrived at four (two for each hemisphere) linear empirical relationships between the reconstructed polar magnetic flux at minimum and the amplitude of the next cycle. The suggested relation correctly reproduced the observed amplitude of Cycle 24 (predicted Version 1 SSN $77\pm 16$ vs. $67$ observed).

Another relevant data set is the polar facular counts recorded in the National Astronomical Observatory of Japan (NAOJ) at Mitaka Observatory during the period 1954–1996 (Li et al 2002). A third series of polar facular data, originating from Kodaikanal Observatory Ca K spectroheliograms, was compiled by Priyal et al (2014). (As these are chromospheric features, the authors prefer to call their index a polar network index.) Major disagreements between these data sets are seen in Fig. 11 indicating that the use of polar facular proxies is still not on very firm grounds. This is further shown by a comparison of the deducted MWO polar flux with polar magnetic fluxes computed from WSO polar field data using the calibration formula derived in Muñoz-Jaramillo et al (2012) (magenta curve in the plot). It is apparent that prior to the calibration period 1996–2006 the reconstructed MWO polar fluxes are systematically higher than the measured WSO fluxes, suggesting problems with these calibrations. In particular, while WSO polar field strengths peak at roughly the same amplitude in cycles 21, 22 and 23, in agreement with the comparable amplitudes of the subsequent sunspot cycles, the different polar facular counts indicate greatly different amplitudes for these minima, and they are also mutually incompatible.

Polar faculae are not the only long-term data base relevant for the study of high-latitude magnetic fields. $\mathrm{H}_{\alpha}$ synoptic charts are available from various observatories from as early as 1870. As $\mathrm{H}_{\alpha}$ filaments and filament channels lie on the magnetic neutral lines, these maps can be used to reconstruct the overall topology, if not the detailed map, of the large-scale solar magnetic field. While in between the neutral lines only the polarity of the field can be considered known, the variation of the global dipole moment may be tolerably well estimated even by fixing the field amplitude at a constant value. As higher order multipoles decay rapidly with distance, a potential field model fitted to the given synoptic map will then yield an acceptable representation of the polar field at the source surface. With this approach Makarov et al (2001) computed the multipole coefficients of the solar magnetic field, reconstructing the polar field strength at the source surface back to 1915. They also introduced the so-called dipole-octupole index (aka $A$-index): the sum of the axial dipole and octupole magnetic moments as a simple measure of the polar field amplitude. The method was later applied by Obridko and Shelting (2008) to predict the amplitude of Cycle 24: their forecast proved to be within 10% of the actual value. In turn, Tlatov (2009) has shown that several indices of the polar magnetic field during the activity minimum, determined from these charts, correlate well with the amplitude of the incipient cycle.

2.3.3 Extending the range of the polar precursor: early

forecasts for Cycle 25

The polar precursor is customarily evaluated at the cycle minimum, offering a prediction over a time span of 3 – 4 years, comparable to the rise time of the next cycle. Table 1, however, shows that using the maximal value of the global dipole moment results in a somewhat lower scatter around the mean relationship. This opens the possibility of making a prediction several years before the minimum.

The forecasting potential of the global dipole may also provide the ground for the findings of Hawkes and Berger (2018) who proposed the “helicity flux” as a cycle precursor. Perhaps more aptly called helicity input rate by the differential rotation, their helicity flux is defined by a weighted hemispheric surface integral of (a functional of) the radial magnetic field, where the weight function is fixed by the differential rotation profile. In such an integral, which is not unlike equation (12), the dipole component will naturally give a dominant contribution as the higher order terms change sign over each hemisphere, which largely cancels their effect. The authors found that the helicity input rate anticipates the sunspot numbers with a time shift of 4.5–6.9 years (depending on the cycle considered). This seems to be in line with the variable time delays between maxima of the global dipole moment and the next cycle maximum (Fig. 10). Their prediction for the amplitude of Cycle 25 is 117 — similar or slightly stronger than Cycle 24.

Similarly, using the brightness temperature of the 17 GHz microwave emission as a proxy for the field strength, Gopalswamy et al (2018) found that the correlation between this proxy and the sunspot number is maximal for time shifts of $\sim 4$–$6$ years (depending on cycle and hemisphere). Their forecast for Cycle 25 is a smoothed SSN of 89 for the S hemisphere and 59 for the N hemisphere (the latter being a lower bound as the proxy had not reached its maximum at the time of publication). These values are again comparable to or slightly higher than the amplitudes for Cycle 24.

In this context it is interesting to note that Makarov et al (1989) and Makarov and Makarova (1996) found that the number of polar faculae observed at Kislovodsk anticipates the next sunspot cycle with a time lag of 5 – 6 years on average in cycles 20–22; even short term annual variations or “surges” of sunspot activity were claimed to be discernible in the polar facular record. An apparently conflicting result was obtained by Li et al (2002), who found using the Mitaka data base that the best autocorrelation results with a time shift of about 4 years only. The discrepancy may perhaps be related to the cycle dependence of these time shifts, as partly different periods were considered.

In theory it is conceivable that successful forecasts might be attempted even earlier. After all, the polar field starts to build up at the time of polar reversal, about 5–6 years before the next minimum. Petrovay et al (2018) explored this possibility in a study of the coronal green lime emission at high latitudes. They tried to correlate various features of the characteristic “rush to the poles” (RTTP) feature of this emission around the polar reversal with properties of the following sunspot cycle and found a significant correlation between the speed of the RTTP and the time from the reversal to the next maximum. From this they predict that Cycle 25 will most likely peak in late 2024. Combining this date with the minimax3 method discussed in Section 2.1 above, the cycle is amplitude is estimated as 130, and the minimum is expected for 2019.

It is worth noting here that, in addition to the start of the increase of the polar field 5–6 years before the minimum, early precursors at high latitudes may be expected also on completely different grounds. The concept of the extended solar cycle implies that small ephemeral bipoles belonging to an upcoming solar cycle appear at high latitudes and start to migrate equatorward years before the first spots of the new cycle are observed. Thus, early signs of the equatorward propagating toroidal flux ring at high latitudes may give hints on the amplitude of an upcoming cycle (cf. also Badalyan et al, 2001). (From a formal point of view this would be then an internal cycle precursor, related to the one based on the Waldmeier rule.)

It is not impossible that some of the very early precursors suggested above, if real, may be partly explained by this effect. For example, Makarov and Makarova (1996) considered all faculae poleward of $50^{\circ}$ latitude. Bona fide polar faculae, seen on Hinode images to be knots of the unipolar field around the poles, are limited to higher latitudes, so the wider sample may consist of a mix of such “real” polar faculae and small bipolar ephemeral active regions. These latter are known to obey an extended butterfly diagram, as confirmed by Tlatov et al (2010): the first bipoles of the new cycle appear at higher latitudes about 4 years after the activity maximum.

In interpreting high-latitude migration patterns it should, however, be taken into account that it is not yet clear how far the wings of the butterfly diagram can actually extended backwards, i.e., to what extent a high latitude equatorward propagating branch is contiguous with the low latitude branch or is an unrelated phomenon (cf. the discussion in Cliver 2014 and Petrie et al 2014).

The high-latitude torsional oscillation pattern is usually considered the most pregnant manifestation of the extended solar cycle. This pattern has been unusually week during Solar Cycle 24, apparently suggesting significant further weakening of solar activity (Howe et al 2013). However, later observations (Howe et al 2018; Komm et al 2018) indicate that the low-latitude equatorward branch of the torsional oscillation is actually stronger in Cycle 24 than it was in Cycle 23, if measured against the same background flow.

In order to obtain a precursor that varies smoothly enough to be useful also between successive minima, Schatten and Pesnell (1993) introduced a new activity index, the “Solar Dynamo Amplitude” (SoDA) index, combining the polar field strength with a traditional activity indicator (the 10.7 cm radio flux F10.7). Around minimum, SoDA is basically proportional to the polar precursor and its value yields the prediction for F10.7 at the next maximum; however, it was constructed so that its 11-year modulation is minimized, so theoretically it should be rather stable, making predictions possible well before the minimum.

It remains to be seen whether SoDA actually improves the predictive skill of the polar precursor, to which it is more or less equivalent in those late phases of the solar cycle when forecasts start to become reliable. Using the SoDA index Pesnell and Schatten (2018) predicted Cycle 25 to peak at $R=134\pm 25$ in 2025.

2.4 Geomagnetic and interplanetary precursors

Relations between the cycle related variations of geomagnetic indices and solar activity were noted long ago. It is, however, important to realize that the overall correlation between geomagnetic indices and solar activity, even after 13-month smoothing, is generally far from perfect. This is due to the fact that the Sun can generate geomagnetic disturbances in two ways:

•

(a) By material ejections (such as CMEs or flare particles) hitting the terrestrial magnetosphere. This effect is obviously well correlated with solar activity, with no time delay, so this contribution to geomagnetic disturbances peaks near, or a few years after, sunspot maximum. (Note that the occurrence of the largest flares and CMEs is known to peak some years after the sunspot maximum – see Figure 16 in Hathaway 2015.)

•

(b) By a variation of the strength of the general interplanetary magnetic field and of solar wind speed. Geomagnetic disturbances may be triggered by the alternation of the Earth’s crossing of interplanetary sector boundaries (slow solar wind regime) and its crossing of high speed solar wind streams while well within a sector. The amplitude of such disturbances will clearly be higher for stronger magnetic fields. The overall strength of the interplanetary magnetic field, in turn, depends mainly on the total flux present in coronal holes, as calculated from potential field source surface models of the coronal magnetic field. At times of low solar activity the dominant contribution to this flux comes from the two extended polar coronal holes, hence, in a simplistic formulation this interplanetary contribution may be considered linked to the polar magnetic fields of the Sun, which in turn is a plausible precursor candidate as we have seen in the previous subsection. As the polar field reverses shortly after sunspot maximum, this second contribution often introduces a characteristic secondary minimum in the cycle variation of geomagnetic indices (cf. Fig. 12, somewhere around the maximum of the curve.

The component (a) of the geomagnetic variations actually follows sunspot activity with a variable time delay. Thus a geomagnetic precursor based on features of the cycle dominated by this component has relatively little practical utility. This would seem to be the case, e.g., with the forecast method first proposed by Ohl (1966), who noticed that the minimum amplitudes of the smoothed geomagnetic aa index are correlated to the amplitude of the next sunspot cycle (see also Du et al 2009, Du 2012). This is clearly just a secondary consequence of the previously discussed correlations between solar indices at cycle minimum (occurring somewhat earlier than the minimum in $aa$, see Fig. 12) and the amplitude of the next maximum. While the polar precursor is more directly relevant and is available sooner, the amplitude of the $aa$ index at minimum is still useful as a proxy as these data have been available since 1868.

An indication that the total geomagnetic activity, resulting from both mechanisms does contain useful information on the expected amplitude of the next solar cycle was given by Thompson (1993), who found that the total number of disturbed days in the geomagnetic field in cycle $n$ is related to the sum of the amplitudes of cycles $n$ and $n+1$ (see also Dabas et al, 2008).

A method for separating component (b) was proposed by Feynman (1982) who correlated the annual aa index with the annual mean sunspot number and found a linear relationship between $R$ and the minimal value of aa for years with such $R$ values. She interpreted this linear relationship as representing the component (a) discussed above, while the amount by which aa in a given year actually exceeds the value predicted by the linear relation would be the contribution of type (b) (the “interplanetary component”). The interplanetary component usually peaks well ahead of the sunspot minimum and the amplitude of the peak seemed to be a good predictor of the next sunspot maximum. However, it is to be noted that the assumption that the “surplus” contribution to aa originates from the interplanetary component only is likely to be erroneous, especially for stronger cycles. It is known that the number of large solar eruptions shows no unique relation to $R$: in particular, for $R$ > 100 their frequency may vary by a factor of 3 (see Figure 18 in Hathaway, 2015), so in some years they may well yield a contribution to aa that greatly exceeds the minimum contribution. A case in point was the “Halloween events” of 2003, that very likely resulted in a large false contribution to the derived “interplanetary” aa index (Hathaway, 2010). As a result, the geomagnetic precursor method based on the separation of the interplanetary component predicted an unusually strong cycle 24 ($R_{m}\sim 150$), in contrast to most other methods, including Ohl’s method and the polar field precursor, which suggested a weaker than average cycle ($R_{m}\sim 80\mbox{\,--\,}90$).

An alternative method for the separation of the interplanetary component, based on the use of the use of the $F_{10.7}$ index to model the variation of the activity-related component, was proposed by Pesnell (2014).

In addition to the problem of neatly separating the interplanetary contribution to geomagnetic disturbances, it is also wrong to assume that this interplanetary contribution is dominated by the effect of polar magnetic fields at all times during the cycle. Indeed, Wang and Sheeley (2009) pointed out that the interplanetary magnetic field amplitude at the Ecliptic is related to the equatorial dipole moment of the Sun that does not survive into the next cycle, so despite its more limited practical use, Ohl’s original method, based on the minima of the aa index is physically better founded, as the polar dipole dominates around the minimum. (Cf. also Ng 2016.) The total amount of open interplanetary flux, more closely linked to polar fields, could still be determined from geomagnetic activity if the interplanetary contribution to it is further split into:

•

(b1) A contribution due to the varying solar wind speed (or to the interplanetary magnetic field strength anticorrelated with it), which in turn reflects the strength of the equatorial dipole.

•

(b2) Another contribution due to the overall interplanetary field strength or open magnetic flux, which ultimately reflects the axial dipole.

Clearly, if the solar wind speed contribution (b1) could also be subtracted, a physically better founded prediction method should result. While in situ spacecraft measurements for the solar wind speed and the interplanetary magnetic field strength do not have the necessary time coverage, Svalgaard and Cliver (2005, 2007) and Rouillard et al (2007) devised a method to reconstruct the variations of both variables from geomagnetic measurements alone. Building on their results, Wang and Sheeley (2009) arrived at a prediction for the maximum amplitude of solar cycle 24 which agreed well with that of Bhatt et al (2009), who applied a modified form of Ohl’s method. Both forecasts proved to be somewhat too high (by $\sim 20$–$25\,$%, a little over $1\sigma$).

The open magnetic flux can also be derived from the extrapolation of solar magnetograms using a potential field source surface model. The magnetograms applied for this purpose may be actual observations or the output from surface flux transport models, using the sunspot distribution (butterfly diagram) and the meridional flow as input. Such models indicate that the observed latitude independence of the interplanetary field strength (“split monopole” structure) is only reproduced if the source surface is far enough ($>10\,R_{\odot}$) and the potential field model is modified to take into account the heliospheric current sheet (current sheet source surface model, Schüssler and Baumann 2006; Jiang et al 2010a). The extrapolations are generally found to agree well with in situ measurements where these are available. A comprehensive review of this topic is given in Lockwood (2013).

2.5 The quest for a precursor of the polar precursor

Just as the suggestion of a polar precursor was based on a qualitative physical understanding of the process generating the strong toroidal magnetic field that gives rise the observed active regions, an extension of the temporal range of our forecasting capability would clearly benefit from a similar qualitative physical understanding of how the strong polar fields prevalent around sunspot minimum are formed.

Magnetograph observations rather clearly indicate that the polar field is built up as a result of the the poleward transport of trailing polarity flux from active regions, while much of the leading polarity flux cancels with its counterpart on the other hemisphere by cross-equatorial diffusion. In the currently widely popular Babcock–Leighton scenario the poleward transport is interpreted as a combination of turbulent diffusion and advection by a poleward meridional flow.171717Note that the real situation may well be more complex than this simple scenario suggests: in numerical simulations of spherical turbulent dynamos latitudinal transport by pumping effects is quite often prevalent (see, e.g., Racine et al 2011, Simard et al 2016, Warnecke et al 2018). This suggests that the buildup of the polar field may be controlled by either of two effects:

(a) variations of the poleward flow speed

(b) variations in the tilt angles and/or latitudinal distribution of bipolar active regions which ultimately determine the net flux imbalance in the meridional direction.

In the following subsections these two influencing factors are considered in turn.

2.5.1 Photospheric flow variations

Considering the effect of meridional flow variations on intercycle variations is a delicate task as such

changes are also associated with the normal course of the solar activity cycle, the overall flow at mid-latitudes being slower before and during maxima and faster during the decay phase. Therefore, it is just the cycle-to-cycle variation in this normal pattern that may be associated with the activity variations between cycles. In this respect it is of interest to note that the poleward flow in the late phases of Cycle 23 seems to have had an excess speed relative to the previous cycle (Hathaway and Rightmire, 2010). If this were a latitude-independent amplitude modulation of the flow, then most flux transport dynamo models (e.g., Belucz et al 2015) would predict a stronger than average polar field at the minimum, contrary to observations. On the other hand, in the surface flux transport model of Wang et al (2009) an increased poleward flow actually results in weaker polar fields, as it lets less leading polarity flux to diffuse across the equator and cancel there. As the analysis by Muñoz-Jaramillo et al (2010) has shown, the discrepancy resulted from the form of the Babcock–Leighton source term in flux transport dynamo models, and it can be remedied by substituting a pair of opposite polarity flux rings representing each individual AR as source term instead of the $\alpha$-term. With this correction, 2D flux transport dynamos and surface flux transport models agree in predicting a weaker polar field for faster meridional flow.

It is known from helioseismology and magnetic correlation tracking that meridional flow speed fluctuations follow a characteristic latitudinal pattern associated with torsional oscillations and the butterfly diagram, consisting of a pair of axisymmetric bands of latitudinal flows converging towards the activity belts, migrating towards the equator, and accompanied by similar high-latitude poleward branches (Snodgrass and Dailey 1996, Chou and Dai 2001, Beck et al 2002, Liang et al 2018, Lin and Chou 2018).

This suggests interpreting the unusual meridional flow speeds observed during cycle 23 as an increased amplitude of this migrating modulation, rather than a change in the large-scale flow speed (Cameron and Schüssler, 2010). In this case, the flows converging on the activity belts tend to inhibit the transport of following polarities to the poles, resulting in a lower polar field (Jiang et al 2010b; note, however, that Švanda et al 2007 find no change in the flux transport in areas with increased flows). It is interesting to note that the torsional oscillation pattern, and thus presumably the associated meridional flow modulation pattern, was shown to be fairly well reproduced by a microquenching mechanism due to magnetic flux emerging in the active belts (Petrovay and Forgács-Dajka, 2002). Alternatively, the modulation pattern may also be thermally induced (Spruit 2003) or it may result from large-scale magnetic field torques (Passos et al 2017, Hazra and Choudhuri 2017, Ruždjak et al 2017). This suggests that stronger cycles may be associated with a stronger modulation pattern, introducing a nonlinearity into the flux transport dynamo model (Jiang et al 2010b, Cameron and Schüssler 2012a, Karak and Choudhuri 2012). The relationship between activity level and flow modulation, however, seems more complex than a simple proportionality (Komm et al 2017). In particular, the modulation signal in the Cycle 24 activity belt seems to be too strong in comparison with the low amplitude of the sunspot cycle.

In addition to a variation in the amplitude of migrating flow modulations, their migration speed may also influence the cycle. Howe et al (2009) pointed out that in the minimum of Cycle 24 the equatorward drift of the torsional oscillation shear belt corresponding to the active latitude of the cycle was slower than in the previous minimum. They suggested that this slowing may explain the belated start of cycle 24.

Under the assumption that meridional flow modulations are the main factor controlling the buildup of the poloidal field from AR sources, Hung et al (2015), Hung et al (2017) suggest an inverse approach to derive flow variations from magnetic data. As, however, we will see in the next subsection, the validity of the underlying assumption is open to question.

In summary: while magnetically induced modulations of the meridional flow and their effect on flux transport may be a potentially important nonlinear feedback mechanism controlling intercycle activity variations, the limited observational record and the apparent complexities of the interplay have as yet not permitted their use as a precursor.

2.5.2 Active region tilts

As luck would have it, soon after the first version of this review was finished, a paper was published that prompted a flurry of activity in a completely new field. In the paper Dasi-Espuig et al (2010) analysed sunspot group catalogues extracted from the Mt.Wilson and Kodaikanal photoheliograms for solar cycles 15–21, focusing on the distribution of tilt angles of the longitudinal axis of bipolar active regions to the azimuthal direction. Area-weighted averages of the tilt angles of sunspot groups were calculated in latitudinal bins, then normalized by the latitude to yield a tilt parameter. Two effects emerged from the analysis:

(1) Tilt quenching (TQ): an anticorrelation between the amplitude of a solar cycle and its mean tilt parameter.

(2) Tilt precursor (TP): the product of the amplitude of a sunspot cycle with its mean tilt parameter turned out to be a good predictor of the amplitude of the subsequent cycle. Muñoz-Jaramillo et al (2013b) further demonstrated that this product is also a good predictor of the polar magnetic flux at minimum (as reconstructed from polar facular counts), suggesting that the predictive potential of this method is based on the role of tilt angles in controlling the amount of net flux transported towards the poles.

The combination of TP with a TQ relationship upon which random fluctuations in the tilt are superposed implies that intercycle variations in solar activity will be controlled by a nonlinear feedback mechanism, into which a stochastic element is incorporated. This realization prompted intense activity in the development of model-based cycle prediction, to be discussed in the following section.

Ironically, in parallel with the major impact of the Dasi-Espuig et al (2010) paper on theoretical work, it was soon subjected to criticism on observational grounds. Ivanov (2012) repeated the analysis, now including the Pulkovo sunspot group catalogue (covering a shorter period, cycles 18–21). The TQ effect in the Mt.Wilson data set was found to depend crucially on the low value of the mean tilt for the anomalous strong cycle 19. In the Kodaikanal data the effect appeared more robustly but it still seems to depend on high tilt values in cycles 15 and 16, for which the Mt.Wilson set yields lower values. The Pulkovo data are consistent with the Kodaikanal series but they only start with Cycle 18, so no definitive conclusion was drawn from them.

In their hemispherically separated analysis of the Mt.Wilson data McClintock and Norton (2013) found that in Cycle 19 a strong suppression of the tilt was only present in the Southern hemisphere. Accordingly, the TQ effect is only seen in the South.

Kitchatinov and Olemskoy (2011) analysed the Pulkovo data set, focusing on the TP effect. Instead of considering average tilts, they evaluated the area-weighted latitude difference between leading and trailing subgroups and averaged this quantity for each cycle in their data set. While data were available for 3 cycles only, they confirmed the good correlation between this predictor and a measure (the dipole–octupole index) of the amplitude of the polar magnetic field in the next minimum.

A fourth tilt database was compiled by Baranyi (2015) based on the Debrecen sunspot catalogues, while a fifth set of tilts was measured by Işık et al (2018) from solar drawings made at Kandilli Observatory for cycles 19–24. Overall, the results for these cycles seem to compare well for the Kodaikanal, Pulkovo, Debrecen and Kandilli data, but as cycles 15 and 16 are only covered by the Kodaikanal data, it is no surprise that the TQ effect is not clearly seen in the shorter Debrecen and Kandilli data sets.

All the previously considered data sets were based on sunspot positions alone, without magnetic polarity information. Tilts of active regions taking into account the magnetic polarities of spots were determined by Li and Ulrich (2012) from Mt.Wilson and MDI magnetograms. McClintock et al (2014) compared these measurements with the Debrecen tilt data, focusing on anti-Hale regions which are the major reason for the discrepancies. The occurrence rate of anti-Hale regions was found to be 8.5%. Tilts of active regions taking into account the magnetic polarities of spots were recently also determined by Tlatova et al (2018) from the archive of solar drawings at Mt.Wilson Observatory: these drawings include polarity information since 1917. Their results concerning cycle dependence of tilts are inconclusive.

In summary, the suspected tilt quenching and tilt precursor effects have opened new directions in sunspot cycle forecasting, but the evidence is still controversial, especially for TQ. These inconclusive results underline the importance of obtaining more data and, especially, longer data sets. Contributions towards this goal like that of Senthamizh Pavai et al (2015b), who determined tilts for cycles 7–10 from sunspot drawings by Schwabe, are therefore highly valuable and may be key for a later clarification of these issues.

In view of the doubts regarding the reality of TQ, attention has recently started to turn towards alternatives. One such alternative, proposed by J. Jiang, might be latitude quenching: the mean latitude of active regions in a given cycle was found to be positively correlated with cycle amplitudes (Jiang et al 2011a). From higher latitudes a lower fraction of leading flux will manage to diffuse across the equator, leaving less trailing flux to contribute to the polar fields. Therefore, the correlation found represents a negative feedback effect —hence the name latitude quenching.

2.5.3 Alternatives to tilt quenching: inflows and latitude quenching

In view of the doubts regarding the reality of TQ, attention has recently started to turn towards alternatives.

The tilts of active regions are usually attributed to the effect of Coriolis force on the rising magnetic flux loops (Schmidt 1968; D’Silva and Choudhuri 1993; Fisher et al 1995; Pevtsov et al 2014;). The tilt angle varies inversely with the initial field strength, making the suggested tilt quenching effect quite plausible: stronger toroidal fields should simply come with weaker tilt. (Note that the tilt also depends on the entropy of the rising loop: a cycle-dependent variation in thermal properties of the layer where toroidal flux is stored may therefore also explain TQ, as shown by Işık 2015.)

An alternative explanation for tilt quenching was put forward by Jiang et al (2010b) and Cameron and Schüssler (2012a) who pointed out that the meridional inflows directed towards the activity belts, discussed in Sect. 2.5.1 will tend to reduce the AR tilts.

On the other hand, the activity-related meridional flow pattern may partly originate from the superposition of more localized circular inflows towards individual active regions (Haber et al 2004, Švanda et al 2008, Löptien et al 2017). It is unclear whether such concentric flows can exert a torque on the flux loops that may reduce their tilt; however this may not even be necessary: simulations by Martin-Belda and Cameron (2017a) suggested that the hindering effect on the inflows on the separation of leading and trailing polarity fluxes is sufficient to significantly reduce the amplitude of the dipolar seed field being built up for the next cycle. This inflow effect might then even completely substitute TQ as the main nonlinearity mechanism controlling intercycle variations

Another alternative to TQ, proposed by J. Jiang, might be latitude quenching: the mean latitude of active regions in a given cycle was found to be positively correlated with cycle amplitudes (Jiang et al 2011a). From higher latitudes a lower fraction of leading flux will manage to diffuse across the equator, leaving less trailing flux to contribute to the polar fields. Therefore, the correlation found represents a negative feedback effect —hence the name latitude quenching.

3 Model-based predictions

3.1 Surface flux transport models

Surface flux transport (SFT) models describe the transport of magnetic flux across the solar surface, modelling it as an advective-diffusive transport:

[TABLE]

where $\theta$ and $\phi$ are colatitude and longitude, $B$ is the radial component of the magnetic flux density, $\eta$ is the turbulent magnetic diffusivity, $v$ is the speed of the meridional flow, and $\Omega(\theta)$ is the differential rotation profile.

Equation (14) can be interpreted as the radial component of the induction equation where the neglected radial terms have been replaced by the source term $S$ describing flux emergence and the heuristic decay term $B/\tau$, supposed to represent vertical diffusion. The latter term is only occasionally used, to improve agreement with observations.

SFT models were first developed in the 1980s for the interpretation of the then newly available synoptic magnetogram record. This “age of enlightenment”, to use the expression of Sheeley’s (2005) historical review, was followed by less intense activity in the field until, from about 2010, a renaissance of SFT modelling ensued. The revival was prompted by the increasing acceptance of the polar precursor as the most reliable physical precursor technique: SFT models offered a way to “predict the predictor”, promising to increase the temporal scope of forecasts.

For detailed discussions of the advance made in SFT modelling we refer the reader to the reviews by Jiang et al (2014b) and Wang (2017). Here we only present brief highlights of the main results from the past decade.

3.1.1 Parameter optimization

For an application of the method, the parameters in equation (14) such as $\eta$, $v(\theta)$ or $\tau$ need to be specified. This problem has reopened the issue of parameter optimization (Lemerle et al 2015, Virtanen et al 2017, Whitbread et al 2017). Choosing the right parameters is a nontrivial task as the outcome depends on the data used for calibration and on the choice of the merit function for the optimization.

Recently Petrovay and Talafha (2019) presented the results of a large-scale systematic study of the parameter space in an SFT model where the source term representing the net effect of tilted flux emergence was chosen to represent a typical, average solar cycle as described by observations, comparing the results with observational constraints on the spatiotemporal variation of the polar magnetic field. It was found that without a significant decay term in the SFT equation (i.e., for $\tau>10$ yr) the global dipole moment reverses too late in the cycle for all flow profiles and parameters, providing independent supporting evidence for the need of a decay term, even in the case of identical cycles. An allowed domain is found to exist for $\tau$ values in the 5–10 yr range for all flow profiles considered. Generally higher values of $\eta$ (500–800) are preferred though some solutions with lower $\eta$ are still allowed.

3.1.2 Nonlinear feedback effects

Nonlinear feedback mechanisms such as TQ (Cameron et al 2010), a variable modulation of the meridional flow in the form of inflow belts flanking the active latitudes (Jiang et al 2010b, Cameron and Schüssler 2012b), concentric inflows around individual AR (Martin-Belda and Cameron 2016, 2017b) or latitude quenching (Jiang et al 2011b) have been considered in many studies. With the right parameterization, feedback was found to allow magnetic field reversals even when a stronger cycle is followed by a considerably weaker cycle. [A difficulty in obtaining reversal in such situations was a main motivation for the introduction of the decay term in equation (14).]

It is a curious circumstance that TQ currently has stronger support from models than from observations. Indeed, Cameron et al (2010) found that a combination of TQ with an SFT model results in polar field strengths that approximately correctly predict the amplitude of the next cycle for Cycles 15/16–21/22. The same holds also if the TQ results in the model as a consequence of meridional inflow belts (Cameron and Schüssler 2012b).

3.1.3 Extending the polar precursor

Sunspot observations are available for several centuries, while proxy indicators of the polar field strength start from the late 19th century only. Using sunspot records to reconstruct the source term in SFT models, the polar field may also be reconstructed for cycles where no polar field proxies are known. Such a programme was carried out by Jiang et al (2011a) who reconstructed the butterfly diagram from the sunspot number record for the period since 1700; properties of the butterfly “wings” were correlated with the amplitudes and lengths of solar cycles based on the GPR sunspot catalogue, then these correlations were used for reconstruction in earlier times. In a subsequent work (Jiang et al, 2011b), the same authors used this butterfly diagram as a source in an SFT model to determine polar fields. The derived polar field values were found to correlate rather well with the amplitude of the subsequent cycle, thereby extending the period of time for which we have evidence for the polar precursor.

A deficiency of the source reconstruction based on sunspot catalogues is the lack of information on magnetic polarities and on the distribution of weaker plage fields. A promising new method based on the use of Ca ii K synoptic maps combined with available sunspot magnetic field measurements was recently successfully tested by Virtanen et al (2019).

3.1.4 The importance of fluctuations: rogue active regions

The nonlinear feedback effects discussed above define an essentially deterministic mechanism of intercycle variations. While on the level of the spatiotemporal distribution of individual active regions there may be numerous different realizations of a sunspot cycle with a given sunspot number profile, one might expect that the variation of at least the statistical average taken over all realizations can be reliably predicted.

Starting from 2013, however, it was gradually realized that the behaviour of an individual realization (like the real Sun) can strongly deviate from the statistical expectations. The magnetic flux of the largest ARs in a solar cycle is comparable to the flux in the polar caps where the polar magnetic field is concentrated around the minimum. The amplitude of the polar field built up during a cycle is therefore highly sensitive to the exact balance of leading vs. trailing polarity flux transported to the poles and cancelled by cross-equator diffusion. Major plumes on the observational magnetic butterfly diagram were identified by Cameron et al (2013) as originating from large cross-equatorial AR where cancellation between the two polarities is minimal due to their advection in opposite directions by the meridional flow diverging on the equator. Such cross-equatorial plumes were incorporated into an SFT model by Cameron et al (2014).

The role of random scatter in active region properties was further investigated by Jiang et al (2014a) who showed that the dipole contribution of a single active region drops quite fast with heliographic latitudes. For those low-latitude active regions with an unusually large contribution to the global dipole Nagy et al (2017) introduced the name “rogue” AR. Representing an active region as a simple bipole with tilt angle $\alpha$ relative to the azimuthal direction, it is straightforward to show from equation (12) that its contribution to the axial dipole is

[TABLE]

where $F$ is the magnetic flux and $d$ is the angular separation of the two polarities. Large values of $F$, $d$ and/or $\alpha$ are therefore conditions for a significant “dynamo effectivity” of an AR, i.e., a significant contribution to the polar field built up in the cycle. However, equation (15) only gives the initial contribution to the dipole moment. For AR further from the equator cross-equatorial cancellation of the leading polarity will be less efficient, fluxes of both polarities will be largely transferred to the pole and their net effect will mostly cancel: this is the reason why the final dipole contributions and therefore the dynamo effectivity of AR also drops quite fast with heliographic latitude.

3.1.5 Explaining the end of the Modern Maximum

Beside the general investigations discussed, a special objective of SFT modelling efforts was to correctly “hindcast” the unusually weak polar fields in the minimum of Cycle 24 that brought the Modern Maximum to an end. Initial efforts (Yeates 2014, Upton and Hathaway 2014a) encountered difficulties in reproducing the polar field, until Jiang et al (2015) were finally able to correctly reproduce the evolution of the polar field by incorporating in their source term individual observed active regions (modelled as idealized bipoles, but with tilt values, fluxes and separations derived from observations). After carefully excluding recurrent ARs from the source term they found that the chief responsibility for the deviation of the polar flux from its expected value lies with a low number of large low-latitude rogue AR with non-Hale or non-Joy orientations.

In a similar research Upton and Hathaway (2014b) focused on the predictability of the evolution of the axial dipole moment. From some selected instant onwards, they substitute actual ARs with the ARs of another solar cycle of similar amplitude, and they find that the dipole evolution can be well predicted over $\sim 3$ years. It should be noted, however, that this study does not cover the period 2003–2006 which seems to have been crucial in the development of an anomalously low dipole moment at the end of Cycle 23. Predictability was also considered by Whitbread et al (2018) who addressed the issue how many AR need to be taken into account to reproduce the dipole moment at the end of a cycle. Ordering active regions by decreasing dipole moment contribution they found that, despite the fact that the polar magnetic flux is comparable to the flux a single large AR, the dipole moment at the end of a cycle cannot be reproduced without accounting for the net contribution from hundreds of active regions during the cycle.

The tilt of active regions is a manifestation of the writhe of the underlying flux loop, and writhe is one form of helicity (cf. Petrovay et al 2006), while nonzero helicity is a condition for free magnetic energy available for eruptions. On this ground Petrovay and Nagy (2018b) tentatively suggested that there may be a large overlap between rogue AR and flare/CME-productive AR. A relevant study was recently undertaken by Jiang et al (2019) who found that flare productivity and dynamo effectivity of ARs are governed by different parameters. Flare productivity primarily depends on the structural complexity of ARs, large flares being much more common $\delta$-spots, while the dipole moment contribution of an AR, which ultimately determines its effect on the dynamo, is determined by the latitudinal separation of polarities. So while there is indeed a large overlap between the flare-productive ARs and “rogue” or exceptional ARs, the two characteristics do not necessarily go hand in hand.

3.1.6 Forecasting Cycle 25

The buildup of the polar dipole moment during the ongoing Cycle 24 has been followed with keen attention. Researchers analysed a number of important episodes (e.g., Petrie and Ettinger 2017). Yeates et al (2015) examined the origin of a prominent poleward surge in the magnetic butterfly diagram in 2010–11 by a combination of analysis of observational data and SFT simulations, concluding that the episode is not expected to have a major imact of the dipole buildup. Sun et al (2015) presented an observational analysis of the polar reversal process in Cycle 24. This is of particular interest owing to the ill-defined nature of polarity reversal in the N hemisphere: the field strength here lingered around zero for well over two years until it finally started to increase towards the end of 2014. (It may be worth noting that this was correctly predicted in the SFT simulation of Upton and Hathaway 2014b.) As a result, the phase shift between the hemispheres has changed sign: while activity peaked first on the N hemisphere in the last few cycles, indications are that in Cycle 25 activity will first peak in the South (cf. also Labonville et al 2019). The phase shift was the consequence of a few surges of opposite polarity that hindered the growth of flux in the N polar region.

Going further than focusing on selected events, a number of researchers attempted to predict Cycle 25 by incorporationg all ARs of Cycle 24 in the source term of an SFT simulation and modelling further evolution under some more or less plausible assumptions.

Allowing a random scatter in AR tilts and also in the time-dependent meridional flow Hathaway and Upton (2016) and Upton and Hathaway (2018) arrived at a prediction similar to or marginally lower than Cycle 24. A similar conclusion (Cycle 25 somewhat weaker than Cycle 24) was reached by Iijima et al (2017) who, based on the plateau-like nature of the dipole moment maximum discussed in Sect. 2.3 above, assume no further contributions to the dipole moment.

Considering 50 different random realizations drawn from a statistical ensemble of ARs Cameron et al (2016) predict that Cycle 25 will be similar or slightly stronger than Cycle 24. In a similar later study with improved technical details Jiang et al (2018) arrived at the prediction that Cycle 25 will peak in the range 93 to 159 (see also Jiang and Cao 2018).

3.2 The solar dynamo: a brief summary of current models

While attempts to predict future solar cycles on the basis of the empirical sunspot number record have a century-old history, predictions based on physical models of solar activity only started one solar cycle ago. The background of this new trend is, however, not some significant improvement in our understanding of the solar dynamo. Rather, it is the availability of increasingly fast new computers that made it possible to fine-tune the parameters of certain dynamo models to reproduce the available sunspot record to a good degree of accuracy and to apply data assimilation methods (such as those used in terrestrial weather prediction) to these models. This is not without perils. On the one hand, the capability of multiparametric models to fit a multitude of observational data does not prove the conceptual correctness of the underlying model. On the other hand, in chaotic or stochastic systems such as the solar dynamo, fitting a model to existing data will not lead to a good prediction beyond a certain time span, the extent of which can only be objectively assessed by “postdiction” tests, i.e., checking the models predictive skill by trying to “predict” previous solar cycles and comparing those predictions to available data. Apparently successful postdiction tests have led some groups to claim a breakthrough in solar cycle prediction owing to the model-based approach (Dikpati and Gilman, 2006; Kitiashvili and Kosovichev, 2008). Yet, as we will see in the following discussion, a closer inspection of these claims raises many questions regarding the role that the reliance on a particular physical dynamo model plays in the success of their predictions.

Extensive summaries of the current standing of solar dynamo theory are given in the reviews by Petrovay (2000), Ossendrijver (2003), Solanki et al (2006), Charbonneau (2010) and Cameron et al (2017). As explained in detail in those reviews, all current models claiming to acceptably represent the solar dynamo are based on the mean-field theory approach wherein a coupled system of partial differential equations governs the evolution of the toroidal and poloidal components of the large-scale magnetic field. Until recently, the large-scale field was assumed to be axially symmetric in practically all models. In some nonlinear models the averaged equation of motion, governing large-scale flows is also coupled into the system.

In the simplest case of homogeneous and isotropic turbulence, where the scale $l$ of turbulence is small compared to the scale $L$ of the mean variables (scale separation hypothesis), the dynamo equations have the form

[TABLE]

Here $\mathbf{B}$ and $\mathbf{U}$ are the large-scale mean magnetic field and flow speed, respectively; $\eta_{T}$ is the magnetic diffusivity (dominated by the turbulent contribution for the highly conductive solar plasma), while $\alpha$ is a parameter related to the non-mirror symmetric character of the magnetized plasma flow.

In the case of axial symmetry the mean flow $\mathbf{U}$ may be split into a meridional circulation $\mathbf{U}_{c}$ and a differential rotation characterized by the angular velocity profile $\Omega_{0}(r,\theta)$:

[TABLE]

where $r$, $\theta$, $\phi$ are spherical coordinates and $\mathbf{e}_{\phi}$ is the azimuthal unit vector. Now introducing the shear

[TABLE]

assuming $\alpha\ll\Omega L$ and ignoring spatial derivatives of $\alpha$ and $\eta_{T}$, Eq. (16) simplifies to the pair

[TABLE]

[TABLE]

where $B$ and $A$ are the toroidal (azimuthal) components of the magnetic field and of the vector potential, respectively, and $\frac{\partial A}{\partial x}$ is to be evaluated in the direction $90^{\circ}$ clockwards of $\vec{\Omega}$ (along the isorotation surface) in the meridional plane. These are the classic $\alpha\Omega$ dynamo equations, including a meridional flow.

In the more mainstream solar dynamo models the strong toroidal field is now generally thought to reside near the bottom of the solar convective zone. Indeed, it is known that a variety of flux transport mechanisms such as pumping (Petrovay, 1994) remove magnetic flux from the solar convective zone on a timescale short compared to the solar cycle. Following earlier simpler numerical experiments, MHD numerical simulations have indeed demonstrated this pumping of large scale magnetic flux from the convective zone into the tachocline below, where it forms strong coherent toroidal fields (Browning et al 2006, Fan and Fang 2014, Warnecke et al 2018). As this layer is also where rotational shear is maximal, it is plausible that the strong toroidal fields are not just stored but also generated here, by the winding up of poloidal field.181818Some doubts concerning this deep-seated storage of toroidal flux recently arose as in at least one MHD simulation (Nelson et al 2013) the toroidal flux was found to be floating in a wreath-like structure in the middle of the convective envelope. The two main groups of dynamo models, interface dynamos and flux transport dynamos, differ mainly in their assumptions about the site and mechanism of the $\alpha$-effect responsible for the generation of a new poloidal field from the toroidal field.

In interface dynamos $\alpha$ is assumed to be concentrated near the bottom of the convective zone, in a region adjacent to the tachocline, so that the dynamo operates as a wave propagating along the interface between these two layers. While these models may be roughly consistent and convincing from the physical point of view, they have only had limited success in reproducing the observed characteristics of the solar cycle, such as the butterfly diagram.

Flux transport dynamos, in contrast, rely on the Babcock–Leighton mechanism for $\alpha$, arising due to the action of the Coriolis force on emerging flux loops, and they assume that the corresponding $\alpha$-effect is concentrated near the surface. They keep this surface layer incommunicado with the tachocline by introducing some arbitrary unphysical assumptions (such as very low diffusivities in the bulk of the convective zone). The poloidal fields generated by this surface $\alpha$-effect are then advected to the poles and there down to the tachocline by the meridional circulation – which, accordingly, has key importance in these models. The equatorward deep return flow of the meridional circulation is assumed to have a significant overlap with the tachocline (another controversial point), and it keeps transporting the toroidal field generated by the rotational shear towards the equator. By the time it reaches lower latitudes, it is amplified sufficiently for the flux emergence process to start, resulting in the formation of active regions and, as a result of the Babcock–Leighton mechanism, in the reconstruction of a poloidal field near the surface with a polarity opposed to that in the previous 11-year cycle. While flux transport models may be questionable from the point of view of their physical consistency, they can be readily fine-tuned to reproduce the observed butterfly diagram quite well. This ready adaptability made flux transport dynamos hugely popular in the research community (Charbonneau 2007, Karak et al 2014). In flux transport dynamos the $\alpha$ term is usually nonlocal (generating poloidal field at the surface out of toroidal field at the bottom of the convective zone).

It should be noted that while the terms “interface dynamo” and “flux transport dynamo” are now very widely used to describe the two main approaches, the more generic terms “advection-dominated” and “diffusion-dominated” would be preferable in several respects. This classification allows for a continuous spectrum of models depending on the numerical ratio of advective and diffusive timescales (for communication between surface and tachocline). In addition, even at the two extremes, classic interface dynamos and circulation-driven dynamos are just particular examples of advection or diffusion dominated systems with different geometrical structures.

3.3 Is model-based cycle prediction feasible?

As it can be seen even from the very brief and sketchy presentation given above, all current solar dynamo models are based on a number of quite arbitrary assumptions and depend on a number of free parameters, the functional form and amplitude of which is far from being well constrained. For this reason, Bushby and Tobias (2007) rightfully say that all current solar dynamo models are only of “an illustrative nature”. This would suggest that as far as solar cycle prediction is concerned, the best we should expect from dynamo models is also an “illustrative” reproduction of a series of solar cycles with the same kind of long-term variations (qualitatively and, in the statistical sense, quantitatively) as seen in solar data. Indeed, Bushby and Tobias (2007) demonstrated that even a minuscule stochastic variation in the parameters of a particular flux transport model can lead to large, unpredictable variations in the cycle amplitudes. And even in the absence of stochastic effects, the chaotic nature of nonlinear dynamo solutions seriously limits the possibilities of prediction, as the authors find in a particular interface dynamo model: even if the very same model is used to reproduce the results of one particular run, the impossibility of setting initial conditions exactly representing the system implies that predictions are impossible even for the next cycle. Somewhat better results are achieved by an alternative method, based on the phase space reconstruction of the attractor of the nonlinear system – this is, however, a purely empirical time series analysis technique for which no knowledge of the detailed underlying physics is needed. (Cf. Sect. 4.3 above.)

Despite these very legitimate doubts regarding the feasibility of model-based prediction of solar cycles, in recent years several groups have claimed to be able to predict upcoming solar cycles on the basis of dynamo models with a high confidence. So let us consider these claims.

3.4 Explicit models

3.4.1 Axisymmetric models

The first attempt at model-based solar cycle prediction was the work of the solar dynamo group in Boulder (Dikpati et al, 2006; Dikpati and Gilman, 2006). Their model is a flux transport dynamo, advection-dominated to the extreme. The strong suppression of diffusive effects is assured by the very low value (less than 20 km${}^{2}/$s) assumed for the turbulent magnetic diffusivity in the bulk of the convective zone. As a result, the poloidal fields generated near the surface by the Babcock–Leighton mechanism are only transported to the tachocline on the very long, decadal time scale of meridional circulation. The strong toroidal flux residing in the low-latitude tachocline, producing solar activity in a given cycle is thus the product of the shear amplification of poloidal fields formed near the surface about 2 – 3 solar cycles earlier, i.e., the model has a “memory” extending to several cycles. The mechanism responsible for cycle-to-cycle variation is assumed to be the stochastic nature of the flux emergence process. In order to represent this variability realistically, the model drops the surface $\alpha$-term completely (a separate, smaller $\alpha$ term is retained in the tachocline); instead, the generation of poloidal field near the surface is represented by a source term, the amplitude of which is based on the sunspot record, while its detailed functional form remains fixed.

Dikpati and Gilman (2006) found that, starting off their calculation by fixing the source term amplitudes of sunspot cycles 12 to 15, they could predict the amplitudes of each subsequent cycle with a reasonable accuracy, provided that the relation between the relative sunspot numbers and the toroidal flux in the tachocline is linear, and that the observed amplitudes of all previous cycles are incorporated in the source term for the prediction of any given cycle. For cycle 24 the model predicted peak smoothed annual relative sunspot numbers of 150 (v1) or more. Elaborating on their model, they proceeded to apply it separately to the northern and southern hemispheres, to find that the model can also be used to correctly forecast the hemispheric asymmetry of solar activity (Dikpati et al, 2007).

Even though, ultimately, the prediction proved to be off by about a factor of 2, the extraordinary claims of this pioneering research prompted a hot debate in the dynamo community. Besides the more general, fundamental doubt regarding the feasibility of model-based predictions (see Sect. 4.2 above), more technical concerns arose. In order to understand the origin of the predictive skill of the Boulder model, Cameron and Schüssler (2007) studied a version of the model simplified to an axially symmetric SFT model, wherein only the equation for the radial field component is solved as a function of time and latitude. The equation includes a source term similar to the one used in the Boulder model. As the toroidal flux does not figure in this simple model, the authors use the cross-equatorial flux $\Phi$ (the amount of magnetic flux crossing the equator in a given time) as a proxy, arguing that this may be more closely linked to the amplitude of the toroidal field in the upcoming cycle than the polar field. They find that $\Phi$ indeed correlates quite well (correlation coefficients $r\sim 0.8\mbox{\,--\,}0.9$, depending on model details) with the amplitude of the next cycle, as long as the form of the latitude dependence of the source term is prescribed and only its amplitude is modulated with the observed sunspot number series (“idealized model”). But surprisingly, the predictive skill of the model is completely lost if the prescribed form of the source function is dropped and the actually observed latitude distribution of sunspots is used instead (“realistic model”). Cameron and Schüssler (2007) interpreted this by pointing out that $\Phi$ is mainly determined by the amount of very low latitude flux emergence, which in turn occurs mainly in the last few years of the cycle in the idealized model, while it has a wider temporal distribution in the realistic model. The conclusion is that the root of the apparently good predictive skill of the truncated model (and, by inference, of the Boulder model it is purported to represent) is actually just the good empirical correlation between late-phase activity and the amplitude of the next cycle, discussed in Sect. 2.1 above. This correlation is implicitly “imported” into the idealized flux transport model by assuming that the late-phase activity is concentrated at low latitudes, and therefore gives rise to cross-equatorial flux which then serves as a seed for the toroidal field in the next cycle. So if Cameron and Schüssler (2007) are correct, the predictive skill of the Boulder model is due to an empirical precursor and is thus ultimately explained by the good old Waldmeier effect (cf. Sect. 1.3.3) In view of the fact that no version of the Boulder model with a modified source function incorporating the realistic latitudinal distribution of sunspots in each cycle was ever presented, this explanation seems to be correct, despite the fact that the effective diffusivity represented by the sink term in the reduced model is significantly higher than in the Boulder model, and consequently, the reduced model will have a more limited memory, cf. Yeates et al (2008).

Another flux transport dynamo code, the Surya code, originally developed by A.  Choudhuri and coworkers in Bangalore, has also been utilized for prediction purposes. The crucial difference between the two models is in the value of the turbulent diffusivity assumed in the convective zone: in the Bangalore model this value is 240 km${}^{2}/$s, 1 – 2 orders of magnitude higher than in the Boulder model, and within the physically plausible range (Chatterjee et al, 2004). As a result of the shorter diffusive timescale, the model has a shorter memory, not exceeding one solar cycle. As a consequence of this relatively rapid diffusive communication between surface and tachocline, the poloidal fields forming near the surface at low latitudes due to the Babcock–Leighton mechanism diffuse down to the tachocline in about the same time as they reach the poles due to the advection by the meridional circulation. In these models, then, polar magnetic fields are not a true physical precursor of the low-latitude toroidal flux, and their correlation is just due to their common source. In the version of the code adapted for cycle prediction (Choudhuri et al, 2007; Jiang et al, 2007), the “surface” poloidal field (i.e., the poloidal field throughout the outer half of the convection zone) is rescaled at each minimum by a factor reflecting the observed amplitude of the Sun’s dipole field. The model showed reasonable predictive skill for the last three cycles for which data are available, and could even tackle hemispheric asymmetry (Goel and Choudhuri, 2009). For Cycle 24, the predicted amplitude was 30 – 35% lower than Cycle 23.

In retrospect, these first attempts at model-based solar cycle prediction are now generally seen as precursor methods in disguise. Cameron and Schüssler (2007) convincingly argued that the apparent predictive (postdictive) skill of the Boulder model was related to the “minimax” family of internal precursors, based on the Waldmeier law combined with the overlap of consecutive cycles (see Sect. 2.1). The predictive skill of the Surya model, on the other hand, is based on the polar precursor: essentially, the polar field amplitude at minimum is set to its observed value, and the model simply mechanically winds up the poloidal field into a toroidal field, ensuring a proportionality.

A more sophisticated approach to dynamo based predictions was spurred by the realization of the importance of nonlinearities like TQ and stochastic effects in SFT models . In particular, the significance of the effect of individual active regions in the buildup of poloidal fields gave an impetus to the development of non-axisymmetric models where individual AR can be treated. Nevertheless, Kitchatinov et al (2018) recently constructed an axisymmetric flux transport dynamo where stochastic effects were still retained as fluctuations of the $\alpha$ parameter. It was found that with a correlation time on the order of a month the model is able to mimick the effects of rogue AR such as intercycle variations and even the triggering of grand minima, as in the model of Nagy et al (2017).

3.4.2 Nonaxisymmetric models

The first dynamo model to explicitly deal with individual AR was constructed by Yeates and Muñoz-Jaramillo (2013). Emerging flux loops in this model were created by imposed upflows with properties calibrated to reproduce observed AR characteristics. Another 3D code capable of dealing with individual AR, the STABLE code was developed in Boulder (Miesch and Dikpati 2014, Miesch and Teweldebirhan 2016).

The third such code is the 2 $\times$ 2D code of the Montreal group (Lemerle et al 2015, Lemerle and Charbonneau 2017). This latter code ingenuously combines an axially symmetric flux transport dynamo code with a 2D SFT code. The dynamo code couples to the SFT code by an emergence function specifying the locations and properties of the randomly created bipolar source regions for the SFT, based on the distribution of the toroidal field. The azimuthally averaged magnetic field resulting from the SFT component, in turn, provides the upper boundary condition for the flux transport model. The model is computationally highly efficient, optimizable, and it was carefully calibrated to reproduce the observed statistics of active regions in Cycle 21. It can be run for hundreds or thousands of cycles, thereby providing a data base for analysis that far exceeds the observational record of solar cycles. While the model’s behaviour also displays some differences relative to the real Sun (esp. in the phase relationship of the poloidal and toroidal fields), intercycle variations comparable to those seen in the real Sun are displayed by the solutions, including Dalton-like minima and grand minima (Lemerle and Charbonneau 2017, Nagy et al 2017).

Investigations in close analogy with this were also undertaken using the STABLE model. Karak and Miesch (2017) considered the effect of tilt quenching and tilt scatter, and found that the scatter can be a main factor behind long-term activity variations, while even a subtle tilt quenching effect is sufficient to limit the growth of the dynamo. Karak and Miesch (2018) found that fluctuations in AR tilt can be responsible for both the onset of and recovery from grand minima. While these results are in general agreement with those obtained with other codes, in another work based on the STABLE model Hazra et al (2017) reported that the dynamo effectivity of individual ARs increases with heliographic latitude, in contradiction to other results and to expectations based on the importance of cross-equatorial diffusion for the removal of leading polarity flux. It may be that diffusion of weak magnetic flux through the surface might be responsible for these results — an effect that, by construction, is avoided in the 2 $\times$ 2D model and which can also be suppressed by introducing downwards pumping of the magnetic field (Karak and Cameron 2016).

These nonaxisymmetric dynamo models have now reached a level where actual physics-based prediction is within reach. A new generation of dynamo-based cycle predictions is heralded by the work of Labonville et al (2019) who used data assimilation in the 2 $\times$ 2D to arrive to a state closely mimicking the current state of the solar dynamo, then run the code further to see its future development. Repeating these runs with many different random realizations of the AR emergences the authors generated an ensemble of future solar dynamo models (Fig. 13). This enabled them to give both a mean forecast and to attribute meaningful errors to their forecast. They predicted a peak $R$ value of $89^{+29}_{-14}$ for cycle 25, the maximum occurring in $2025.3^{+0.89}_{-1.05}$. The forecast is also broken down to hemispheres: activity in the North is predicted to peak 6 months later but at a 20 % higher level than in the South.

3.5 Truncated models

The “illustrative” nature of current solar dynamo models is nowhere more clearly on display than in truncated or reduced models where some or all of the detailed spatial structure of the system is completely disregarded, and only temporal variations are explicitly considered. This is sometimes rationalised as a truncation or spatial integration of the equations of a more realistic inhomogenous system; in other cases, no such rationalisation is provided, representing the solar dynamo by an infinite, homogeneous or periodic turbulent medium where the amplitude of the periodic large-scale magnetic field varies with time only.

In the present subsection we deal with models that do keep one spatial variable (typically, the latitude), so growing wave solutions are still possible — these models, then, are still dynamos even though their spatial structure is not in a good correspondence with that of the solar dynamo.

This approach in fact goes back to the classic migratory dynamo model of Parker (1955) who radially truncated (i.e., integrated) his equations to simplify the problem. Parker seems to have been the first to employ a heuristic relaxation term of the form $-B_{r}/\tau_{d}$ in the poloidal field equation to represent the effect of radial diffusion; here, $\tau_{d}=d^{2}/\eta_{T}$ is the diffusive timescale across the thickness $d$ of the convective zone. His model was generalized by Moss et al (2008) and Usoskin et al (2009b) to the case when the $\alpha$-effect includes an additive stochastic noise, and nonlinear saturation of the dynamo is achieved by $\alpha$-quenching. These authors do not make an attempt to predict solar activity with their model but they can reasonably well reproduce some features of the very long term solar activity record, as seen from cosmogenic isotope studies.

The other classic reduced dynamo model is that of Leighton (1969), the first mathematical formulation of the flux transport dynamo concept. An updated version of this dynamo, recently developed by Cameron and Schüssler (2017a) is a very promising, simple and versatile tool to reproduce many of the observed features of the solar dynamo with simple parameterizations. While it has not yet been utilized for cycle prediction, it has already proved to be useful in understanding long-term activity variations such as the dominant periods in hemispheric asymmetry (Schüssler and Cameron 2018).

Another radially truncated model, this time formulated in a Cartesian system, is that of Kitiashvili and Kosovichev (2009). In this model stochastic effects are not considered and, in addition to using an $\alpha$-quenching recipe, further nonlinearity is introduced by coupling in the Kleeorin–Ruzmaikin equation (Zeldovich et al, 1984) governing the evolution of magnetic helicity, which in the hydromagnetic case contributes to $\alpha$. Converting the toroidal field strength to relative sunspot number using the Bracewell transform, Eq. (3), the solutions reproduce the asymmetric profile of the sunspot number cycle. For sufficiently high dynamo numbers the solutions become chaotic, cycle amplitudes show an irregular variation. Cycle amplitudes and minimum–maximum time delays are found to be related in a way reminiscent of the Waldmeier relation.

Building on these results, Kitiashvili and Kosovichev (2008) attempted to predict solar cycles using a data assimilation method. The approach used was the so-called Ensemble Kalman Filter method. Applying the model for a “postdiction” of the last 8 solar cycles yielded astonishingly good results, considering the truncated and arbitrary nature of the model and the fundamental obstacles in the way of reliable prediction discussed above. The question may arise whether the actual physics of the model considered has any significant role in this prediction, or we are dealing with something like the phase space reconstruction approach discussed in Sect. 4.3 where basically any model with an attractor that looks reasonably similar to that of the actual solar dynamo would do. Either way, the method is remarkable. Its prediction for cycle 24 essentially proved to be correct. An early prediction for Cycle 25 (Kitiashvili 2016) yielded a maximum occurring in 2023/24 at a level $R=90\pm 15$ (v2 values), with the cycle starting in 2019/20. It is noteworthy that a forecast considering observational data up to the previous minimum (2008) already yielded quite similar results for the amplitude, although the timing of the maximum would have been expected earlier. Understanding the origin of this impressive apparent predictive skill certainly deserves more in-depth research.

3.6 The Sun as an oscillator

An even more radical simplification of the solar dynamo problem ignores any spatial dependence in the solutions completely, concentrating on the time dependence only. Spatial derivatives appearing in Eqs. (17) and (18) are estimated as $\nabla\sim 1/L$ and the resulting terms $U_{c}/L$ and $\eta_{T}/L^{2}$ as $1/\tau$ where $\tau$ is a characteristic time scale. This results in the pair

[TABLE]

[TABLE]

which can be combined to yield

[TABLE]

where $D=\alpha\Omega\tau^{2}/L$ is the dynamo number. For $D<1$, Eq. (21) clearly describes a damped linear oscillator. For $D>1$, solutions have a non-oscillatory character. The system described by Eq. (21), then, is not only not a true dynamo (missing the spatial dependence) but it does not even display growing oscillatory solutions that would be the closest counterpart of dynamo-like behaviour in such a system. Nevertheless, there are a number of ways to extend the oscillator model to allow for persistent oscillatory solutions, i.e., to turn it into a relaxation oscillator:

(1) The most straightforward approach is to add a forcing term $+\sin(\omega_{0}t)$ to the r.h.s. of Equations (21). Damping would cause the system to relax to the driving period $2\pi/\omega_{0}$ if there were no stochastic disturbances to this equilibrium. Hiremath (2006) fitted the parameters of the forced and damped oscillator model to each observed solar cycle individually; then in a later work (Hiremath, 2008) he applied linear regression to the resulting series to provide a forecast (see Sect. 4.1 above).

(2) Another trick is to account for the $\pi/2$ phase difference between poloidal and toroidal field components in a dynamo wave by introducing a phase factor $i$ into the first term on the r.h.s. of Eq. (20). This can also be given a more formal derivation as equations of this form result from the substitution of solutions of the form $A\propto e^{ikx}$, $B\propto e^{i(kx+\pi/2)}$ into the 1D dynamo equations. This route, combined with a nonlinearity due to magnetic modulation of differential rotation described by a coupled third equation, was taken by Weiss et al (1984). Their model displayed chaotic behaviour with intermittent episodes of low activity similar to grand minima. The generic normal form model recently introduced by Cameron and Schüssler (2017b) is also loosely related to this approach.

(3) Wilmot-Smith et al (2006) showed that another case where dynamo-like behaviour can be found in an equation like (21) is if the missing effects of finite communication time between parts of a spatially extended system are reintroduced by using a time delay $\Delta t$, evaluating the first term on the r.h.s. at time $t-\Delta t$ to get the value for the l.h.s. at time $t$. This time-delay approach has been further developed by Hazra et al (2014) and Turner and Ladde (2018).

(4) Yet another possibility is to introduce a nonlinearity into the model by assuming $D=D_{0}[1-f(B)]$ where $f(B=0)=0$ and $f\geq 0$ everywhere. (Note that any arbitrary form of $\alpha$- or $\Omega$-quenching can be cast in the above form by series expansion.) The governing equation then becomes one of a nonlinear oscillator:

[TABLE]

In the most commonly assumed quenching mechanisms the leading term in $f(B)$ is quadratic; in this case Eq. (22) describes a Duffing oscillator (Kanamaru, 2008). For large positive dynamo numbers, $D_{0}\gg 1$, then, the large nonlinear term dominates for high values of $B$, its negative sign imposing oscillatory behaviour; yet the origin is a repeller so the oscillation will never be damped out. The Duffing oscillator was first considered in the solar context by Paluš and Novotná (1999).

More generally, a complete spatial truncation of the nonlinear $\alpha\Omega$ dynamo equations with a dimensional analysis of some of the terms (see review by Lopes et al 2014) can be shown to give rise to a nonlinear oscillator equation of the form

[TABLE]

where $B$ is the amplitude of the toroidal magnetic field, and the parameters $\mu$, $\xi$ and $\gamma$ may be expressed by the dynamo parameters (dynamo number, meridional flow amplitude, nonlinearity parameters). This is a combination of the van der Pol and Duffing oscillators, the two most widely studied nonlinear oscillator problems.

Under certain conditions on the parameters, it is reduced to a van der Pol oscillator (Adomian, 1989; Mininni et al, 2001; Kanamaru, 2007):

[TABLE]

with $\mu>0$. From this form it is evident that the problem is equivalent to that of an oscillator with a damping that increases with amplitude; in fact, for small amplitudes the damping is negative, i.e., the oscillation is self-excited.

These simple nonlinear oscillators were among the first physical systems where chaotic behaviour was detected (when a periodic forcing was added). Yet, curiously, they first emerged in the solar context precisely as an alternative to chaotic behaviour. Considering the mapping of the solar cycle in the differential phase space $\{B,dB/dt\}$, Mininni et al (2000) got the impression that, rather than showing signs of a strange attractor, the SSN series is adequately modelled by a van der Pol oscillator with stochastic fluctuations. This concept was further developed by Lopes and Passos (2009) who fitted the parameters of the oscillator to each individual sunspot cycle. Subsequently, this parameter fitting was also exploited for cycle prediction purposes (Passos 2012).

The parameter $\mu$ is related to the meridional flow speed and the fit indicates that a slower meridional flow may have been responsible for the Dalton minimum. This was also corroborated in an explicit dynamo model (the Surya code) — however, as we discussed in Sect. 2.5.1, this result of flux transport dynamo models is spurious and the actual effect of a slower meridional flow is likely to be opposite to that suggested by the van der Pol oscillator model.

In an alternative approach to the problem, Nagovitsyn (1997) attempted to constrain the properties of the solar oscillator from its amplitude–frequency diagram, suggesting a Duffing oscillator driven at two secular periods. While his empirical reconstruction of the amplitude–frequency plot may be subject to many uncertainties, the basic idea is certainly noteworthy.

A further twist to the oscillator representation of the solar cycle is to consider a system of two coupled nonlinear oscillators (Kuramoto model). This approach has been taken in a series of papers by Blanter et al. (2014, 2016). The two variables were taken to represent the sunspot numbers and the geomagnetic $aa$ index, considered to be proxies for the toriodal and poloidal field amplitudes, respectively. In a later work (Blanter et al 2017) a coupled oscillator model was employed to account for the time variation of hemispheric asymmetries of solar activity, as known from observations (Norton and Gallagher 2010; Norton et al 2014).

In summary: despite its simplicity, the oscillator representation of the solar cycle is a relatively new development in dynamo theory, and its obvious potential for forecasting purposes is yet to be fully exploited.

4 Extrapolation methods

In contrast to precursor methods, extrapolation methods only use the time series of sunspot numbers (or whichever solar activity indicator is considered) but they generally rely on more than one previous point to identify trends that can be used to extrapolate the data into the future. They are therefore also known as time series analysis or, for historic reasons, regression methods.

A cornerstone of time series analysis is the assumption that the time series is homogeneous, i.e., the mathematical regularities underlying its variations are the same at any point of time. This implies that a forecast for, say, three years ahead has equal chance of success in the rising or decaying phase of the sunspot cycle, across the maximum or, in particular, across the minimum. In this case, distinguishing intracycle and intercycle memory effects, as we did in Sects. 1.3.2 and 2, would be meaningless. This concept of solar activity variations as a continuous process stands in contrast to that underlying precursor methods, where solar cycles are thought of as individual units lasting essentially from minimum to minimum, correlations within a cycle being considerably stronger than from one cycle to the next. While, as we have seen, there is significant empirical evidence for the latter view, the possibility of time homogeneity cannot be discarded out of hand. Firstly, if we consider the time series of global parameters (e.g., amplitudes) of cycles, homogeneity may indeed be assumed fairly safely. This approach has rarely been used for the directly observed solar cycles as their number is probably too low for meaningful inferences — but the long data sets from cosmogenic radionuclides are excellent candidates for time series analysis.

In addition, there may be good reasons to consider the option of homogeneity of solar activity data even on the scale of the solar cycle. Indeed, in dynamo models the solar magnetic field simply oscillates between (weak) poloidal and (strong) toroidal configuration: there is nothing inherently special about either of the two, i.e., there is no a priori reason to attribute a special significance to solar minimum. While at first glance the butterfly diagram suggests that starting a new cycle at the minimum is the only meaningful way to do it, there may be equally good arguments for starting a new cycle at the time of polar reversal. And even though SFT and dynamo models strongly suggest that spatial information regarding, e.g., the latitudinal distributions of sunspots may well be essential for cycle prediction, some studies point to a possibility to reconstruct this spatial information from time series alone (Jiang et al 2011a; Mandal et al 2017).

There is, therefore, plenty of motivation to try and apply standard methods of time series analysis to sunspot data.

Indeed, as the sunspot number series is a uniquely homogeneous and long data set, collected over centuries and generated in what has long been perceived to be a fairly carefully controlled manner, it has become a favorite testbed of time series analysis methods and is routinely used in textbooks and monographs for illustration purposes (Udny Yule, 1927; Box and Jenkins, 2008; Wei, 2005; Tong, 1993). This section will summarize the various approaches, proceeding, by and large, from the simplest towards the most complex.

4.1 Linear regression

Linear (auto)regression means representing the value of a time series at time $t$ by a linear combination of values at times $t-\Delta t$, $t-2\Delta t$, $\dots$, $t-p\Delta t$. Admitting some random error $\epsilon_{n}$, the value of $R$ in point $n$ is

[TABLE]

where $p$ is the order of the autoregression and the $c_{i}$’s are weight parameters. A further twist on the model admits a propagation of errors from the previous $q$ points:

[TABLE]

This is known as the ARMA (AutoRegressive Moving Average) model.

Linear regression techniques have been widely used for solar activity prediction during the course of an ongoing cycle. Their application for cycle-to-cycle prediction has been less common and successful (Lomb and Andersen, 1980; Box and Jenkins, 2008; Wei, 2005).

Brajša et al (2009) applied an ARMA model to the series of annual values of $R$. A successful fit was found for $p=6$, $q=6$. Using this fit, the next solar maximum was predicted to take place around 2012.0 with an amplitude $90\pm 27$, and the following minimum occurring in 2017.

Instead of applying an autoregression model directly to SSN data, Hiremath (2008) applied it to a forced and damped harmonic oscillator model claimed to well represent the SSN series. This resulted in a predicted amplitude of $110\pm 10$ for Solar Cycle 24, with the cycle starting in mid-2008 and lasting 9.34 years.

4.2 Spectral methods

“…the use of any mathematical algorithm to derive hidden periodicities from the data always entails the question as to whether the resulting cycles are not introduced either by the particular numerical method used or by the time interval analyzed.”

(de Meyer, 1981)

Spectral analysis of the sunspot number record is used for prediction under the assumption that the main reason of variability in the solar cycle is a long-term modulation due to one or more periods.

The usual approach to the problem is the purely formal one of representing the sunspot record with the superposition of eigenfunctions forming an orthogonal basis. From a technical point of view, spectral methods are a complicated form of linear regression. The analysis can be performed by any of the widely used means of harmonic analysis:

(1) Least squares (LS) frequency analysis (sometimes called “Lomb–Scargle periodogram”) consists in finding by trial and error the best fitting sine curve to the data using the least squares method, subtracting it (“prewhitening”), then repeating the procedure until the residuals become indistinguishable from white noise. The first serious attempt at sunspot cycle prediction, due to Kimura (1913), belonged to this group. The analysis resulted in a large number of peaks with dubious physical significance. The prediction given for the upcoming cycle 15 failed, the forecasted amplitude being $\sim$ 60 while the cycle actually peaked at 105 (SSN v1 values). However, it is interesting to note that Kimura correctly predicted the long term strengthening of solar activity during the first half of the 20th century! LS frequency analysis on sunspot data was also performed by Lomb and Andersen (1980), with similar results for the spectrum.

(2) Fourier analysis is probably the most commonly used method of spectral decomposition in science. It has been applied to sunspot data from the beginning of the 20th century (Turner, 1913c, b; Michelson, 1913). Vitinsky (1973) judges Fourier-based forecasts even less reliable than LS periodogram methods. Indeed, for instance Cole (1973) predicted cycle 21 to have a peak amplitude of 60, while the real value proved to be nearly twice that.

(3) The maximum entropy method (MEM) relies on the Wiener–Khinchin theorem that the power spectrum is the Fourier transform of the autocorrelation function. Calculating the autocorrelation of a time series for $M\ll N$ points and extrapolating it further in time in a particular way to ensure maximal entropy can yield a spectrum that extends to arbitrarily low frequencies despite the shortness of the data segment considered, and also has the property of being able to reproduce sharp spectral features (if such are present in the data in the first place). A good description of the method is given by Ables (1974), accompanied with some propaganda for it — see Press et al (1992) for a more balanced account of its pros and cons. The use of MEM for sunspot number prediction was pioneered by Currie (1973). Using maximum entropy method combined with multiple regression analysis (MRA) to estimate the amplitudes and phases, Kane (2007) arrived at a prediction of 80 to 101 for the maximum amplitude of cycle 24 (v1 values). It should be noted that the same method yielded a prediction (Kane, 1999) for cycle 23 that was far off the mark.

(4) Singular spectrum analysis (SSA) is a relatively novel method for the orthogonal decomposition of a time series. While in the methods discussed above the base was fixed (the trigonometric functions), SSA allows for the identification of a set of othogonal eigenfunctions that are most suitable for the problem. This is done by a principal component analysis of the covariance matrix $r_{ik}=\langle R_{i}R_{i+k}\rangle$. SSA was first applied to the sunspot record by Rangarajan (1998) who only used this method for pre-filtering before the application of MEM. Loskutov et al (2001) who also give a good description of the method, already made a prediction for cycle 24: a peak amplitude of 117 (v1 value).The forecast was later corrected slightly downwards to 106 (Kuzanyan et al, 2008).

The dismal performance of spectral predictions with the methods (1) – (3) indicates that the sunpot number series cannot be well represented by the superposition of a limited number of fixed periodic components. Instead,

• •

the periods may be time dependent,

• •

the system may be quasiperiodic, with a significant finite width of the periodic peaks (esp. the 11-year peak),

• •

there may be non-periodic (i.e., chaotic or stochastic) components in the behaviour of the system, manifested as a continuous background in the spectrum.

In practice, all three effects suggested above may play some part. The first mentioned effect, time dependence, may in fact be studied within the framework of spectral analysis. MEM and SSA are intrinsically capable of detecting or representing time dependence in the spectrum, while LS and Fourier analysis can study time dependence by sliding an appropriate data window across the period covered by observations. If the window is Gaussian with a width proportional to the frequency we arrive at the popular wavelet analysis. This method was applied to the sunspot number series by Ochadlick et al (1993), Vigouroux and Delachie (1994), Frick et al (1997), Fligge et al (1999), and Li et al (2005) who could confirm the existence and slight variation of the 11-year cycle and the Gleissberg-cycle. Recently, Kolláth and Oláh (2009) called attention to a variety of other generalized time dependent spectral analysis methods, of which the pseudo-Wigner transform yields especially clear details (see Fig. 15). The time varying character of the basic periods makes it difficult to use these results for prediction purposes but they are able to shed some light on the variation as well as the presistent or intermittent nature of the periods determining solar activity.

In summary, it is fair to say that forecasts based on harmonic analysis are notoriously unreliable. The secular variation of the basic periods, obeying as yet unknown rules, would render harmonic analysis practically useless for the prediction of solar cycles even if solar activity could indeed be described by a superposition of periodic functions. Although they may be potentially useful for very long term prediction (on centennial scales)191919But even here care is needed not to read more into the data than they contain, as discussed by Usoskin (2018), when it comes to cycle-to-cycle forecasts the best we can hope from spectral studies is apparently an indirect contribution, by constraining dynamo models with the inambiguously detected periodicities.

In what remains from this subsection, we briefly review what these apparently physically real periods are and what impact they may have on solar cycle prediction.

4.2.1 The 11-year cycle and its harmonics

As an example of the period spectrum obtained by these methods, in Fig. 14 we present the FFT based power spectrum estimate of the smoothed sunspot number record. Three main features are immediately noticed:

• •

The dominant 11-year peak, with its sidelobes and its 5.5-year harmonic.

• •

The 22-year subharmonic, representing the even–odd rule.

• •

The significant power present at periods longer than 50 years, associated with the Gleissberg cycle.

The dominant peak in the power spectrum is at $\sim$ 11 years. Significant power is also present at the first harmonic of this period, at 5.5 years. This is hardly surprising as the sunspot number cycles, as presented in Fig. 3, have a markedly asymmetrical profile. It is a characteristic of Fourier decomposition that in any periodic series of cycles where the profiles of individual cycles are non-sinusoidal, all harmonics of the base period will appear in the spectrum.

Indeed, were it not for the 13-month smoothing, higher harmonics could also be expected to appear in the power spectrum. It has been proposed (Krivova and Solanki, 2002) that these harmonics are detected in the sunspot record and that they may be related to the periodicities of $\sim$ 1.3 years intermittently observed in solar wind speed (Richardson et al, 1994; Paularena et al, 1995; Szabo et al, 1995; Mursula and Zieger, 2000; Lockwood, 2001) and in the internal rotation velocity of the Sun (Howe, 2009, Sect. 10.1). An analoguous intermittent 2.5 year variation in the solar neutrino flux (Shirai, 2004) may also belong to this group of phenomena. It may be worth noting that, from the other end of the period spectrum, the 154-day Rieger period in solar flare occurrence (Rieger et al, 1984; Bai and Cliver, 1990) has also been tentatively linked to the 1.3-year periodicity. Unusually strong excitation of such high harmonics of the Schwabe cycle may possibly be explained by excitation due to unstable Rossby waves in the tachocline (Zaqarashvili et al, 2010).

The 11-year peak in the power spectrum has substantial width, related to the rather wide variation in cycle lengths in the range 9 – 13 years. Yet Fig. 14 seems to suggest the presence of a well detached second peak in the spectrum at a period of $\sim$ 14 years. The presence of a distinct peak at the first harmonic and even at the subharmonic of this period seems to support its reality. Indeed, peaks at around 14 and 7 years were already found by other researchers (e.g., Kimura, 1913; Currie, 1973) who suggested that these may be real secondary periods of sunspot activity.

The situation is, however, more prosaic. Constraining the time interval considered to data more recent than 1850, from which time the sunspot number series is considered to be more reliable, the 14.5-year secondary peak and its harmonics completely disappear. On the other hand, the power spectrum for the years 1783 – 1835 indicates that the appearance of the 14.5-year secondary peak in the complete series is almost entirely due to the strong predominance of this period (and its harmonic) in that interval. This interval covers the unusually long cycle 4 and the Dalton minimum, consisting of three consecutive unusually weak cycles, when the “normal” 11-year mode of operation was completely suppressed.

As pointed out by Petrovay (2010), this probably does not imply that the Sun was operating in a different mode during the Dalton minimum, the cycle length being 14.5 years instead of the usual 11 years. Instead, the effect may be explained by the well known inverse correlation between cycle length and amplitude, which in turn is the consequence of the strong inverse correlation between rise rate and cycle amplitude (Waldmeier effect), combined with a much weaker or nonexistent correlation between decay rate and amplitude (see Sect. 1.3.3). The cycles around the Dalton minimum, then, seem to lie at the low amplitude (or long period) end of a continuum representing the well known cycle length–amplitude relation, ultimately explained by the Waldmeier effect.

A major consequence of this is that the detailed distribution of peaks varies significantly depending on the interval of time considered. Indeed, Kolláth and Oláh (2009) recently applied time dependent harmonic analysis to the sunspot number series and found that the dominant periods have shown systematic secular changes during the past 300 years (Fig. 15). For instance, the basic period seems to have shortened from 11 years to 10 years between 1850 and 1950, with some moderate increase in the last 50 years. (This is consistent with the known anticorrelation between cycle length and amplitude, cf. Sect. 1.3.3.)

4.2.2 The even–odd (a.k.a. Gnevyshev–Ohl) rule

A cursory look at Fig. 3 shows that solar cycles often follow an alternating pattern of higher and lower maxima. In this apparent pattern, already noticed by the early observers (e.g, Turner, 1913a), odd cycles have been typically stronger than even cycles in the last two centuries.

This even–odd rule can be given two interpretations: a “weak” one of a general tendency of alternation between even and odd cycles in amplitude, or a “strong” one of a specific numerical relation between the amplitudes of consecutive cycles.

Let us first consider the rule in its weak interpretation. At first sight the rule admits many exceptions, but the amplitude of solar cycles depends on the particular measuring method used. Exceptions from the even–odd alternation rule become less common if a long term trend (calculated by applying a 12221 or 121 filter, see Sect. 1.3.1) is subtracted from the data (Charbonneau, 2001), or if integrated cycle amplitudes (sums of annual mean sunspot numbers during the cycle) are used (Gnevyshev and Ohl, 1948).

In fact, as evident from, e.g., the work of Mursula et al (2001) where cycle amplitudes are based on group sunspot numbers and the amplitude of a cycle is defined as the sum of the annual GSN value over the course of the cycle, the odd–even alternation may be considered as strictly valid with only four exceptions:

• •

In the pairs 7 – 8 and 17 – 18, odd cycles are followed by stronger even cycles at the end of Dalton minimum and at the beginning of the Modern Maximum. These exceptions could be made to disappear by the subtraction of the long term trend as suggested by Charbonneau (2001).

• •

The pair 22 – 23 represents another apparent break of the weak even–odd rule which is not easily explained away, even though the relative difference is smaller if the Kislovodsk sunspot number series is used (Nagovitsyn et al, 2009). The possibility is obviously there that the subtraction of the long term trend may resolve the problem but we have no way to tell in the near future.

• •

Prior to cycle 5, the phase of the alternation was opposite, even cycles being stronger than odd cycles. As cycle 4 is known to have been anomalously long anyway (the so-called “phase catastrophe” in the solar cycle, Vitinsky et al 1986) and its decaying phase is not well covered by observations (Vaquero, 2007), this gave rise to the suggestion of a “lost solar cycle” between cycles 4 and 5 (Usoskin et al 2001). This cycle, however, would have been even more anomalous than cycle 4 and despite intensive searches in historic data the evidence is still not quite conclusive (Krivova et al 2002; see, however, Usoskin et al 2009a).

The issue whether the even–odd rule can go through phase jumps or not is important with respect to its possible origin. One plausible possibility is that the alternation is due to the superposition of a steady primordial magnetic field component on the oscillatory magnetic field generated by the dynamo (Levy and Boyer, 1982). In this case, any phase jump in the Gnevyshev–Ohl rule should imply a phase jump in Hale’s polarity rules, too. Alternatively, persistent even–odd alternation may also arise in nonlinear dynamos as a period–2 limit cycle (Durney, 2000); with a stochastic forcing occasional phase jumps are also possible (Charbonneau, 2001; Charbonneau et al, 2007).

While we have no information on this from the 18th century phase jump, we can be certain that there was no such phase jump in polarities in the last two decades, even though the even–odd rule seems to have been broken again. It will be interesting to see when (and if) the even–odd rule settles in again, whether it will have done so with a phase jump or not. For instance, if cycle 25 will again exceed cycle 24 it would seem that no phase jump occurred and both theoretical options are still open. But if cycle 25 will represent a further weakening from cycle 24, followed by a stronger cycle 26, a phase jump will have occurred, which may exclude the primordial field origin of the rule if Hale’s polarity rules remain unchanged.

Let us now discuss the stronger interpretation of the even–odd rule. In the first quantitative study of the relative amplitudes of consecutive cycles, Gnevyshev and Ohl (1948) found a rather tight correlation between the time integrated amplitudes of even and subsequent odd cycles, while the correlation between odd cycles and subsequent even cycles was found to be much less strong. This gave rise to the notion that solar cycles come in “two-packs” as even–odd pairs. Nagovitsyn et al (2009) confirmed this puzzling finding on the basis of data covering the whole period of telescopic observations (and renumbering cycles before 1790 in accordance with the lost cycle hypothesis); they also argue that cycle pair 22 – 23 does not deviate strongly from the even–odd correlation curve so it should not be considered a “real” exception to the even–odd rule. Javaraiah (2012) analyzed the validity of the rule considering large and small sunspot groups separately, and found that while for large groups the rule holds with a few exceptions, for small groups a ‘reverse G-O rule’ holds where odd numbered cycles are consistently stronger than the preceding, rather than the following even numbered cycle.

Shortly after its formulation by Gnevyshev and Ohl (1948), the (strong) even–odd rule was used by Kopecký (1950) to successfully predict the unusually strong cycle 19. This made this rule quite popular for forecast purposes. However, forecasts based on the even–odd rule completely failed for cycle 23, overpredicting the amplitude by

 50% (see review by Li et al, 2001). Taken together with the implausibility of the suggested two-pack system, this shows that it is probably wiser to take the position that “extraordinary claims need extraordinary evidence” — which is yet to be provided in the case of the “strong” even–odd rule.

Finally, in the context of the even–odd rule, it is also worth mentioning the three-cycle regularity proposed by Ahluwalia (1998). Even though the evidence presented for the alleged triadic pattern is not overwhelming, this method resulted in one of the few successful predictions for the amplitude of cycle 23.

4.2.3 The Gleissberg cycle

Besides the changes in the length of the 11-year cycle related to the amplitude–cycle length correlation, even more significant are the variations in the period of the so-called Gleissberg cycle (Gleissberg, 1939). This “cycle”, corresponding to the 60 – 120 year “plateau” in Fig. 14 was actually first noticed by Wolf, who placed it in the range 55 – 80 years (see Richard, 2004, for a discussion of the history of the studies of the Gleissberg cycle). Researchers in the middle of the 20th century characterized it as an 80 – 100 year variation. Figure 15 explains why so widely differing periods were found in different studies: the period has in fact shown a secular increase in the past 300 years, from about 50 years in the early 18th century, to a current value exceeding 140 years. This increased length of the Gleissberg cycle also agrees with the results of Forgács-Dajka and Borkovits (2007).

The detection of $\sim$ 100 year periods in a data set of 300 years is of course always questionable, especially if the period is even claimed to be varying. However, the very clear and, most importantly, nearly linear secular trend seen in Fig. 15 argues convincingly for the reality of the period in question. This clear appearance of the period is due to the carefully optimized choice of the kernel function in the time–frequency analysis, a method resulting in a so-called pseudo-Wigner distribution (PWD). In addition, in their study Kolláth and Oláh (2009) presented a conscientious test of the reliability of their methods, effectively proving that the most salient features in their PWD are not artefacts. (The method was subsequently also applied to stellar activity, Oláh et al 2009.) This is the most compelling evidence for the reality of the Gleissberg cycle since 1750 yet presented. Further evidence was more recently presented by Le Mouël et al (2017) using singular spectrum analysis.

4.2.4 Supersecular cycles

For the 210-year Suess (known also as de Vries) cycle, McCracken and Beer (2008) presented further evidence for the temporally intermittent nature of this marked peak in the spectrum of solar proxies. The Suess cycle seems to have a role in regulating the recurrence rate of grand minima. Grand minima, in turn, only seem to occur during < 1 kiloyear intervals (“Spörer events”) around the minima of the $\sim$ 2400-year Hallstatt cycle.

For further discussion of long term variations in solar activity we refer the reader to the reviews by Beer et al (2006) and Usoskin (2017).

4.3 Nonlinear methods

“…every complicated question has a simple answer which is wrong. Analyzing a time series with a nonlinear approach is definitely a complicated problem. Simple answers have been repeatedly offered in the literature, quoting numerical values for attractor dimensions for any conceivable system.”

(Hegger et al, 1999)

The nonlinearities in the dynamo equations readily give rise to chaotic behaviour of the solutions. The long term behaviour of solar activity, with phenomena like grand minima and grand maxima, is also suggestive of a chaotic system. While chaotic systems are inherently unpredictable on long enough time scales, their deterministic nature does admit forecast within a limited range. It is therefore natural to explore this possibility from the point of view of solar cycle prediction.

4.3.1 Attractor analysis and phase space reconstruction: the

pros …

Assuming that the previous $(M-1)$ values of the sunspot number do in some way determine the current expected value, our problem becomes restricted to an $M$-dimensional phase space, the dimensions being the current value and the $(M-1)$ previous values. With a time series of length $N$, we have $N-M+1$ points fixed in the phase space, consecutive points being connected by a line. This phase space trajectory is a sampling of the attractor of the physical system underlying the solar cycle (with some random noise added to it). The attractor represents a mapping in phase space which maps each point into the one the system occupies in the next time step: if this mapping is known to a good degree of precision, it can be used to extend the trajectory towards the future.

For the mapping to be known, $M$ needs to be high enough to avoid self-crossings in the phase space trajectory (otherwise the mapping is not unique) but low enough so that the trajectory still yields a good sampling of the attractor. The lowest integer dimension satisfying these conditions is the embedding dimension $D$ of the attractor (which may have a fractal dimension itself).

Once the attractor has been identified, its mathematical description may be done in three ways.

(1) Parametric fitting of the attractor mapping in phase space. The simplest method is the piecewise linear fit suggested by Farmer and Sidorowich (1987) and applied in several solar prediction attempts, e.g., Kurths and Ruzmaikin (1990). Using a method belonging to this group, Kilcik et al (2009) gave a correct prediction for Cycle 24. Alternatively, a global nonlinear fit can also be used: this is the method applied by Serre and Nesme-Ribes (2000) as the first step in their global flow reconstruction (GFR) approach.

(2) Nonparametric fitting. The simplest nonparametric fit is to find the closest known attractor point to ours (in the $(M-1)$-dimensional subspace excluding the last value) and then using this for a prediction, as done by Jensen (1993). (This resulted in so large random forecast errors that it is practically unsuitable for prediction.) A more refined approach is simplex projection analysis, recently applied by Singh and Bhargawa (2017) for the problem of solar cycle prediction. (See also Sarp et al 2018.) A most remarkable extension of these methods was presented by Covas (2017) who, instead of focusing on the time series of SSN only, considered the problem of extending the whole spatiotemporal data set of sunspot positions (butterfly diagram) into the future.

Neural networks, discussed in more detail in Sect. 4.3.4 below, are a much more sophisticated nonparametric fitting device.

(3) Indirectly, one may try to find a set of differential equations describing a system that gives rise to an attractor with properties similar to the observed. In this case there is no guarantee that the derived equations will be unique, as an alternative, completely different set may also give rise to a very similar attractor. This arbitrariness of the choice is not necessarily a problem from the point of view of prediction as it is only the mapping (the attractor structure) that matters. Such phase space reconstruction by a set of governing equations was performed, e.g., by Serre and Nesme-Ribes (2000) or Aguirre et al (2008). On the other hand, instead of putting up with any arbitrary set of equations correctly reproducing the phase space, one might make an effort to find a set with a structure reasonably similar to the dynamo equations so they can be given a meaningful physical interpretation. Methods following this latter approach were discussed in Sects. 3.5 and 3.6.

4.3.2 …the cons …

Finding the embedding dimension and the attractor structure is not a trivial task, as shown by the widely diverging results different researchers arrived at. One way to find the correct embedding dimension is the false nearest neighbours method (Kennel et al, 1992), essentially designed to identify self-crossings in the phase space trajectory, in which case the dimension $M$ needs to be increased. But self-crossings are to some extent inevitable, due to the stochastic component superimposed on the deterministic skeleton of the system.

As a result, the determination of the minimal necessary embedding dimension is usually done indirectly. One indirect method fairly popular in the solar community is the approach proposed by Sugihara and May (1990) where the correct dimension is basically figured out on the basis of how successfully the model, fit to the first part of the data set, can “predict” the second part (using a piecewise linear mapping).

Another widely used approach, due to Grassberger and Procaccia (1983), starts by determining the correlation dimension of the attractor, by simply counting how the number of neighbours in an embedding space of dimension $M\gg 1$ increases with the distance from a point. If the attractor is a lower dimensional manifold in the embedding space and it is sufficiently densely sampled by our data then the logarithmic steepness $d$ of this function should be constant over a considerable stretch of the curve: this is the correlation dimension $d$. Now, we can increase $M$ gradually and see at what value $d$ saturates: that value determines the attractor dimension, while the value of $M$ where saturation is reached yields the embedding dimension.

The first nonlinear time series studies of solar activity indicators suggested a time series spacing of 2 – 5 years, an attractor dimension $\sim$ 2 – 3 and an embedding dimension of 3 – 4 (Kurths and Ruzmaikin, 1990; Gizzatullina et al, 1990). Other researchers, however, were unable to confirm these results, either reporting very different values or not finding any evidence for a low dimensional attractor at all (Calvo et al, 1995; Price et al, 1992; Carbonell et al, 1994; Kilcik et al, 2009; Hanslmeier and Brajša, 2010). In particular, I would like to call attention to the paper by Jensen (1993), which, according to ADS and WoS, has received a grand total of zero citations up to 2010, yet it displays an exemplary no-nonsense approach to the problem of sunspot number prediction by nonlinear time series methods. Unlike so many other researchers, the author of that paper does not enforce a linear fit on the logarithmic correlation integral curve (his Figure 4); instead, he demonstrates on a simple example that the actual curve can be perfectly well reproduced by a simple stochastic process.

These contradictory results obviously do not imply that the mechanism generating solar activity is not chaotic. For a reliable determination a long time series is desirable to ensure a sufficiently large number of neighbours in a phase space volume small enough compared to the global scale of the attractor. Solar data sets (even the cosmogenic radionuclide proxies extending over millennia but providing only a decadal sampling) are typically too short and sparse for this. In addition, clearly distinguishing between the phase space fingerprints of chaotic and stochastic processes is an unsolved problem of nonlinear dynamics which is not unique to solar physics. A number of methods have been suggested to identify chaos unambiguously in a time series but none of them has been generally accepted and this topic is currently a subject of ongoing research — see, e.g., the work of Freitas et al (2009) which demonstrates that the method of “noise titration”, somewhat akin to the Sugihara–May algorithm, is uncapable of distinguishing superimposed coloured noise from intrinsically chaotic systems.

4.3.3 …and the upshot

Starting from the 1980s, many researchers jumped on the chaos bandwagon, applying nonlinear time series methods designed for the study of chaotic systems to a wide variety of empirical data, including solar activity parameters. From the 1990s, however, especially after the publication of the influential book by Kantz and Schreiber (1997), it was increasingly realized that the applicability of these nonlinear algorithms does not in itself prove the predominantly chaotic nature of the system considered. In particular, stochastic noise superposed on a simple, regular, deterministic skeleton can also give rise to phase space characteristics that are hard to tell from low dimensional chaos, especially if strong smoothing is applied to the data. As a result, the pendulum has swung in the opposite direction and currently the prevailing view is that there is no clear cut evidence for chaos in solar activity data (Panchev and Tsekov, 2007).

One might take the position that any forecast based on attractor analysis is only as good as the underlying assumption of a chaotic system is: if that assumption is unverifiable from the data, prediction attempts are pointless. This, however, is probably a too hasty judgment. The potentially most useful product of phase space reconstruction attempts is the inferences they allow regarding the nature of the underlying physical system (chaotic or not), even offering a chance to constrain the form of the dynamo equations relevant for the Sun. As discussed in the previous section, such truncated models may be used for forecast directly, or alternatively, the insight they yield into the mechanisms of the dynamo may be used to construct more sophisticated dynamo models.

4.3.4 Neural networks

Neural networks are algorithms built up from a large number of small interconnected units (“neurons” or “threshold logic units”), each of which is only capable of performing a simple nonlinear operation on an input signal, essentially described by a step function or its generalized (rounded) version, a sigmoid function. To identify the optimal values of thresholds and weights parameterizing the sigmoid functions of each neuron, an algorithm called “back propagation rule” is employed which minimizes (with or without human guidance) the error between the predicted and observed values in a process called “training” of the network. Once the network has been correctly trained, it is capable of further predictions.

The point is that any arbitrary multidimensional nonlinear mapping may be approximated by a combination of stepfunctions to a good degree — so, as mentioned in Sect. 4.3.1 above, the neural network can be used to find the nonlinear mapping corresponding to the attractor of the given time series.

More detailed introductions to the method are given by Blais and Mertz (2001), Conway (1998), and by Calvo et al (1995); the latter authors were also the first to apply a neural network for sunspot number prediction. Unfortunately, despite their claim of being able to “predict” (i.e., postdict) some earlier cycles correctly, their prediction for Cycle 23 was off by a wide margin (predicted peak amplitude of 166 [v1] as opposed to 121 observed). One of the neural network forecasts for Cycle 24 (Maris and Oncica, 2006) was equally far off, while another one (Uwamahoro et al, 2009) yielded a more conservative value. A prediction for Cycle 25 based on a version of the neural networks approach was given by Attia et al (2013), resulting in a cycle amplitude slightly below that of cycle 24.

5 Summary evaluation

The performance of various forecast methods in Cycles 21 – 23 was discussed by Li et al (2001) and Kane (2001). Predictions for Cycle 24 were presented in Petrovay (2010) (Table 1), Pesnell (2008) and Pesnell (2012); the experiences gained in this cycle were discussed in Pesnell (2016).

Precursor methods generally stand out with their internally consistent forecasts which for Cycles 21 and 22 proved to be correct. For Cycle 23, these methods were still internally consistent in their prediction, mostly scattering in a narrow range between202020Version 1 sunspot numbers are used throughout this section, unless otherwise indicated. 150 and 170; however, the cycle amplitude proved to be considerably lower ($R_{\mathrm{max}}=121$). It should be noted, however, that one precursor based prediction, that of Schatten et al (1996) was significantly lower than the rest ($138\pm 30$) and within $0.6\,\sigma$ of the actual value. For Cycle 24 most precursor methods again consistently indicated a lower-than-average cycle amplitude in the range 70 – 100, except Feynman’s geomagnetic precursor method which mistakenly resulted in a very high value of 150. (The likely reasons were discussed in Sect. 2.4 above.) The closest hit at the actual peak value of 67 [v1] or 116 [v2] was produced by the Minimax3 and the polar field precursor methods (Petrovay 2010; Svalgaard et al 2005).

Indeed, the polar precursor method of Schatten and Sofia (1987) and Schatten et al (1996), has consistently proven its skill in all cycles. As discussed in Sect. 2.3, this method is essentially based on the polar magnetic field strength as precursor.

Model based methods are a new development that have only had limited occasion to prove their skill. For Cycle 24 only three conceptually different such predictions were made, all of which were based on dynamo models. The pioneering attempt by Dikpati and Gilman (2006) proved to be way too high (see Sect. 3.4.1 for a discussion of the possible reasons). The flux transport dynamo based predictions of Choudhuri et al (2007) and Jiang et al (2007) were close hits; however, as already mentioned, these employed a technique (adjusting the dipole moment at the minimum to observations) which renders them essentially a polar field precursor method in disguise. Another correct model-based forecast was given by Kitiashvili and Kosovichev (2008); the good performance of this dynamo model, seemingly rather far removed from physical reality, still needs to be understood and it may possibly be equivalent to a phase space reconstruction method, as in item (3) of Sect. 4.3.1.

Extrapolation methods as a whole have shown a much less impressive performance. Overall, the statistical distribution of maximum amplitude values predicted by “real” forecasts made using these methods (i.e., forecasts made at or before the minimum epoch) for any given cycle does not seem to significantly differ from the long term climatological average of the solar cycle quoted in Sect. 1.3 above. It would of course be a hasty judgement to dismiss each of the widely differing individual approaches comprised in this class simply due to the poor overall performance of the group. In particular, some novel methods introduced in the last decades, such as SSA, phase space reconstruction or neural networks have hardly had a chance to debut, so their further performance will be worth monitoring in upcoming cycles.

The effect of the sunspot number revision on solar cycle prediction methods is limited to numerical corrections of minor importance.

6 Forecasts for Cycle 25 [NEW!]

Table 2 presents a collection of forecasts for the amplitude of Cycle 25, without claiming completeness. The objective was to include one or two representative forecasts from each category. As the time of the minimum starting Cycle 25 is yet to be established, all these forecast qualify as “early”, in the sense that the most well-established methods, relying on precursor values evaluated at the time of minimum, cannot yet be applied.

It is clear also from this table that the issue of cycle prediction is less contentious for Cycle 25 than it was for Cycle 24. The overwhelming majority of forecasts agree that the amplitude of Cycle 25 is most likely to lie within $\pm 20\%$ of Cycle 24, i.e., no major change in the level of solar activity is expected. The remaining controversy mostly concerns where in this range the cycle will peak. Dynamo based predictions indicate that Cycle 25 will peak at somewhat lower values than Cycle 24, while precursor techniques and SFT modelling suggest a cycle amplitude comparable to or slightly higher than the previous cycle. Two recent neural network based forecasts yield a weak cycle peaking quite early.

Following the development of the sunspot number during the next few years will be most interesting in the light of these predictions, and it may shed more light on the strong and weak points in our understanding of the roots of solar activity variations.

Acknowledgements.

This work was supported by the Hungarian National Research, Development and Innovation Fund (grant no. NKFI K-128384) and by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 739500. Some figures (or their input data) have been provided by J. Jiang, A. Kilcik and P. Charbonneau, as noted in the captions. Source for International Sunspot Numbers: WDC-SILSO, Royal Observatory of Belgium, Brussels. Wilcox Solar Observatory data used in this study was obtained via the web site wso.stanford.edu, courtesy of J. T. Hoeksema. The Wilcox Solar Observatory is currently supported by NASA.