© The Authors 2026

1. Introduction

The surface flux transport (SFT) model has successfully reproduced the magnetic flux patterns associated with the Sun’s 22-year magnetic cycle (Yeates et al. 2023). Incorporating meridional flow, differential rotation, and turbulent diffusion, the model, first formulated by Leighton (1964), provides an essential framework for interpreting polar field formation and for predicting the amplitude of future solar cycles through empirical correlations (Muñoz-Jaramillo et al. 2013; Jiang et al. 2018).

Key SFT parameters include the meridional flow speed, surface diffusivity, differential rotation, and the decay time associated with radial flux loss. The emergence of bipolar magnetic regions (BMRs), represented through a source term following Cameron & Schüssler (2007) and implemented as in Petrovay & Talafha (2019), governs the buildup and reversal of the axial dipole moment (Yeates et al. 2023). Because these parameters interact nonlinearly and are not independently constrained, their optimisation requires a systematic exploration of parameter space (Petrovay & Talafha 2019; Whitbread et al. 2017). Uncertainties in BMR emergence properties further complicate this process, as the relationship between sunspot groups and their magnetic flux is not always well constrained (Yeo et al. 2021).

Several limitations of classical SFT formulations motivate efforts to improve optimisations. They include the reliance on uniform diffusivity and steady, axisymmetric flow profiles, the assumption of linear parameter dependences, and the neglect of radial diffusion in some models (Yeates et al. 2023). Such simplifications can introduce errors into predictions of polar field strength and dipole evolution (Petrovay et al. 2020). Moreover, model sensitivity to meridional flow, tilt properties, and decay time can significantly affect the timing of polar field reversal and the ultimate dipole amplitude (Jiang et al. 2013). In parameterised SFT formulations based on statistically averaged source terms, the absence of a decay term can lead to unrealistically persistent dipole fields and delayed reversals (Petrovay & Talafha 2019). However, recent data-assimilative simulations that incorporate observed active-region emergence directly have demonstrated that the observed polar-field evolution can be reproduced without invoking an explicit radial diffusion term (Yeates et al. 2025; Wang et al. 2025), highlighting the model-dependent nature of this parameterisation.

Nonlinear feedbacks related to BMR emergence provide an important mechanism for regulating solar cycle amplitude. Observations indicate that both the mean tilt angle and the average emergence latitude vary with cycle strength, resulting in tilt quenching (TQ) and latitude quenching (LQ). These effects reduce the efficiency of dipole buildup in strong cycles and have been examined using SFT and dynamo models (Jiang 2020; Talafha et al. 2022). In coupled dynamo frameworks, such as the 2 × 2D model of Lemerle et al. (2015), LQ has been shown to play a significant role, comparable to that of TQ, with historical observations (1923–1985) providing stronger evidence of LQ (Yeates et al. 2025). The relative importance of these effects also depends on the ratio of meridional flow speed to diffusivity (Talafha et al. 2022).

Previous studies have investigated the nonlinear regulation of the solar cycle within SFT and dynamo frameworks using different approaches. Cameron et al. (2010) demonstrated that cycle-dependent variations in the Joy’s law tilt can effectively limit dipole growth without invoking an explicit decay term. In contrast, Talafha et al. (2022) introduced observable nonlinearities in the form of TQ and LQ and examined their impact on dipole moment modulation. Optimisation studies based on genetic algorithms, such as Lemerle & Charbonneau (2017), have identified preferred transport regimes that reproduce key features of the solar cycle, but they did not explicitly address how observable nonlinear feedbacks reshape the admissible parameter space. More recently, algebraic treatments have shown that the combined action of TQ and LQ naturally produces a saturation (‘ceiling’) of dipole growth with increasing cycle amplitude (Talafha et al. 2025).

Building on the optimisation framework of Petrovay & Talafha (2019), the present work incorporates both TQ and LQ into the SFT source term, within a parameter-space optimisation that also includes a finite flux-decay timescale. We quantified how these nonlinear feedbacks reshape the admissible domains of meridional flow speed, surface diffusivity, and decay time, and we assessed how their combined action constrains polar field formation within the present parameterised modelling framework. Section 2 describes the model setup and the implementation of the quenching prescriptions. The optimisation results are presented in Sect. 3. The interpretation of the results is discussed in Sect. 4, and the conclusions are summarised in Sect. 5.

2. Methodology

