© M. Nagy et al., Published by EDP Sciences 2020

1 Introduction

Dynamo models based on the Babcock–Leighton mechanism of dipole regeneration by the decay of bipolar magnetic active regions have become a leading explanatory framework for the solar magnetic activity cycle, well on par with mean-field models based on turbulent induction (see Petrovay & Christensen, 2010; Charbonneau, 2014, for recent reviews). The operation of the Babcock–Leighton mechanism is now very well charted by magnetographic observations (Jiang et al., 2014b; Hathaway & Upton, 2016; Upton & Hathaway, 2018), and the surface magnetic flux so generated is in principle sufficient to account for the magnetic flux emerging as active regions during the following activity cycle (Cameron & Schüssler, 2015). The primary remaining uncertainty, both observationally and theoretically, lies with the ability of the polar field so generated to be subducted into the solar interior, where differential rotation shear induces a toroidal component and presumably closes the dynamo loop. Out-of-ecliptic observations of the solar polar regions, such as those to be provided by the ongoing Solar Orbiter mission and proposed Solaris mission, can hopefully provide the information needed to address this key question.

For a bipolar magnetic region of unsigned magnetic flux Φ emerging at latitude λ, the contribution δD to the surface axial dipole moment is given by

$$\delta D=\frac{3d\cos \lambda }{4\pi R^2}\mathrm{\Phi } \sin \alpha ,$$(1)

where d is the angular separation between the leading and trailing magnetic poles, and α is the tilt angle between the line joining the two poles and the East–West line. Observationally, this angle is known to increase with increasing latitude, so that the leading polarity (with respect to the direction of rotation) is on average closer to the equator than its trailing counterpart. This pattern is known as Joy’s Law, and, at any given latitude, is characterized by a large scatter about the mean value (see Pevtsov et al., 2014, and references therein). This scatter is considered a primary source of stochastic fluctuations in Babcock–Leighton dynamos (see, e.g., Cameron & Schüssler, 2017; Karak & Miesch, 2017; Lemerle & Charbonneau, 2017; Nagy et al., 2017; Bhowmik & Nandy, 2018).

The ultimate contribution to the surface dipole, however, also depends on surface flux transport processes. In particular and somewhat counterintuitively, everything else being equal bipolar magnetic regions emerging close to the equator contribute more to the ultimate dipole (Jiang et al., 2014a; Whitbread et al., 2018), despite the observed increase of the tilt angle α with latitude. This happens because, ultimately, cross-equatorial diffusive cancellation of leading polarity flux is what allows accumulation of trailing polarity flux at the poles (DeVore et al., 1984; Cameron et al., 2013). In addition to this purely kinematic effect, magnetically-mediated variations of surface flows can also alter magnetic flux transport to the pole, and thus alter the dynamo efficiency of decaying bipolar active regions.

Numerical models of dynamos based on the Babcock–Leighton mechanism of poloidal field regeneration have amply demonstrated the ability for exponential growth in the linear regime (see, e.g., Fig. 6 in Lemerle & Charbonneau, 2017, also Kumar et al., 2019), so that a nonlinear backreaction mechanism is required to achieve (statistically) stable cycle amplitudes. Potentially important physical mechanisms range from magnetic backreaction on internal differential rotation and internal meridional flows (e.g. Rempel, 2006; Karak & Choudhuri, 2012; Passos et al., 2012), to various processes mediating surface flux transport (Hathaway & Rightmire, 2010; Jiang et al., 2010), up to the alteration of the emergence latitudes and characteristics of bipolar magnetic regions (Fan, 2009, and references therein). Of these various possible amplitude-limiting nonlinearities, that most often incorporated (or invoked) in extant Babcock–Leighton models of the solar cycle are the so-called tilt quenching, namely the reduction of the tilt angle α in equation (1) with increasing field strength in sunspot-generating toroidal flux ropes, as demonstrated by thin flux tube simulations (Fan, 2009).

