© ESO 2022

1. Introduction

The origin of intercycle variations in the series of 11-year solar activity cycles and the possibilities to predict them have been the focus of intense research over the past decade (Petrovay 2020; Charbonneau 2020; Nandy 2021). In a physical system displaying periodic or quasiperiodic behaviour, intercycle variations must be due to either inherent nonlinearities or external forcing. A plausible example is the way in which the strictly sinusoidal, periodic solutions of the simple linear oscillator equation turn into less regular quasiperiodic variations when nonlinear or stochastic forcing terms are added. Indeed, under some simplifying assumptions, the equations describing a solar flux transport dynamo can be truncated to a nonlinear oscillator equation (Lopes et al. 2014).

Here, “external” forcing refers to the physical model under consideration and not necessarily to the object of study, namely: in mean-field models of the solar dynamo, stochastic variation triggered by small-scale disturbances that are not explicitly modelled is the most plausible source of such forcing. Possibilities for making solar cycle predictions are determined by the nature and interplay of these nonlinear and stochastic effects.

Based on its observed characteristics, the solar dynamo is likely to be an αω-dynamo, wherein strong toroidal magnetic fields are generated by the windup of a weaker poloidal (north–south) field by differential rotation. The strong toroidal field buried below the surface then occasionally and locally emerges to the surface and protrudes into the atmosphere in the form of east–west oriented magnetic flux loops that are observed as solar active regions. Besides their dominant azimuthal magnetic field component, these loops, presumably under the action of Coriolis force, also develop a poloidal field component which is typically oriented opposite to the preexisting global poloidal field. The evolution of the poloidal magnetic field on the solar surface (where it is effectively radial) can be followed on solar full-disk magnetograms and found to be well described by phenomenological surface flux transport (SFT) models (Jiang et al. 2014). In such SFT models, the evolution of the radial magnetic field at the surface is described by an advective-diffusive transport equation, with an added source term representing emerging active regions and an optional sink term to describe the effects of radial diffusion. Advection is attributed to the poleward meridional flow, while diffusion is due to supergranular motions.

At solar minimum, the field is found to be strongly concentrated at the Sun’s poles, so the dominant contribution to the solar dipole moment, namely:

(1)

(where λ: latitude, t: time, B: azimuthally averaged field strength) comes from the polar caps. It is this polar field that serves as seed for the windup into the toroidal field of the next solar cycle, so for intercycle variations it is sufficient to consider the azimuthally averaged SFT equation. This is expressed as:

(2)

where R is the solar radius, η is the supergranular diffusivity, U is the meridional flow speed, τ is the time scale of decay due to radial diffusion, and s is the source representing flux emergence.

As Eq. (2) is linear, the amplitude of the dipole moment, D_i, built up by the end of cycle i − 1 is uniquely set by the inhomogeneous term, s. The windup of the poloidal field, D_i, in turn, is expected to be linearly related to the peak amplitude of the toroidal field in cycle, i, and, hence, to the amplitude of the next activity cycle – this is indeed confirmed by both observations and dynamo models (Muñoz-Jaramillo et al. 2013, Kumar et al. 2021). The only step where nonlinearities and stochastic noise can enter this process is thus via flux emergence, linking the subsurface toroidal field to the poloidal source s(λ, t) at the surface. As s(λ, t) is directly accessible to observations, this implies that the nonlinearities involved may, in principle, be empirically constrained.

Indeed, Dasi-Espuig et al. (2010) reported that the average tilt of the axis of bipolar active regions relative to the azimuthal direction is anticorrelated with regard to the cycle amplitude. Starting with the pioneering work of Cameron et al. (2010), this tilt quenching (TQ) effect has become widely used in solar dynamo modeling, especially as the product of the mean tilt with the cycle amplitude has been shown to be a good predictor of the amplitude of the next cycle. The observational evidence, however, is still inconclusive, and the numerical form of the effect remains insufficiently constrained (see, however, Jiao et al. 2021 for more robust recent evidence and also Jha et al. 2020 for the effect’s dependence on field strength).

