Characterizing the Solar Cycle Variability Using Nonlinear Time Series Analysis at Different Amounts of Dynamo Supercriticality: Solar Dynamo is Not Highly Supercritical

The solar dynamo is essentially a cyclic process in which the toroidal component of the magnetic field is converted into the poloidal one and vice versa. This cyclic loop is disturbed by some nonlinear and stochastic processes mainly operating in the toroidal to poloidal part. Hence, the memory of the polar field decreases in every cycle. On the other hand, the dynamo efficiency and, thus, the supercriticality of the dynamo decreases with the Sun’s age. Previous studies have shown that the memory of the polar magnetic field decreases with the increase of supercriticality of the dynamo. In this study, we employ popular techniques of time series analysis, namely, compute Higuchi’s fractal dimension, Hurst exponent, and Multi-Fractal Detrended Fluctuation Analysis to the amplitude of the solar magnetic cycle obtained from dynamo models operating at near-critical and supercritical regimes. We show that the magnetic field in the near-critical regime is governed by strong memory, less stochasticity, intermittency, and breakdown of self-similarity. On the contrary, the magnetic field in the supercritical region has less memory, strong stochasticity, and shows a good amount of self-similarity. Finally, applying the same time series analysis techniques in the reconstructed sunspot data of 85 cycles and comparing their results with that from models, we conclude that the solar dynamo is possibly operating near the critical regime and not too much supercritical regime. Thus the Sun may not be too far from the critical dynamo transition.


INTRODUCTION
The solar cycle is not regular; the amplitude of the cycle has a considerable amount of variation.This is best seen in the observed sunspot number plot for the last 400 years and its proxy, such as the concentration of the cosmogenic isotopes ( 14 C and 10 Be) data for the last several thousands of years (Usoskin 2017;Biswas et al. 2023).Another prominent feature of the cycle irregularity is the Gnevyshev-Ohl/Even-Odd rule which is an alternating pattern of strong and weak cycles (Gnevyshev & Ohl 1948;Hathaway 2015).While in some cycles, the amplitude changed drastically from one cycle to the next (e.g., Cycles 6, 20), there is a reasonably smooth variation as seen by some envelopes in the amplitude, e.g., Gleissberg cycle (Gleissberg 1939).This indicates that although the solar cycle is irregular, there is some (temporal) memory in the underlying system.
We have good support that an αΩ type dynamo model, operating in the solar convection zone (SCZ), is responsible for causing the solar magnetic cycle (Karak et al. 2014;Cameron & Schüssler 2015;Charbonneau 2020).In this type of dynamo the poloidal field (which is observed to become strongest near the solar minimum) gives rise to the toroidal field and thus the sunspots for the next cycle.Hence, there is an unavoidable memory of about 5 years in the solar dynamo.This is indeed observed in the observations because there is a strong cor-relation between the polar field (or its proxy) at the solar minimum (or the peak field/proxy) and the amplitude of the next sunspot cycle (Jiang et al. 2007;Wang & Sheeley 2009;Kitchatinov & Olemskoy 2011;Priyal et al. 2014;Kumar et al. 2021b).In fact, this is true in any αΩ type dynamo model as long as the poloidal field gives rise to the toroidal field (Charbonneau & Barlet 2011).The polar precursor method of the solar cycle prediction is indeed based on this idea (Schatten et al. 1978;Choudhuri et al. 2007;Kumar et al. 2021b;Bhowmik & Nandy 2018;Kumar et al. 2022;Biswas et al. 2023).Now, is the memory of the polar field limited to the next cycle only, or is it propagated to multiple following cycles?At first glance, one may think that the memory should be propagated to multiple cycles as the solar dynamo is just the oscillation between the two components of the magnetic field: poloidal and toroidal.However, if we carefully analyze the dynamo chain, then we find that the toroidal to poloidal part of the solar dynamo is the one that involves some randomness arising due to the distributions of the BMR properties, primarily due to scatters around Joy's law and randomness in the BMR emergences (Jiang et al. 2014;Karak & Miesch 2017, 2018;Nagy et al. 2017;Karak 2020).Hence, in every cycle during the generation of the poloidal field, the memory of the polar field is degraded.The toroidal to poloidal component of the solar dynamo also involves some nonlinearities, which at least include the flux loss due to magnetic buoyancy in the formation of BMR (Biswas et al. 2022), latitude quenching (Jiang 2020;Karak 2020), and tilt quenching (Jha et al. 2020).The nonlinearity plays an important role in determining the memory of the polar field when the dynamo becomes supercritical.Kumar et al. (2021a) have shown that if the dynamo operates near the critical dynamo transition, then the dynamo tends to be linear and the polar field of a cycle strongly correlates to the toroidal field of multiple following cycles.On the other hand, this correlation is limited to one cycle if the dynamo operates in a highly supercritical region.They have made this conclusion by performing various types of dynamo simulations at different parameters, namely diffusivity, meridional circulation, and nonlinearity and in all models, they found that the correlation between the polar field and the toroidal field of the subsequent cycles is consistently shortened from multiple cycles to one cycle as the supercriticality of the dynamo is increased.They have further shown that this degradation of memory is independent of the amount of diffusion and the advection of the magnetic field, as suggested by Yeates et al. (2008).The diffusion and advection, of course, determine the memory within a cycle and a little bit beyond one cycle (Charbonneau & Dikpati 2000;Jiang et al. 2007).However, when these are kept constant, whether the memory of the polar field is propagated to multiple cycles or not is determined by the amount of the supercriticality.
While in the previous work (Kumar et al. 2021a), the memory was measured just by measuring the correlation between the peaks of the polar fields and the peaks of the following cycle toroidal fields.In the present work, we shall apply some well-known techniques of nonlinear time series analyses.In the time series analyses, two popular quantities, namely the Higuchi's dimension (D) and Hurst exponent (H) are generally used to determine the complexity of the system.D is a fractal dimension that determines the geometrical structure at multiple scales.On the other hand, H helps to identify the presence of long-term memory in a time series (Oliver & Ballester 1998;Sánchez Granero et al. 2008;Souza Corrêa et al. 2017;Das et al. 2022).
In the present study, we shall compute D and H of the cycles produced in the different regimes of the dynamo to measure the memory of the cycles.Using these measures, we shall independently demonstrate that the memory of the solar cycle beyond one cycle is indeed determined by the supercriticality of the solar dynamo.We shall also compare the results with those from the observed sunspot data and comment on the supercriticality of the solar dynamo.Knowing the amount of supercriticality will tell how far our Sun is from the critical dynamo transition (below which the dynamo action ceased).Answering this is essential because, with age, the Sun's dynamo efficiency decreases (as the rotation rate and thus the generation of poloidal field reduces), and at some point in its life, the Sun will stop producing its (large-scale) dynamo action.Incidentally, some studies hint that our Sun is not too far from the critical dynamo transition (Rengarajan 1984;Metcalfe et al. 2016;Kitchatinov & Nepomnyashchikh 2017;Vashishth et al. 2023).