Observationally, the signature of these various nonlinearities is very weak at best; the amplitude of solar torsional oscillations amount to a mere 1–2% of differential rotation, and decreases very rapidly with depth (see Gizon, 2004; Howe, 2009, and references therein). The mean poleward meridional flow varies only by 25% between the maximum and minimum phases of the solar cycle or from one cycle to another (Komm et al., 1993; Hathaway & Rightmire, 2010), at least part of that variation being associated with surface inflows towards large active regions (Gizon, 2004; Cameron & Schüssler, 2010; Jiang et al., 2010). These inflows are believed to be driven by horizontal photospheric temperature gradients resulting from localized cooling in active regions (Spruit, 2003; Cameron & Schüssler, 2010, see also Fig. 2 in Martin-Belda & Cameron 2017), though see Petrovay & Forgács-Dajka (2002) for an alternative mechanism. They can reach many tens of percent of the mean axisymmetric poleward meridional flow at mid-latitude, but remain spatially confined to active region belts. Helioseismological measurements over activity cycle 23 indicate that their amplitude drops very rapidly with depth below the photosphere (see Fig. 1 in González Hernández et al., 2010). As a consequence they do not strongly alter the global equator-to-pole advection time, even though they do have a significant impact on surface magnetic flux evolution and the buildup of the surface dipole (Jiang et al., 2010; Martin-Belda & Cameron, 2017).

Tilt quenching is marginally detected in white light and magnetographic data, with the effect more pronounced in some activity cycles than in others and dependent to a significant extent on the data and analysis technique used (see, e.g., Dasi-Espuig et al., 2010; McClintock & Norton, 2013; Pevtsov et al., 2014; Tlatova et al., 2018; Jha et al., 2020). Nonetheless, any one of these weak effects may be the manifestation of nonlinear saturation of the dynamo if the latter only operates in the moderately supercritical regime, which may well be the case (on this point see Cameron & Schüssler, 2015, and references therein). In particular, independent numerical implementations of Babcock–Leighton dynamo models indicate that even mild levels of tilt quenching can readily saturate the cycle amplitude (cf. Karak & Miesch, 2017; Lemerle & Charbonneau, 2017).

The purpose of this paper is to examine the potential role of plasma inflows into active region belts as a nonlinear amplitude-limiting mechanism for the solar magnetic cycle, extending in this respect the calculations of Cameron & Schüssler (2012). More specifically, we adopt the model formulation of Jiang et al. (2010) but instead of imposing fixed latitudinal profiles and inflow speeds, we make the width and speed of these inflows proportional to emerging magnetic flux, thus turning them into a bona fide nonlinear magnetic backreaction mechanism. Working with the Babcock–Leighton 2 × 2D dynamo model of Lemerle & Charbonneau (2017), we study solar cycle simulations including tilt quenching, various level of stochastic fluctuations in tilt angles, in the presence or absence of magnetic flux-dependent inflows of varying magnitude. We focus on the ability – or lack thereof – such inflows may have to stabilize the cycle amplitude in the mildly supercritical regime, in the presence of stochastic perturbations. We also examine to what degree flux-dependent inflows can reduce or enhance the impact of so-called “rogue” active regions (Nagy et al., 2017), and thus decrease or increase the probability of shutting off the dynamo and entering Maunder-Minimum-like phases of suppressed cyclic activity.

The remainder of the paper is structured as follows. In Section 2 we summarize the most salient features of the Lemerle & Charbonneau (2017) Babcock–Leighton solar cycle model used in what follows, and describe our reformulation of the parametrization of inflows introduced by Jiang et al. (2010), so as to include dependencies on the flux of emerging active regions (Sect. 3). Results of our solar cycle simulations are then presented, first with regards to the stabilization effects of inflows on cycle amplitude, duration and cross-hemispheric coupling (Sect. 4), and then on mitigating the shutoff of the dynamo (Sect. 5) in response to stochastic forcing.

2 The 2 × 2D Babcock–Leighton solar cycle model

The Babcock–Leighton dynamo model on which this study is based has been described in detail in a number of recent publications (see Lemerle et al., 2015; Lemerle & Charbonneau, 2017; Nagy et al., 2017), so that only a synopsis of its most salient features is presented here.

Defined in the usual spherical polar coordinates (r, θ, ϕ), the model couples a conventional two-dimensional surface flux transport (SFT) simulation for the photospheric magnetic field B _r (r = R _⊙, θ, ϕ, t)${\widehat{\mathbf{e}}}_r$ to a conventional axisymmetric (∂/∂ϕ = 0) flux transport dynamo (FTD) describing the coupled evolution of the internal toroidal magnetic field component B _ϕ (r, θ, t)$\widehat{\mathbf{e}}$ _ϕ and poloidal component, expressed via a toroidal vector potential A _ϕ (r, θ, t)$\widehat{\mathbf{e}}$ _ϕ .

