# Can irregularities of solar proxies help understand quasi-biennial solar variations?

A. Shapoval, J. L. Le Mouël, M. Shnirman, V. Courtillot

► To cite this version:

A. Shapoval, J. L. Le Mouël, M. Shnirman, V. Courtillot. Can irregularities of solar proxies help understand quasi-biennial solar variations?. Nonlinear Processes in Geophysics, 2014, 21, pp.797-813. 10.5194/npg-21-797-2014. insu-03581114

HAL Id: insu-03581114

https://insu.hal.science/insu-03581114

Submitted on 19 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Can irregularities of solar proxies help understand quasi-biennial solar variations?

A. Shapoval<fnref n="1" />,<fnref n="3" />,<fnref n="4" />, J. L. Le Mouël<fnref n="2" />, M. Shnirman<fnref n="1" />,<fnref n="2" />, and V. Courtillot<fnref n="2" />

*Correspondence to:* V. Courtillot (courtil@ipgp.fr)

Received: 22 January 2014 – Published in Nonlin. Processes Geophys. Discuss.: 13 March 2014

Revised: 10 June 2014 – Accepted: 22 June 2014 – Published: 1 August 2014

**Abstract.** We define, calculate and analyze irregularity indices $\lambda_{ISSN}$ of daily series of the International Sunspot Number ISSN as a function of increasing smoothing from $N = 162$ to 648 days. The irregularity indices $\lambda$ are computed within 4-year sliding windows, with embedding dimensions $m = 1$ and $2$. $\lambda_{ISSN}$ displays Schwabe cycles with ~5.5-year variations ("half Schwabe variations" HSV). The mean of $\lambda_{ISSN}$ undergoes a downward step and the amplitude of its variations strongly decreases around 1930. We observe changes in the ratio $R$ of the mean amplitude of $\lambda$ peaks at solar cycle minima with respect to peaks at solar maxima as a function of date, embedding dimension and, importantly, smoothing parameter $N$. We identify two distinct regimes, called Q1 and Q2, defined mainly by the evolution of $R$ as a function of $N$: Q1, with increasing HSV behavior and $R$ value as $N$ is increased, occurs before 1915–1930; and Q2, with decreasing HSV behavior and $R$ value as $N$ is increased, occurs after ~1975. We attempt to account for these observations with an autoregressive (order 1) model with Poissonian noise and a mean modulated by two sine waves of periods $T_1$ and $T_2$ ($T_1$ = 11 years, and intermediate $T_2$ is tuned to mimic quasi-biennial oscillations QBO). The model can generate both Q1 and Q2 regimes. When $m = 1$, HSV appears in the absence of $T_2$ variations. When $m = 2$, Q1 occurs when $T_2$ variations are present, whereas Q2 occurs when $T_2$ variations are suppressed. We propose that the HSV behavior of the irregularity index of ISSN may be linked to the presence of strong QBO before 1915–1930, a transition and their disappearance around 1975, corresponding to a change in regime of solar activity.

## 1 Introduction

Regular and irregular features of solar activity reflect the behavior of the solar dynamo. Their spectrum contains low-frequency “cycles”, from decadal to centennial scales, whose durations and amplitudes vary with time, and a higher frequency spectrum with much stronger irregularities, notably in the 1–3-year pseudoperiod range. The case of quasi-biennial oscillations (QBO) has been widely discussed in the recent literature (e.g., McIntosh et al., 1992; Lawrence et al., 2008; Mursula et al., 2003; Rouillard and Lockwood, 2004). The range of 1–3-year quasi-periodicities has been studied in a number of time series, using different techniques such as power spectral analysis (Rouillard and Lockwood, 2004; Valdes-Galicia et al., 1996), wavelet analysis (Kudela et al., 2002; Mursula et al., 2003), empirical mode decomposition (Vecchio et al., 2010), or the successive approximation technique (Mavromichalaki et al., 2003). All techniques confirm the reality of these quasi-periodicities, with time-varying amplitude and “frequency”.

Several papers discuss variations with periods close to 27 days (related to the Sun’s rotation as seen from Earth). For instance, in an earlier paper (Le Mouël et al., 2007), we considered the series of the International Sunspot Number (ISSN; SIDC team, 2005) and magnetic *aa* index: we computed their energy for periods of around 27 days and found that this energy roughly followed the initial time series it was computed from. More detailed analysis revealed a significant increase of energy approximately two decades prior to the increase in solar activity that occurred in the 1930s. Other papers deal with the long-term evolution of short-term variations of

<footnote n="1">¹IEPT RAS, Profsoyuznaya str. 84/32, 117 997 Moscow, Russia</footnote>

<footnote n="2">²IPGP, 1 rue Jussieu, 75005 Paris, France</footnote>

<footnote n="3">³Financial University, Leningradsky pr. 49, 125 167 Moscow, Russia</footnote>

<footnote n="4">⁴National Research University Higher School of Economics, 20 Myasnitskaya Ulitsa, 101 000 Moscow, Russia</footnote>
different time series by standard wavelet analysis (Lawrence et al., 2008) or using some modification of Kolmogorov entropy (Blanter et al., 2005, 2006). These papers reveal the existence of different regimes in the long-term evolution of the high-frequency part of the spectrum (estimated locally in time).

In a previous paper (Shapoval et al., 2013), we introduced the irregularity index of a given time series as the convergence (or divergence) rate of nearby points in a certain phase space, under a “one-step” translation. In the case of low-dimensional dynamical systems, the irregularity index corresponds to the maximal Lyapunov exponent (e.g., Bergé et al., 1984). Lyapunov exponents characterize the convergence (resp. divergence) rate of infinitesimally close trajectories of a dynamical system to (resp. from) its attractor in phase space. There is a link between the magnitude of the Lyapunov exponent and the regularity of the process: the larger the exponent, the stronger the irregularities. In contrast to the maximal Lyapunov exponent, the irregularity index can be computed for shorter time series with a significant random component.

In Shapoval et al. (2013), we explored variations of the irregularity index $\lambda_m(t)$ of the daily ISSN series as a function of time for intermediate values (4–6) of the embedding dimension $m$: $\lambda_m(t)$ generally attains strong main maxima at ISSN minima, has secondary maxima at ISSN maxima and minima at the time of the descending and ascending phases of the Schwabe cycles. Such a pattern of “half-Schwabe cycles”, with a large amplitude of $\lambda$ main maxima, remained stable between 1850 and 1915, then changed to a new pattern (with significantly smaller maxima) that remained stable from 1935 to 2005. We interpreted this pattern change as an indication of a “hidden” change in the regime of solar activity, the years 1915–1935 being a transitional interval. We could reproduce the observed behavior of $\lambda$ with a synthetic signal, consisting of an autoregressive process of order 1 with Poisson noise, modulated by an 11-year sine. The switch between the two regimes was obtained by a change in autocorrelation, itself linked to the lifetime of sunspots. In a second paper (Shapoval et al., 2014), we found additional evidence of the two regimes of the irregularity index using embedding dimensions from 3 up to 32. During the first regime R1, from 1850 to 1915, $\lambda$ values were larger than during the second regime R2. The difference is most remarkably seen at the minima of the Schwabe cycles. The value of $\lambda$ at the recent minimum between cycles 23 and 24 was found to be as large as the largest value of $\lambda$ prior to 1915, and much larger than values between 1915 and 2000. This could signal a return of solar activity to regime R1. In Shapoval et al. (2014), we established that the two regimes of $\lambda$ were stable with respect to the parameters used in the computation and to detrending (“decycling”) of the Schwabe cycles.

In Shapoval et al. (2013, 2014), we studied the two regimes of the irregularity index with embedding dimension $m$ between 3 and 16. In the present paper, we concentrate

on the smallest values of $m$ (1 and 2). However our analysis cannot be performed at many solar minima because the distances between nearby points in the phase space contain too many zeros. Therefore, we first preprocess the data by smoothing them over $N=162, 324$ or $648$ successive values (these numbers are chosen as multiples of 27 to suppress the influence of solar rotation on the times series).

Several authors have suggested that observed solar (magnetic) time series are generated by an (as yet) unknown low-dimensional dynamical system (see Zhang, 1996, and Sello, 2001, for a review and original results). Attempts to reconstruct the dynamical system and to use it to predict the future behavior of the time series have led, according to the reports of their authors, to reasonable medium-term predictions of the Schwabe cycle. The efficiency of different predictions is out of the scope of this paper. Pesnell’s review (2012) of the prediction of ongoing cycle 24 together with Love and Rigler’s (2012) and Choudhuri and Karak’s (2012) finding of random walk properties exhibited by some cycle-to-cycle characteristics constitute a useful introduction to the subject. The horizon of the predictions based on chaotic models is linked to the estimates of Lyapunov exponents (Bershadskii, 2009; Zhang, 1996; Sello, 2001). In these studies, Lyapunov exponents are focused on the low-frequency part of the data spectrum, and the dynamical system is reconstructed based on at least decades of observation. In the present paper, we use the irregularity index with embedding dimensions $m=1$ and 2 to characterize higher-frequency variations of ISSN in the period range of the QBO.