MODEL DATA
We have used the data from three Babcock-Leighton type dynamo models, namely Models I, II, and Time Delay dynamo, operating them in the near critical and supercritical regimes.Models I and II are essentially flux transport dynamo models built using the Surya code (Nandy & Choudhuri 2002;Chatterjee et al. 2004;Choudhuri 2018)  η 0 = 2×10 12 cm 2 s −1 and for the meridional circulation: v 0 = 15 m s −1 ), and advection dominates in Model II (same diffusivity parameters but v 0 = 26 m s −1 ).Yeates et al. (2008) named these two models as diffusion and advection-dominated models, respectively.To operate these models in near-critical and supercritical regimes, we take α0 = α 0 /α crit 0 = 2 and 4, respectively (where α crit 0 is the minimum α 0 needed to obtain dynamo transition).For the time delay dynamo model, we follow the one presented in Wilmot-Smith et al. (2006).The time delay model is also of Babcock-Leighton type in which the delays involved in communicating the toroidal and poloidal fields to their respective source regions are captured by suitable time delays in their sources in the differential equations (Sec.6.4 of Karak 2023).The details of all the models are presented in Kumar et al. (2021a).For operating the time delay dynamo model in near critical and supercritical, we took α0 = 1 and 3, respectively.We note that the time delay model quickly goes to supercritical regime (e.g., see Figure 5 of Kumar et al. 2021a) and thus, to operate the model in near-critical and supercritical regimes, we have taken the value of α0 as 1 and 3 instead of 2 and 4 as taken in models I and II.From each model, we first compute the absolute value of the toroidal flux at low latitude at the base of the convection zone as a function of time (toroidal field in case of time delay model).Then, the peak values of the cycles of the toroidal flux/field are taken as the time series for our analyses.A representative example of the time series for Model I at critical regime is shown in Figure 1.The number of data points in each time series lies between 1800 and 2800.
Let us comment on the solar cycles from our models I and II.So far, no dynamo model is completely realistic (Charbonneau 2020;Karak 2023); our present ones are not the exception.Our dynamo models are axisymmetric and kinematic.The magnetic buoyancy and the Babcock-Leighton process are parameterised in simplified manners (Nandy & Choudhuri 2002;Biswas et al. 2022).However, the models still reproduce some basic features of the solar cycles, namely, the regular polarity reversals, a strong correlation between the polar field at the minimum and the amplitude of the next cycle toroidal field, a strong correlation between the rise rate and the amplitude of the cycle (Waldmeier effect), longterm modulation of the amplitude, grand minima and maxima.To highlight some of these features, in Figures 9 and 10 of Appendix, we show the variation of the toroidal field and the butterfly diagram from Model I at near-critical and supercritical regimes.We can already see from these figures that the long-term modulation of the cycle in the two regimes is different.The near-critical model produces a strong long-term modulation and frequent grand minima-like events.Our aim is to quantify this change of solar cycle memory more rigorously and to identify at what parameter regime of dynamo supercriticality, the model reproduces the observation best.

