© The Authors 2025

1. Introduction

The sunspot cycle – characterized by an approximately 11-year quasi-periodic rise and fall in solar activity – is a striking manifestation of the Sun’s magnetic behavior (Charbonneau 2020; Usoskin 2023). However, the strength of the Sun’s activity cycle is not uniform, and varies from one cycle to another, resulting in a variable forcing of the heliosphere that seamlessly bridges the physical phenomena originating in our host star’s interior to their manifestation on the solar surface, outer atmosphere, and planetary impacts (Nandy et al. 2018; Daglis et al. 2021; Dash et al. 2023; Nandy et al. 2023). Observational evidence has uncovered that over the past century, the Sun has exhibited a prolonged period of unusually high activity known as the Modern Maximum (Solanki et al. 2004). The strongest cycle during this phase – in fact, over the past four centuries since direct sunspot cycle observations began with Galileo Galilei and his contemporaries – was sunspot cycle 19 peaking around 1957; this was followed by an unexpectedly weak sunspot cycle 20. Understanding and reconstructing this extreme solar activity fluctuation over cycles 19-20 through physics-based models has remained an outstanding challenge (Cameron et al. 2010; Jiang et al. 2013; Bhowmik & Nandy 2018; Bhowmik 2019; Virtanen et al. 2022; Pal & Nandy 2024; Yeates et al. 2025).

It is understood that nonlinearities inherent in the magnetohydrodynamic solar dynamo mechanism and stochastic perturbations in the dynamo source terms can lead to amplitude fluctuations from one solar cycle to another (Charbonneau & Dikpati 2000; Saha et al. 2025). The Sun’s large-scale dipolar field–the poloidal component (of which the polar field is a proxy)–is inducted by solar differential rotation within the Sun’s convection zone to produce the toroidal magnetic field component of the following sunspot cycle (Parker 1955). Magnetic (Lorentz) feedback of strong toroidal flux tubes on the Sun’s differential rotation can, in principle, result in amplitude modulation; however, observations show that inter-cycle variations (known as torsional oscillations) in the solar differential rotation which acts as the source of the Sun’s toroidal field – is quite small (≤5%) (Mahajan et al. 2024).

Strong toroidal fields are unstable to magnetic buoyancy and rise up to emerge through the solar surface, giving rise to bipolar sunspot pairs, which are observed to be systematically tilted relative to the local latitude (due to the action of the Coriolis force on rising flux tubes). Observations (Dasi-Espuig et al. 2010; Muñoz-Jaramillo et al. 2012), theoretical considerations (Cameron & Schüssler 2015), as well as data-driven physical models of the long-term evolution of solar magnetic fields–such as solar surface flux transport (SFT) models and dynamo models (Bhowmik & Nandy 2018)–have suggested that the dispersal of the magnetic flux of these tilted bipolar sunspots (mediated by plasma flows) is the primary mechanism for recreation of the Sun’s large-scale dipolar field; the latter mechanism is known as the Babcock-Leighton (BL) mechanism (Babcock 1961; Leighton 1969). As magnetic flux tubes rise through the solar convection zone, they are buffeted by vigorous turbulence, generating an observed scatter in the tilt angles around the mean (Joy’s law tilt) expected from the Coriolis force (Cheung & Isobe 2014). The amplitude of this scatter is observed to be much more than the mean tilt angle, which acts as a significant source of perturbation on the poloidal field source (Jiang et al. 2014; Nagy et al. 2017). Anti-Hale active regions (ARs) – bipolar sunspot pairs whose polarity orientation does not conform to the conventional solar cycle trend – are an additional source of perturbation (Nagy et al. 2017; Pal et al. 2023). Other mechanisms, such as tilt quenching and latitudinal quenching, are also known to act as amplitude-limiting mechanisms for the BL poloidal field generation mechanism (Jiao et al. 2021; Yeates et al. 2025).

Although the interplay of these nonlinearities and stochastic forcing is theorized to drive solar activity fluctuations, capturing the extreme variation between cycles 19 and 20 has remained elusive. Here we address the outstanding question of whether solar cycles 19 and 20 can be accurately reconstructed using known physics and physical models of solar magnetic field evolution.

2. Results and discussions

