Bart: Modeling the historical sunspot record from planetary periods

Talkshop contributor ‘Bart’ has used a powerful mode of analysis of the sunspot record to reveal periodicities and also devised an elegant algorithm which can produce a time series very similar to the historical sunspot record. Bart does not explicitly link the two periods found to planetary periods as I have in the title of this thread to make it easy for Google to find us, but the two periods are very close to twice the orbital period of Jupiter and the synodic period of Jupiter and Saturn, which when averaged are close to the Hale cycle period. I got Bart’s permission to make an article, and Bart told me in comments here, that adjusting the two values to the exact planetary periods would make little difference to the result.

Modeling the historical sunspot record

I observed four significant peaks in the PSD of the SSN process at 10 years, 11.8 years, 10.8 years, and 131 years (unfortunately, the one at 131 years is obscured in the PSD I showed because I had to make the usual tradeoff between bias (resolution) and variance in my plot, and I chose to present a smoother version, the better to resolve the higher frequency peaks).

I made the assumption that the Sun spot number is a proxy for the magnitude of the process which is occurring. I redid the estimate using the magnitude squared because the autocorrelation function for the square of a Gaussian process is well known, as I related here:

Preliminary PSD analysis informs me that the solar cycle is governed by two quasi-periodic processes with periods of roughly T1 = 20 and T2 = 23.6 years.

The physical basis of the solar cycle was elucidated in the early twentieth century by George Ellery Hale and collaborators, who in 1908 showed that sunspots were strongly magnetized (this was the first detection of magnetic fields outside the Earth), and in 1919 went on to show that the magnetic polarity of sunspot pairs:

“Hale’s observations revealed that the solar cycle is a magnetic cyclewith an average duration of 22 years. However, because very nearly all manifestations of the solar cycle are insensitive to magnetic polarity, it remains common usage to speak of the “11-year solar cycle”.’

Is always the same in a given solar hemisphere throughout a given sunspot cycle;
Is opposite across hemispheres throughout a cycle;
Reverses itself in both hemispheres from one sunspot cycle to the next.

Figure 1

As per Papoulis (for me, 2nd edition, page 233), if a Gaussian process with autocorrelation function R(t) is squared, the resulting autocorrelation function, which I will call Q(t), is

Q(t) = R(0)^2 + 2*R(t)^2

The PSD of this is the Fourier transform of a constant R(0)^2 plus the scaled convolution in the frequency domain of the Fourier transform of R(t) (PSD of the underlying process) with itself (i.e., multiplication in the time domain results in convolution in the frequency domain, and vice versa).

The sunspot count appears to reflect the energy of these combined processes at around 20 and 23.6 years, which necessarily has apparent periods of 0.5*T1, 0.5*T2, T1*T2/(T2+T1), and T1*T2/(T2-T1) years, or 10 years, 11.8 years, 10.8 years, and 131 years. This latter appears as a quasi-beat period in the data. I say “quasi-” because these are not rigidly defined periods of steady state sinusoids, but mean periods of random excitation of resonance phenomena.

It is pretty clear that the model for the process governing Sun spot occurrence is the correct one, even if the parameterization is somewhat statistically uncertain (and even if some parameters may be randomly or deterministically varying slowly and/or narrowly in time, as well as the precise frequency distribution of noise energy, though we really only care about that within a narrow band around the resonances).

The processes can be characterised by the following simple model:

Figure 2

I started with zero initial conditions, ran it for the time period indicated, and re-looped around a couple more times to make sure initial transients were settled out to get a stationary process. I implemented it in MATLAB Simulinkusing the continuous time simulation blocks exactly as in my diagram. For the random inputs, I used the Band Limited White Noise with a bandwidth of 5 years^-1 (sample time of 0.1 years). Doing this just happened to create the results shown in a couple of runs.

Here is an example of the output generated:

Figure 3

If I had just four random periods, I could not generally fit them in this way. I have a mapping of a 2d plane to 4d space, and the appropriate 4 parameters just happen to lie on the projected 2d submanifold described by the set of points (T1/2,T2/2,T2*T1/(T2+T1),T2*T1/(T2-T1)) which has measure zero in that 4d space.

Measure zero means the probability of the line up of these frequency lines by chance is ZERO. There is no doubt about it. Zero. Zilch. None.

I could devise a model for the amplitude of a sinusoidal process “z” driven by noise “e” as