Section 2 recalls the definition of the irregularity index and previous attempts to use them in trying to characterize the solar dynamo. Section 3 illustrates further applications of the irregularity index to the Wolf number ISSN and also to the geomagnetic index aa, with results on the evolution of its higher frequency content. A simple autoregressive model is next constructed in Sect. 4, in order to try and reproduce some of the observed properties of the irregularity index and, in particular the appearance (depending on the fundamental parameters of the irregularity index and of model parameters) of half-Schwabe cycle peaks. The discussion and conclusion are given in Sect. 5.

## 2 Basic tools

This section recalls the definition of classical Lyapunov exponents and of the irregularity index first introduced in Shapoval et al. (2013). Further remarks that may be useful to better appreciate the characteristics of the method are given in the Appendix.
## 2.1 Theoretical background

### 2.1.1 Definition

Lyapunov exponents are well defined for dynamical systems. Let $F$ map a $m$ dimensional Euclidian space $\Omega$ into itself. The Lyapunov exponent $\lambda$ measures the rate of exponential convergence or divergence of initially close points in a phase space under the map $F$:

$$ \|\mathbf{J}\varepsilon\| \sim \|\varepsilon\| e^{\lambda}, \quad \varepsilon \in \Omega, $$

where $\mathbf{J}$ is the linear part (Jacobian matrix) of $F$, $\| \cdot \|$ is the norm in the phase space, and $\| \varepsilon \|$ is small.

Formally, we define the trajectory $U_0, U_1, U_2, \dots,$ as

$$ U_1 = F(U_0), \quad U_2 = F(U_1), \quad \dots, $$

for an arbitrary point $U_0$ of the phase space. The small distance $\varepsilon_n$ in the neighborhood of $U_n$ becomes $\varepsilon_{n+1} = J(U_n)\varepsilon_n$ under the map $F$. Thus,

$$ \varepsilon_{n+1} = J_n \varepsilon_0, \quad J_n = J(U_n) J(U_{n-1}) \dots J(U_0). $$

The limit,

$$ \lambda = \lim_{n \to \infty} \lim_{\varepsilon_0 \to 0} \frac{1}{n} \log \left( \frac{\| J_n \varepsilon_0 \|}{\| \varepsilon_0 \|} \right), \quad (1) $$

is the Lyapunov exponent. For so-called ergodic systems, the limit (Eq. 1) is the same for almost any initial point $U_0$ (Oseledets, 1968; Eckmann and Ruelle, 1985).

### 2.1.2 Reconstruction of a dynamical system

Sometimes, when an underlying dynamical system does exist but is not known, a Lyapunov exponent can still be computed for a time series $u_1, u_2, \dots, u_L$, allowing one to reconstruct the key features of the dynamical system, i.e., the embedding dimension $m$ and the map $F$ (Wolf et al., 1985; Rostenstein et al., 1993; Kantz, 1994). The vectors of the phase space are supposed to be

$$ \begin{aligned} U_1 &= (u_1, u_{T+1}, \dots, u_{(m-1)T+1}), \\ U_2 &= (u_2, u_{T+2}, \dots, u_{(m-1)T+2}), \end{aligned} \quad (2) $$

and so on, where $T$ is a delay. The Lyapunov exponent is computed for the map $F$ defined on the set $\{U_i\}$ by

$$ F(U_n) = U_{n+1}. \quad (3) $$

A corollary of the fundamental Takens theorem (Takens, 1981) underlies this computation: let the time series be a projection of the orbit of a dynamical system that lies on its attractor $A$ and be dense on it. Then the Lyapunov exponents of the attractor $A$ and of the set $\{U_i\}$ are the same for an arbitrary delay $T$.

### 2.1.3 Standard computational technique

In practice, the time series under study are always finite and noisy. Values of the delay $T$ and of the embedding dimension $m$ must first be selected in order to estimate the Lyapunov exponent. The delay is frequently taken to be the time of the first minimum of the autocorrelation function of the series, or that of its mutual information (Fraser and Swinney, 1986). The embedding dimension is chosen to be the minimal value $m$ such that the map $F$ transforms a neighborhood of each point $U_i$ defined in Eq. (2) into a neighborhood of $F(U_i)$.

Given $T$ and $m$, the Lyapunov exponent is to be inferred from the quantities

$$ \log \left( \frac{\| J(U - V) \|}{\| U - V \|} \right) \quad (4) $$

for sufficiently close points $U, V$ in the phase space $\Omega$. Algorithms introduced by Rostenstein et al. (1993) and Kantz (1994) have been used with success in recent analyses of solar time series (Macek et al., 2006; Li and Li, 2007). In order to circumvent the rather slow computation of the Jacobian, Ding and Li (2007) use the initial nonlinear map $F$ rather than its linearization $J$ when computing the ratio (Eq. 4).

## 2.2 Definition of the irregularity index

Based on the standard technique described in Sects. 2.1.2 and 2.1.3, this section introduces a straightforward definition of the quantity computed in the paper. In order to determine the irregularity index, we relax the requirement that close points in the phase space must be remote along the time axis, contrary to what is done for the Lyapunov exponent.

### 2.2.1 Phase space

We consider a sliding window of $L$ values $u_1, u_2, \dots, u_L$, where $u_i$ is the $i$th daily value of a given index, counted within the window. Given the embedding dimension $m$ and delay $T$, define the vectors $U_i$ in the phase space by Eq. (2).

### 2.2.2 The map

Let $F$ be the displacement along the orbits given by Eq. (3).

### 2.2.3 Nearest neighbors

For each $U_i$, find the nearest point $U_j$ which does not coincide with $U_i$. Specifically, take $j = \Psi(i)$ such that $\text{dist}(U_i, U_j) = \min_{U_i \neq U_j} \text{dist}(U_i, U_j)$, and $l = 1, 2, \dots, L$; the distance between two vectors $U_i$ and $U_j$ is the square root of the sum of the squares of the differences between each vector coordinate: $(\sum_{k=1}^{m} (u_{(k-1)T+j} - u_{(k-1)T+i})^2)^{1/2}$.
### 2.2.4 Space-close points

We next build the sequence $Θ$ of the distances corresponding to the different<fnref n="1" /> pairs $(U_i, U_{\Psi(i)})$, where $i$ goes from 1 to $L$. Let $\tilde{L} = |Θ|$ be the number of these pairs and $d^*$ be the left $α$ quantile ($α \in [0, 1]$) of $Θ$; in other words, the pairs $(U_i, U_{\Psi(i)})$, and $i \in \{1, ..., L\}$ are ordered according to the distance between the two elements of each pair, so that the ordered sequence is $\{U_{i_k}, U_{\Psi(i_k)}\}$ $k=1, ..., \tilde{L}$, where $\text{dist}(U_{i_k}, U_{\Psi(i_k)}) \le \text{dist}(U_{i_k+1}, U_{\Psi(i_k+1)})$. $P$ is defined as the first $α$ fraction of the ordered pairs, i.e., $P = \{(U_{i_k}, U_{\Psi(i_k)}) : k \le \tilde{L}α\}$.

### 2.2.5 Small distances for Eq. (2)

We enlarge $P$ to the set $\tilde{P}$ by adding the pairs displaced along the orbit of each $(U_{i_k}, U_{\Psi(i_k)}) \in P$ until the distance between the elements in each pair becomes large enough (see below) or the end of the window is reached. Formally,

$$ \tilde{P} = \bigcup_{k=1}^{\tilde{L}} \left\{ (U_{i_k}, U_{\Psi(i_k)}), (F U_{i_k}, F U_{\Psi(i_k)}), \dots, (F^l U_{i_k}, F^l U_{\Psi(i_k)}) \right\} $$

such that (i) $\text{dist}(F^{l'} U_{i_k}, F^{l'} U_{\Psi(i_k)}) \le d^*$, $l' = 0, 1, ..., l$, and (ii) either $l + i_k = L$ or $l + \Psi(i_k) = L$ (the end of the window is reached) or $\text{dist}(F^{l+1} U_{i_k}, F^{l+1} U_{\Psi(i_k)}) > d^*$ (the distance is large enough).

### 2.2.6 Irregularity index

For each pair $(U, V) \in \tilde{P}$, we compute $\log[\text{dist}(FU, FV)/\text{dist}(U, V)]$ and define the irregularity index $λ$ as the median of these numbers. The computed irregularity index is assigned to the middle of the sliding window of length $L$. Let $g$ be the lag between two successive sliding windows (we use $g = L/8$; i.e., 6 months); we construct the new time series:

$$ \lambda_{L/2}, \lambda_{L/2+g}, \lambda_{L/2+2g}, \dots, $$

consisting of the irregularity indices found for sliding windows $[1, L]$, $[g+1, g+L]$, $[2g+1, 2g+L]$, and so forth. This new time series is considered as an additional solar index series.

## 2.3 Some specifics of the irregularity index

Many papers have aimed at reconstructing the dynamical system underlying long time series, such as the daily Wolf numbers (Spiegel and Wolf, 1987; Lawrence et al., 1995), the monthly Wolf numbers (Ruzmaikin et al., 1992; Price et al., 1992), and also some yearly series (Ostryakov and

Usoskin, 1990). These series are either sufficiently long (tens of thousands of points for daily Wolf numbers) or smooth (since at least the 27-day variations are averaged). The embedding dimension for these systems is generally taken to be at least 7, and the delay $T$ is of the order of months (see e.g., Greenkorn, 2009, for a summary table). However, smaller embedding dimensions ($4 < m < 9$) have been used by Greenkorn (2009) for daily data over one Schwabe cycle. The latter paper shows that the Lyapunov exponent is only weakly sensitive to the value of delay $T$ as long as it remains small (a few days). The orbit corresponding to the time series mentioned above now and then returns to the same regions in the phase space. Usually, only points that are close in the phase space but far from each other on the time axis are used to estimate the Lyapunov exponent (Rosenstein et al., 1993).

In this paper, on the contrary, we do not set any limit to the distance in time of points that are close in the phase space, and the exact definition of points being “close” is adapted to the data being studied, using the $α$ quantile of the smallest distances, as explained above. When the embedding dimension $m=1$, the smallest positive distance between two points in the phase space is 1 because of the integer nature of ISSN. We find that at the minima of the ISSN series many points lying at distance 1 are mapped along the corresponding trajectories to points lying at exactly the same distance, so that the values of the irregularity index computed at the signal minima can be inadequate (many values of the ratio – Eq. 4 – are zero). The transition from $m=1$ to $m=2$ changes the properties of the exponent. Whereas we studied embedding dimensions $m$ from 4 to 6 in Shapoval et al. (2013) and up to 32 in Shapoval et al. (2014), we concentrate here on the cases $m=1$ and 2 that are the simplest, and because they shed light on the occurrence of the quasi-biennial variations, at the focus of the present paper.

We first smooth the ISSN daily series and investigate the properties of the irregularity index computed for different smoothing, with delay $T=1$. We compute the irregularity index within a 4-year sliding window; the choice of this 4-year length is a compromise between two opposite requirements: first, the window must be sufficiently large to obtain a stable determination of the irregularity index; second, it should be shorter than the Schwabe cycles. We have checked that the values of the irregularity index calculated as explained are only weakly sensitive to changes of window length inside a 3–5-year interval.

## 3 Data analysis

Solar activity is estimated in the paper with the Wolf (ISSN, sunspot) numbers, involving the number of groups and the number of spots in each particular group. The number of groups reflects the emerging magnetic field and is an indicator of activity. The number of spots within a group depends on the magnetic field as such and also on the interaction

<footnote n="1">¹If $U_i$ is the nearest neighbor of $U_j$ and $U_j$ is the nearest neighbor of $U_i$ then the distance $\text{dist}(U_i, U_j)$ is considered only once.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The irregularity index $\lambda$ (blue) computed in a 4-year sliding window for the ISSN averaged over 162 (top panel), 324 (middle panel), and 648 (bottom panel) days; embedding dimension $m=1$. The Wolf numbers averaged first over $N=162$, 324, and 648 days, and then over 4 years are shown in red, together with solar cycle number. Dashed black vertical lines are located at the maxima of ISSN.</figcaption></figure>

between the magnetic and velocity fields. In this paper, we mainly study sunspot numbers, but we also present some preliminary results on group sunspot numbers (Hoyt and Schatten, 1998) in order to check that results are not affected by the way in which ISSN is determined (including possible data heterogeneities).

The Wolf (sunspot) numbers (ISSN) are defined as $K$ ($10G+s$), where $G$ is the number of sunspot groups, $s$ is the number of individual spots, and $K$ is a factor that is relative to the observer. Daily data series of ISSN are available from 1849 onwards (SIDC-team, 2005). We now apply our algorithm to the ISSN series $w(t)$. The time series is strongly affected by the ~27-day signal connected to solar rotation and reflected axis asymmetry of solar activity (Bartels, 1934; Kitchatinov and Olemskoy, 2005; Howe 2009; Le Mouël et al., 2007). Therefore the series is first smoothed in order to reduce the influence of the Sun's rotation. Namely, it is averaged over multiples of 27 (days), i.e., $N=162 (27 \times 6)$, 324 (27 × 12), and 648 (27 × 24) days. This results in new series

$$ISSN_{(t)} = \sum_{k=t-[N/2]}^{t-[N/2]+N} w(k), \text{ where } [x] \text{ is the integer part of } x.$$

### 3.1 Case $m=1$

The evolution of the irregularity index $\lambda$ computed with $m=1$ and the three values of $N$ given above is shown (in blue) in Fig. 1, together with the original ISSN series smoothed over the same 4-year window (in red). With

162-day averaging (Fig. 1a), there is a clear one-to-one correspondence between Schwabe cycles of ISSN and $\lambda$. The maxima of the cycles coincide precisely with each other in time. The $\lambda$ “Schwabe cycles” exhibit asymmetry: the rising segments are shorter and steeper than the decreasing ones. There is some structure in the decreasing segments, sometimes in the form of a secondary maximum; minima in the irregularity index cycles occur later than minima in the ISSN series.

When the data are smoothed over larger windows, oscillations with a period close to 5.5 years, i.e., half the period of the Schwabe cycle, appear (Fig. 1b, c). We call these “half-Schwabe variations” (HSV; see Shapoval et al., 2013, 2014). In the following, we use HSV to refer to the presence of irregularity maxima at solar minima (thus generating a 5.5-year quasi-periodicity), since the irregularity maxima at solar maxima are almost always present. We also sometimes refer to HSV to refer to the amplitude or amplitude changes of the irregularity peaks and their ratios (see below). In Fig. 1c, both maxima and minima of the ISSN Schwabe cycles correspond to maxima of HSV. Averaging over 324 days leads to an intermediate behavior of the irregularity index (Fig. 1b): secondary peaks at solar cycle minima appear clearly in the 1870s and 1880s and after 1950, but some are not or hardly visible at 1865, 1900 or 1975.

In order to provide a more quantitative measure of HSV behavior, we determine the ratio $R$ of the amplitude of $\lambda$ oscillations near maxima of ISSN ($\Delta S_{max}$) to that near minima
of ISSN ($\delta S_{min}$; $R = \delta S_{min} / \Delta S_{max}$ in Fig. 2; the method is introduced in Shapoval et al. (2013) and further explained in the present paper). Figure 2 presents a schematic “Schwabe cycle” (smoothed artificial signal in red; actually somewhat more than one full period) and its irregularity index $\lambda$ in blue. $\lambda$ attains main maxima $\lambda_{Smax}$ at the maxima of the (smoothed) original signal, secondary maxima $\lambda_{Smin}$ at the minima of the signal, and its minima $\lambda_{mid}$ on the descending and ascending phases of the signal (subscripts in this notation correspond to the signal not to $\lambda$ itself). Three local minima occur in Fig. 2 (because somewhat more than one cycle is represented) $\lambda_{mid}^1$, $\lambda_{mid}^2$, and $\lambda_{mid}^3$. $\lambda_{mid}$ is defined as their mean. Let $\Delta S_{max} = \lambda_{Smax} - \lambda_{mid}$ and $\delta S_{min} = \lambda_{Smin} - \lambda_{mid}$. $\Delta S_{max}$ and $\delta S_{min}$ measure the amplitude of $\lambda$ peaks at signal maxima (minima) when they exist. Finally, $R = \delta S_{min} / \Delta S_{max}$ measures “HSV performance”, such that a decrease of R accompanies a clearer appearance of HSV.

The quantities $\lambda_{Smax}$, $\lambda_{Smin}$, and $\lambda_{mid}$ can be determined even if the HSV structure is subdued. $\lambda_{Smin}$ is taken to be the maximal value of $\lambda$ in a neighborhood of the minimum of the smoothed original signal. We extend the construction of R to several solar cycles. In such a case the quantities $\lambda_{Smax}$, $\lambda_{Smin}$, and $\lambda_{mid}$ are obtained by averaging the corresponding quantities for all cycles included in the time window of interest.

We see in Fig. 1 that $\delta S_{min}$ increases significantly (from ~0.1 to ~0.7) when N is increased from 162 to 648 (HSV is actually hardly visible when N = 162 and almost as strong as peaks at solar maxima when N = 648). This is mainly due to the large drop of $\lambda$ on descending and ascending phases as smoothing is increased. In Fig. 1b, we note differences in behavior between the period before and after ~1930 (for instance lower overall mean value of $\lambda$ and larger amplitude of HSV and Schwabe cycles after 1930).