Jiang (2020) recently called attention to another nonlinear modulation mechanism: latitude quenching (LQ). This is based on the emprirical finding based on an analysis of a long sunspot record, which shows that the mean latitude where active regions emerge in a given phase of the solar cycle is correlated with the cycle amplitude. From higher latitudes, a lower fraction of leading flux can manage become diffused across the equator, leaving less trailing flux to contribute to the polar fields. Therefore, the correlation found here represents a negative feedback effect. Assuming a linear dependence for both the mean tilt and the mean latitude on cycle amplitude, based on one particular SFT setup, Jiang (2020) found that TQ and LQ yield comparable contributions to the overall nonlinearity in the process of the regeneration of the poloidal field from the poloidal source. The net result showed that the net dipole moment change during a cycle tends to be saturated for stronger cycles.

Our objective in the present work is to further explore the respective roles played by TQ and LQ in the solar dynamo. In Sect. 2, we extend the analysis to other SFT model setups and study how the respective importance of TQ and LQ depend on model parameters. In Sect. 3, we further consider a case where the form of TQ is explicitly assumed to be nonlinear, namely: the 2×2D dynamo model. Our conclusions are summarized in Sect. 4.

2. Quadratic nonlinearities in a surface flux transport model

2.1. Model

For the study of how different nonlinearities compare in SFT models with different parameter combinations, our model setup is a generalization of the approach of Petrovay & Talafha (2019) (hereafter Paper I).

In the prior study, our intention was to model “typical” or “average” solar cycles, so that our source function would not consist of individual ARs but a smooth distribution (interpreted as an ensemble average) representing the probability distribution of the emergence of leading and trailing polarities on the solar surface. By nature, this source is thus axially symmetric, so our whole SFT model is reduced to one dimension, as described by Eq. (2).

For the meridional flow, we consider a sinusoidal profile with a dead zone around the poles,

(3)

with λ₀ = 75°. This profile has been used in many other studies, including Jiang (2020), who took the values u₀ = 11 m s⁻¹, η = 250 km² s⁻¹, and τ = ∞.

Our source term is a generalization of the source used in Paper I. The source is a smooth distribution representing an ensemble average or probability distribution of the emergence of leading and trailing polarities on the solar surface (Cameron & Schüssler 2007, Muñoz-Jaramillo et al. 2010), represented as a pair of rings of opposite magnetic polarities:

(4)

where k = ±1 is a factor that depends on the sign of the toroidal field and A_n is the amplitude for cycle n. In the numerical implementation of this source profile, care was taken to ensure zero net flux on the spherical surface by reducing the amplitude of the equatorward member of each pair by an appropriate sphericity factor f_sph.

s₁ is the time profile of activity in a cycle while s₂ characterizes the latitude dependence of activity at a given cycle phase. In contrast to Paper I where a series of identical cycles was considered, here we allow for intercycle variations, so A_n will be different for each cycle.

The time profile of solar activity in a typical cycle was determined by Hathaway et al. (1994) from the average of many cycles as:

(5)

with a = 0.00185, b = 48.7, c = 0.71, where t_c is the time since the last cycle minimum.

The latitudinal profile s₂[λ;λ₀(t),δλ] is a Gaussian migrating equatorward during the course of a cycle:

(6)

with the constant fixed as . Following the empirical results of Jiang et al. (2011), the standard deviation is given by:

(7)

where P = 11 year is the cycle period. We note that alternate fitting formulae were determined by Lemerle et al. (2015) on the basis of a single cycle.

We obtain the latitudinal separation 2Δλ of the rings as a consequence of Joy’s law, while the mean latitude λ₀ was again empirically determined by Jiang et al. (2011). These parameters, in turn, are now dependent on cycle amplitude introducing two nonlinearities into the problem.

Tilt Quenching (TQ): Joy’s law varies depending on the amplitude of the cycle as

(8)

Here, A₀ is an arbitrary reference value corresponding to a “typical” or average cycle amplitude, as set by the nonlinear limitation of the dynamo. Equation (8), thus, measures the departure of the feedback Δλ = f(A_n) from its value in this reference state; its assumed linear form should be valid for a restricted range in Δλ and is not incompatible with the rather lax observational constraints,– cf. Fig. 1b in Jiang (2020). The nonlinearity parameter b_joy is another free parameter to be explored; our reference value for it will be 0.15¹.

Latitude Quenching (LQ): following Jiang et al. (2011) for the mean latitude λ₀ of activity at a given phase of cycle i we have:

(9)

with the mean latitude in cycle i given by:

(10)

Here, the nonlinearity parameter is fairly tightly constrained by the empirical results of Jiang et al. (2011) as b_lat ≃ 2.4.