A systematic exploration of the SFT model’s parameter space was performed by Petrovay & Talafha (2019) to optimise its ability to reproduce the solar polar magnetic field. The parameters vary across a wide range of values, and the resulting magnetic field evolution is evaluated against a set of observational constraints, including the timing of polar field reversals, the relative amplitude of the polar field, and the latitudinal structure of the polar field distribution. This approach enables the identification of parameter combinations that produce physically realistic and observationally consistent outcomes.

2.1. SFT model

The evolution of the large-scale radial magnetic field on the solar surface is described by the SFT model through an advective-diffusive transport equation. In its general form, the governing equation is given by

$$\begin{array}{cc}\hfill \frac{\partial B}{\partial t}& =\frac{1}{Rcos\lambda }\frac{\partial }{\partial \lambda }(Bucos\lambda )\hfill \\ \hfill & +\frac{\eta }{R^{2}cos\lambda }\frac{\partial }{\partial \lambda }(cos\lambda \frac{\partial B}{\partial \lambda })-\frac{B}{\tau }+S(\lambda ,t),\hfill \end{array}$$(1)

where B(λ, t) is the radial magnetic field as a function of heliographic latitude λ and time t, u(λ) is the meridional flow profile, η is the surface diffusivity, τ is the decay timescale, and S(λ, t) represents the source term accounting for flux emergence. The model assumes axial symmetry and a radial field approximation, reducing the problem to one dimension in latitude. This formulation captures the essential transport mechanisms: advection by meridional flow characterised by u(λ), carries flux polewards, while differential rotation introduces shear in the longitudinal direction. In the axisymmetric SFT model, longitudinal variations are either averaged out or treated as a background process. Diffusion represented by η, accounts for the dispersal of magnetic flux due to super-granular motions. This process helps to smooth out small-scale magnetic features, leading to a more uniform large-scale field distribution. The diffusion coefficient is assumed to be uniform over the solar surface in this work. A decay term $-\frac{B}{\tau }$ is included as a phenomenological representation of vertical diffusion or other three-dimensional effects not explicitly resolved in the model. This term prevents the indefinite accumulation of the global dipole moment by mimicking the loss processes of magnetic flux. The source term, S(λ, t), represents the emergence of new magnetic flux, typically through bipolar active regions. It is designed to capture the statistical behaviour of flux emergence over an average solar cycle.

2.2. Parameter optimisation framework

To determine the optimal combination of SFT model parameters, we adopted the approach detailed in Petrovay & Talafha (2019), which involves a systematic, large-scale exploration of the three-dimensional parameter space defined by the meridional flow amplitude (u₀), the magnetic diffusivity (η), and the decay timescale (τ). Rather than relying on algorithmic optimisation techniques such as genetic algorithms, this method evaluates the model’s performance across a dense grid of parameter values. This approach involves computing many SFT model realisations across a structured grid of parameter combinations. Each model evolves over a hundred solar cycles until a steady cyclic pattern is established. The resulting magnetic field configurations are then evaluated against key observationally motivated constraints, including the timing of polar-field reversals relative to sunspot minimum (T_rev), the polar-field minimum-to-extremum amplitude ratio (R_pol), defined as R_pol = B_pol(t_min)/B_pol, ext, where B_pol(t_min) is the polar-field amplitude at cycle minimum and B_pol, ext is the extremal (peak absolute) polar-field value during the cycle, the latitude of the polar cap boundary (λ_cap). In addition, we considered constraints on the axial dipole moment evolution, namely the dipole reversal timing (D_rev) and the axial dipole moment minimum-to-extremum amplitude ratio (D_min/ext), defined as D_min/ext = D(t_min)/D_ext. The first three constraints are used to define the primary admissible parameter domains, while the dipole-based constraints provide complementary diagnostics of the global magnetic-field evolution. This separation allows us to distinguish between constraints directly linked to observable polar-field properties and those probing the global dipole evolution, thereby providing a more comprehensive assessment of the SFT parameter space. Models that satisfy all observational criteria are classified as admissible, and the corresponding parameter sets are identified as optimal within the context of the selected meridional flow profile. This approach allows for a transparent delineation of allowed and excluded regions in parameter space and offers clear physical insight into the role of each parameter in shaping solar magnetic field evolution.

To ensure that the optimisation results are not biased by the peculiarities of any specific cycle, the source term S(λ, t) is constructed to represent a statistically averaged solar cycle. It is modulated as a smooth, axisymmetric distribution that captures the net effect of tilted BMR emergence. It consists of pairs of Gaussian flux rings of opposite polarity, whose latitudinal positions and separations evolve throughout the cycle according to empirical fits derived from observational studies. The amplitude and temporal profile of the source are based on the cycle-averaged sunspot emergence rate, ensuring that the resulting field evolution reflects typical solar behaviour rather than cycle-specific anomalies (the exact formulation is given in Petrovay & Talafha 2019).