## 3.2 Case m = 2

The above computations are repeated for m = 2 (Fig. 3). Different epochs can readily be distinguished. First, both the amplitude of variations and actual values of $\lambda$ change around 1930, in a much more visible way than in the case m = 1, confirming that 1930 is a time of first-order regime change. This is particularly clear in Fig. 3b where $\lambda$ drops from a mean value of about 0.3 to 0.2. This implies that the ISSN series becomes less irregular after 1930 (see discussion). In Fig. 3a (162-day averaging), HSV is clearly visible in cycles 21–23, i.e., from 1975 onwards; it is still visible from about 1915 to 1975 but is barely recognizable prior to 1915. 1915 and 1975 therefore appear as possible second-order regime changes or at least singularities. Increased smoothing strengthens HSV behavior. In Fig. 3b, with 324-day averaging, HSV is seen with varying shapes and amplitudes from 1867 to 1930 and from 1945 to 2005, with a gap at the times of cycle minima 16–17 and 17–18. In Fig. 3c, with 648-day averaging, HSV is quite clearly present from 1867 to 1915, but it is subdued from 1975 to 2005. We further note that some $\lambda$ maxima at

solar minima are enhanced by increasing smoothing prior to 1915 (cycles 12–13 and 13–14), whereas others remain similar or are reduced after 1975 (cycles 21–22 and 22–23).

Using the same notations as introduced for the case m = 1, we see (Table 1a) that indeed $\lambda_{Smin}$ is larger for the period 1867–1915 than for 1975–2008, both decreasing with increasing N, the former from 0.37 to 0.22 and the latter from 0.30 to 0.10, when N increases from 162 to 648. We see in Fig. 3 that HSV is better marked when N is increased for the period 1867–1915 and when N is decreased for the period 1975–2008. This can be expressed by the evolution of the ratio $R = \delta S_{min} / \Delta S_{max}$ of the mean amplitudes of HSV peaks at solar minima vs. solar maxima as a function of N. When N increases from 162 to 648, this ratio increases from 0.67 to 0.79 for the period 1867–1915 but decreases from 0.71 to 0.34 for 1975–2008 (Table 1a).

In summary, the $\lambda$ curves shown in Figs. 1 and 3 allow us to distinguish different epochs. In Fig. 1, we see the very strong appearance of HSV with an amplification of $\lambda$ peaks at solar cycle minima, and an indication of a change in behavior of $\delta S_{min}$ values as a function of N before and after ~1930. In Fig. 3, we also see that $\lambda$ decreases in both mean value and amplitude of variations in the 1930s. The R ratio increases with N for the period 1867–1915, but it decreases in 1975–2008. So, evolution of the irregularity index reveals a first-order singular date ~1930. Finer analysis of HSV properties (R ratio evolution as a function of smoothing) reveals second-order singular dates around 1915 and 1975.

## 4 A model

We have already made a first attempt at constructing a model that would embody HSV behavior of irregularity in solar activity as observed in real data in Shapoval et al. (2013). We extend here the method, with slightly different choices of

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Construction of $R = \delta S_{min} / \Delta S_{max}$ (see text); smoothed artificial signal (red) and its irregularity index (blue); main maxima $\lambda_{Smax}$, secondary maxima $\lambda_{Smin}$, and local minima $\lambda_{mid}^i$. Black lines show construction of $\Delta S_{max}$ and $\delta S_{min}$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 3. Blue curves: the irregularity index $\lambda$ computed in 4-year sliding windows for the Wolf numbers averaged over $N$ days; red curves: the Wolf numbers averaged over $N$ days and then over 1461 days (4 years), where $N$ is 162 (top panel), 324 (center panel), and 648 (bottom panel) days; $m=2$. Dashed black vertical lines are located at times of possible regime change of $\lambda$.</figcaption></figure>

relevant parameters, to concentrate on what could be related to the QBO.

## 4.1 Definition of the model

Consider a first-order autoregressive AR(1) process $x(t)$:

$$x(t) = ax(t-1) + \eta(t),$$

where the random variable $\eta(t)$ is Poissonian with mean $\mu(t)$, $P\{\eta = \eta\} = e^{-\mu}\mu^n/n!$. The mean $\mu(t)$ is modulated by the sum of two periodic functions (as opposed to only one in Shapoval et al., 2013) with periods $T_1$ and $T_2$, $T_1 > T_2$. The longer period $T_1$ is set to 11 years (11 × 365 days), corresponding to the Schwabe cycle. We choose the shorter (which we call intermediate) period $T_2$ in the range from 1 to 3 years, such that it includes the QBO. Therefore,

$$\mu(t) = h\left(-\cos\frac{2\pi t}{T_1} - k\cos\frac{2\pi t}{T_2} + c\right), \quad (5)$$

where $h > 0$, $0 < k < 1$; $c > 1 + k$ is a vertical shift. The synthetic signal $w(t)$ is defined by

$$w(t) = [Mx(t)]. \quad (6)$$

$M$ is set to 10 to mimic the factor in the definition of the group sunspot number. This yields the same order of magnitude for synthetic and observed ISSN values.

The model is a function of variable $t$ (time, in days) and depends on five adjustable parameters: $a, h, T_2, k$ and $c$. The

value of parameter $a$ in the auto-regressive process determines the correlation of the data. Modeling the sunspot series by an autoregressive AR(1) model connects $a$ to the lifetime of sunspots (Blanter et al., 2005). $h$ controls the smoothness of the signal. The factor $k < 1$ controls the relative amplitude of the $T_2$ vs. $T_1$ modulations. The vertical shift $c$ controls the ratio of the maximum to the minimum of $\lambda$. Figure 4 shows an example of a realization of the model with parameters given in the legend.

## 4.2 Modeling results

### 4.2.1 Appearance of HSV with data smoothing (case $m=1$)