METHODS
Nonlinear time series analysis techniques have been employed in various diverse fields, starting from the prediction in the stock market to understand the dynamics of various complex systems (Cervantes et al. 2013;Salcedo-Sanz et al. 2022).In astronomical and solar data, these have been used to identify the existence of low-dimensional chaotic or stochastic signature in the underlying processes (Mundt et al. 1991;Carbonell et al. 1994;Qin 1998;Hanslmeier & Brajša 2010;Karak et al. 2010), and the persistence of the memory in the solar cycle (Maddanu & Proietti 2022).Recently, Das et al. (2022) identified some memory in the solar cycle asymmetry data using different nonlinear dynamics parameters like correlation dimension and fractal dimension techniques.In the present work, we will exploit the nonlinear time series analysis to investigate the memory of the solar cycles in the dynamo models.

Higuchi's dimension
Fractal dimension D, used to characterize nonlinear time series, is computed using the method given in Higuchi (1988).Here we only briefly discuss it.
We start with a time series X(i) containing N observations that have been sampled at regular intervals.Thus, X(i) : X(1), X(2), X(3), ..., X(N ). (1) From X(i), we construct a new time series where, 1 ≤ m ≤ s and n = N −m s denotes the greatest integer less than or equal to N −m s .Associated with each X m s , we calculate the length of the curve L m (s) as where k = (N − 1)/ns.The average length of the curve ⟨L(s)⟩ is obtained by taking the mean of L m (s) for 1 ≤ m ≤ s.For s min < s < s max , if we obtain ⟨L(s)⟩ ∝ s −D , the time series is a fractal of dimension D in that range of s.We compute ⟨L(s)⟩ for 2 ≤ s ≤ 32 and 2 ≤ s ≤ 256 for the data obtained from critical and supercritical dynamo regimes, respectively.D is computed from the slope of the double logarithmic plot of ⟨L(s)⟩ vs s.
The value of D is a fraction with 1 < D < 2. A value close to 2 denotes a space-filling curve, while a value close to 1 is a straight line.

Hurst Exponent
We follow the R/S method introduced by Mandelbrot & Wallis (1969) to find the Hurst exponent H.
We have a time series, X(i), i = 1, 2, ..., N , whose Hurst exponent we want to compute.Now, choose a temporal window s, with s t < s < N .Here, s t is the Theiler window.Now, we make (N − s + 1) subsets of the series X(i) as follows: x t0 (s) : X(t 0 ), X(t 0 +1), X(t 0 +2), ..., X(t 0 +s−1), ( 4) where, t 0 = 1, 2, ..., N − s + 1.Now, average of the subset x t0 (s) is given by Now, the standard deviation of x t0 (s) corresponding to window s is given by We define the set of cumulative deviations of x t0 (s) from the mean as where, i = 1, 2, ..., s.Thus, the Range of y i (t 0 , s) is and the rescaled range measure R/S is given by For a particular s, we have 1 ≤ t 0 ≤ N − s + 1 and thus, the corresponding rescaled range is This value of the rescaled range is found to vary as where, k is a constant and H is the Hurst exponent.A plot is made for log(R/S) vs log(s) for 1 ≤ s ≤ N and the linear portion of the graph is fitted to obtain the Hurst exponent.A time series obtained from a white noise process yields H = 0.5.When H > 0.5 for a time series, the series is said to be persistent or has a longterm memory.An increase in the value of the time series at a particular step is more probable to be followed by its increase in the next step (rather than a decrease).When H < 0.5 for a time series, it is said to be antipersistent or has a short-term memory.An increase in the value of the time series at a particular step is more probable to be followed by its decrease in the next step.