The FTD simulation drives emergence of bipolar magnetic regions in the SFT simulation, via an emergence function determined by the magnitude and spatial distribution of the interior magnetic field, and setting the probability of bipolar magnetic region emergence per unit time. Physical characteristics of active regions – flux, tilt angle, pole separation, etc. – are drawn from statistical distributions constructed from analysis of magnetographic data by Wang & Sheeley (1989). The SFT simulation then tracks the spatiotemporal evolution of the photospheric field, and attendant buildup of the global dipole. The corresponding zonally-averaged radial surface field is then fed to the FTD via its outer boundary condition at r = R _⊙, closing the dynamo loop. This so-called “2 × 2D” model, so dubbed because it runs concurrently two distinct but coupled two-dimensional simulations, has the advantage of being much faster than full three-dimensional Babcock–Leighton models (see, e.g., Yeates & Muñoz-Jaramillo, 2013; Karak & Miesch, 2017; Kumar et al., 2019), while still incorporating a detailed longitude–latitude representation of the model photosphere. Bhowmik & Nandy (2018) present a conceptually similar dynamo model implementation.

The emergence function is a key component of this dynamo model. It includes a lower field strength threshold below which the emergence probability of bipolar magnetic regions rapidly vanishes. This implies that the resulting dynamo is not self-excited, i.e., it cannot amplify an arbitrarily weak seed magnetic field. The various free parameters defining the emergence function, as well as other defining numerical parameters of the model not tightly constrained by observations, are formally optimized to produce a best-fit to the sunspot butterfly diagram of activity cycle 21 (1976–1986); see Section 3.2 in Lemerle & Charbonneau (2017).

A reduction of the tilt angle α in equation (1) with increasing internal magnetic field strength is the only amplitude-limiting nonlinearity included in the model (Lemerle & Charbonneau, 2017, Sect. 4.2). The following simple and largely ad hoc algebraic formula is used:

$${\alpha }_q=\frac{\alpha \left(\theta \right)}{1+(B_t/B_q{)}^2},$$(2)where B _t is the internal toroidal field strength at the latitude where emergence is triggered, and B _q is a parameter setting the field strength at which significant tilt quenching begins. In all simulations described in what follows we retain the value B _q  = 5 × 10² G used in Lemerle & Charbonneau (2017). Note that this expression is meant to capture the reduction of the tilt angle developing in a flux ropes buoyantly rising through the solar convection zone, as predicted by thin flux tubes simulations for tubes of increasing internal field strength (see Fan, 2009, and references therein), rather than variations driven by surface processes affecting the structure of bipolar active regions after emergence.

3 Implementation of surface inflows into active regions

For the implementation of the meridional inflow we opted to use the axisymmetric parametric form introduced by Jiang et al. (2010). In their implementation a sinusoidal latitudinal inflow profile is used with fixed width and amplitude during the course of their cycles. In their model the only time-dependence results from having the centroid of the inflow belts follow the midpoint of the active latitudes. In order to make the inflow a true nonlinear feedback mechanism, as advocated by Cameron & Schüssler (2010), we make the width, amplitude and latitudinal position of the parametrized inflows depend on the total magnetic flux having emerged in the past interval of 35 days (corresponding to the internal dynamo timestep in the Lemerle & Charbonneau, 2017 model). We do so separately in each solar hemisphere, so that the impact of inflows on hemispheric asymmetry can also be investigated.

The global axisymmetric surface flow is parameterized as:

$$v(t,\lambda )=v\left(\lambda \right)+\Delta v(t,F,\mathrm{\lambda }),$$(3)where λ is the latitude, and v(λ) represents the steady, global poleward surface meridional flow. The quantity Δv(t, F, λ) is the time-dependent surface inflow, defined as:

$$\Delta v\left(t,F,\lambda \right)=v_0\left(t,F\right)\cdot \sin \left(\frac{\lambda -{\lambda }_c\left(t,F\right)}{\Delta \lambda \left(t,F\right)}\right), \mathrm{for}:-180\mathrm{°}\le \frac{\lambda -{\lambda }_c\left(t,F\right)}{\Delta \lambda \left(t,F\right)}<180\mathrm{°}.$$(3a)

In these expressions the flux-weighted quantities are calculated as follows:

$$v_0(t,F)=v_{00}\cdot \arctan \left(\frac{\sum \mathrm{}{F}_i}{F_{00}}\right)$$(3b)

$${\lambda }_c(t,F)=\frac{\sum \mathrm{}({\lambda }_i\cdot F_i)}{\sum \mathrm{}{F}_i}$$(3c)