Already with $m=1$ and without intermediate period modulation $T_2$ ($k=0$), HSV appears in response to increasing data smoothing (Fig. 5). For $N=162$ days, $\lambda$ shows no HSV maxima at solar minima. When $N$ is increased to 648, all $\lambda$ values decrease, but their overall structure changes markedly. $\lambda$ remains approximately the same at solar minima, decreases slightly at solar maxima and falls dramatically in intermediate intervals (corresponding to the ascending and descending phases of the “solar” cycle). As a result, $\lambda$ now peaks sharply not only at solar maxima but also at solar minima: this behavior is indeed reminiscent of that observed for ISSN (compare Figs. 5 and 1a, c), as looked for when building the model.
<table><caption>**Table 1.** Calculation of *R* ratio values (see text). *N* is the number of days over which the ISSN data are smoothed to remove in particular the effect of solar rotation; $\lambda_{S_{min}}$ is the mean value of the maxima of the irregularity index at the times of solar minima over the period range indicated under the heading “years” (by default 1867–2008 when not indicated); $\lambda_{S_{max}}$ is the mean value of the maxima of the irregularity index at the times of solar maxima over the same period range; $\lambda_{mid}$ is the mean value of the minima of the irregularity index as indicated in Fig. 2 in the text; $\delta_{S_{min}} = \lambda_{S_{min}} - \lambda_{mid}$; $\Delta S_{max} = \lambda_{S_{max}} - \lambda_{mid}$; and $R = \delta_{S_{min}} / \Delta S_{max}$. Each part of the table corresponds to a figure as indicated.</caption><thead><tr><th>N</th><th>λ<sub>S<sub>min</sub></sub></th><th>λ<sub>S<sub>max</sub></sub></th><th>λ<sub>mid</sub></th><th>δ<sub>S<sub>min</sub></sub></th><th>δ<sub>S<sub>max</sub></sub></th><th>R</th><th>years</th></tr></thead><tbody><tr><td colspan="8" style="text-align: left;">(a) (corresponding to Fig. 3)</td></tr><tr><td>162</td><td>0.37</td><td>0.43</td><td>0.23</td><td>0.14</td><td>0.20</td><td>0.67</td><td>1867–1910</td></tr><tr><td>324</td><td>0.31</td><td>0.35</td><td>0.15</td><td>0.16</td><td>0.19</td><td>0.83</td><td></td></tr><tr><td>648</td><td>0.22</td><td>0.25</td><td>0.10</td><td>0.12</td><td>0.15</td><td>0.79</td><td></td></tr><tr><td>162</td><td>0.30</td><td>0.34</td><td>0.18</td><td>0.12</td><td>0.17</td><td>0.71</td><td>1970–2008</td></tr><tr><td>324</td><td>0.21</td><td>0.26</td><td>0.12</td><td>0.09</td><td>0.14</td><td>0.64</td><td></td></tr><tr><td>648</td><td>0.10</td><td>0.18</td><td>0.06</td><td>0.04</td><td>0.12</td><td>0.34</td><td></td></tr><tr><td colspan="8" style="text-align: left;">(b) (corresponding to Fig. 6)</td></tr><tr><td colspan="8" style="text-align: left;">Fig. 6a (T = 450, a = 0.8, h = 0.4, c = 1.7, k = 0.35)</td></tr><tr><td>162</td><td>0.34</td><td>0.48</td><td>0.30</td><td>0.03</td><td>0.17</td><td>0.19</td><td></td></tr><tr><td>648</td><td>0.21</td><td>0.32</td><td>0.11</td><td>0.10</td><td>0.21</td><td>0.47</td><td></td></tr><tr><td colspan="8" style="text-align: left;">Fig. 6b (T = 610, a = 0.8, h = 0.4, c = 1.7, k = 0.35)</td></tr><tr><td>162</td><td>0.37</td><td>0.47</td><td>0.31</td><td>0.06</td><td>0.16</td><td>0.39</td><td></td></tr><tr><td>648</td><td>0.24</td><td>0.35</td><td>0.10</td><td>0.14</td><td>0.25</td><td>0.56</td><td></td></tr><tr><td colspan="8" style="text-align: left;">Fig. 6c (T = 700, a = 0.8, h = 0.4, c = 1.7, k = 0.35)</td></tr><tr><td>162</td><td>0.35</td><td>0.46</td><td>0.30</td><td>0.05</td><td>0.15</td><td>0.31</td><td></td></tr><tr><td>648</td><td>0.24</td><td>0.31</td><td>0.10</td><td>0.14</td><td>0.21</td><td>0.70</td><td></td></tr><tr><td colspan="8" style="text-align: left;">Fig. 6d (T = 800, a = 0.8, h = 0.4, c = 1.7, k = 0.35)</td></tr><tr><td>162</td><td>0.36</td><td>0.51</td><td>0.32</td><td>0.05</td><td>0.19</td><td>0.25</td><td></td></tr><tr><td>648</td><td>0.18</td><td>0.34</td><td>0.10</td><td>0.08</td><td>0.23</td><td>0.35</td><td></td></tr><tr><td colspan="8" style="text-align: left;">(c) (corresponding to Fig. 7)</td></tr><tr><td colspan="8" style="text-align: left;">Fig. 7b (a = 0.8, h = 0.4, k = 0, c = 1.2)</td></tr><tr><td>162</td><td>0.14</td><td>0.54</td><td>0.09</td><td>0.05</td><td>0.45</td><td>0.12</td><td></td></tr><tr><td>648</td><td>0.09</td><td>0.30</td><td>0.05</td><td>0.04</td><td>0.25</td><td>0.16</td><td></td></tr><tr><td colspan="8" style="text-align: left;">Fig. 7c (a = 0.8, h = 0.4, k = 0, c = 1.6)</td></tr><tr><td>162</td><td>0.43</td><td>0.54</td><td>0.32</td><td>0.12</td><td>0.22</td><td>0.53</td><td></td></tr><tr><td>648</td><td>0.20</td><td>0.33</td><td>0.09</td><td>0.11</td><td>0.24</td><td>0.46</td><td></td></tr><tr><td colspan="8" style="text-align: left;">Fig. 7d (a = 0.8, h = 0.4, k = 0, c = 1.7)</td></tr><tr><td>162</td><td>0.48</td><td>0.53</td><td>0.33</td><td>0.15</td><td>0.20</td><td>0.77</td><td></td></tr><tr><td>648</td><td>0.27</td><td>0.30</td><td>0.10</td><td>0.17</td><td>0.20</td><td>0.83</td><td></td></tr><tr><td colspan="8" style="text-align: left;">(d) (corresponding to Fig. 8)</td></tr><tr><td colspan="8" style="text-align: left;">(a = 0.8, h = 0.8, k = 0, c = 1.7)</td></tr><tr><td>162</td><td>0.40</td><td>0.46</td><td>0.28</td><td>0.12</td><td>0.17</td><td>0.69</td><td></td></tr><tr><td>648</td><td>0.15</td><td>0.26</td><td>0.07</td><td>0.08</td><td>0.18</td><td>0.43</td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** A realization of the AR(1) process introduced in Sect. 3.1, shown prior to smoothing by *N* days and over the 4-year interval over which it will next be averaged (see text). Model parameters: *a* = 0.8, *h* = 0.4, *k* = 0.35, *c* = 1.7, *T*₁ = 11 years, and *T*₂ = 700 days.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The irregularity index computed in 4-year sliding windows for synthetic data averaged over 162 (solid gray) and 648 (dashed black) days (*m* = 1). Model parameters: *a* = 0.8, *h* = 0.4, *k* = 0.0, and *c* = 1.7. Red dashed vertical lines at solar sunspot minima.</figcaption></figure>

### 4.2.2 Increase in HSV behavior as a function of data smoothing (case *m* = 2)

The behavior of the irregularity index for *m* = 2 is significantly richer. We set intermediate period variations at a relatively strong level (*k* = 0.35; Fig. 6): the irregularity index exhibits HSV that increases as smoothing is increased (see Fig. 6a–c when averaging interval goes from 162 days – solid gray lines – to 648 days – dashed black lines). In the case when *T*₂ = 610, *N* = 648, strong HSV peaks are always present except in 1 out of 10 possible occurrences (at *y* = 11; Fig. 6b). On the contrary, for the same *T*₂ but with *N* = 162, HSV peaks are quite subdued (yet generally visible), though again only 1 out of 10 is missing (at *y* = 41; Fig. 6b). In that case, increasing data smoothing in the model results in amplifying HSV behavior: *R* increases from 0.39 to 0.56
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The irregularity index computed in 4-year sliding windows for synthetic data averaged over 162 (solid gray) and 648 (dashed black) days ($m=2$). Model parameters: $a=0.8$, $h=0.4$, $c=1.7$, $k=0.35$, and intermediate period variation $T_2$ set at 450 (a), 610 (b), 700 (c), 800 (d) days. Red dashed vertical lines at solar sunspot minima.</figcaption></figure>

(Table 1b). HSV increases (that is *R* increases) with smoothening when $T_2$ is in the interval [450, 700] (Fig. 6a-c). The effect slowly disappears when $T_2$ reaches 800 (Fig. 6d).

Figure 7 illustrates the effect of changing the value of the “vertical shift” *c* when there is no intermediate period modulation ($k=0$). HSV behavior becomes increasingly significant as *c* is increased: *R* grows from ~0.1 to 0.8 as *c* increases from 1.2 to 1.7 (Table 1c).

Comparing Fig. 6c (where $k=0.35$) to Fig. 7d (where $k=0$, $a=0.8$, $h=0.4$, and $c=1.7$ in both cases), we see that HSV behavior is more visible in the case of a smaller *N* (162) when *k* is smaller (*R* is then respectively 0.31 vs. 0.77; Table 1c). When $k=0$, there is little or no HSV increase (*R* increases from 0.77 to 0.83; Table 1c; Fig. 7b-d) whereas with $k=0.35$ it grows significantly (from 0.31 to 0.70; Table 1b). Therefore, *k* is an important factor controlling HSV behavior.

### 4.2.3 Decrease of HSV as a function of data smoothing (case $m=2$)

When intermediate period ($T_2$) variations are suppressed and parameter $h$ (that controls the smoothness of the signal) is increased (Fig. 8), we find another regime in which HSV decreases with increasing smoothing. When *N* increases

from 162 to 648, the *R* ratio decreases from 0.69 to 0.43 (Table 1d).

### 4.2.4 A direct comparison of model with observations

In Fig. 9, we model the data of cycles 21–23 (which have similar durations, to allow a comparison with a model where $T_1=11$ years). Parameter *k* is set to zero and *a* to 0.9. This choice of *a* reflects the increase of the lifetime of sunspots found by Blanter et al. (2005). We can now directly compare the irregularity index computed for the synthetic and actual signals.

Since the model contains a random ingredient and the computation of the irregularity index is sensitive to particular realizations, two of them are shown in Fig. 9. The irregularity index for the model follows rather precisely that for ISSN for *N* = 162 (Fig. 9, middle row panels, except at the minimum between cycles 22 and 23 in the realization on the right side). The quality of the fit is somewhat less for *N* = 648 (Fig. 9, bottom row panels). Nevertheless, the model realizations globally reproduce the *λ* pattern of the real data quite faithfully.