To incorporate nonlinear effects into the source term of the SFT model, we followed the formulation described by Talafha et al. (2022), where the emergence properties of active regions are modulated by the solar cycle amplitude. Two observable nonlinearities are introduced: TQ and LQ. In the case of TQ, the latitudinal separation Δλ of the BMRs (representing Joy’s law) is reduced for stronger cycles, following the relation

$$\begin{array}{c}\hfill \mathrm{\Delta }\lambda =1^{°}.5sin{\lambda }_0(1-b_{\mathrm{joy}}\frac{A_n-A_0}{A_0}),\end{array}$$(2)

where A_n is the amplitude of the nth cycle, A₀ is the reference amplitude for an average cycle, and b_joy is the nonlinearity parameter controlling the strength of TQ with a reference value of 0.15.

For LQ, the mean latitude (λ₀) of flux emergence is increased with cycle amplitude according to the empirical formula

$$\begin{array}{c}\hfill {\lambda }_n={14}^{°}.6+b_{\mathrm{lat}}\frac{A_n-A_0}{A_0},\end{array}$$(3)

where b_lat ≈ 2.4 is derived from observations (Jiang et al. 2011). This modified mean latitude feeds into a time-dependent latitudinal profile,

$$\begin{array}{c}\hfill {\lambda }_0(t;n)=[26.4-34.2\left(\frac{t}{P}\right)+16.1{\left(\frac{t}{P}\right)}^2]\left(\frac{{\lambda }_n}{{14}^{°}.6}\right),\end{array}$$(4)

where P is the cycle period. Both effects alter the spatiotemporal structure of the source term S(λ, t), leading to cycle-dependent variations in flux emergence characteristics.

To evaluate the realism of each SFT model realisation, we compared its output against a set of well-established observational constraints that characterise the solar polar magnetic field and dipole moment during a typical cycle. Specifically, T_rev must fall within observed ranges, typically around 3.4–4.3 years after the minimum, and the field amplitude ratios must reflect the gradual decline of polar field strength after its peak. Additionally, λ_cap is constrained by the requirement that the field drops to half its polar value between latitudes 65° and 75°. These criteria, summarised in Table 1 of Petrovay & Talafha (2019), provide quantitative benchmarks for model validation. Parameter combinations that produce magnetic field evolutions satisfying all of these constraints are classified as admissible.

The numerical implementation of the SFT model follows a one-dimensional explicit time-stepping scheme, designed to ensure magnetic flux conservation throughout the simulation. The model is run on a latitudinal grid with a resolution of 0.5°, which is sufficient to accurately capture the structure of the polar cap field and meet the constraints on the topknot width. A fixed timestep of 6 hours is used, ensuring numerical stability across all parameter combinations tested. Each simulation is initialised with a dipolar field configuration and evolved over 100 solar cycles, allowing transient effects to decay and a quasi-steady cyclic behaviour to emerge. Once the model reaches steady oscillations, the results from the final cycle are evaluated against the observational constraints.

3. Results

To investigate the impact of LQ and TQ on the admissible parameter space, we first examined the case with no quenching applied (linear case), focusing on how the accepted regions are constrained by the three constraints mentioned earlier. For these tests, the decay timescale was fixed at τ = 8 yr. Figure 1 shows the results for the linear case: the admissible domain is concentrated at η values between 600 and 650 km² s⁻¹ when u₀ lies between 18 and 19 m s⁻¹, with an additional smaller region appearing around η ≈ 250 and u₀ between 14 and 15 m s⁻¹. These baseline regions provide a useful reference for identifying how the parameter space is modified once quenching mechanisms are introduced in subsequent simulations.

Fig. 1. Admissible domains in the (u0, η) parameter space for the case without TQ or LQ, with a decay timescale of τ = 8 yr. The coloured maps indicate the parameter domains that satisfy the three model constraints: the polar-field reversal timing (Trev), the polar-field minimum-to-extremum amplitude ratio (Rpol), and the polar cap boundary latitude (λcap). The numerical annotations denote the ±1σ intervals associated with each constraint, following the optimisation framework of Petrovay & Talafha (2019). The grey-shaded region shows the intersection of these constraints, corresponding to parameter combinations that simultaneously satisfy all three criteria. This case serves as the reference configuration for comparison with the quenching models shown in Figs. 2–5.