$$\Delta \lambda (t,F)=\frac{1}{2}\sqrt{\frac{1}{n}\sum \mathrm{}{\left({\lambda }_i-{\lambda }_c(t,F)\right)}^2.}$$(3d)

The parameters involved are the unsigned flux F _i of the ith emerging BMR, its centroid latitude of emergence (λ _i ), and the number n of BMRs emerged within an internal dynamo timestep, all compiled separately for each hemisphere. The other two parameters are the average flux for all emergences, F ₀₀ = 4.99 × 10²¹ Mx, used for normalization, and the amplitude of the inflow, v ₀₀ = −500 (or −750) cm s⁻¹, set to cover the range of observed values. The arctangent dependency in equation (3b) prevents unrealistically high inflow speeds in cases of rogue BMR emergences (Nagy et al., 2017). All flux-dependent quantities are recomputed at the end of an internal dynamo timestep (Δt = 35 d) and kept constant during the subsequent time step. This amounts to setting a flux-dependent “lifetime” for inflows also at 35 days.

Figure 1 shows a representative case of a simulation using the 2 × 2D solar cycle model of Lemerle & Charbonneau (2017), with inflows introduced following the procedure just described, and inflow amplitude v ₀₀ = −500 cm s⁻¹. Panel (A) shows a time series of pseudo-SSN (solid line), panel (B) a butterfly diagram of BMR emergence, panel (C) a time-latitude magnetogram of the surface radial magnetic field, as returned by the SFT module, and panel (D) shows the total latitudinal velocity, as given by equation (3). This simulation segment is chosen because it exhibits an interval of strongly suppressed activity, during which BMR emergences cease completely for two cycles in the Northern hemisphere. Note how during this time interval, despite activity being present only in the Southern hemisphere, a well-defined surface dipole is sustained nonetheless and continues to undergo regular polarity reversals, in response to magnetic flux diffusing from the Southern to the Northern hemisphere, as clearly visible on panel (C). Emergence of BMRs eventually resumes in the Northern hemisphere, but a significant hemispheric asymmetry in both amplitude and phasing still persists for three activity cycles, from t = 2262–2290 yr in simulation time.

Fig. 1 A segment of solar cycle simulation including surface inflows, with amplitude parameter v 00 = −500 cm s−1. Panel (A): pseudo sunspot number (pSSN) timeseries (solid line); panel (B): surface density of BMR emergences, the simulation’s equivalent to the sunspot butterfly diagram; panel (C): time-latitude diagram of the zonally-averaged surface radial magnetic field; panel (D) time-latitude diagram of the surface latitudinal flow speed, including the contribution from inflows into active region belts. The dashed line on panel (A) is the pSSN time series for a simulation where inflows were artificially turned off during the second cycle (2213 ≤ t ≤ 2225 yr) and resuming afterwards, but otherwise identical.

Not surprisingly, the spatiotemporal pattern of latitudinal flow perturbations due to active regions follows closely the density of BMR emergences (cf. panels (B) and (D)). Note however, on panel (D), the brief surges of transequatorial flow taking place in the descending phase of the second cycle (2212–2223 yr), reflecting the emergence of high-flux BMRs very close to the equator. These turn out to have a profound impact in the unfolding of subsequent magnetic cycles.

The dashed line on panel (A) of Figure 1 shows the result of a numerical experiment whereby the inflows were artificially zeroed out altogether throughout the second complete activity cycle (2213–2224 yr). In contrast to the simulation including inflows, cyclic activity now perdures, showing significant fluctuations in amplitude but still sustaining its normal mean amplitude throughout the remainder of the simulation segment plotted. This simple experiment demonstrates quite clearly that the inflows can have a large impact on the long-term cyclic behavior, by all appearances much more important than suggested by earlier investigation of the inflow phenomenon using only surface flux transport simulations (as in, e.g. Jiang et al., 2010; Cameron & Schüssler, 2012; Martin-Belda & Cameron, 2017). As discussed in what follows, this larger impact ultimately results from the workings of the internal solar dynamo process, as captured by the Lemerle & Charbonneau (2017) solar cycle model.

4 Impact of inflows on cycle properties