We first reconstructed century-scale solar polar field variations using an SFT model. We prepared the synthetic input of ARs using observed sunspot emergence statistics derived from the RGO/USAF/NOAA Data Centre (1874-2023) database–which provides details on the area, emergence time, and location of sunspots. This information was used to drive the SPhoTraM model (a brief description of which is available in Appendix A.1). Past work has shown that long-term simulations using this database fail to reproduce variations across all solar cycles, including extreme solar cycles (Bhowmik & Nandy 2018; Pal & Nandy 2024). This limitation arises because this database lacks two crucial pieces of information necessary for precise modeling. The first is the unavailability of the tilt angle information of ARs, which significantly influences polar field generation. Although models often assume that tilt angles follow the mean Joy’s law tilt expected from Coriolis force, scatter around the mean tilt is crucial for capturing polar field variations (Jiang et al. 2014; Nagy et al. 2017). The second limitation is the lack of observational constraints on the polarity orientation of sunspots, which results in an inability to account for the impact of anomalous ARs (anti-Hale configurations). These anomalous sunspots can significantly reduce the polar field and the overall solar cycle strength (Pal et al. 2023). To overcome these limitations, we incorporate these two properties of ARs into our SPhoTraM simulations.

Our simulation starts in 1902, at the beginning of solar cycle 14, initializing it with a dipolar magnetic field configuration. First, we generated the ensembles of sunspot input sources for a single solar cycle. Each ensemble introduced random scatter into the standard empirically calculated tilt angles (for technical details, see Appendix A.1). Additionally, we distributed approximately 8%–8.4% of the total ARs of that cycle as anti-Hale regions–which we note is consistent with the observed range gleaned from the most recent solar cycles for which polarity information exists (Li 2018). Accounting for tilt scatter and anti-Hale active region distribution in this manner, we generated multiple input sources for solar cycles 14 to 23 and performed Monte Carlo (MC) SPhoTraM simulations. From these ensembles, we chose the optimal solution for each solar cycle, which minimized the difference between the simulated and observed polar field at cycle minima (for detailed methodology of optimizing, see Appendix A.1).

Fig. 1. Simulated and observed large-scale magnetic fields of the Sun. Top panel: Variations in the optimized simulated polar flux (solid curves) alongside the observed polar flux obtained from the MWO polar faculae database (dashed curves) for solar cycles 15–23. The red and blue curves represent the northern and southern hemispheres, respectively. Additionally, the black and gray curves in the same panel depict the total yearly and monthly averaged unsigned sunspot flux, derived from the RGO/USAF/NOAA database. Bottom panel: Time-latitudinal distribution of the radial magnetic field (Br) based on the optimized simulation. Here, yellow and blue shades respectively indicate magnetic fields of positive and negative polarity.

Fig. 2. Simulated solar cycle reconstruction versus observations. The violet curve represents the total yearly averaged unsigned sunspot flux from solar cycle 16 to 24, derived from the USAF/RGO/NOAA database. The orange curve indicates the total unsigned flux simulated using the dynamo model (refer A.2), driven by the poloidal field generated by the SPhoTraM model. A Pearson correlation coefficient of 0.93, calculated with a 99% confidence level, highlights the strong agreement between the observed and simulated flux strengths at cycle maxima.

Figure 1 (top panel) shows the optimized simulated polar flux variations for the northern and southern hemispheres (solid red and blue, respectively) compared with the observed polar flux (dashed curves). This figure demonstrates that by accounting for tilt angle scatter and anti-Hale sunspots, we are able to successfully reproduce the polar field evolution of solar cycles 18 and 19 (which serve as sources for sunspot cycles 19 and 20, respectively), along with other cycles over a century-scale. The butterfly diagram (Figure 1 bottom panel) depicts the time-latitude evolution of the longitudinally averaged radial magnetic field (B_r). The diagram shows the surface field dynamics associated with polar field reversal and build-up in our optimized simulation. Additionally, the magnetic flux surges toward the poles, representing both positive and negative polarity contributions, including perturbations due to anomalous ARs, are clearly discernible. These surges play a crucial role in determining polar field amplitude.