Figure 2 illustrates the case where only TQ is applied, which highlights the isolated effect of TQ on the admissible parameter space. In contrast, Fig. 3 shows the influence of LQ alone, where the best-fit parameter values are concentrated around u₀ = 17–19 m s⁻¹ and η ≈ 650 km² s⁻¹. When the two nonlinearities are included simultaneously, Fig. 4, the admissible domain becomes even more restricted, with allowed values clustering near u₀ = 18–19 m s⁻¹ and η = 600–650 km² s⁻¹.

Fig. 2. Same as Fig. 1 but for the case with TQ included (bjoy = 0.15). Compared to the no-quenching case, the effect of TQ is modest, leading to only a slight contraction of the admissible domains.

Fig. 3. Same as Fig. 1 but for the case with LQ included (blat = 2.4). In contrast to the TQ case (Fig. 2), LQ has a stronger impact, narrowing the admissible parameter ranges for both meridional flow speed and diffusivity.

Fig. 4. Same as Fig. 1 but for the case with both TQ (bjoy = 0.15) and LQ (blat = 2.4) included. Compared to the LQ and TQ cases (Figs. 2 and 3), the combined effect of TQ and LQ significantly restricts the admissible domains, shifting them towards narrower bands in u0 and η.

Figure 5 presents the case with both quenching mechanisms and a longer decay timescale of τ = ∞. Compared to the short-decay case in Fig. 4 for τ = 8 yr, the admissible grey-shaded region is considerably larger, demonstrating how the absence of effective decay broadens the parameter space.

Fig. 5. Same as Fig. 1 but for the case with both TQ (bjoy = 0.15) and LQ (blat = 2.4) included, with a long decay timescale of τ = ∞. Compared to the short decay case (Fig. 4), the admissible domain expands considerably, illustrating that the absence of effective flux decay relaxes the parameter constraints, although such models tend to produce unrealistically delayed dipole reversals.

In the absence of quenching, the optimal values for η lie between 300 and 350 km² s⁻¹ and around 650 km² s⁻¹, with u₀ concentrated between 16 and 17 m s⁻¹. When only TQ is included, the best values remain nearly identical to the linear case, which shows that TQ has only a modest effect on the parameter space. By contrast, LQ produces a broader admissible range, with η spanning 300–350 and 600–650 km² s⁻¹, and u₀ extending from 16 to 18 m s⁻¹. When the two quenching mechanisms are applied simultaneously, the admissible domain narrows again, with η clustering near 350 and 650 km² s⁻¹, while u₀ is restricted to the lower range of 14–16 m s⁻¹.

The inclusion of quenching mechanisms has a clear physical impact on the admissible domains of the u₀ − η parameter space. When no quenching is applied, the allowed regions are relatively broad, with multiple disjoint islands of admissibility. Introducing TQ alone produces only a modest contraction of these regions, indicating that the suppression of the Joy’s law tilt does not drastically alter the large-scale transport balance. In contrast, LQ exerts a stronger influence. By systematically shifting flux emergence towards higher latitudes in stronger cycles, it enhances the poleward cancellation of flux and thereby narrows the admissible windows of both u₀ and η. When the two nonlinearities act together, the admissible domains shrink further and shift towards lower meridional flow amplitudes, resulting in a more selective parameter space that reflects the ‘ceiling effect’ reported in recent algebraic approaches (Talafha 2025). This ceiling arises because additional flux input no longer translates into stronger polar fields once quenching limits the efficiency of dipole build-up.

The T_rev constraint shows only a weak dependence on the inclusion of TQ or LQ and remains primarily governed by the balance between meridional flow speed and the effective flux-decay timescale. By contrast, the R_pol is strongly affected by nonlinear feedbacks, with latitude quenching producing the most pronounced reduction of the admissible parameter ranges. The λ_cap exhibits comparatively little sensitivity to source-term modulation across all model configurations, reflecting its predominantly geometrical character. When all three constraints are applied simultaneously, the intersection of admissible ranges contracts substantially under combined TQ and LQ, yielding a markedly narrower region of parameter space than in the unquenched case.