For the SFT model parameters, our model grid will be a subset of the grid described in Paper I. For τ, we only consider two values here: a decay time scale of eight years, comparable to the cycle length and supported by a number of studies (see Paper I and references therein), and the reference case with effectively no decay (τ = ∞ years). For u₀ = 10 m s⁻¹, the full set of η values in the grid is considered; and similarly for η = 500 km² s⁻¹, the full set of u₀ values is studied (as listed in Table 2). Furthermore, for the study of the effects of individual nonlinearities, for each parameter combination in the SFT model four cases will be considered, as listed in Table 1.

For each combination of the parameters u₀, η, τ, b_lat, and b_joy we run the SFT code solving Eq. (2) for 1000 solar cycles to produce a statistically meaningful sample of cycles.

2.2. Results

In order to study how the nonlinearities affect the dipole moment built up during cycles of different amplitudes, here we consider the purely stochastic case where individual cycle amplitudes are drawn from a statistical distribution without any regard for previous solar activity. The factor k = ±1 in Eq. (4) is assumed to alternate between even and odd cycles. The form of the statistical distribution chosen is not important for our current purpose but, for concreteness, fluctuations are represented as multiplicative noise applied to the amplitude of the poloidal field source. This noise then follows lognormal statistics: A_n = A₀ × 10^G, where G(0, σ) is a gaussian random variable of zero mean and standard deviation of σ = 0.13. The coefficient for an average cycle is set to A₀ = 0.001 e^7/τ to ensure that the resulting dipolar moments roughly agree with observed values (in Gauss).

We note that the latitude-integrated amplitude of our source term s(λ, t) must be linearly related to the observed value of the smoothed international sunspot number:

The proportionality coefficient, K, was arbitrarily tuned to closely reproduce the observed mean value of S for each value of τ.

Introducing the notation r = e^−11/τ (with τ in years), the net contribution of cycle n to the change in dipolar moment is given by ΔD_n = D_n + 1 − rD_n. Figures 1 and 2 display |ΔD| against the cycle amplitude S_n for some example runs.

Fig. 1. Net contribution of a cycle to the solar dipolar moment vs. cycle amplitude for the linear case and for the cases with tilt quenching and/or latitude quenching, with quadratic fits. Parameter values for this run were u0 = 5 m s−1, η = 250 km2 s−1 and τ = 8 years.

Fig. 2. Net contribution of a cycle to the solar dipolar moment vs. cycle amplitude for the linear case and for the cases with tilt quenching and/or latitude quenching, with quadratic fits. Parameter values for this run were u0 = 20 m s−1, η = 600 km2 s−1 and τ = ∞ years.

The linear case (noQ) is obviously well suited to a linear fit, while all three nonlinear cases are seen to be well represented by quadratic fits forced to intersect the linear fit at the mean value of S_n.

The respective importance of TQ and LQ is seen to differ for these cases. For a quantitative characterization of the influence of nonlinearities we use the deviation of |ΔD_n| from the linear case at an arbitrary selected value of S_n (roughly twice the mean), as illustrated by the solid vertical lines in Figs. 1 and 2. To give an example, in the TQ case, the deviation is dev_TQ = 7.2281−6.1647 = 1.0634 for the model in Fig. 1. The relative importance of latitude quenching versus tilt quenching is then characterized by the ratio of these deviations for the LQ and TQ cases, dev_LQ/dev_TQ. Comparing these values to several parameter combinations reveals that the relative importance of the two effects is primarily determined by the ratio u₀/η. This is borne out in Fig. 3.

Fig. 3. Relative importance of LQ vs. TQ, plotted against the parameter u0/η for the parameter combinations listed in Table 2. Separate quadratic fits for the case τ = ∞ (orange) and τ = 8 years (blue) are shown.

2.3. Interpretation

Over a short interval of time [t, t + dt], the smooth source distribution assumed, Eq. (4) is reduced to a bipole with a flux of

(11)

in each polarity.

The ultimate contribution of such a bipole to the solar axial dipole moment after several years has been analytically derived by Petrovay et al. (2020) as:

(12)

where δD₁ is the initial contribution of the bipole to the global dipole moment, t_n + 1 is the time of the minimum of the next cycle, while the asymptotic dipole contribution factor f_∞ is given by:

(13)

Here, a is a constant and defining Δ_u as the divergence of the meridional flow at the equator and the dynamo effectivity range λ_R is a universal function of Δ_u/η for plausible choices of the meridional flow profile.

For the source considered here:

(14)