In the next stage, we employed a two-dimensional (axisymmetric) kinematic dynamo model (for model details see Appendix A.2) to reconstruct the solar cycle time series over a century covering the Modern Maximum phase. In this approach, we assimilated the longitudinally averaged surface magnetic field from the optimized SFT simulation as the source of the poloidal field into the dynamo model at the end of each solar cycle (see Appendix A.2 for a more detailed methodology). Our solar dynamo model integrates a buoyancy algorithm to simulate sunspot emergence as eruptions of the toroidal field when it exceeds a specified threshold magnetic field strength. These toroidal field eruptions were then used as a proxy for the total sunspot flux that has erupted during the cycle. The simulated magnetic flux from these sunspot eruptions is compared with the observed unsigned sunspot flux, which represents the total unsigned flux derived from sunspot emergence data in the RGO/USAF/NOAA database. The result is depicted in Figure 2. Our results show that the SPhoTraM-generated poloidal field successfully reproduces the cycle strength of solar cycles 19 and 20, while capturing the overall trend of other cycles over the century scale. This outcome highlights the importance of properly accounting for tilt angle scatter and anomalous regions in governing extreme variations in solar activity.

To deconstruct the physics of extreme solar variability over sunspot cycles 19-20, we studied the mean flux-weighted tilt angle distribution of sunspots in the simulation runs that successfully reproduce the polar flux for solar cycles 18 and 19 (which subsequently act as seeds of sunspot cycles 19 and 20, respectively). If no stochasticity is included in tilt-angle variations, the values of flux-weighted tilt coefficient, i.e., Φ>T_A remain nearly identical for cycles 18 and 19, as shown in the top panel of Figure 3. However, upon incorporating tilt-angle fluctuations, Φ>T_A increases substantially for cycle 18, while it decreases for cycle 19. Moreover, the correlation between Φ>T_A of the nth solar cycle and the observed average polar flux at the end of the same cycle improves significantly when stochasticity is incorporated into the sunspot emergence statistics (compare the magenta squares in the top and bottom panels of Figure 3). Similarly, we observe a strong positive correlation (r = 0.88, with 99.9% confidence) between Φ>T_A of the nth cycle and the dipole moment at the end of the same cycle, as generated by the optimized SFT simulation (see the red circles in the bottom panel of Figure 3). This result suggests that solar cycle 18, with its higher Φ>T_A, leads to a stronger dipole moment at the end of the cycle, which seeds the extreme sunspot cycle 19. This is because when the tilt angles of sunspot pairs are high, the latitudinal separation between opposite polarities is more efficiently amplified by plasma flows, avoiding intra-active region flux cancellation and allowing transport of higher amount of flux efficiently toward the poles. A stronger polar field enhances the seed poloidal field for the following cycle, which explains the exceptionally strong solar cycle 19. The highest value of Φ>T_A further indicates that cycle 18 had ARs with very high tilt-angle scatter – or, equivalently, the highest degree of stochasticity – accounting for its exceptionally high amplitude.

Fig. 3. Statistical correlation analysis. Top panel: Correlation between the flux-weighted tilt coefficient multiplied by total flux (Φ>TA) for the nth solar cycle and the observed average polar flux at the end of the same cycle. This correlation is estimated using the sunspot database without incorporating stochasticity. Bottom panel: Correlation between the Φ>TA and two quantities at the end of the same solar cycle (n): (i) the dipole moment simulated by SPhoTraM simulation (left x-axis, red circles), and (ii) the observed average polar flux (right x-axis, magenta squares). These correlations are estimated after incorporating stochasticity in the sunspot emergence statistics. The gray dashed lines in both panels represent the best-fit linear regressions.

On the other hand, the sudden drop in Φ>T_A observed for solar cycle 19 offers a contrasting and compelling explanation for the weak amplitude of solar cycle 20. When the tilt angle of a sunspot pair is small, a larger portion of magnetic flux cancels internally within the active region itself, thereby reducing the efficiency of poleward flux transport. This diminished transport of magnetic flux weakens the buildup of the polar field, which in turn contributes to a weaker subsequent solar cycle – as observed in cycle 20. Additionally, the Φ>T_A parameter accounts for the influence of anomalous sunspots, including anti-Hale ARs, which possess reversed polarities relative to standard sunspots. Being a stronger cycle, solar cycle 19 exhibits a higher number of such anomalous ARs compared to cycle 18, as shown in Figure A.1 in Appendix A.1. These anomalous ARs contribute a greater amount of opposite-polarity flux toward the poles, ultimately diminishing the net polar field strength by the end of cycle 19. Consequently, anti-Hale ARs play a crucial role in generating extreme solar cycle variability, where a high-amplitude cycle culminates in a weak polar field, leading to a subdued following cycle.

3. Summary and conclusions