The role of the decay term is equally important. With a short decay timescale (τ = 8 yr), corresponding to strong radial diffusion, the model requires tightly constrained values of η and u₀ to satisfy the observational benchmarks, and the combined effect of quenching mechanisms sharply limits the admissible domains. In contrast, when the decay term is effectively absent (τ = ∞), the admissible parameter space expands considerably, even under both TQ and LQ. Physically, this reflects that without significant flux loss, the system can sustain a broader range of transport conditions while still meeting the polar field constraints. In parameterised SFT models based on statistically averaged source terms, the absence of a decay term can lead to unrealistically persistent dipole fields and delayed reversals. However, recent data-assimilative SFT simulations that incorporate observed active-region emergence directly (e.g. Yeates et al. 2025; Wang et al. 2025) have demonstrated that the observed polar-field evolution can be reproduced without invoking an explicit radial diffusion term.

It is therefore important to emphasise that the flux-decay timescale τ introduced in the present optimisation framework should not be interpreted as a literal physical radial diffusivity. Instead, it serves as a phenomenological parameter that represents unresolved three-dimensional processes, such as turbulent radial mixing or flux submergence below the photosphere, which are not explicitly captured by surface-only transport models. In SFT simulations that assimilate observed active regions, such loss processes can be implicitly accounted for through the source term itself, potentially reducing or eliminating the need for an explicit decay term (Yeates et al. 2025; Wang et al. 2025; Luo et al. 2025). The optimisation results presented here should therefore be understood as conditional on the adopted parameterised source formulation.

This model dependence is further highlighted by recent developments in alternative numerical realisations of the SFT equation. In particular, Athalathil et al. (2026) employed a physics-informed neural network framework to solve the SFT equation in a mesh-free setting, incorporating nonlinear TQ and LQ without invoking an explicit radial diffusion or decay term. Despite the absence of such a loss term, their simulations reproduced the observed cycle-to-cycle modulation of the polar field, suggesting that the need for a decay term may depend on the numerical implementation of nonlinear feedback mechanisms. This supports the interpretation that the admissible parameter domains identified in the present optimisation study are conditional on the adopted parameterised source formulation and may not represent a unique physically required regime.

Figures 6 and 7 illustrate how the admissible domains vary with the decay timescale, expressed as log₁₀(τ), in the presence of both nonlinear quenching mechanisms. In Fig. 6, log₁₀(τ) is plotted against diffusivity η, showing that longer decay timescales (weaker decay) broaden the admissible range, while shorter τ values confine it to narrower intervals. Figure 7 presents the complementary case, where log₁₀(τ) is plotted against the meridional flow amplitude u₀. For a fixed diffusivity, the admissible u₀ values shift towards lower ranges as τ decreases, indicating that stronger flux decay can be compensated by slower poleward transport. Together, these results confirm that, when quenching mechanisms are active, τ, η, and u₀ become tightly coupled, defining a narrow surface of admissible solutions. This contrasts with the broader parameter domains identified by Petrovay & Talafha (2019), where the absence of nonlinear quenching allowed a wider range of (u₀, η) combinations for τ ≃ 5–10 yr. The present results therefore extend their framework by showing that the inclusion of TQ and LQ further constrains the physically viable region, particularly for shorter decay timescales.

Fig. 6. Dependence of the admissible parameter domain on the decay timescale, expressed as log10(τ), plotted against the surface diffusivity (η) for the case including both LQ and TQ (blat = 2.4, bjoy = 0.15). These panels illustrate the dipole-based constraints on the global magnetic-field evolution. Left: Axial dipole moment reversal time (Drev) as a function of η and log10(τ). Right: Normalised axial dipole moment minimum-to-extremum amplitude ratio (Dmin/ext). The colour-coded contours show how the decay timescale modulates both the timing and amplitude of the global dipole moment: shorter decay timescales (lower τ) restrict the admissible domain to lower diffusivities, whereas longer decay timescales broaden the range of permissible η.

Fig. 7. Dependence of the admissible parameter domain on the decay timescale, expressed as log10(τ), plotted against the meridional flow amplitude (u0) for the case including both LQ and TQ (blat = 2.4, bjoy = 0.15), at a fixed diffusivity of η = 450 km2 s−1. These panels illustrate the dipole-based constraints on the global magnetic-field evolution. Left: Axial dipole moment reversal time (Drev) as a function of u0 and log10(τ). Right: Normalised axial dipole moment minimum-to-extremum amplitude ratio (Dmin/ext). The colour-coded contours show that, for a fixed diffusivity, shorter decay timescales (lower τ) restrict the admissible range to slower meridional flows, whereas longer decay timescales allow higher u0.