Past investigations of the impact of inflows into active region belts carried out using surface flux transport simulations (e.g., Jiang et al., 2010; Martin-Belda & Cameron, 2017, and references therein) have identified two somewhat distinct effects altering the buildup of global dipole in the course of a cycle: (1) enhancing the mutual cancellation of leading and trailing polarity fluxes of a BMR, whatever its emerging latitude may be; and (2) driving transequatorial advection of magnetic flux, when BMRs emerge close to the equator (an effect visible on Fig. 1D). The first is a deterministic effect leading to a reduction of the global dipole; the second can either enhance or suppress dipole buildup, depending if the BMRs emerging close to the equator deviate strongly or not from Joy’s Law. This second effect of inflow is thus intimately tied to the stochasticity of BMR emergences.

In order to disentangle the impact of inflows from the strong stochastic fluctuating behavior inherent to the Lemerle & Charbonneau (2017) solar cycle model, we carried out two distinct sets of simulations, which in what follows are dubbed:

• Full stochasticity: this is the original formulation of the Lemerle & Charbonneau (2017) model, in which the flux, centroid latitude, pole separation and tilt angle of emerging BMRs are all extracted from distributions constructed from observations.

• Reduced stochasticity: as previous, except for pole separation and tilt angle, which are fixed at the means of their respective statistical distributions; these two factors having been previously identified as the major contributions to stochastic fluctuations developing in dynamo models of the Babcock–Leighton variety.

We carried out several sets of 300 statistically independent simulations, all with the optimal parameter values determined in Lemerle & Charbonneau (2017), except for the dynamo number being raised slightly to K = 1.0. We ran sets in full and reduced stochasticity modes, as described previously. In addition to a no-inflow reference set, we considered solar-like flow speed parameter v ₀₀ = −500 and −750 cm s⁻¹. In each simulation we delineate individual activity cycles (i.e., half magnetic cycles), and measure (1) their peak amplitude, measured in “pseudo-SSN” (labeled pSSN), a simple count of the number of emergence event per unit time; (2) the cycle duration, based on the pSSN time series; (3) the dipole strength at the end of each cycle. As detailed in what follows, from these primary cycle characteristics we also calculate other secondary quantities of interest, for example the correlation coefficient between pSSN cycle amplitude and dipole strength, and measures of hemispheric asymmetry (as defined in Nagy et al., 2019, and summarized in Table 1 below).

Figure 2 and Tables 1 and 2 summarize our simulations results, partitioned according to stochasticity level. It is clear from Figure 2 that the probability density histograms for cycle amplitude and dipole strength are more strongly affected by the level of stochasticity than they are by the presence and magnitude of inflows. The distributions do shift to slightly smaller values as the inflow magnitude parameter v ₀₀ is increased, as evidenced by the corresponding shift of the distributions means. Tables 1 and 2 indicate that inflows v ₀₀ = −750 cm s⁻¹ lead to a 17–19% drop in cycle amplitude, and 10–11% decrease in dipole moment. The level of stochasticity impacts mostly the scatter of the cycle properties about their mean values, which allows to distinguish the impact of stochasticity from that of surface inflows.

Fig. 2 Frequency histograms for dipole moments on the left panels, and cycle amplitudes on the right. The top row shows the plots for the reduced stochasticity case (Table 1), and the bottom row for full stochasticity (Table 2). The vertical dashed lines indicate the mean values of the different distributions.

Introducing inflows leads to a small but noticeable decrease in the scatter of global cycle properties about their mean values. One may note also another interesting inflow-related change in the histograms of cycle amplitudes, easier to notice in the full stochasticity case (Fig. 2D). The sparse but extended tail of very high amplitude cycles present in the set of no-inflow simulations (blue) disappear upon the introduction of inflows. No such change is observed in the low tail of the distributions. With regards to cycle amplitude, this is arguably the only significant instance of a differential variation in the impact of surface inflows (on this point see also Cameron & Schüssler, 2012). For the range of parameters explored here, neither the cycle duration, nor the correlation between the peak dipole moment and amplitude of the subsequent cycle, change significantly upon introducing surface inflows. It is also important to note that surface inflows do not shift the time delay between the simulated cycle maxima and the peak dipole moment epochs, even if the dipole is defined in terms of the polar cap flux. This indicates that the presence of inflows does not significantly slow down the flux transport towards the poles.

Another potential impact of inflows concerns the hemispheric asymmetry in cycle characteristics. To quantify such effects we use the asymmetry parameters described in Nagy et al. (2019), namely: the polar cap flux asymmetry (ΔΦ), the pSSN amplitude asymmetry (ΔpSSN), and the timelag (ΔT) between the cycle onset in each hemisphere. Figure 3 shows correlation plots between these quantities, for the set of fully stochastic simulations. The corresponding linear correlation coefficient are given in the last two lines of Table 1, and in Table 2 for the reduced stochasticity simulations. The only noteworthy trend emerging from these Tables is an increase of the correlation between cycle onset lag and polar cap flux asymmetry as v ₀₀ increases, more pronounced in the fully stochastic simulations.