Combining the last three equations and integrating for the full duration of the cycle yields the net change in the solar dipole moment:

(15)

with

(16)

This shows that |ΔD_n| and, hence, the nonlinear effects in it, only depend on model parameters through τ and the dynamo effectivity range λ_R. For a fixed meridional flow profile, as used in this work, Δ_u ∼ u₀, hence λ_R is uniquely related to u₀/η, explaining the dependence on the latter variable, as shown in Fig. 3. On the other hand, usinq Eqs. (26) and (27) of Petrovay et al. (2020), λ_R can be derived for each of our models, allowing us to plot the underlying relation of the nonlinearity ratios with λ_R. This plot is displayed in Fig. 4.

Fig. 4. Relative importance of LQ vs. TQ, plotted against the dynamo effectivity range, λR, for the parameter combinations listed in Table 2. Separate fits of the form are shown for the case τ = ∞ (orange) and τ = 8 years (blue). Large coloured symbols mark the loci of SFT and dynamo models studied in recent literature: Bhowmik & Nandy (2018) (), Jiang & Cao (2018) (), Lemerle & Charbonneau (2017) (), Whitbread et al. (2018) (+).

In order to understand the form of this relation we recall from Eqs. (8)–(10) that

(17)

Upon substituting the Eqs. (8)–(10) into (15), the deviations from the linear case correspond to the terms containing b_joy and b_lat. We will consider these terms in the limit τ → ∞ and, relying on the mean value theorem for integrals, time integrals will be substituted by the product of the cycle period and the value of the integrand for an appropriately defined average. These weighted averages are not expected to differ from the unweighted average by more than a few degrees, so they will be substituted here simply by for an approximation.

In the TQ case, using Eq. (8), the term containing b_joy is immediately found to scale as

(18)

The LQ case is slightly more complex. In the limit τ → ∞, the dependence on the latitude of activity enters Eq. (15) through the combination . A Taylor expansion yields:

(19)

Substituting here and plugging the expression back into Eq. (15), the term involving b_lat is found to scale as:

(20)

Dividing Eq. (20) by (18) yields

(21)

This expectation is indeed confirmed by the fits in Fig. 4.

2.4. Discussion

It is to be noted in Fig. 4 that for high values of λ_R, the value of dev_LQ/dev_TQ becomes slightly negative. This is due to the circumstance that latitudinal quenching can affect ΔD in two ways:

(i) by modulating the fraction of active regions above λ_R that do not contribute to ΔD, which results in a negative feedback; (ii) by modulating the mean tilt angle of active regions in accordance with Joy’s law, which is a positive feedback effect.

For the first effect dominates. However, when λ_R significantly exceeds the fraction of ARs above λ_R becomes negligible, so only the second, positive feedback effect remains. As a result, dev_LQ changes sign.

A slight dependence on τ is also noticed. This can be attributed to the fact that for shorter values of τ the dipole moment contribution of ARs appearing early in the cycle will have more time to decay until the next minimum, hence, these ARs will have a lower relative contribution to |ΔD|. But the ARs in the early part of the cycle are the ones at the highest latitudes, that is, they are most affected by LQ via effect (i) above. Hence, for shorter values of τ, the relative importance of LQ is expected to be slightly suppressed. This suppression may be expected to vanish towards higher values of λ_R, in agreement with Fig. 4.

The larger, coloured symbols in Fig. 4 mark the positions of various dynamo and SFT models discussed in recent publications. It is apparent that the case where LQ results in a positive feedback is not strictly of academic interest: at least one SFT model discussed recently, namely, that of Whitbread et al. (2018) is in this parameter range. The parameter choice for this model was motivated by an overall optimization of observed versus simulated magnetic supersynoptic maps (Whitbread et al. 2017). At the other end of the scale, the low diffusivity, medium flow speed models of Bhowmik & Nandy (2018) and Jiang & Cao (2018) are in a range where LQ is manifest as a negative feedback giving the dominating contribution to the nonlinearity, although the effect of TQ remains non-negligible.

Conversely, the locus of the SFT component of the 2×2D model of Lemerle & Charbonneau (2017) in Fig. 4 indicates that if the TQ and LQ effects were formulated as in our model, the contributions of TQ should dominate but the negative feedback due to LQ would also contribute quite significantly to the nonlinearity. We take a closer look at this model in the following section.