Figure 7, which examines the case at η = 450 km² s⁻¹ as a representative cut, demonstrates that admissibility collapses to a narrow band in u₀ once nonlinear feedbacks are included. For fixed diffusivity, only a limited range of meridional flow speeds can counterbalance the suppression effects of TQ and LQ while maintaining agreement with observational constraints. This nonlinear coupling is consistent with, but more restrictive than, the results of Petrovay & Talafha (2019), who found that models with τ ≳ 10 yr and no quenching reversed the dipole too late in the cycle. In contrast, acceptable solutions existed for τ ∼ 5–10 yr with higher η values (500–800 km² s⁻¹). In the quenched models presented here, the same τ range remains optimal, yet the admissible (u₀, η) combinations are significantly reduced, forming narrow ridges rather than extended ‘islands’ in parameter space. This demonstrates that nonlinear quenching mechanisms reinforce the decay constraint identified by Petrovay & Talafha (2019), yielding a more self-limiting and physically consistent configuration for solar cycle modulation.

In this context, the decay timescale τ plays a dual role that merits clarification. Physically, it is commonly interpreted as a parametrisation of vertical diffusion or other three-dimensional processes not explicitly captured in SFT models, such as flux submergence or radial turbulent mixing. Within the present optimisation framework, however, τ should be regarded primarily as a phenomenological control that limits the long-term memory of the surface field and prevents the accumulation of unrealistically persistent dipole moments. The preferred range τ ≃ 8–10 yr therefore reflects the effective timescale required for the SFT model to reproduce observed polar-field reversal timing in the presence of nonlinear quenching. Importantly, the inclusion of TQ and LQ does not eliminate the need for a finite decay term; instead, the two effects act together to further restrict the admissible parameter space and enforce a physically self-consistent, self-limiting regime of polar field evolution.

4. Discussion

The admissible domains of η and u₀ listed in Table 1 can be directly compared with the earlier parameter-space study of Petrovay & Talafha (2019). In their work, the inclusion of a decay term was found to be essential for reproducing realistic polar field reversals: models without decay produced excessively late reversals, even if the amplitude constraints could be satisfied. Our results are consistent with this conclusion by showing that the admissible parameter space is much broader when τ → ∞ as seen in Fig. 5, whereas it becomes sharply restricted for τ = 8 yr (Figs. 1–4).

The present optimisation considers both nonlinearities: TQ alone introduces only minor modifications compared to the linear case, while LQ exerts a stronger influence by narrowing the admissible domains of both η and u₀. When both quenching mechanisms are included, the admissible domains converge towards the narrow bands already identified by Petrovay & Talafha (2019), suggesting that the combination of finite decay and nonlinear feedbacks naturally drives the system towards restricted parameter regimes. This provides a physical interpretation for the ‘islands’ of admissibility previously mapped, which can be understood as the combined outcome of decay-driven flux loss and quenching-driven flux suppression.

The quenching parameters adopted in this study, b_joy = 0.15 for TQ and b_lat = 2.4 for LQ, are fixed at values motivated by observational and modelling studies (Jiang et al. 2011; Talafha et al. 2022). While a systematic exploration of the sensitivity to these parameters is beyond the scope of the present work, we note that moderate variations in their amplitudes primarily affect the quantitative extent of the admissible domains rather than their qualitative structure. In particular, the contraction of the admissible parameter space, the stronger influence of LQ relative to TQ, and the emergence of a saturation (‘ceiling’) effect in dipole amplification are robust features that persist for physically reasonable choices of the quenching strengths. The fixed-parameter approach adopted here therefore suffices to assess how observable nonlinear feedbacks reshape the SFT optimisation landscape.

A direct comparison can be made with the SFT simulations of Cameron et al. (2010), who modelled solar cycles 15–21 by incorporating the observed active-region properties into the source term. In addition to the cycle-to-cycle variations of the Joy’s law tilt angles, their use of observed emergence latitudes implicitly introduced latitude-dependent modulation of the source term. Their model successfully reproduced the timing of the polar-field reversals and the amplitude of the open flux without requiring any additional decay term, demonstrating that the empirically observed anti-correlation between cycle strength and mean tilt angle, together with cycle-dependent emergence latitudes, can provide an effective self-regulating mechanism. In contrast, the present optimisation study introduces both TQ and LQ as explicit nonlinear functions in the parameterised source term, together with an adjustable flux-decay timescale τ. The results confirm the central conclusion of Cameron et al. (2010), that nonlinear suppression of tilt plays a critical role in limiting dipole growth, but further show that, within this parameterised framework, the inclusion of latitude-dependent quenching and flux decay substantially tightens the admissible (u₀, η, τ) domains. Whereas Cameron et al. (2010) achieved stable reversals over a broad parameter range using observationally constrained emergence properties, our results indicate that when both quenching mechanisms and decay are accounted for in a statistically averaged source formulation, the acceptable solutions cluster along narrow ridges in parameter space.