Fig. 3 Correlation between hemispheric asymmetry measures in fully stochastic simulations with and without inflows. These quantities are the polar cap flux asymmetry (ΔΦ), the pSSN amplitude asymmetry (ΔpSSN), and the timelag (ΔT) between the cycle onsets in each hemisphere (for details see Nagy et al., 2019). In both panels dashed black line shows the ±1σ range around the linear fit to the asymmetry measures in the no-inflow simulations. Orange dotted (green dot dashed) line shows the same for inflow amplitudes v 00 = −500 (−750) cm s−1. Black crosses indicate cycles outside of the ±1σ range for simulation runs without inflows. Orange asterisks (green triangles) flag the same cycles, but with v 00 = −500 (−750) cm s−1. With inflows included, a large fraction of cycles formerly outside the ±1 standard deviation range are now back within that range (see text).

In order to further investigate the origin of these variations (or lack thereof), we selected the asymmetric cycles in fully stochastic simulations, lying outside of the ±1σ range in the correlation plots of Figure 3, as labeled by “+”. For these cycles we run further tests with inflows turned on with v ₀₀ = −500 cm s⁻¹. According to the results, plotted as orange asterisks in Figure 3, the introduction of surface inflows brings 49% (28 out of 57) of these strongly amplitude-asymmetric cycles back within the ±1σ range, and brings the hemispheric onset lag within the ±1σ range for 79% (23 out of 29) of cycles. With inflow amplitude raised to v ₀₀ = −750 cm s⁻¹ (green triangles), these percent reductions become 43% (25 out of 58) for amplitude asymmetry, and 55% (16 out of 29) for onset lag. The presence of inflows clearly increases significantly the level of cross-hemispheric coupling in our dynamo simulations.

To sum up the results presented in this section: our simulations indicate that the introduction of magnetic flux-dependent inflows into activity belts has a modest but significant stabilizing effect on cycle amplitude, as compared to no-inflow simulations. This stabilizing effect is also observed and more prominent in various measures of hemispheric asymmetry, most notably perhaps in the hemispheric lag in cycle onset. Part of the stabilizing effects results from the suppression of most “outlier” cycles, characterized by amplitude largely in excess of the simulation mean, when inflows are turned on.

5 Impact of inflows on the onset of Grand Minima

As discussed in Section 4.3 of Lemerle & Charbonneau (2017) and Section 4.1 of Nagy et al. (2017), the occasional emergence of a large “rogue” active region can shut down cyclic activity. This happens when the “rogue” reduces the buildup of the dipole to an extent such that an internal toroidal magnetic field strong enough to exceed the emergence threshold can no longer be produced by differential rotation shear on the cycle timescale. The presence of this threshold also implies that the 2 × 2D solar cycle model of Lemerle & Charbonneau (2017) is not self-excited, so that once the cycle shuts off a secondary dynamo process, such as turbulent induction, is required to restart the cycle. The viability of such a Grand Minima scenario has been already explored in various types of dynamo models, and for a variey of “kickstart” inductive mechanisms (see, e.g. Schmitt et al., 1996; Ossendrijver, 2000; Charbonneau et al., 2004; Hazra et al., 2014; Passos et al., 2014), including in the present dynamo model (Ölçek et al., 2019).

Inflows into active regions can potentially impact this probability of shutdown; as documented in Section 4 (viz. Tables 1 and 2), including inflows tends to produce on average weaker dipole moments at the end of the cycle, which then make it easier for the dynamo to shut down in response to a large fluctuation. On the other hand, inflows can also reduce the “anomalous” transequatorial magnetic flux associated with the emergence, close to the equator, of a large active regions strongly deviating from Joy’s Law; this should reduce the negative impact of such a “rogue” emergence on dipole buildup, and therefore help to sustain normal cyclic activity.

In our dynamo model at least, the first of these effects turns out to dominate over the second. Two ensembles of 100 simulations were assembled, operating in full stochasticity mode and using inflow amplitude parameters v ₀₀ = −500 and −750 cm s⁻¹, together with a third “reference” 100-member ensemble, with inflows turned off. Each simulation in each ensemble was run over a timespan of 100 mean cycle duration. For each simulation in each ensemble, we monitored whether a simulation underwent shutdown, and if so, after how many cycles. We then constructed the cumulative probability density histograms plotted on Figure 4, where each bin gives the fraction of simulation within each ensemble still running at normal amplitude after a number cycles given by the bin range.