An explicit comparison (Fig. 9) of the irregularity index for the model and ISSN time series is possible for cycles 21–23 (1975–2005) because *λ* constructed with ISSN exhibits a smooth, quasi-cyclic and regular behavior, as is the case for
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The irregularity index computed in 4-year sliding windows for synthetic data averaged over 162 (solid gray) and 648 (dashed black) days (m=2). Model parameters are: a=0.8, h=0.4, k=0, and c=1.05 (a), 1.2 (b), 1.6 (c) and 1.7 (d). Red dashed vertical lines at solar sunspot minima.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** The irregularity index computed in 4-year sliding windows for synthetic data averaged over 162 (solid gray) and 648 (dashed black) days (m=2). Model parameters: a=0.8, h=0.8, k=0, and c=1.7. Red dashed vertical lines at solar sunspot minima.</figcaption></figure>

model realizations. The regime observed in 1867–1915 does not display such regular cycles of $\lambda$ and therefore does not allow such an easy comparison.

**5 Summary, discussion and conclusion**

The evolution of the daily values of sunspot number ISSN from 1850 to 2005 has been studied in this paper, using tools

from dynamical systems. Some interesting results are obtained with the irregularity index, a new method introduced in Shapoval et al. (2013, 2014). The method computes the rate of divergence of close trajectories in the phase space under a one-step translation mapping. This index is akin to the maximal Lyapunov exponent for time series calculated for low-dimensional dynamical systems, but is applicable to short time series with a random component. We have computed the irregularity index $\lambda$ of ISSN for embedding dimensions $m=1$ and $2$ within a 4-year sliding window, after first averaging the data over $N=162, 324$ and $648$ days (multiples of the solar rotation period). The irregularity index for $N=162$ follows the Schwabe cycle (Fig. 1), with sharp, high peaks at solar cycle maxima. But when $N$ becomes large enough, it also exhibits sharp maxima at solar cycle minima (see also Shapoval et al., 2013), resulting in 5.5-year time variations, i.e., half the period of the Schwabe cycles (Fig. 1 middle and bottom panels; Fig. 3); we call them *half-Schwabe variations* (HSV).

The mean level of the irregularity index for ISSN undergoes a downward step around 1930 (particularly clear with embedding dimension $m=2$, as seen in Fig. 3). This can be linked to the observation by Bershadskii (2008) that a change in the fractal properties of ISSN took place at that time. For a given time period, HSV can be characterized by the mean differences $\delta_{S_{min}}$ of $\lambda$ values at the times of Schwabe
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top row panels: two synthetic signals (blue; see text) and ISSN (red) averaged over 4 years. The two columns of the figure are relative to different synthetic signals. The irregularity index ($m=2$) for model (blue) and ISSN (red) series are shown averaged over 162 (middle row) and 648 days (bottom row). $a=0.9$, $h=0.4$, and $k=0$.</figcaption></figure>

cycle minima ($\lambda_{S_{min}}$) and the mean of $\lambda$ minima ($\lambda_{mid}$) at the middle times of the descending and ascending phases of the Schwabe cycles; we use the ratio $R$ of $\delta_{S_{min}}$ over its equivalent $\Delta_{S_{max}}$ taken at the times of Schwabe cycle maxima. A first regime (denoted by Q1) is characterized by $R$ increasing with $N$, and a second one (Q2) has $R$ decreasing with $N$. For ISSN, with $m=2$, Q1 is observed until ~1915, whereas Q2 appears after ~1975; the main transition may be around 1930 (Fig. 3).

HSV as such may not be regarded as a result with great importance. Our functional $\lambda$ can attain its extrema on both ascending and descending phases. If such is the case, HSV appears because of a certain similarity between these ascending and descending phases. That is why we cannot yet discuss the physics underlying the essence of HSV. However, we are entitled to look for simple time series that would display the properties observed for ISSN.

A synthetic signal, generated by a simple autoregressive model of order 1 exhibits many of the above-mentioned properties of the irregularity index of ISSN. The random part of this synthetic signal is taken to be Poissonian. Its mean is modulated by the sum of two periodic functions with periods $T_1 = 11$ years and $T_2 < T_1$, the latter being tunable in an interval that can range from months to years. The introduction of intermediate oscillations ($T_2$) allows one to reproduce both the Q1 and Q2 regimes (Table 1). When the embedding dimension $m$ is 1, HSV (5.5-year pseudoperiod) oscillations appear even if $T_2$ variations are absent ($k=0$; Fig. 5). When $m$ is equal to 2, the behavior of the irregularity index

series becomes richer: regime Q1, in which HSV behavior increases with smoothing $N$, is observed for larger values of $k$ (Fig. 6), whereas regime Q2, in which HSV decreases with $N$, is obtained when intermediate period ($T_2$) variations are absent ($k=0$) and parameter $h$ (which controls the smoothness of the signal) is increased (Fig. 8). We conclude that high frequency components of ISSN have much in common with an AR(1) process. The presence and then disappearance of ~1-2-year ($T_2$) oscillations seem to be required to produce a transition between regimes Q1 and Q2 (when $m=2$). We propose that these oscillations may be linked to the QBO, the second most powerful solar variation after the 11-year cycle (e.g., Ivanov et al., 2002).

At first order, the observed change in the mean level of $\lambda$ around 1930 found in this paper could mark a shift of solar activity to a new regime (a transition of the solar dynamical system to a new state). This regime change is also marked by a (second order) change in the way HSV amplitude varies as the data is increasingly smoothed. These observed features can be reproduced by the model: the $R$ ratio increases with increasing $N$ (Q1) prior to 1915 and decreases after 1975 (Q2). Although several model parameters interact to promote one or the other regime, the most important one appears to be parameter $k$, which reflects the presence or absence of intermediate $T_2$ variations in the process. The shift of the irregularity index of ISSN from regime Q1 to Q2 may be due to the decrease or even disappearance of QBO.

In contrast to a standard statistical analysis, one cannot introduce a reasonable null hypothesis in the present study. Instead, we have checked the stability of the observed phenomena with respect to the parameters that control the irregularity index and we have tested the significance of our conclusions with the auto-regressive model.