A meaningful comparison can also be made with the optimisation results of Lemerle & Charbonneau (2017), who employed a genetic algorithm to calibrate a two-dimensional flux-transport dynamo model against observed solar cycle features. Their study identified an optimal regime characterised by a diffusivity of η ≃ 450–600 km² s⁻¹ and a meridional flow amplitude of u₀ ≃ 12–18 m s⁻¹, which produced the best agreement with the timing and amplitude of the solar dipole reversals. The parameter ranges obtained in the present study are broadly consistent with those results, particularly for models including moderate decay timescales (τ ∼ 8–10 yr). However, when nonlinear TQ and LQ mechanisms are introduced, the admissible domains shrink markedly and align along narrow ridges in the (u₀, η) plane, rather than the broader basins of attraction reported by Lemerle & Charbonneau (2017). This indicates that while the two optimisation approaches converge on similar transport amplitudes, the inclusion of nonlinear feedback enforces tighter coupling between advection, diffusion, and decay processes. In this sense, the present model extends the empirical optimisation of Lemerle & Charbonneau (2017) by providing a physically constrained explanation for the restricted range of solutions that reproduce observed solar cycle variability.

The inclusion of explicit nonlinear quenching mechanisms provides a physically grounded explanation for the self-regulation of the solar cycle amplitude. Tilt and latitude quenching together reproduce the observed saturation of the axial dipole moment, linking flux-transport parameters to the nonlinear backreaction of magnetic activity on the source term. The resulting narrow admissible domains suggest that the solar cycle operates near a marginally stable regime, where modest variations in surface flow or diffusivity can produce significant modulation in dipole strength. This reinforces the interpretation of the solar dynamo as a weakly nonlinear system, consistent with the empirical correlations between tilt angle, emergence latitude, and cycle strength reported in observational studies.

From the perspective of solar cycle predictability, these results imply that the dynamo operates close to a saturation threshold. When quenching feedbacks are active, small perturbations in surface transport or emergence statistics can lead to disproportionately large changes in polar field buildup, thereby limiting the predictive horizon of purely kinematic models. The quenching mechanisms thus provide a natural explanation for the stochastic component of solar cycle variability: while nonlinear feedbacks constrain the mean cycle amplitude, the detailed timing and asymmetry of reversals may remain sensitive to fluctuations in active region emergence.

We also note that the adopted cycle profile follows the analytic formulation of a quadratic fit derived by Jiang et al. (2011) from many observed solar cycles. Although this profile provides a useful statistically averaged description of solar cycle emergence, observed solar cycles are not strictly self-similar. In particular, Jiang et al. (2018) show that the rising phase is cycle-dependent, whereas the declining phase is comparatively less variable. The use of a simplified cycle profile may therefore influence the inferred admissible parameter ranges, and exploring the impact of more realistic cycle-dependent source profiles will be an important direction for future work.

5. Conclusion

This study extends the SFT optimisation of Petrovay & Talafha (2019) by incorporating explicit nonlinear quenching mechanisms, TQ and LQ, together with a tunable flux-decay term. The inclusion of this feedback significantly refines the admissible parameter space of (u₀, η, τ), transforming the broad islands of solutions identified in earlier studies into narrow, physically consistent ridges. The results demonstrate that LQ exerts a stronger influence than TQ, producing a pronounced ceiling effect that limits the amplification of the axial dipole moment in strong cycles. The optimal decay timescale is found to lie near τ ≃ 8–10 yr, consistent with the physically realistic range inferred from previous optimisation studies.

The combined action of TQ, LQ, and flux decay establishes a tightly coupled balance between advection, diffusion, and flux loss. This coupling naturally leads to a self-limiting, weakly nonlinear regime in which the solar cycle amplitude is saturated by feedback from the surface field itself. The admissible domains revealed by the optimisation suggest that the Sun operates close to this saturation threshold, where small perturbations in surface transport or emergence properties can lead to substantial modulation in polar field strength. Consequently, these nonlinear feedbacks inherently limit the predictability of the solar cycle: while the mean amplitude and timing remain constrained by the quenching-regulated balance, stochastic fluctuations in active region emergence or flow perturbations can still drive cycle-to-cycle variability.