Fig. 4 Cumulative probability density histograms showing the fraction of ensemble members still running at normal cycle amplitude after N cycles. These are constructed from 100-members ensembles of simulations with inflows amplitude v 00 = 0, −500 and −750 cm s−1 (blue, orange and yellow as labeled). All simulations in each ensemble use a dynamo number K = 1.0 and full stochasticity in the properties of emerging BMRs.

The three cumulative histograms differ markedly, solutions including inflows showing a much greater probability of shutdown than in their absence, the more so the larger the inflow amplitude v ₀₀. The effect is quite pronounced: here 43% of simulations without inflows are still operating after 100 cycles, while only 4% of simulations including inflows at v ₀₀ = −500 cm s⁻¹ are still doing so, and all simulations at v ₀₀ = −750 cm s⁻¹ have stopped before 90 cycles. In all three cases the approximately exponential falloff of the cumulative histograms is consistent with the shutdown being triggered by a random process, here the emergence of one or more “rogue” active regions close to the equator (see Sect. 4.1 in Nagy et al., 2017, for discussion of this effect). No simulation operating in reduced stochasticity mode underwent shutdown before the 100th cycle, whether at v ₀₀ = 0, −500 or −750 cm s⁻¹. Qualitatively similar results are obtained upon varying the dynamo number within the range yielding solar-like butterfly diagrams.

The fact that surface inflows towards activity belts stabilize global cycle measures (viz. Sect. 4) but enhance the probability of dynamo shutdown may appear contradictory, but it is not. The former is primarily a surface flux transport effect, and can be readily captured by SFT simulations (such as, e.g., Jiang et al., 2010; Martin-Belda & Cameron, 2017); the latter is a dynamo effect associated with the field strength threshold below which the internal toroidal field can no longer generate emerging BMRs showing the systematic pattern of E–W tilt required to regenerate the surface dipole. Thus, the strong nonlinearity associated with this lower threshold acts as an amplifier of the small reduction of the surface dipole moment caused by the inflows.

6 Conclusions

In this paper we presented sets of solar cycle simulations using the Lemerle & Charbonneau (2017) Babcock–Leighton dynamo model, including an axisymmetric formulation for surface inflows into activity belts. Inspired by the surface flux transport simulations presented in Jiang et al. (2010), our results improved on the realism of these simulations in two ways: (1) the spatiotemporal profile of inflows is set by the magnetic flux emerging in the SFT module of the simulation, thus yielding a truly nonlinear backreaction mechanism; and (2) our simulations use a full dynamo model, allowing to trace in an internally consistent manner the nonlinear feedback of inflows on the development of subsequent cyclic activity.

Our simulations indicate that the introduction of axisymmetric inflows into active regions belts reduces the amplitude of the surface dipole ultimately generated at the end of each cycle, the more so the stronger those inflows; and consequently the amplitude of subsequent cycles is also reduced. We find a reduction in cycle amplitude, as measured by our pseudo-SSN, by 12% for inflow amplitude v ₀₀ = 500 cm s⁻¹, approaching 20% at v ₀₀ = 750 cm s⁻¹. This is comparable to the reduction estimated by Jiang et al. (2010) on the basis of their SFT simulation, but significantly smaller than the 30% reported in Martin-Belda & Cameron (2017). The inflows also produce a modest but noticeable decrease in the scatter of global cycle characteristics about their mean values. This is in part, but not exclusively, due to the fact that the introduction of inflows with significant amplitudes tends to suppress “outlier” cycles of very high amplitude (i.e., more than twice the simulation mean). Another significant and related impact of inflows is a reduction in the degree of hemispheric amplitude asymmetry and asynchrony in cycle onset, as determined from our pseudo-SSN time series. Overall, magnetic flux-dependent inflows into active region belts thus have a significant stabilizing effect on global cyclic behavior.

Dynamos of the Babcock–Leighton variety are not self-excited, in the sense that they cannot amplify an arbitrarily weak seed magnetic field; thin flux tube simulations indicate that a minimal magnetisation of a few tens of kG is necessary for toroidal flux ropes to survive their buoyant rise through the turbulent solar convection zone and emerge with a systematic pattern of tilt resembling Joy’s Law (see Fan, 2009, and references therein). The solar cycle model used in this paper includes such a threshold, and its presence turns out to act as an amplifier for the effect of inflows. Even though inflows only reduce the mean cycle amplitude by a few tens of percent compared to simulations without inflows, the results presented in Section 5 demonstrate that their presence greatly increases the probability of the cycle shutting down following an unfavorable sequence of BMR emergences. Our results, taken at face value, thus indicate that surface inflows into active regions belts can have a strong impact on supra-cycle timescales.