Solar cycle 19 stands out as exceptionally strong, yet the polar field build-up at its end is surprisingly weak, leading to a weaker subsequent solar cycle 20. Such variations in solar cycles can arise from multiple factors, including nonlinear mechanisms and stochastic processes involved in magnetic field generation and transport. The lack of long-term observations makes it challenging to constrain these processes. For this study, we employed a novel ensemble-run methodology applied to a coupled, data-driven solar surface flux transport model and solar dynamo model to reconstruct the past ten solar cycles. We specifically focused on understanding the extreme variation from sunspot cycles 19 to 20.

Our results reinforce the idea that the key factor driving cycle-to-cycle variations on timescales ranging from decades to centuries is the random scatter in sunspot tilt angles. More importantly, our simulations show that reasonable fluctuations in the tilt angle of ARs – within the observed range of variabilities – are able to recover the significant drop in polar field amplitude from cycle 18 to cycle 19 that act as seeds for the historically strong sunspot cycle 19, and the much weaker cycle 20, respectively. Our results imply that no exotic new physics needs to be invoked to explain extreme fluctuations observed during the Modern Maximum, or over such ‘supradecadal envelops’ of unusual solar activity. Our data-driven, observationally constrained physics-based simulations lend further credence to the emerging understanding that stochastic forcing – and not nonlinear quenching – is the primary driver of solar variability.

Building upon our findings, we can surmise that such random (stochastic) fluctuations, manifest in solar ARs with highly anomalous tilt (and sometimes with large flux), may indeed result in a catastrophic reduction in the value of the polar field – thus precipitating a Maunder-like grand minimum (alluded to in Nagy et al. (2017)). This appears to be a distinct possibility based on our work, which recovers the extreme fall in amplitude from sunspot cycle 19 to 20. While Nagy et al. (2017) suggest a single rogue or extreme anomalous AR may achieve this, a number of anomalous ARs may achieve an analogous effect in precipitating grand minima episodes. Which of these scenarios is more probable is an intriguing question; for a relevant study see Nagy et al. (2020). Our work motivates further investigations into these distinct possibilities.

Acknowledgments

This research was conducted at the Center of Excellence in Space Sciences India (CESSI), supported by IISER Kolkata, the Ministry of Education, Government of India. We acknowledge Prantika Bhowmik, Paul Charbonneau, and Kristóf Petrovay for useful discussions.

References

Appendix A: Numerical model description

A.1. Data-driven SPhoTraM simulation

A.1.1. Model equation

The Solar Photospheric-Flux Transport Model (SPhoTraM) is a newly developed two-dimensional numerical model designed to simulate the solar surface magnetic flux transport processes. It solves the radial component of the magnetic induction equation in the spatial domain in the presence of diffusion η (mimicking the effect of super-granular flows) and transport profiles like meridional circulation v(θ) and differential rotation ω(θ). SPhoTraM employs various finite difference schemes combined with flux limiter algorithms to numerically solve the magnetic induction equation with improved accuracy. In general, the surface flux transport model is the replication of the well-known BL mechanism, which can be expressed by the following mathematical equation:

$$\begin{array}{c}\hfill \frac{\partial B_r}{\partial t}=-\omega (\theta )\frac{\partial B_r}{\partial \varphi }-\frac{1}{{\mathrm{R}}_{\odot }sin\theta }\frac{\partial }{\partial \theta }(v(\theta )B_{r}sin\theta )\\ \hfill +\frac{\eta }{{\mathrm{R}}_{\odot }^2}[\frac{1}{sin\theta }\frac{\partial }{\partial \theta }(sin\theta \frac{\partial B_r}{\partial \theta })+\frac{1}{{sin}^2\theta }\frac{{\partial }^2{B}_r}{\partial {\varphi }^2}]+S(\theta ,\varphi ,t).\end{array}$$(A.1)

In our model, the advective terms, v(θ) follow an axisymmetric profile, peaking at mid-latitudes with a velocity of 15 m/s and going to zero beyond ±75°, while ω(θ) is derived from helioseismic observations, with the profile following Bhowmik & Nandy (2018), Pal et al. (2023). We set η to 250 Km²/s. These models are driven by incorporating as input the emergence profile and statistics of bipolar sunspot pairs. In our formulation, this is encapsulated in the source term S(θ, ϕ, t), which represents the emergence of magnetic flux. This term incorporates contributions from newly formed sunspot regions based on observational or theoretical input. Each emergence is modeled as an idealized bipolar magnetic region, following a prescribed mathematical profile consistent with previous studies (Jiang et al. 2011; Bhowmik & Nandy 2018; Pal et al. 2023).