As a caveat regarding the validity of these conclusions, we recall that the results presented in this section are based on the assumption of one particular form of the time profile of the source function, Eq. (5), with constant parameters and a strictly fixed periodicity of 11 years. As long as deviations from it are stochastic, this source may still be considered representative as an ensemble average. However, there have also been indications of systematic deviations.

It is well known that the length of solar cycles is inversely related to their amplitude: this effect implies an inverse correlation of the parameter b in Eq. (5) with a (or, equivalently, S_max). It is, however, straightforward to see that a variation in b will not affect our results. If c is kept constant, as suggested by Hathaway et al. (1994), b is the only time scale appearing in s₁, so, with S_max fixed, its variation will simply manifest as a self-similar temporal expansion or contraction of the profile (5). As a result, cycle lengths will then vary and so will the area below the curve, namely, the total amount of emerging flux, and, hence, the amplitude of the dipole moment source. As a result, the dipole moment built up by the end of the cycle will linearly scale with . This scaling factor, however, will be the same for the TQ, LQ, and subsequent cases, so the ratio between their deviations, plotted in Fig. 4, remains unaffected.

Furthermore, Jiang et al. (2018) noted that actual cycle profiles tend to show a systematic upward deviation from the profile (5) in their late phases. Such an effect would indeed had an impact on our results somewhat. As in these late phases the activity is concentrated near the equator, at latitudes well below λ_R, the extra contribution to the dipole moment due to this excess is not expected to be subject to LQ, while it is still affected by TQ. Hence, dev_LQ/dev_TQ, plotted in Fig. 4, is expected to be somewhat lower with this correction. This relative difference, in turn, may be expected to vanish towards higher values of λ_R, where LQ becomes insignificant anyway.

3. Nonlinear effects in the 2×2D dynamo model

In this section, we consider the role of nonlinearities in a full dynamo model, namely, the 2×2D dynamo model (Lemerle et al. 2015 and Lemerle & Charbonneau 2017). This model, also used for solar cycle forecast (Labonville et al. 2019), couples a 2D surface SFT model with an axisymmetric model of the dynamo operating in the convective zone. The azimuthal average of the poloidal magnetic field resulting from the SFT component is used as upper boundary condition in the subsurface component, while the source term in the SFT model is a set of bipolar magnetic regions (BMRs) introduced instantaneously with a probability related to the amount of toroidal flux in a layer near the base of the convective zone, as a simplified representation of the flux emergence process. The model has a large number of parameters optimized for best reproduction of the observed characteristics of solar cycle 21.

3.1. Nonlinearity vs. stochastic scatter in the 2×2D dynamo

Figure 5 presents the equivalent of Figs. 1 and 2 for 5091 cycles simulated in the 2×2D model. A marked nonlinearity is clearly seen: the net contribution of a cycle to the dipole moment increases slower than in the linear case with cycle amplitude, representing a negative feedback, similarly to the case of the pure SFT model.

Fig. 5. Contribution of all active regions in a cycle to the net change in the solar dipole moment during dynamo cycles vs. the cycle amplitude in the 2×2D dynamo model. Solid line: Mean; dashed curves: ±1σ range.

In the dynamo model, the dipole moment built up during a cycle feeds back into the cycle as the amplitude of the toroidal field generated by the windup of the poloidal fields (and, hence, the amplitude of the upcoming cycle) increases with the amplitude of the poloidal field. This is borne out in Fig. 6. We note that the relation is not exactly linear and there is a slight threshold effect at work, which is presumably related to the requirement that the toroidal field must exceed a minimal field strength for emergence.

Fig. 6. Solar dipole moment at cycle minimum vs. the amplitude of the subsequent cycle in the 2×2D dynamo model.

The combined result of the two feedbacks presented in Figs. 5 and 6 is a nonlinear coupling between the amplitudes of subsequent cycles as shown in Fig. 7. The convex shape of the median curve implies that on average, very weak cycles will be followed by stronger cycles and very strong cycles will be followed by weak cycles. However, there is a very large amount of scatter present in the plot. A comparison with the scatter in Figs. 5 and 6 shows that the main contribution to the scatter comes from Fig. 5. As ΔD is uniquely determined by the source term in the SFT equation representing flux emergence, the main cause of random scatter around the mean nonlinear relation is identified as the stochastic nature of the flux emergence process, namely, the scatter in the time, location, and properties of individual active regions (BMRs in the 2×2D model) statistically following the evolution of the underlying toroidal field.