All results reported upon in this paper include a parameterized tilt-quenching, the only upper amplitude-limiting mechanism included in the Lemerle & Charbonneau (2017) solar cycle model used as a starting point of our analyses. We also carried out simulation with tilt quenching turned off, to assess the ability of inflows into activity belts to act as a cycle amplitude saturation mechanism. Is is important to distinguish here the gradual decrease of the tilt angle occuring after emergence under the action of the surface inflow (González Hernández et al., 2010), from the tilt quenching embodied in equation (2), which reflects a field-strength-dependent physical process taking place prior to emergence, as toroidal flux ropes rise buoyantly across the solar convection zone (Fan, 2009; Dasi-Espuig et al., 2010). In the range of model parameters yielding solar-like butterfly diagrams, inflows in themselves could not achieve this amplitude-limiting effect, even though for the dynamo number used in all our simulations the dynamo operates close to criticality. This is very much a model-dependent result, tied in particular to our adopted prescription for the inflow parameterization as a function of magnetic flux in emerging BMRs (Eqs. (3)–(3d) above), and clearly requires further investigation. Likewise, the 2D latitude-longitude representation of the solar photosphere also make it possible to treat active region inflows as non-axisymmetric flow patterns, instead of the longitudinally-averaged representation used in this paper. The simulation results reported upon here motivate the further development of the inflow model along this direction.

Acknowledgments

This research was supported by the Hungarian National Research, Development and Innovation Office (NKFI grant no. K128384) of Hungary, by the Natural Sciences and Engineering Research Council of Canada (Discovery Grant Program), by the Fonds de Recherche de Recherche du Québec – Nature et Technologie (Programme de recherche collégiale), and by Campus Mundi Scholarship Program (Tempus Public Foundation). The collaboration of the authors was facilitated by support from the International Space Science Institute in ISSI Team 474. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

Fig. 1 A segment of solar cycle simulation including surface inflows, with amplitude parameter v 00 = −500 cm s−1. Panel (A): pseudo sunspot number (pSSN) timeseries (solid line); panel (B): surface density of BMR emergences, the simulation’s equivalent to the sunspot butterfly diagram; panel (C): time-latitude diagram of the zonally-averaged surface radial magnetic field; panel (D) time-latitude diagram of the surface latitudinal flow speed, including the contribution from inflows into active region belts. The dashed line on panel (A) is the pSSN time series for a simulation where inflows were artificially turned off during the second cycle (2213 ≤ t ≤ 2225 yr) and resuming afterwards, but otherwise identical.

In the text

	Fig. 2 Frequency histograms for dipole moments on the left panels, and cycle amplitudes on the right. The top row shows the plots for the reduced stochasticity case (Table 1), and the bottom row for full stochasticity (Table 2). The vertical dashed lines indicate the mean values of the different distributions.
In the text

Fig. 3 Correlation between hemispheric asymmetry measures in fully stochastic simulations with and without inflows. These quantities are the polar cap flux asymmetry (ΔΦ), the pSSN amplitude asymmetry (ΔpSSN), and the timelag (ΔT) between the cycle onsets in each hemisphere (for details see Nagy et al., 2019). In both panels dashed black line shows the ±1σ range around the linear fit to the asymmetry measures in the no-inflow simulations. Orange dotted (green dot dashed) line shows the same for inflow amplitudes v 00 = −500 (−750) cm s−1. Black crosses indicate cycles outside of the ±1σ range for simulation runs without inflows. Orange asterisks (green triangles) flag the same cycles, but with v 00 = −500 (−750) cm s−1. With inflows included, a large fraction of cycles formerly outside the ±1 standard deviation range are now back within that range (see text).

In the text

Fig. 4 Cumulative probability density histograms showing the fraction of ensemble members still running at normal cycle amplitude after N cycles. These are constructed from 100-members ensembles of simulations with inflows amplitude v 00 = 0, −500 and −750 cm s−1 (blue, orange and yellow as labeled). All simulations in each ensemble use a dynamo number K = 1.0 and full stochasticity in the properties of emerging BMRs.

In the text