Multifractal Analysis
Many nonlinear systems, including solar activity, are characterized by intermittent phenomena and in such scenarios, a single Hurst exponent is not sufficient to capture the essential characteristics of the system.The Multifractal Detrended Fluctuation Analysis (MF-DFA) helps reveal the complexity or multifractal structure of the time series by characterizing the amplitude fluctuations in the data.The MF-DFA method has been developed from Detrended Fluctuation Analysis (DFA) by Kantelhardt et al. (2002).This method has been applied to understand the underlying dynamics that lead to the flux variability of quasars (Belete et al. 2018) and the nature of solar flare activity (McAteer et al. 2007;Sen 2007).
We have a time series, X(i), i = 1, 2, ..., N , and we want to understand its multifractal structure.Let ⟨X⟩ be the mean of the time series X(i).We compute the profile as Next, we choose a segment size s and divide the new series Y (i) into N s = [N/s] segments, where [k] denotes the integer less than or equal to any real number k.
Starting from the first element of the time series Y (i), we obtain N s segments, each of size s and will be left with N (mod s) elements at the last, which do not belong to any segment.To account for the remaining part of the series, we repeat the same process by starting from the end of the time series, obtaining N s segments, and leaving N (mod s) elements at the beginning.Hence, in total, we obtain 2N s segments.For each of these 2N s segments of length s, we compute the first-order polynomial fit y ν (i), 1 ≤ i ≤ s using the least square method.Next, the variance for each of the segments is computed as for ν = 1, 2, ..., N s , and as We now compute the average variance of each of the segments as , q ̸ = 0 (15) We are interested in knowing how the average variance F q (s) scales with the size of the segments s for different values of q.Suppose F q (s) ∼ s h (q)  (17) and we are interested to find the value of the exponent h(q) known as the generalized Hurst exponent.For very small segment size s ≈ 10, there are systematic deviations from the scaling behavior, while for large segment size s > N/4, the method is statistically unreliable (Kantelhardt et al. 2002;Ihlen 2012).Thus, we have chosen the segment size s, such that 15 ≤ s ≤ N/4.For q = 2, the method corresponds to the method of Detrended Fluctuation analysis (DFA) and the generalized Hurst exponent h(q = 2) corresponds to the Hurst exponent H.
For any monofractal time series, the generalized Hurst exponent is independent of q (or varies weakly with q).While a strong q dependence of h(q), indicates the time series is multifractal.Moreover, if a homogeneous scaling behavior of F q (s) (in Equation ( 17)) is obtained, the h(q) values for q < 0 are usually larger than the values when q > 0.
This generalized Hurst exponent h(q) is related to the scaling exponent τ (q) of the standard partition function based multifractal formalism, by Difference between the slope of τ (q) at q < 0 and q > 0 indicates the strength of multifractality in the data.Another way to see the multifractal structure of the data is to use the singularity spectrum given by where, α = τ ′ = dτ dq .Here, α denotes the singularity strength and f (α) denotes the dimension of the subset of the time series that is characterized by α.
Once we have the multifractal spectrum, we can note down the minimum (maximum) value of the singularity strength α min (α max ).The width of the multifractal spectrum is given by ∆α = α max − α min .A broad spectrum, indicated by large values of ∆α, represents a higher degree of multifractality in the data (Ihlen 2012).In the monofractal limit, the spectrum reduces to a single point and ∆α goes to zero.However, for real-world signals which always have a finite length, the multifractal spectrum always has a small (non-zero) width.
The value of α at which f (α) assumes a maximum value is denoted by α 0 .The symmetry in the shape of the spectrum is given by the quantity A = (α max − α 0 )/(α 0 − α min ).For a right-skewed (left-skewed) spectrum, we have A > 1 (A < 1).When A = 1, the multifractal spectrum is symmetric.A right-skewed (leftskewed) multifractal spectrum arises due to left (right) truncation.A left (right) truncation of the multifractal spectrum occurs due to the leveling of the q th -order Hurst exponents for q > 0 (q < 0) (Ihlen 2012).A left (right) truncation indicates that the multifractal structure is sensitive to small-scale temporal fluctuations with small (large) amplitudes.Moreover, a left (right) truncation shows a higher abundance of small (large) amplitude fluctuations in the time series.