The homogeneity of the ISSN-series is a long-debated question. Svalgaard (2010, 2012) points to an abrupt increase of ISSN in ~1945 and argues that this increase is caused by changes in the measurement rules. The NASA website (http://solarscience.msfc.nasa.gov/greenwch.shtml) also notes that the sunspot series is not uniform; abrupt changes occurred in 1941-1942 (sunspot numbers) and 1976-1977 (sunspot areas, not used in our paper). However, our conclusions about regime changes are not seriously affected by such events, because we use ratios (Eq. 4). Moreover, the date of the ~1930 singularity is remote from 1941 to 1942 (or 1945).

In order to see whether the observed behavior of the irregularity index of ISSN could be affected by such data problems, we have computed the irregularity index for another proxy influenced by solar activity but derived completely independently, namely the geomagnetic index *aa* (available at http://isgi.latmos.ipsl.fr/source/indices/aa/). With $m=1$, we computed the irregularity index for *aa* as such (without any prior averaging over multiples of 27 days). The values of $\lambda(t)$ (Fig. 10) exhibit both a noticeable increase in mean level (from about 1.75 to 2) and a decrease in range (from 0.5 to 0.2) in the 1930s. The sign of the change in mean value
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Red curve: the daily *aa* averaged over 4 years. Blue curve: the irregularity index computed in a 4-year sliding window with *m* = 1. Vertical dashed lines are at the maxima of Wolf numbers (averaged over 4 years). The number of each Schwabe cycle is indicated.</figcaption></figure>

is opposite to that found for the irregularity index of ISSN (Fig. 10 vs. Fig. 3), but the same singularity in solar behavior could be at the origin of both.

Although this paper focuses mainly on changes of the ir-regularity index with smoothing, we also describe briefly changes of $\lambda$ with time, in order to provide further evi-dence of the robustness of the technique. Our previous pa-per (Shapoval et al., 2013) examines the time evolution of $\lambda$ computed for ISSN with a 4-year sliding window and differ-ent embedding dimensions and finds a change of regime in approximately 1915–1940. The same computation has been repeated for the Hoyt and Schatten group sunspot numbers (GSN, Hoyt and Schatten, 1998). Despite differences in in-homogeneities and potential problems with the two series, the main results are quite similar, excluding the possibility of an artefact due to the choice of an imperfect time series. The irregularity index of GSN exhibits two different regimes with a clear transition in the period 1915–1940 (details and figures in Appendix B). This strengthens the result obtained for ISSN and published in Shapoval et al. (2013) and further supports our approach, as used in the present paper.

It would be good to find some physical evidence to sup-port our hypothesis that the change of regime can be linked to QBO rarefactions. There has been a significant amount of research on oscillations in the 1–2-year period range in cosmic rays (Valdes-Galicia et al., 1996; Kudela et al., 2002; Rouillard and Lockwood, 2004) and in midlatitude coronal holes area (McIntosh et al., 1992). Obridko and Shelting (2007) give a brief review of these works, together with new results (see also e.g., Ivanov et al., 2002). Inter-mediate variations do not seem to have been reported up to now for ISSN, but further interesting observations have been made for *aa* (Lockwood, 2001; Mursula et al., 2003). Using an extended *aa* index over 160 years, Mursula et al. (2003) have found that the power of “mid-term quasi-periodicities”

(identical to QBO) is larger at periods alternating between 1.3 and 1.6 years; maxima of 1.3-year oscillations occur at the maxima of Schwabe cycles 18 and 22, while the 1.6-year oscillations peak at the maxima of cycles 16 and 21. The spectral power of *aa* is high during periods of high solar activity and would reflect the strength of the solar dynamo. Sudden disappearance of power is considered as a precursor for long-term decreases in solar activity. Our observations of a post-1975 decrease of *R* ratio with smoothing and our mod-eling of this regime by using an AR(1) process without $T_2$ variations are in line with the work of Mursula et al. (2003).

At least two different mechanisms could generate QBO. On one hand, Ivanov et al. (2002) show that the QBO of solar magnetic fields are mainly revealed in their large-scale component; they argue that QBO actually reflect variations in the equatorial dipole (and to a lesser extent quadrupole); for these authors, QBO sources are located near the base of the convection zone and remain invariable. Vecchio et al. (2012), using magnetic synoptic maps from 1976 to 2003, propose that QBO are fundamental modes associated with poleward magnetic flux migration from low to high latitudes (part of meridional circulation) during the maximum and descend-ing phases of the solar cycle. A strong link between QBO and the solar dynamo is inferred from these and other works. Time variations of QBO might therefore provide information on changes in meridional flow. On the other hand, nonlinear-ity of the solar dynamo itself could be the source of QBO. Using a nonlinear Babcock–Leighton model, Charbonneau et al. (2007) support the hypothesis that the nonlinear com-ponent of the solar dynamo prevails over the stochastic one. Mayr and Schatten (2012) argue that the strong nonlinearity in the Charbonneau et al. (2007) equations could generate QBO without any time-dependent solar excitation.

The irregularity index of ISSN computed in this paper may provide a measure of the irregular behavior of the solar dy-namo. Duhau and de Jager (2008) propose that ISSN may be used as a proxy of the toroidal component of the Sun’s magnetic field and *aa* of the poloidal component. The irregu-larity index of *aa* as such presents a change in the 1930s, with a sign opposite to that for ISSN. We could therefore interpret our observations of changes in regime of the irregularity in-dices of ISSN and *aa* as indicating respectively a decrease in the irregular character of the toroidal field and an increase in the irregularity of the poloidal field in the 1930s, date of the advent of a grand maximum period in solar activity. Our analysis also suggests that another change may have started around 1975, as witnessed by decreasing HSV as a function of smoothing (see Fig. 3a–c in that order). This may have heralded the 2005 change found in our complementary stud-ies of the irregularity index (Shapoval et al., 2013, Figs. 2–5; Shapoval et al., 2014, Figs. 1–3).

The irregularity index method is promising but still not a fully understood tool. It appears to be able to uncover singu-lar phenomena and solar activity changes that cannot easily be seen by other means, but the tool depends on a number
of parameters, particularly the embedding dimension and
changes of behavior as the embedding dimension is changed.
We note that the regime changes R1 and R2 uncovered by
Shapoval et al. (2013, 2014) are not identical to the regime
changes Q1 and Q2 found in the present paper. The R1/R2
regimes are marked by different levels of the irregularity in-
dex (computed with embedding dimensions from 4 to 32).
In the present paper, when embedding dimension *m* is 2, we
also find some evidence of the R1/R2 regimes (Fig. 3a, b).
But we introduce an additional tool, the analysis of the ir-
regularity index as a function of data smoothing (*N*), and
this is what reveals the Q1/Q2 regimes. Although R1 and
Q1, ending around 1915–1930, could correspond to the same
regime, R2 (starting after 1930 and possibly ending in 2008)
and Q2 (emerging clearly only after 1975) do not coincide.
The physical nature of these singularities and the differences
in their timing and behavior remain to be deciphered.
## Appendix A: Examples of computation of the irregularity index

In order to illustrate some aspects of the computations involved in this paper, we generate synthetic data with the formula

$$u(t) = \sin \frac{2\pi t}{20} + \eta(t), \quad t = 1, 2, \dots, 19,$$

where $\eta(t)$ is a random variable uniformly distributed over [-0.05, 0.05]. The embedding space is one-dimensional ($m=1$).

Figure A1 exhibits 19 points of the sample. For each point, we first find the nearest neighbor. In the present case, only 13 pairs of nearest neighbors are different, so that the set of distances consists of 13 values. Let $\alpha = 0.25$. Since the integer part of $0.25 \times 13$ is 3, the $\alpha$ quantile $d^*$ is the third distance from the lowest one. Distances less than or equal to $d^*$ are “small”, according to our definition. For each pair of nearest neighbors, the map $F$ defined in Eq. (3) that moves points along their trajectory is applied, the distance between the corresponding images being small. In the example above, the pair at times 6 and 19 (red points in Fig. A1) possesses the smallest distance. A move to the right is impossible because 19 is the largest time in the window. Thus, this pair does not contribute to the computation of the exponent.

The points at times 4 and 11 form the next pair (blue filled circles). The distance between $u(5)$ and $u(12)$ (blue empty circles), which is $|u(12) - u(5)|$, is larger than $d^*$; therefore, this pair generates a single quantity log $[|u(12) - u(5)|/|u(11) - u(4)|]$ as a candidate irregularity index.

The last pair under consideration is $[u(1), u(11)]$. Since the value $u(11)$ appears twice in the pairs we consider, the corresponding point in the graph is marked first by a blue circle and second by a green circle (the green circle is smaller). Although the values $u(2)$ and $u(12)$ are very close, the distance between $F(u(1))=u(2)$ and $F(u(11))=u(12)$ is larger than the critical distance. Thus, the quantity log $[|u(12) - u(2)|/|u(11) - u(1)|]$ becomes the second candidate to the irregularity index. The irregularity index is chosen as the median of the candidate values.

## Appendix B: Regime change of $\lambda$ of GSN and ISSN

The irregularity index $\lambda$ of both ISSN and GSN is computed with a 4-year sliding window, 8-day delay, and embedding dimensions 2–32. According to Figs. B1 and B2, $\lambda$ of both the series exhibits two different patterns before ~1915 and after ~1940 with a transition during 1915–1940. The patterns for ISSN differ by the values of the irregularity index. High values of the irregularity index, most markedly seen at the minimum of cycles 14 and 15, underlie the pattern prior to 1915. The second pattern continues with the minimum

<figure><img src="image_14.png" /><figcaption>Figure A1. Computation of $\lambda$. Synthetic signal vs. time (see text).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure B1. The irregularity index $\lambda$ computed for ISSN within 4-year sliding windows; the embedding dimension $m$ is indicated; vertical lines are at solar cycle minima.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

**Figure B2.** The irregularity index $\lambda$ computed for GSN (Hoyt and Schatten, 1998) within 4-year sliding windows; the embedding dimension $m$ is indicated; vertical lines are at solar cycle minima.

of cycles 23 and 24 when $\lambda$ achieves another remarkable maximum.

The irregularity index of GSN also exhibits four high maxima at the minima of cycles 11–15. We cannot check the existence of the cycle 23–24 peak since recent data, following the Hoyt and Schatten technique, are not available as open sources. This confirmation of the regime change of solar activity between 1915 and 1940 vindicates our approach and results.
*Acknowledgements.* IPGP provided support to A. Shapoval and M. Shnirman during their visit to the institute. A. Shapoval was partially supported by RFBR grants 14-01-00346 and 14-01-00773. IPGP contribution NS 3503.

Edited by: G. Lapenta

Reviewed by: L. Svalgaard and two anonymous referees

## References

Bartels, J.: Twenty-Seven Day Recurrences in Terrestrial-Magnetic and Solar Activity, 1923–1933, Terr. Magnet. Atmos. Elect., 39, 201–202, 1934.

Bergé, P., Pomeau, Y., and Vidal, C.: L'Ordre dans le Chaos, Hermann, Paris, France, 353 pp., 1984.

Bershadskii, A.: New dynamics of the Sun convection zone and global warming, arXiv:0805.2108v1 [astro-ph.SR], 2008.

Bershadskii, A.: Chaotic mean wind in turbulent thermal convection and long-term correlations in solar activity, arXiv:0908.4008v4 [astro-ph.SR], 2009.

Blanter, E. M., Shnirman, M. G., and Le Mouël, J.-L.: Solar variability: Evolution of correlation properties, J. Atmos. Solar-Terr. Phys. 67, 521–534, 2005.

Blanter, E. M., Le Mouël, J.-L., Perrier, F., and Shnirman, M. G.: Short-term correlation of solar activity and sunspot: evidence of lifetime increase, Solar Phys. 237, 329–350, 2006.

Charbonneau, P., Beaubien, G., and St-Jean, C.: Fluctuations in Babcok-Leighton dynamos: II. Revisiting the Gnevyshev-Ohl rule, Astrophys. J., 658, 657–662, 2007.

Choudhuri, A. R. and Karak, B. B.: Origin of Grand Minima in Sunspot Cycles, Phys. Rev. Lett., 109, 171103–171106, 2012.

Ding, R. and Li, J.: Nonlinear finite-time Lyapunov exponent and predictability, Phys. Lett. A 364, 396–400, 2007.

Duhau, S. and de Jager, C.: The Solar Dynamo and its Phase Transitions during the Last Millenium, Solar Phys. 250, 1–15, 2008.

Eckmann, J.-P. and Ruelle, D.: Ergodic theory of chaos and strange attractors, Rev. Modern Phys., 57, 617–656, 1985.

Fraser, A. M. and Swinney, H. L.: Independent coordinates for strange attractors from mutual information, Phys. Rev. A, 33, 1134–1140, 1986.

Greenkorn, R. A.: 2009, Analysis of Sunspot Activity Cycles, Solar Phys., 255, 301–323, doi:10.1007/s11207-009-9331-z, 2009.

Howe, R.: Solar interior rotation and its variation, Living Rev. Solar Phys., 6, 1–91, 2009.

Hoyt, D. V. and Schatten, K. H.: Group Sunspot Numbers: A new Solar Activity Reconstruction, Solar Phys., 181, 491–512, 1998.

Ivanov, E. V., Obrídko, V. N., and Shelting, B. D.: Quasi-biennial oscillations of the solar magnetic fields, in: Solar variability: from core to outer frontiers, The 10th European Solar Physics Meeting, 9–14 September 2002, Prague, Czech Republic, edited by: Wilson, A., ESA SP-506 v. 2, ESA Publications Division, Noordwijk, 847–850, 2002.

Kantz, H.: A robust method to estimate the maximal Lyapunov exponent of a time series, Phys. Lett. A, 185, 77–87, 1994.

Kitchatinov, L. L. and Olemskoi, S. V.: Active longitudes of the sun: The rotation period and statistical significance, Astron. Lett., 31, 280–284, 2005.

Kudela, K., Rybak, J., Antalova, A., and Storini, M.: Time evolution of low-frequency periodicities in cosmic ray intensity, Solar Phys. 205, 165–175, 2002.

Lawrence, J. K., Cadavid, A. C., and Ruzmaikin, A. A.: Turbulent and chaotic dynamics underlying solar magnetic variability, Astron. Astrophys., 455, 366–375, 1995.

Lawrence, J. K., Cadavid, A. C., and Ruzmaikin, A. A.: Rotational Quasi-Periodicities and the Sun-Heliosphere Connection, Solar Phys., 252, 179–193, 2008.

Le Mouël, J.-L., Shnirman, M. G., and Blanter, E. M.: The 27-Day Signal in Sunspot Number Series and the Solar Dynamo, Solar Phys., 246, 295–307, 2007.

Li, Q.-X. and Li, K.-J.: 2007, Low dimensional chaos from the group sunspot numbers, Chinese J. Astron. Astrophys., 7, 435–440, 2007.

Lockwood, M.: Long-term variations in the magnetic fields of the sun and the heliosphere: their origin, effects and implications, J. Geophys. Res., 106, 16021–16038, 2001.

Love, J. J. and Rigler, E. J.: Sunspot random walk and 22-year variation, Geophys. Res. Lett., 39, L10103–L10108, 2012.

Macek, W. M., Bruno, R., and Consolini, G.: Testing for multifractality of the slow solar wind, Adv. Space Res. 37, 461–466, 2006.

Mavromichalaki, H., Preka-Papadema, P., Petropoulos, B., Tsagouri, I., Georgakopoulos, S., and Polygiannakis, J.: Low- and high-frequency spectral behavior of cosmic-ray intensity for the period 1953–1996, Ann. Geophys., 21, 1681–1689, doi:10.5194/angeo-21-1681-2003, 2003.

Mayr, H. G. and Schatten, K. H.: Nonlinear oscillators in space physics, J. Atmos. Sol.-Terr. Phys., 74, 44–50, 2012.

McIntosh, P. S., Thompson, R. J., and Willock, E. C.: A 600-day periodicity in solar coronal holes, Nature, 360, 322–324, 1992.

Mursula, K., Zieger, B., and Vilppola, J. H.: Mid-term quasi-periodicities in geomagnetic activity during the last 15 solar cycles: connection to solar dynamo strength, Solar Phys., 212, 201–207, 2003.

Obrídko, V. N. and Shelting, B. D.: Occurrence of the 1.3-year periodicity in the large-scale solar magnetic field for 8 solar cycles, Adv. Space Res., 40, 1006–1014, 2007.

Oseledets, V. I.: A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems, Trans. Moscow Math. Soc., 19, 197–231, 1968.

Ostryakov, V. N. and Usoskin, I. G.: On the dimension of solar attractor, Solar Phys., 127, 405–412, 1990.

Pesnell, W. D.: Solar Cycle Predictions (Invited Review), Solar Phys., 281, 507–532, 2012.

Price, C. P., Prichard, D., and Hogenson, E. A.: Do the sunspot numbers form a chaotic set?, J. Geophys. Res., 97, 19113–19120, 1992.

Rosenstein, M. T., Collings, J. J., and De Luca, C. J.: A practical method for calculating largest Lyapunov exponents from small data sets, Physica D, 65, 117–134, 1993.

Rouillard, A. and Lockwood, M.: Oscillations in the open solar magnetic flux with a period of 1.68 years: imprint on galactic cosmic rays and implications for heliospheric shielding, Ann. Geophys., 22, 4381–4395, doi:10.5194/angeo-22-4381-2004, 2004.
Ruzmaikin, A., Feynman, J., Kosacheva, V.: On Long-Term Dynamics of the Solar Cycle, in: The solar cycle; Proceedings of the National Solar Observatory/Sacramento Peak 12th Summer Workshop 27, edited by: Harvey, K. L., ASP Conference Series, San Francisco, 547–556, 1992.

Sello, S.: Solar cycle forecasting: A nonlinear dynamics approach, Astron. Astrophys., 377, 312–320, 2001.

Shapoval, A., Le Mouël, J.-L., Courtillot, V., and Shnirman, M.: Two regimes in the regularity of sunspot numbers, Astrophys. J., 779, 108–116, 2013.

Shapoval, A., Le Mouël, J.-L., Courtillot, V., and Shnirman, M.: submitted, Is a sudden increase of irregularity of sunspot numbers a precursor of a return to low solar activity?, J. Geophys. Res., in review, 2014.

SIDC-team: World Data Center for the Sunspot Index, Royal Observatory of Belgium, Monthly Report on the International Sunspot Number, online catalogue of the sunspot index: http://www sidc.be/sunspot-data/, 1850–2005, 2005.

Spiegel, E. A. and Wolf, A.: Chaos and the Solar Cycle, in: Chaotic Phenomena in Astrophysics, vol. 497, edited by: Buchier, J.-R. and Eichhorn, H., Ann. N.Y. Acad. Sci., New York, 55–60, 1987.

Svalgaard, L.: Updating the Historical Sunspot Record, arXiv:1003.4666 [astro-ph.SR], 2010.

Svalgaard, L.: How well do we know the sunspot number? Comparative Magnetic Minima: Characterizing quiet times in the Sun and Stars, Proceedings of the International Astronomical Union, IAU Symposium 286, Mendoza, Argentina, 27–33, 2012.

Takens, F.: Detecting strange attractors in turbulence, edited by: Rand, D. A., and Young, L. S., Springer, Berlin, 366–381, 1981.

Valdes-Galicia, J. F., Perez-Enrizuez, R., and Otaola, J. A.: The cosmic ray 1.68 year variation: a clue to understand the nature of the solar cycle?, Solar Phys., 167, 409–417, 1996.

Vecchio, A., Laurenza, M., Carbone, V., and Storini, M.: Quasi-biennial modulation of solar neutrino flux and solar and galactic comtic rays by solar cyclic activity, Astrophys. J. Lett., 709, L1–L5, 2010.

Vecchio, A., Laurenza, M., Meduri, D., Carbone, V., and Storini, M.: The dynamics of the solar magnetic field: polarity reversals, butterfly diagrams, and quasi-biennial oscillations, Astrophys. J., 749, 27–36, doi:10.1088/0004-637X/749/1/27, 2012.

Wolf, A., Swift, J. B., Swinney, H. L., and Vastano, J. A.: Determining Lyapunov Exponent from a Time Series, Physica D, 16, 285–317, 1985.

Zhang, Q.: A nonlinear prediction of the smoothed monthly sunspot numbers, Astron. Astrophys. 310, 646–650, 1996.