Fig. 7. Cycle amplitude vs. the amplitude of the subsequent cycle in the 2×2D dynamo model. Solid line: Mean; dashed curves: ±1σ range.

3.2. Latitude quenching in the 2×2D dynamo model

The question arises as to whether the nonlinearity plotted in Fig. 5 is a consequence of TQ, LQ, or both. Apart from a threshold effect which was not found to be strong or even necessarily needed in the optimized model, the only nonlinearity explicitly built into the 2×2D model is TQ. The tilt angle (i.e., initial dipole contribution δD₁) of a BMR is assumed to scale with

(22)

where B_q is the quenching field amplitude.

As the latitudinal distribution of BMRs is fed into the SFT component of the 2×2D model from the mean-field component, LQ cannot be imposed. Nevertheless, this does not mean that LQ is not present. Indeed, as demonstrated in Fig. 8, the 2×2D model displays a very marked latitude quenching effect, with a tight correlation between the mean emergence latitude of BMRs and cycle amplitude. This LQ is even stronger than witnessed in the Sun, varying by ∼10° as varies between 0.1 and 2. In contrast, according to Eq. (10), only varies by about 5° (from 12° to 17°) for the same range in S_n.

Fig. 8. Latitude quenching in the 2×2D dynamo model. Mean latitude of BMRs in a cycle vs. cycle amplitude.

The origin of this LQ effect is elucidated in Fig. 9. Three subsamples of all simulated cycles are considered: cycles stronger or weaker than the mean amplitude by at least 1 standard deviation are displayed in red or blue, respectively, while black marks cycles not farther from the mean amplitude than 0.1 standard deviation. The plots clearly demonstrate that tracks of latitudinal migration of these three subsets are systematically shifted relative to each other. The alternative possibility that all cycles follow the same track with LQ resulting from the different temporal variation of the activity level can thus be discarded.

Fig. 9. Butterfly tracks in the 2×2D dynamo model sorted by cycle amplitude: (a) tracks of the active latitude vs. cycle phase; (b) tracks of the active latitude vs. time; Red, black, and blue curves represent flux-weighted mean latitudes for subsamples of strong, average, weak cycles, respectively, as described in the text. Shaded areas mark the standard deviations of the actual data around the mean curves.

Cameron & Schüssler (2016) considered the time evolution of mean latitude and width of the activity belts in observed solar cycles, interpreting the decay phase of cycles as a result of the cancellation of oppositely oriented toroidal flux bundles across the equator by turbulent diffusion. They pointed out that this leads to a certain universality of the characteristics of the decay phases in different solar cycles. Our finding that the tracks of migrating toroidal flux bundles systematically differ for strong and weak cycles may at first sight seem to contradict this interpretation. However, a closer look at the properties of decay phases of cycles simulated in the 2×2D model confirms the universality posited by Cameron & Schüssler (2016). In Fig. 10a, we present the evolution of the (pseudo) sunspot number simulated in the models against the latitude of activity belts, rather than against time, while Fig. 10b presents the evolution of the width of the latitudinal distribution of BMRs in a similar way (see Figs. 3 and 4 in Cameron & Schüssler 2016 for the observational equivalents of these plots.) It is apparent that the curves for the strong, average and weak subsamples are nicely aligned towards their low-latitude (i.e., late-phase) ends, confirming the universality proposition.

Fig. 10. Universality of the cycle decay phase in the 2×2D dynamo model: (a) Pseudo sunspot number against flux-weighted mean latitude of BMRs. Curves show flux-weighted mean values; (b) Width 2δλ of the latitudinal distribution of BMRs vs. their flux-weighted mean latitude. Colour coding is the same as in Fig. 9 Shaded areas mark the standard deviations of actual data around the mean curves. Solid lines correspond to a 0.5° binning, dotted lines do to a 2° binning.

4. Conclusion

In this work, we explore the relative importance and the determining factors behind the LQ and TQ effects. The degree of nonlinearity induced by TQ, LQ, and their combination was systematically probed in a grid of surface flux transport (SFT) models. The relative importance of LQ versus TQ has been found to correlate with the ratio u₀/η in the SFT model grid, where u₀ is the meridional flow amplitude and η is diffusivity. An analytical interpretation of this result has been given, further showing that the main underlying parameter is the dynamo effectivity range λ_R, which, in turn, is determined by the ratio of equatorial flow divergence to diffusivity. The relative importance of LQ versus TQ was found to scale as .