xdot = -y
ydot = w^2 * x + e
z = sqrt((w*x)^2+y^2)

You could also model it directly as a phase “phi” and z using a transformation of variables so as to obtain a pair of differential equations

phidot = f1(z,phi,e)
zdot = f2(z,phi,e)

for appropriate functions f1() and f2(). The periodicity of z would be twice that of x and y, so you would tell me the frequency is 2w and I would say it was w. Both models are valid. Both are right. But, the linear model is much easier to work with.

All these numbers are preliminary, rough estimates, mind you. There are many other steps which would be needed to nail them down precisely.  Once completed  a Kalman Filter run backwards and forwards on the data set will properly initialize the states, and a projection forward with error bounds from the propagated covariance can easily be generated. (but, it takes time, and I have a job and a life).

Postscript:

Now, an interesting aside I noticed on another thread is:
Here’s another possibility I happened on the other day, though. Let’s suppose that the Sun’s instantaneous average temperature is dominated by the sum of a bias and two sinusoidally varying processes

Tsun = A + B*cos(w1*t + phi1) + C*cos(w2*t + phi2)

where w1 = 2*pi/T1 is the radial frequency of a process with period T1, and similarly for w2. A, B, and C are constants of proportionality. The phases phi1 and phi2 determine the relative phase of the cyclical processes.

Let us assume that Sun Spot numbers are proportional to the magnitude of the cyclical part:

Sspot ~ | B*cos(w1*t + phi1) + C*cos(w2*t + phi2) |

The square of this process is

Sspot^2 = (B*cos(w1*t + phi1))^2 + (C*cos(w2*t + phi2))^2 +
2*B*C*cos(w1*t + phi1)*cos(w2*t + phi2)

or, using trig identities

Sspot^2 = (B^2+C^2)/2 + (B^2/2)*cos(2*(w1*t+phi1)) + (C^2/2)*cos(2*(w2*t+phi2)) +
B*C*( cos((w1+w2)*t+phi1+phi2) + cos((w1-w2)*t+phi1-phi2) )

Thus, we have frequencies of 2*w1, 2*w2, w1+w2, and w1-w2, which is to say, periods of T1/2, T2/2, T1*T2/(T2+T1), T1*T2/(T2-T1). Spectral analysis shows that this is essentially what we have, with T1 = ~20 years, and T2 = ~23.6 years. This begets periods in the Sun Spot data of T1/2 = 10 years, T2/2 = 11.8 years, T1*T2/(T2+T1) = 10.8 years, and T1*T2/(T2-T1) = 131 years.

Now, radiative heat transfer goes as temperature to the fourth power, so Earth temperatures should be modulated at all the combinations above taken one step further, and including the bias offset of temperature not reflected in the Sun Spot data. So, we have frequencies of 4*w1, 4*w2, 2*(w1+w2), 2*(w1-w2), 3*w1+w2, w1-w2, 3*w1-w2, w1+w2, w1+3*w2, -w1+3*w2, 2*w1, 2*w2, all modulo sign. The periods are then T1/4, T2/4, 0.5*T1*T2/(T2+T1), 0.5*T1*T2/(T2-T1), T1*T2/(3*T2+T1), T1*T2/(T2-T1), T1*T2/(3*T2-T1), T1*T2/(T2+T1), T1*T2/(T2+3*T1), T1*T2/(3*T1-T2), T1/2, T2/2, T1, and T2. These work out to 5, 5.9, 5.4, 65.6, 5.2, 131, 9.3, 10.8, 5.6, 13, 10, 11.8, 20, and 23.6 years. Some of these may cancel each other out depending on the phases and amplitudes. And, some may be attenuated or amplified by the Earth’s actual thermal response, particularly the shorter period (higher frequency) ones ought to tend to be attenuated.

Most of these periods appear to describe the locations of major and minor peaks in the 20th century HADCRUT data within the resolution of the PSD. In particular, there is the ~65 year one, and what may be an average of the 20 and 23.6 year ones, which would resolve into two peaks with more data.
Thus, it is possible that the entire shooting match is the result of just two guns at 20 and 23.6 years, but more extensive analysis is needed.

========================================================================

How about it solar physics institutions? Don’t you think it would be worth negotiating a contract with Bart to continue this promising work so we can get a potentially useful prediction rather than the Hathaway ‘prediction’ which changes monthly?