From model data
To explore how the nature of the magnetic cycle changes with the supercriticality of the solar dynamo, we shall now look into the results from the nonlinear time series analysis methods as discussed above.
For Model I, we find D = 1.657 and D = 1.956 for the time series in the near-critical (Figure 2a) and supercritical (Figure 2b) regimes, respectively.Again, for Model II, the time series in the near-critical and supercritical regimes produce D = 1.501 (Figure 2c) and D = 1.972 (Figure 2d), respectively.Finally, for the Delay Dynamo Model, we obtain D = 1.616 and D = 1.964 when the dynamo operates in the critical (Figure 2e) and supercritical (Figure 2f) regime, respectively.
A periodic time series will be fractal dimension D = 1, and a highly stochastic time series will be fractal dimension D = 2 (Bhatt et al. 2017).In the critical regime, a value of the fractal dimension D close to 1.5 shows that the system may be chaotic and irregular.However, in the supercritical regime of the dynamo, the fractal dimension is close to 2, and thus, the process must be stochastic.
Moreover, we note that when the dynamo operates near the critical regime, the value of the fractal dimension is obtained at a timescales of about 2-32 solar cycles.As we look into larger timescales (beyond 32 solar cycles), the slope changes, which indicates that some different fractal dimension is needed to characterize the time series at these larger timescales.These are typical of structures that are not self-similar.However, when the dynamo operates in the supercritical regime, a single fractal dimension is obtained in a wide range of timescales, ranging from 2 to 256 solar cycles.This is indicative of the self-similar nature of the time series obtained in the supercritical regime but not in the critical regime.
A time series having a fractal dimension between 1.0 to 2.0 will have a memory effect (Souza Corrêa et al. 2017).The range of fractal dimensions for the time series are varying from 1.50 to 1.97 as we go from nearcritical to supercritical regimes of the dynamo.It gives a clue about the persistence of memory.Finally, we computed the Hurst exponent to look for the persistence in the solar cycle memory.
For Model I, we see that the Hurst exponent in the near-critical regime is 0.963 for the range s = 11 − 50 and in the supercritical regime is 0.747 for the range s = 11 − 25 (Figure 3a-b) which is > 0.5, i.e., irregular dynamo time series contains memory.We also see that in the near-critical regime, the value of the Hurst exponent is greater than in the supercritical regime.Moreover, the number of solar cycles over which the memory persists is greater when the dynamo operates in the near-critical regime (11-50) compared to when it operates in the supercritical regime (11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25).
For Model II, the Hurst exponent in the near-critical regime is 0.999 for the range s = 11 − 80 and in the supercritical regime is 0.745 for the range s = 11 − 25 (Figure 3c-d Thus, all the dynamo models give similar results for the Hurst exponents.We observe that a dynamo operating in the near-critical regime has a stronger memory that persists over a higher number of solar cycles when compared to the memory in a dynamo operating in the supercritical regime.This supports our previous conclusion made based on the correlations between the poloidal and toroidal fields (Kumar et al. 2021a).
The results from the MF-DFA method for Model I are shown in Figure 4.The generalized Hurst exponent depends on the value of q, which shows that both the time series have some degree of multifractality.However, for the near-critical regime, the dependence of h(q) on q is stronger than in the super-critical case, implying the multifractal structure is pronounced in the near-critical regime.For q = 2, the generalized Hurst exponent h(q) = 0.890 and 0.603 respectively in the near-critical and supercritical regimes.This shows a large persistence of the time series in the near-critical regime compared to the supercritical regime.The result is qualitatively similar to what has been obtained in the R/S method.But, the values of the Hurst exponent may be slightly smaller since we are looking at much larger timescales in the MF-DFA method.The magnitude of the change in slope of the graph of τ (q) versus q at q = 0 indicates the multifractal structure in the time series.The change in slope in the near-critical (0.652) and supercritical (0.372) regime can also be seen from the jump in the value of dτ dq around q = 0 in Figure 4c.Finally, we obtain the multifractal spectrum (f (α) vs α) with left (right) truncation when the dynamo operates in the near-critical (supercritical) regime.The widths of the spectrum (∆α) in the near-critical (0.770) and supercritical (0.521) regimes also tell us that the multifractality is higher in the near-critical regime when compared to the supercritical regime.
Figure 5 shows the results from the MF-DFA method when applied to the time series obtained from Model II.Again, we notice a stronger dependence of the generalized Hurst exponent h(q) on q in the near-critical regime.For q = 2, the generalized Hurst exponent h(q) = 1.24 in the near-critical regime and h(q) = 0.659 in the supercritical regime.The value of h(q = 2) > 1 implies that the non-stationarity in the data could not be successfully removed by detrending in the near-critical regime (Bryce & Sprague 2012;Ceballos & Largo 2018).However, the value of h(q) in the supercritical regime suggests little or no memory in this regime.The change in slope of τ (q) vs q around q = 0 is higher in the near-critical regime (0.987) compared to the supercritical regime (0.317).The multifractal spectrum (f (α) vs α) is again obtained with left (right) truncation when the dynamo operates in the near-critical (supercritical) regime.The width of the spectrum (∆α) also tells us that the multifractality is higher in the near-critical regime (1.092) when compared to the supercritical regime (0.432). Figure 6 shows the results of the MF-DFA method from the Delay Dynamo Model.As obtained earlier, the generalized Hurst exponent h(q) strongly depends on q in the near-critical regime, implying a higher degree of multifractality.We also find that h(q = 2) = 1.07 in the near-critical regime and h(q = 2) = 0.552 in the supercritical regime.Again, a value of h(q = 2) greater than unity suggests that the non-stationarity or trend in the data from the near-critical regime could not be successfully removed.However, a value slightly larger than unity may be interpreted as high persistence in the data, similar to that obtained by the R/S method (0.983).In the supercritical regime, a value close to 0.5 indicates little or no memory.The slope of τ (q) vs q changes more sharply around q = 0 in the near-critical regime (0.699) compared to the supercritical regime (0.096); also, see Table 1 for the key parameters of MF-DFA methods from all models for comparison.The width of the spectrum (∆α) confirms that the multifractality is higher in the near-critical regime (0.874) when compared to the supercritical regime (0.166).The multifractal spectrum is obtained with a left truncation in the near-critical regime, whereas it is mostly symmetric in the supercritical regime.
The unequal abundance of the small and large amplitude variations gives rise to the multifractal structures in the magnetic field data obtained from the solar dynamo models.A higher multifractality in the nearcritical regime is due to the relatively larger inequality in the abundance of small and large amplitude variations in this regime compared to the supercritical case.Moreover, in the near-critical regime, the multifractal spectrum is obtained with left truncation, implying the multifractal nature of the time series arises due to the higher abundance of local fluctuations with small-scale amplitude variations.In the supercritical regime, the width of multifractal spectrum is narrow and is obtained with occasional right truncation.This suggests that large amplitude fluctuations may be more common but not by a large margin compared to the small amplitude fluctuations.Thus, the time series in this regime exhibits close to monofractal behavior.

From reconstructed sunspot number
Regular homogeneous data of sunspot number are available only for about the last 30 solar cycles (Usoskin 2023), which is too small to study the long-term behavior of the solar activity.However, recently, the annual solar activity series of the last millennium has been reconstructed from 14 C data by Usoskin et al. (2021).They have also isolated the individual solar cycles and the corresponding cycle-averaged sunspot number, ⟨SN⟩, during the period from 971 to 1899.This ⟨SN⟩ is a measure of the cycle strength.The data was found to contain 85 solar cycles.We construct a time series of these 85 data of ⟨SN⟩ and study the nonlinear characteristics of this time series.We note that as the sunspot data we use here is "reconstructed", which involves a sequence of model steps, its quality is compromised.However, as we are using the cycle-averaged SN in our analyses, the effect of noise is somewhat reduced.Furthermore, Usoskin (2023) showed that this reconstructed data closely resembles the available direct observations of solar cycles and reproduces the popular Waldmeier effect (Waldmeier 1935;Karak & Choudhuri 2011).Given these, we hope that this indirectly observed data carries some genuine features of the long-term behaviour of solar activity in the past, which can be used to compare  with the results from dynamo models operating in two different regimes to get a hint of the operation of the solar dynamo.Table 2 tabulates the results of Higuchi's dimension and Hurst exponent applied to the ⟨SN⟩ time series.The Higuchi's dimension and Hurst R/S for the ⟨SN⟩ time series are obtained to be 1.737 and 0.857; see Figure 7. Due to the restricted size of the time series, we were limited to a smaller window range for the Hurst exponent compared to the results obtained from the Model data.To account for the error in the ⟨SN⟩ time series, we compute 100 resampled data sets with the same size as the ⟨SN⟩ time series.An ensemble of 100 data points from a Gaussian distribution is produced, considering the ⟨SN⟩ as the mean and the corresponding errors as the standard deviation.Finally, we compute the Higuchi's dimension and Hurst exponent for the 100 resampled time series and report their means and standard deviations in Table 2.The Higuchi's dimension and the Hurst exponent are found to be 1.749 ± 0.015 and 0.843 ± 0.062, respectively.This suggests that the data is not completely stochastic but may be irregular and chaotic.In addition, there is some memory in the data.However, the results do not precisely locate the mode of operation of the solar dynamo-it only implies that the solar dynamo is operating somewhere between the near-critical and highly supercritical.
Figure 8 shows the multifractal characteristic of the ⟨SN⟩ time series.The time window (segment size s in Section 3) is chosen to be 11-21 cycles.Though systematic deviations may arise at small time windows (11-14 cycles) in the MF-DFA method, we checked there were no qualitative differences in the results when we chose the time window to be 15-21 cycles.The generalized Hurst exponent h(q) strongly depends on q and changes from 2.6 to 0.6 as q varies from −3 to 3. The variation is mostly when q < 0. This shows a very high degree of multifractality in the data.At q = 2, the Hurst exponent is 0.734, indicating there is some persistence in the time series.The slope of the τ (q) vs q graph changes around q = 0 by a large amount (2.19).This change in the slope is also seen in Figure 8c.The multifractal spectrum (f (α) vs α) is obtained with a left truncation and the width of the spectrum (∆α) is quite large (3.010).The spectrum is also found with a left truncation (A = 6.167).We repeat the MF-DFA analysis on the 100 resampled data sets.The spectrum does not follow the expected inverted parabolic shape for 13 data sets and is found to produce f (α) > 1.This may arise when the homogeneous scaling, as assumed in Equation ( 17), does not hold in the considered time window.From the remaining data sets, we find that the width of the multifractal spectrum is ∆α = 2.609±0.633and the asymmetry parameter A = 5.480 ± 1.989.
The large width of the spectrum suggests the need for a large number of exponents to characterize the data and, thus, the breakdown of self-similarity.The left truncation of the multifractal spectrum suggests the abundance of local fluctuations with small-scale amplitude variations.We can also say that the multifractal structure of the data is sensitive to small-scale temporal variations with small-amplitude.The sensitivity of the multifractal structure to the local fluctuations with small amplitude variations is also seen for the time series obtained from the dynamo models when they operate near the near-critical regime (left truncation of the f (α) vs α curve in Figures 4d, 5d and 6d).This suggests that the Sun operates close to the near-critical regime and not in the highly supercritical one.

DISCUSSION AND CONCLUSION
Using nonlinear time series analysis techniques, we have studied the behaviour of long-term variation of the magnetic cycles in near-critical and supercritical regimes of the solar dynamo.For this, we have considered the peak values of the toroidal flux from three different dynamo models.We find that Higuchi's fractal dimension (D) is close to 1.5 and 2, respectively, for the nearcritical and supercritical regimes of the dynamo.On the other hand, the Hurst exponent (H) is near 1 and 0.74 for these two regimes.Values of D and H suggest that the magnetic cycle in the supercritical regime is governed by less memory and a more stochastic process.In other words, when dynamo operates in the supercritical regime, the persistence nature of the cycle is weak, and the stochastic nature is dominant.By observing the long window in H (or in the scaling relation of R/S vs τ as seen in Figure 3 and Table 1), we conclude that the memory of the magnetic cycle in the near-critical regime is very long.This result is congruous with the previous expectation based on the linear correlation between the peaks of the polar field and that of the subsequent cycles (Kumar et al. 2021a).Moreover, we see evidence of selfsimilarity in the time series obtained in the supercritical regime.However, this self-similarity breaks down when we approach the near-critical regime, as we see a single value of the fractal dimension is not good enough for a range of timescales.
The Multi-Fractal Detrended Fluctuation Analysis (MF-DFA) tells us that the breakdown of self-similarity (or intermittency) in the near-critical regime arises due to the large abundance of small-amplitude fluctuations in the time series, with occasional large-amplitude variations.However, the time series in the supercritical regime has little or no disparity in the abundance of small and large amplitude variations.This gives rise to a weak multifractal or close to monofractal behavior in the supercritical regime.
We can understand the above results from the dynamo model in the following way.In the near-critical regime, the effect of nonlinearity is weak, and the dynamo number and, thus, the growth rate is small.In this regime, the dynamo takes a long time to grow the magnetic cycle if it falls to a low value.This produces long-term modulation and intermittent behaviour of the magnetic cycle.On the other hand, in the supercritical regime, the dynamo number and, thus, the growth rate and the nonlinearity are high.These cause the magnetic field to grow rapidly (due to the high growth rate) when it has fallen to a low value or decrease rapidly (due to strong nonlinearity) when the field has enhanced to a high value.These break the long-term modulation, making the cycle more stochastic and tends to give a self-similar pattern.
Further, we compute the nature of the average solar cycle data obtained from annual solar activity series from the last millennium that has been reconstructed from 14 C data (Usoskin et al. 2021).Given the considerable uncertainty in reconstructed data and the limited number of cycles, the computed values of D and H suggest that the solar dynamo is possibly not operating in a highly supercritical regime.The multifractal analyses show a lack of self-similarity and a high multifractal nature in the time series, arising due to the abundance of local fluctuations with small amplitude variation.These suggest that the solar activity is due to an underlying dynamo process that operates in the near-critical regime or, at least not, in the highly supercritical one.This conclusion is in agreement with the previous independent investigations (Rengarajan 1984;Metcalfe et al. 2016;Kitchatinov & Nepomnyashchikh 2017;Vashishth et al. 2023).In this section, we briefly highlight the salient features of the solar cycles obtained from Models I at near critical (Figure 9) and super-critical regimes (Figure 10).As seen from these two figures, the model shows irregular cycles with long-term modulation of the amplitude (top panels), solar-like oscillation of the magnetic field (middle panels), and the regular reversal of the polarity and the equatorward migration of the toroidal field at the base of the SCZ (bottom panels).We note that the cycle period is quite long (∼ 22 years) because of the weak meridional flow in this case.However, model II produces a reasonable cycle duration of about 12 years.

Figure 1 .
Figure 1.Time series of the peak values of the cycles of the toroidal flux obtained from Model I, operating in the (a) nearcritical and (b) supercritical regimes.

Figure 2 .
Figure 2. Variation of length ⟨L(s)⟩ with time interval (s) for toroidal flux data obtained from dynamo models in critical (left) and supercritical (right) regimes.

Figure 3 .
Figure 3. Variation of R/S with time interval (s) for toroidal flux data obtained from different models in the nearcritical (left) and supercritical (right) regimes.

Figure 4 .
Figure 4. Results from the MF-DFA method obtained in the near-critical and supercritical regimes of Model I.

Figure 5 .
Figure 5. Same as Figure 4 but for Model II.

Figure 6 .
Figure 6.Same as Figure 4 but for the delay dynamo model.

Figure 7 .
Figure 7. Variation of length ⟨L(s)⟩ (left) and R/S (right) with time interval (s) for the reconstructed SN during the last millennium (Usoskin et al. 2021).

Figure 8 .
Figure 8. Results from the MF-DFA method obtained from the reconstructed SN.

Figure 9 .
Figure 9. Results from Model I operating at near-critical regime.(a) The temporal variation of the toroidal flux (in code unit) obtained within a radial extent of 0.67R⊙ to 0.72R⊙ and latitude extent of 15 • to 45 • .(b) Highlighting the cycles marked by two vertical black dashed lines in (a).(c) The time-latitude distribution (butterfly diagram) of the toroidal field at 0.72R⊙ during the time interval marked in (b) by vertical black lines.

APPENDIXA.
ADDITIONAL INFORMATION FOR THE MODEL SOLAR CYCLES

Figure 10 .
Figure 10.Same as Figure 9, but from the super-critical regime.

Table 1 .
Summary of the results of Higuchi's dimension, Hurst exponent and the MF-DFA method for different model time series in near-critical and supercritical regimes.

Table 2 .
Summary of results obtained for Higuchi's dimension, Hurst exponent and the MF-DFA method from the reconstructed sunspot number.