As displayed in Fig. 4, for various dynamo and SFT models considered in the literature the contribution of LQ to the amplitude modulation covers a broad range from being insignificant to being the dominant form of feedback. On the other hand, the contribution of a TQ effect (with the usually assumed amplitude) is never negligible.

The role of TQ and LQ was also explored in the 2×2D dynamo model optimized to reproduce the statistical behaviour of real solar cycles. We demonstrated the presence of latitude quenching in the 2×2D dynamo. The locus of the SFT component of the model in Fig. 4 suggests that LQ contributes to the nonlinear modulation by an amount comparable to TQ. This LQ effect is present in the model despite the lack of any modulation in the meridional flow. We note that the meridional inflow module of the 2×2D model developed by Nagy et al. 2020 was not used for these runs. This agrees with the results of Karak (2020) who showed the presence of latitude quenching in the STABLE dynamo model with a steady meridional flow. As LQ is rather well constrained on the observational side (Jiang 2020), its presence and character may potentially be used in the future as a test or merit function in the process of fine-tuning of dynamo models to solar observations.

Acknowledgments

This research was supported by the Hungarian National Research, Development and Innovation Fund (grants no. NKFI K-128384 and TKP2021-NKTA-64), by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 955620, and by the Fonds de Recherche du Québec – Nature et Technologie (Programme de recherche collégiale). The collaboration of the authors was facilitated by support from the International Space Science Institute in ISSI Team 474.

References

All Tables

All Figures

	Fig. 1. Net contribution of a cycle to the solar dipolar moment vs. cycle amplitude for the linear case and for the cases with tilt quenching and/or latitude quenching, with quadratic fits. Parameter values for this run were u0 = 5 m s−1, η = 250 km2 s−1 and τ = 8 years.
In the text

	Fig. 2. Net contribution of a cycle to the solar dipolar moment vs. cycle amplitude for the linear case and for the cases with tilt quenching and/or latitude quenching, with quadratic fits. Parameter values for this run were u0 = 20 m s−1, η = 600 km2 s−1 and τ = ∞ years.
In the text

	Fig. 3. Relative importance of LQ vs. TQ, plotted against the parameter u0/η for the parameter combinations listed in Table 2. Separate quadratic fits for the case τ = ∞ (orange) and τ = 8 years (blue) are shown.
In the text

Fig. 4. Relative importance of LQ vs. TQ, plotted against the dynamo effectivity range, λR, for the parameter combinations listed in Table 2. Separate fits of the form are shown for the case τ = ∞ (orange) and τ = 8 years (blue). Large coloured symbols mark the loci of SFT and dynamo models studied in recent literature: Bhowmik & Nandy (2018) (), Jiang & Cao (2018) (), Lemerle & Charbonneau (2017) (), Whitbread et al. (2018) (+).

In the text

	Fig. 5. Contribution of all active regions in a cycle to the net change in the solar dipole moment during dynamo cycles vs. the cycle amplitude in the 2×2D dynamo model. Solid line: Mean; dashed curves: ±1σ range.
In the text

	Fig. 6. Solar dipole moment at cycle minimum vs. the amplitude of the subsequent cycle in the 2×2D dynamo model.
In the text

	Fig. 7. Cycle amplitude vs. the amplitude of the subsequent cycle in the 2×2D dynamo model. Solid line: Mean; dashed curves: ±1σ range.
In the text

	Fig. 8. Latitude quenching in the 2×2D dynamo model. Mean latitude of BMRs in a cycle vs. cycle amplitude.
In the text

Fig. 9. Butterfly tracks in the 2×2D dynamo model sorted by cycle amplitude: (a) tracks of the active latitude vs. cycle phase; (b) tracks of the active latitude vs. time; Red, black, and blue curves represent flux-weighted mean latitudes for subsamples of strong, average, weak cycles, respectively, as described in the text. Shaded areas mark the standard deviations of the actual data around the mean curves.

In the text

Fig. 10. Universality of the cycle decay phase in the 2×2D dynamo model: (a) Pseudo sunspot number against flux-weighted mean latitude of BMRs. Curves show flux-weighted mean values; (b) Width 2δλ of the latitudinal distribution of BMRs vs. their flux-weighted mean latitude. Colour coding is the same as in Fig. 9 Shaded areas mark the standard deviations of actual data around the mean curves. Solid lines correspond to a 0.5° binning, dotted lines do to a 2° binning.

In the text