The results further suggest a physical connection between latitude-dependent quenching and the inflows towards active regions observed in Doppler and helioseismic measurements. As demonstrated by Talafha et al. (2025), these inflows act as an additional nonlinear regulator that suppresses flux emergence at low latitudes, effectively reinforcing the LQ mechanism included in this work. Such surface inflows, together with global quenching effects, provide a unified picture of how magnetic activity modulates its own transport environment, linking the SFT description at the surface to the deeper dynamo processes operating below.

In summary, the incorporation of nonlinear quenching yields a more constrained and physically self-consistent optimisation of SFT parameters. The emerging picture is that of a dynamo system operating in a marginally stable regime, where nonlinear feedbacks–acting through both surface and subsurface processes–maintain solar cycle amplitudes within a narrow range while allowing for moderate, stochastic variability. Future extensions of this framework will include data assimilation of WSO magnetograms and the explicit treatment of cycle-dependent inflows, enabling the direct observational calibration of the quenching parameters and providing improved constraints for predictive dynamo modelling.

Acknowledgments

This research was funded by the College of Graduate Studies and supported by the Research Institute of Sciences and Engineering (RISE) at the University of Sharjah. This research acknowledges the use of synoptic magnetic field data provided by the Wilcox Solar Observatory (WSO).

References

All Tables

All Figures

Fig. 1. Admissible domains in the (u0, η) parameter space for the case without TQ or LQ, with a decay timescale of τ = 8 yr. The coloured maps indicate the parameter domains that satisfy the three model constraints: the polar-field reversal timing (Trev), the polar-field minimum-to-extremum amplitude ratio (Rpol), and the polar cap boundary latitude (λcap). The numerical annotations denote the ±1σ intervals associated with each constraint, following the optimisation framework of Petrovay & Talafha (2019). The grey-shaded region shows the intersection of these constraints, corresponding to parameter combinations that simultaneously satisfy all three criteria. This case serves as the reference configuration for comparison with the quenching models shown in Figs. 2–5.

In the text

	Fig. 2. Same as Fig. 1 but for the case with TQ included (bjoy = 0.15). Compared to the no-quenching case, the effect of TQ is modest, leading to only a slight contraction of the admissible domains.
In the text

	Fig. 3. Same as Fig. 1 but for the case with LQ included (blat = 2.4). In contrast to the TQ case (Fig. 2), LQ has a stronger impact, narrowing the admissible parameter ranges for both meridional flow speed and diffusivity.
In the text

	Fig. 4. Same as Fig. 1 but for the case with both TQ (bjoy = 0.15) and LQ (blat = 2.4) included. Compared to the LQ and TQ cases (Figs. 2 and 3), the combined effect of TQ and LQ significantly restricts the admissible domains, shifting them towards narrower bands in u0 and η.
In the text

	Fig. 5. Same as Fig. 1 but for the case with both TQ (bjoy = 0.15) and LQ (blat = 2.4) included, with a long decay timescale of τ = ∞. Compared to the short decay case (Fig. 4), the admissible domain expands considerably, illustrating that the absence of effective flux decay relaxes the parameter constraints, although such models tend to produce unrealistically delayed dipole reversals.
In the text

Fig. 6. Dependence of the admissible parameter domain on the decay timescale, expressed as log10(τ), plotted against the surface diffusivity (η) for the case including both LQ and TQ (blat = 2.4, bjoy = 0.15). These panels illustrate the dipole-based constraints on the global magnetic-field evolution. Left: Axial dipole moment reversal time (Drev) as a function of η and log10(τ). Right: Normalised axial dipole moment minimum-to-extremum amplitude ratio (Dmin/ext). The colour-coded contours show how the decay timescale modulates both the timing and amplitude of the global dipole moment: shorter decay timescales (lower τ) restrict the admissible domain to lower diffusivities, whereas longer decay timescales broaden the range of permissible η.

In the text

Fig. 7. Dependence of the admissible parameter domain on the decay timescale, expressed as log10(τ), plotted against the meridional flow amplitude (u0) for the case including both LQ and TQ (blat = 2.4, bjoy = 0.15), at a fixed diffusivity of η = 450 km2 s−1. These panels illustrate the dipole-based constraints on the global magnetic-field evolution. Left: Axial dipole moment reversal time (Drev) as a function of u0 and log10(τ). Right: Normalised axial dipole moment minimum-to-extremum amplitude ratio (Dmin/ext). The colour-coded contours show that, for a fixed diffusivity, shorter decay timescales (lower τ) restrict the admissible range to slower meridional flows, whereas longer decay timescales allow higher u0.

In the text