A.1.2. Model source term

To generate the source term, we require information on the emergence time, heliographic latitudinal and longitudinal positions, and area of active regions for each solar cycle, all of which are taken from the observational database provided by RGO/USAF/NOAA Data Centre (1874-2023). We record the emergence information for each sunspot when it reaches its maximum area. After 1976, the sunspot records transitioned from the RGO database to the USAF/NOAA database. To ensure consistency between the two datasets, we apply a cross-calibration by multiplying a correction factor of 1.4 to all areas less than 206 microhemispheres for active regions appearing post-1976 (Bhowmik & Nandy 2018). Other quantities, such as the magnetic flux associated with a sunspot, tilt angle, sunspot group radius, and the separation between the two polarities of a sunspot pair, are derived from the observational data using empirical and mathematical relationships. The magnetic flux (Φ>) of each active region is calculated using the relation: Φ(A) = 3.5 × 10¹⁹ × A, where A is the sunspot area in microhemispheres, and Φ> is the total magnetic flux in Maxwells, which is assumed to be equally distributed in the two polarities. We estimate the radius of the sunspot from the area information and consider that the radial separation between two polarities of a sunspot pair is proportional to the radius of the sunspot. We compute the standard tilt angle using the square root relation, $\gamma =\mathrm{C}{\mathrm{T}}_{\mathrm{n}}\sqrt{\lambda }$ to incorporate Joy’s tilt law in a sunspot pair. Here, γ is the tilt angle, λ is the latitudinal position of the centroid of the AR and T_n is the tilt coefficient for n^th solar cycle. The variation in tilt angles with solar cycle strength (known as tilt quenching) is modeled through T_n = 1.73  −  0.0035 S_n, as derived from Jiang et al. (2011), where S_n denotes the solar cycle strength. The constant factor C is set to 0.7, accounting for the reduction in tilt angles caused by near-surface localized inflows (Jiang et al. 2014).

A.1.3. Introducing stochasticity in the source term

To introduce tilt scatter, we adopt a methodology commonly used in previous studies (Jiang et al. 2014; Jiao et al. 2021). These studies show that the tilt angles of the sunspots deviate from Joy’s Law in a manner consistent with a Gaussian distribution. This distribution has a zero mean value, indicating no systematic bias, while its standard deviation (σ) depends on the sunspot area (A) according to the empirical relation: σ = −11 × log(A)+35 (For a detailed illustration, see Appendix in Lemerle et al. (2015)). For each sunspot group, we randomly sample a tilt scatter value (ϵ) from this Gaussian distribution (having zero mean and σ standard deviation), which is determined by that sunspot group’s area obtained from the observed database. The resulting tilt angle is then given by $\gamma =\mathrm{C}{\mathrm{T}}_{\mathrm{n}}\sqrt{\lambda }+ϵ$. In our simulations, incorporation of random tilt-angle scatter gives rise to anomalous sunspots, including anti-Joy regions (especially those which carry opposite tilt orientations to standard Joy ARs) and ARs with unusually large tilt angles. We also consider that the anti-Hale ARs are distributed randomly all over the cycle phase, emergence latitude, and the two hemispheres. These anti-Hale sunspot distributions are generated from independent random realizations following the methodology described by Pal et al. (2023). In Figure A.1, such anomalous AR distributions (Latitudinal and flux distribution) for solar cycles 18 and 19 are shown.

To quantify the contribution from stochasticity, we calculate the flux-weighted tilt coefficient, T_A = ($\frac{{\displaystyle {\sum }_{\mathrm{i}}}{\gamma }_{\mathrm{i}}\times {\varphi }_{\mathrm{i}}}{{\displaystyle {\sum }_{\mathrm{i}}}{\varphi }_{\mathrm{i}}|{\lambda }_{\mathrm{i}}|}$), where γ_i, ϕ_i, and |λ_i| represent the tilt angle, flux content, and the absolute latitude of the i^th sunspot. When this coefficient is multiplied by the signed total flux emergence during that cycle this generates ΦT_A – which provides physical insight into the average contribution of tilted sunspots to the resulting polar flux at the end of that cycle (Dasi-Espuig et al. 2010; Jiao et al. 2021).

Fig. A.1. Anomalous Active Region distribution. Top panel: Latitudinal distribution of anomalous ARs (both anti-Hale and anti-Joy regions) for solar cycle 18 (filled histogram) and solar cycle 19 (solid histogram). Bottom panel: Flux distribution of anomalous sunspots, with red filled and solid histograms representing the flux throughout the sunspot cycles 18 and 19, respectively. For clarity and better visualization, the x-axis in this panel is on a logarithmic scale.

A.1.4. Estimation of large-scale magnetic fields

We calculate the axial dipole moment and polar flux using the photospheric map generated by the SPhoTraM simulations utilizing the following relations.

$$\begin{array}{c}\hfill \mathrm{DM}(t)=\frac{3}{4\pi R_{\odot }^2}{\displaystyle {\int }_{\varphi =0}^{2\pi }}{\displaystyle {\int }_{\lambda =-\pi /2}^{\pi /2}}{B}_r(\lambda ,\varphi ,t)sin\lambda cos\lambda d\lambda d\varphi ,\end{array}$$(A.2)

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{N/S}(t)={\displaystyle {\int }_{\varphi =0}^{2\pi }}{\displaystyle {\int }_{{\lambda }_{N/S}}^{}}{R}_{\odot }^2{B}_r({\lambda }_{N/S},\varphi ,t)cos{\lambda }_{N/S}d{\lambda }_{N/S}d\varphi \end{array}$$(A.3)

Here λ and ϕ represent latitude and longitude, respectively. R_⊙ is the solar radius. Φ>_N/S(t) denotes the polar flux in the northern and southern hemispheres, respectively, and DM(t) is the global axial dipole moment. The polar field is calculated based on the surface magnetic field distribution only in the polar cap region (±70° to ±90°), whereas the axial dipole moment corresponds to the entire photospheric magnetic field.

A.1.5. Optimization of the polar flux

The observational magnetograms and polar field measurements are only available from 1974 onward. To bridge this gap, we use polar flux estimates derived from Mount Wilson Observatory calibrated polar faculae data, which extend back to 1906 (Muñoz-Jaramillo et al. 2012). We first generate multiple ensembles of AR input source terms for each solar cycle for the Monte Carlo (MC) SPhoTraM simulation. In each ensemble of this MC simulation, our methodology automatically introduces stochastic variation in the tilt angles and in the distribution of anti-Hale regions in terms of flux content, emergence locations, and phase of appearance. Then we compare the SFT simulated polar flux generated from each ensemble run with observed polar flux values and select a set which satisfies both of the following criteria: (1) the simulated polar flux reversal timing must be within ±1 year of the observed reversal time and (2) the simulated polar flux at the end of each solar cycle must lie within a ±20% range of the observed polar flux value. For each sunspot cycle, we identify the optimal solution from the reduced ensembles by selecting the one that best matches the observed polar flux at the cycle minimum. This is how we optimize the simulated polar flux for our study.

Fig. A.2. Vector potential from solar dynamo model and SPhoTraM and the correction function. Left panel: Surface vector potential at the minima of solar cycle 16 (beginning of cycle 17) from SPhoTraM simulation (APh denoted by red curve) and from dynamo run (AD denoted by black dashed line). Right panel: Latitudinal correction function α(θ).

A.2. Dynamo simulation

A.2.1. Model equations

To generate solar cycles, we utilize the two-dimensional kinematic axisymmetric dynamo model that solves the toroidal (B) and poloidal component (B_P) of the magnetic induction equation in the presence of source terms for sustaining the dynamo mechanism (Chatterjee et al. 2004). Inducing a toroidal field from the poloidal field, followed by the regeneration of the poloidal field from the toroidal field, lies at the heart of the dynamo cycle. In this axisymmetric model, the temporal evolution of the vector potential for the poloidal component A^D(r, θ) and the toroidal component of the magnetic field B(r, θ) are described by the following equations:

$$\begin{array}{cc}\hfill \frac{\partial A^D}{\partial t}+\frac{1}{s}[{\mathbf{v}}_p·\mathrm{\nabla }(s{A}^D)]& ={\eta }_p({\mathrm{\nabla }}^2-\frac{1}{s^2})A^D+\alpha B\hfill \\ \hfill \frac{\partial B}{\partial t}+s[{\mathbf{v}}_p·\mathrm{\nabla }\left(\frac{B}{s}\right)]+(\mathrm{\nabla }·{\mathbf{v}}_p)B& ={\eta }_t({\mathrm{\nabla }}^2-\frac{1}{s^2})B\hfill \\ \hfill +s([\mathrm{\nabla }\times (A^D\widehat{}\}{e}_{\varphi })]·\mathrm{\nabla }\mathrm{\Omega })& +\frac{1}{s}\frac{\partial (sB)}{\partial r}\frac{\partial {\eta }_t}{\partial r}\hfill \end{array}$$(A.4)

Here, s = r sin θ. In these equations, v_p represents the poloidal velocity in the meridional plane (or meridional circulation) that causes advection and distortions in the magnetic field. We utilize a single-cell flow in each hemisphere, threading the convection zone. The terms involving η_P and η_t correspond to the diffusion of poloidal and toroidal magnetic fields, respectively.

Source terms in a dynamo model are crucial to compensate for the dissipation of magnetic fields within the convection zone. The internal rotational shear, as determined from helioseismic observations, acts as a source for generating toroidal fields. We adopt the analytical advective and diffusive flow profiles from Passos et al. (2014) for this study. In our simulation, the generation of the poloidal field from the toroidal component is represented by two distinct mechanisms, which are thought to operate within the Sun’s interior. The first one is the BL Mechanism introduced confined in the near-surface layers. The second mechanism, the mean-field α-effect due to helical turbulent convection, operates ona relatively weaker toroidal field in the bulk of the solar convection zone. These two terms are incorporated in the equation as a single source term α. The mathematical profile with appropriate quenching of these source terms is included following Passos et al. (2014).

A.2.2. SFT maps coupled with dynamo framework

We follow the methodology described in Bhowmik & Nandy (2018) to couple the dynamo model with SFT simulations. We emphasize that our model does not incorporate any intrinsic amplitude fluctuation of the α mechanism over time. Instead, the variability in the poloidal field source term is introduced exclusively by integrating the SPhoTraM-generated surface magnetic map at every cycle minimum. In each ensemble of MC simulations, our methodology automatically introduces variation in the tilt angles and in the distribution of anti-Hale regions in terms of flux content, emergence locations, and phase of appearance, which produces stochasticity.

Our SPhoTraM model generates the radial magnetic field B_r^Ph(R_⊙, θ, t_m) which is related to the magnetic vector potential A^Ph(R_⊙, θ, t_m) through the relation B = ∇ × A, which reduces to:

$$\begin{array}{c}\hfill B_r^{\mathit{Ph}}(R_{\odot },\theta ,t_m)=\frac{1}{R_{\odot }sin\theta }\frac{\partial }{\partial \theta }[sin\theta A^{\mathit{Ph}}(R_{\odot },\theta ,t_m)]\end{array}$$(A.5)

From this equation, we compute the surface vector potential A^Ph by integrated B_r^Ph(R_⊙, θ, t_m) separately for the two hemispheres. The governing equations are:

$$\begin{array}{c}\hfill A^{\mathit{Ph}}(R_{\odot },\theta ,t_m)sin\theta =\{\begin{array}{c}For0<\theta <\pi /2,\hfill \\ {\displaystyle {\int }_0^{\theta }}{B}_r(R_{\odot },{\theta }^{\prime },t_m)sin{\theta }^{\prime }d{\theta }^{\prime }\hfill \\ \hfill \\ For\pi /2<\theta <\pi ,\hfill \\ {\displaystyle {\int }_{\theta }^{\pi }}{B}_r(R_{\odot },{\theta }^{\prime },t_m)sin{\theta }^{\prime }d{\theta }^{\prime }\hfill \end{array}\end{array}$$(A.6)

It should be noted that here B_r^Ph(R_⊙, θ, t_m) is calculated by longitudinally averaging the surface magnetic field during solar minima (t = t_m).

After calculating the vector potential from SPhoTraM run, denoted as A^Ph, it is fed into the dynamo model after proper calibration, which replaces the dynamo-generated vector potential A^D (see left panel of Figure A.2 for an example). To calibrate A^D at solar minima against A^Ph, we scale their amplitudes on the solar surface by a constant factor (k), which is determined at the minimum of solar cycle 16 and kept fixed throughout the simulation. Thus we use A^D = k×A^Ph.

After this amplitude calibration, although the amplitudes of A^D and A^Ph are aligned, their latitudinal distributions, specifically A^D(R_⊙, θ, t_m)sin θ and A^Ph(R ⊙ ,θ, t_m)sin θ still differ significantly. To reconcile this, A^D at each solar minimum is further corrected by a latitudinal function α(θ), such that the product α(θ)×A^D(R ⊙ ,θ, t_m)sinθ matches A^Ph(R ⊙ ,θ, t_m)sin θ at the solar surface (see right panel of Figure A.2 for an example). We assume that the correction to the dynamo’s poloidal field – arising from the BL mechanism – is confined to the radial range between 0.8R⊙ and R⊙. This assimilation is implemented sequentially at each cycle minimum and forces the dynamo generation of the toroidal component of the magnetic field, which acts as the source of the sunspot cycle. To elaborate, at each solar minimum, the dynamo simulation is paused, and A^D(r, θ, t_m)sin θ is multiplied by kα(θ) at all grid points within this mentioned radial layer (remember k is fixed but α(θ) changes for each cycle). Once A^D has been updated to reflect the SPhoTraM-generated vector potential A^Ph, the dynamo simulation resumes and continues its evolution. It is noted that the simulated magnetic flux does not reach the observed low value at minima as depicted in Figure 2. This is a general issue in diverse dynamo models due to cycle overlap (Hazra & Nandy 2016; Kumar et al. 2019). This can be by fine-tuning the meridional flow from cycle to cycle, which we do not consider here (Hazra & Choudhuri 2017). We keep all the other plasma parameters in the dynamo model fixed for each solar cycle simulation.

This iterative approach enables us to produce a sunspot number time series driven by the SPhoTraM model, thereby providing a physics-based reconstruction of solar cycles over an extended period, including the extreme variations.

All Figures

Fig. 1. Simulated and observed large-scale magnetic fields of the Sun. Top panel: Variations in the optimized simulated polar flux (solid curves) alongside the observed polar flux obtained from the MWO polar faculae database (dashed curves) for solar cycles 15–23. The red and blue curves represent the northern and southern hemispheres, respectively. Additionally, the black and gray curves in the same panel depict the total yearly and monthly averaged unsigned sunspot flux, derived from the RGO/USAF/NOAA database. Bottom panel: Time-latitudinal distribution of the radial magnetic field (Br) based on the optimized simulation. Here, yellow and blue shades respectively indicate magnetic fields of positive and negative polarity.

In the text

Fig. 2. Simulated solar cycle reconstruction versus observations. The violet curve represents the total yearly averaged unsigned sunspot flux from solar cycle 16 to 24, derived from the USAF/RGO/NOAA database. The orange curve indicates the total unsigned flux simulated using the dynamo model (refer A.2), driven by the poloidal field generated by the SPhoTraM model. A Pearson correlation coefficient of 0.93, calculated with a 99% confidence level, highlights the strong agreement between the observed and simulated flux strengths at cycle maxima.

In the text

Fig. 3. Statistical correlation analysis. Top panel: Correlation between the flux-weighted tilt coefficient multiplied by total flux (Φ>TA) for the nth solar cycle and the observed average polar flux at the end of the same cycle. This correlation is estimated using the sunspot database without incorporating stochasticity. Bottom panel: Correlation between the Φ>TA and two quantities at the end of the same solar cycle (n): (i) the dipole moment simulated by SPhoTraM simulation (left x-axis, red circles), and (ii) the observed average polar flux (right x-axis, magenta squares). These correlations are estimated after incorporating stochasticity in the sunspot emergence statistics. The gray dashed lines in both panels represent the best-fit linear regressions.

In the text

Fig. A.1. Anomalous Active Region distribution. Top panel: Latitudinal distribution of anomalous ARs (both anti-Hale and anti-Joy regions) for solar cycle 18 (filled histogram) and solar cycle 19 (solid histogram). Bottom panel: Flux distribution of anomalous sunspots, with red filled and solid histograms representing the flux throughout the sunspot cycles 18 and 19, respectively. For clarity and better visualization, the x-axis in this panel is on a logarithmic scale.

In the text

	Fig. A.2. Vector potential from solar dynamo model and SPhoTraM and the correction function. Left panel: Surface vector potential at the minima of solar cycle 16 (beginning of cycle 17) from SPhoTraM simulation (APh denoted by red curve) and from dynamo run (AD denoted by black dashed line). Right panel: Latitudinal correction function α(θ).
In the text