© ESO 2018

1. Introduction

The magnetic solar cycle is attributed to a dynamo process in which motions of a conductive medium lead to the continuous induction of magnetic fields. These motions may arise from convection, which becomes anisotropic in the presence of rotation and stratification and is approximately described by the α-effect (Krause & Rädler 1980). Alternatively, they may be caused by the rise of magnetic structures to the solar surface, again in the presence of rotation and stratification, which leads to the so-called Babcock-Leighton effect (Leighton 1969). While there has been no derivation of this effect from first principles so far, the synergy with the meridional circulation in the convection zone may lead to a periodic magnetic field that satisfactorily explains the solar cycle.

Babcock-Leighton dynamos have been shown to qualitatively reproduce the solar butterfly diagram only if the turbulent magnetic diffusivity in the bulk of the convection zone is lower than 10¹² cm² s⁻¹. This regime is often referred to as the advective regime, since the meridional circulation then determines the cycle period. In the diffusive regime, the cycle period varies with the turbulent magnetic diffusivity and can be reduced to the observed cycle length, but the propagation of the dynamo wave then follows the Parker–Yoshimura rule and is poleward at all latitudes (Yoshimura 1975).

A variety of Babcock-Leighton dynamos has been published in the past 25 years, exemplarily by Choudhuri et al. (1995), Dikpati & Charbonneau (1999), Küker et al. (2001), Chatterjee et al. (2004), Guerrero & deGouveiaDallPino(2007), and Sanchez et al. (2014) for kinematic models with stationary flows and dynamo effect, by Nandy & Choudhuri (2001) for a model with toroidal-field loss by buoyancy, and by Kitchatinov & Olemskoy (2011) for a nonlocal α-effect similar to the Babcock-Leighton effect. A nonkinematic simulation showing a solar-like butterfly diagram has been reported by Rempel (2006); here, the turbulent magnetic diffusivity rises from 10¹⁰ cm² s⁻¹ to 10¹² cm² s⁻¹ in the convection zone. Beyond the Sun, Jouve et al. (2010a), for example, studied the cycle period dependence of Babcock-Leighton dynamos on the stellar rotation rate. The variability of the solar cycle including grand minima has been addressed with Babcock-Leighton dynamos for example by Karak (2010), who varied the meridional circulation, Olemskoy & Kitchatinov (2013), who varied the strength of their nonlocal source term, and Inceoglu et al. (2017), who varied both the generation of the flows and the Babcock-Leighton effect in a nonkinematic setup. In all these attempts, however, the effect of the toroidal magnetic field in the interior of the convection zone on the poloidal-field generation at the surface is instantaneous. We note that this list of papers on the Babcock-Leighton type dynamo is far from complete.

We here address the problem that the Babcock-Leighton effect is not only nonlocal in space, but also in time. Following the pioneering work by Jouve et al. (2010b), we take a further step toward a fully constrained Babcock-Leighton dynamo. Here we use the results of global numerical simulations of flux-tube rise to improve and even constrain a Babcock-Leighton dynamo model.

The global numerical simulations by Fournier et al. (2017) have shown that the rise time of magnetic flux tubes is independent of the magnetic diffusion. We suggest treating this independence with a nonlocality in space and time by designing a Babcock-Leighton dynamo model that has a time delay in the source term.

Nonlocalities have been shown to generate long-term variability in the amplitude and in the period of the magnetic cycle for a wide variety of models. A remarkable early attempt has been made by Yoshimura (1978), who used a 2D setup with source terms nonlocal in both space and time, delivering cyclic magnetic fields interrupted by “low-activity” periods that lasted for a few cycles. A sequence of papers on a 0D model, that is, without any spatial dependence but with a time delay, was published by Wilmot-Smith et al. (2006), and Hazra et al. (2014) (including stochastic variations in the Babcock-Leighton effect), and Tripathi et al. (2018), who found sub-critical dynamo action. This indicates what we show here in a 2D spherical shell with realistic solar differential rotation. It is also in line with the result by Rheinhardt & Brandenburg (2012), who studied the memory effect in a turbulent-α dynamo and found a threshold for dynamo excitation that is lower than without a memory effect. Nonlocalities in turbulent-α dynamos have been studied in various other works that we do not review here. Closest to our approach is the work by Jouve et al. (2010b), who made the time delay magnetic-field dependent. However, the authors always found that the time correlations required to obtain solar-like variability were too long to agree with the model assumptions.

Here we study such nonlocal models with time delays depending on the magnetic field strength, where the dependence is derived from the flux-rise simulations by Fournier et al. (2017). Sect. 2 describes the equations we used, and Sect. 4 shows the general results of using a delay in the induction equation. We discuss the implications of the results in Sect. 6.

2. Model and reference setup

We derived our model from Jouve et al. (2008), and like in Jouve et al. (2010b), we introduced a delay into the Babcock-Leighton source term of a mean-field dynamo. We here present the results of a large parameter study of more than 2000 simulations, which are partly constrained through the results of global numerical simulations.

The model we used here is simplified: we assumed a constant turbulent magnetic diffusivity, η_t = const, in the stellar interior. Since the turbulent diffusivity is a measure of the turbulence intensity, a constant value implies that we can neglect radial turbulent pumping (Rädler 1968). This choice was made to prevent our setup from being polluted with weakly constrained parameters, such as an arbitrary profile of the magnetic diffusivity. Observations of decaying active regions suggest that the turbulent magnetic diffusivity, η_t, at the surface is on the order of 10¹² cm² s⁻¹, while the mixing length theory provides an upper limit of 10¹⁴cm² s⁻¹ at the surface and 10¹³ cm² s⁻¹ at the bottom of the convection zone. The large discrepancy between the various available estimates demonstrates the limit of our current knowledge. Additionally, a constant turbulent magnetic diffusivity has the appreciable side effect that the Reynolds number is small everywhere, allowing coarser grids and shorter computation times, which is adequate for a large parameter study.

Lengths and time were normalized with the stellar radius, R_⋆, and the turbulent magnetic diffusion time, ${\tau }_{\mathrm{diff}}=R_{\star }^2/{\eta }_{\mathrm{t}}$. The free choice of the unit of the magnetic flux density B was made by setting the equipartition value with convective motions u_rms at r = 0.71, $B_{\mathrm{eq}}=u_{\mathrm{rms}}\sqrt{{\mu }_0\rho }$ to unity, where ρ and μ₀ are the gas density and the permeability constant, respectively (see below). The resulting dimensionless set of equations can be written in spherical coordinates (r, θ, ϕ) as

$$\begin{array}{ll}{\partial }_t{B}_{\varphi } =\hfill & \left({\nabla }^2-\frac{1}{{\varpi }^2}\right)B_{\varphi }\hfill \\ \hfill & - Re\varpi u_{\mathrm{P}}\cdot \nabla \left(\frac{B_{\varphi }}{\varpi }\right) - Re(B_{\varphi }\nabla )\cdot u_{\mathrm{P}}\hfill \\ \hfill & + C_{\Omega }\varpi \left[\nabla \times (\varpi A_{\varphi }{e}_{\varphi })\right]\cdot \nabla (\Omega ),\hfill \\ {\partial }_t{A}_{\varphi }=\hfill & \left({\nabla }^2-\frac{1}{{\varpi }^2}\right)A_{\varphi }\hfill \\ \hfill & - Re\frac{u_{\mathrm{P}}}{\varpi }\cdot \nabla (\varpi A_{\varphi }) + C_{S}S,\hfill \end{array}$$(1)

with B_ϕ, u_P, ϖ, A_ϕ, Ω and S, being the azimuthal magnetic flux density, the meridional circulation profile, the cylindrical distance to the rotation axis, the vector potential of the poloidal magnetic field, the angular velocity, and the source term for the poloidal field, respectively.

This system is controlled by three dimensionless parameters, the Reynolds numbers for the rotation and the meridional flow, C_Ω and Re, and the dynamo number C_S,

$$Re=\frac{u_0{R}_{\star }}{{\eta }_{\mathrm{t}}}; C_{\Omega }=\frac{{\Omega }_0{R}_{\star }^2}{{\eta }_{\mathrm{t}}}; C_S=\frac{S_{\mathrm{}0}{R}_{\star }}{B_{\mathrm{eq}}{\eta }_{\mathrm{t}}},$$(2)

with u₀, Ω₀, and S₀ being the maximum meridional velocity, angular velocity, and Babcock-Leighton effect, respectively. The profiles of u_P and Ω are normalized by u₀ and Ω₀, respectively.

The loose constraint on the amplitude and the profile of the meridional circulation gives some freedom for the choice of Re. The original model used a canonical profile aimed at representing the general characteristics of the solar flows. However, the profile remained arbitrary, and the surface shear layer was missing. We suggest to take advantage of the results of Küker & Rüdiger (2011): the authors provide consistent profiles of differential rotation and meridional flows based on the Λ-effect theory (Rüdiger 1989). This theoretical result provides a solar-like profile including a surface shear layer, with a single free coefficient whose value is set to fit the helioseismic observations. The differential rotation together with the meridional circulation profiles and amplitudes are illustrated in Fig. 1.

Fig. 1. Solar-like differential rotation and meridional circulation obtained from a mean-field hydrodynamic simulation (Küker & Rüdiger 2011).

In Babcock-Leighton dynamos, the source term S is based on phenomenological arguments. Observations suggest that the reversed polarity of the poloidal field emerging at the surface leads to a reversal of the large-scale polar magnetic field. The source term S is an attempt to represent the physical processes that generate the poloidal field from the deeply seated toroidal field.

In the traditional frame of the Babcock-Leighton effect, buoyancy transports toroidal magnetic flux to the surface for sufficiently strong magnetic fields. During its rise through the convection zone, a poloidal field can be generated under the action of the Coriolis effect, which then provides the necessary field at the surface for the Babcock-Leighton effect to work. Global simulations were able to show that buoyant magnetic structures locally quench the magnetic turbulent diffusivity, η_t, allowing them to remain coherent throughout their rise (Cattaneo & Hughes 1988). Additionally, Fan et al. (1994) showed that the solar meridional flow does not significantly affect the rise of magnetic flux tubes. Therefore, the buoyant transport of magnetic flux is independent of η_t and of the meridional flow, u_P, and exclusively depends on the magnetic pressure, that is, on B.

In the model presented here, we treat the transport independence of η_t as a memory effect in the mean-field equations. The source of the poloidal field at the surface is correlated with the toroidal field at the bottom of the convection zone at an earlier time. The considered time correlation, τ_delay, represents the time required by a buoyant magnetic structure to rise from the bottom of the convection zone to the surface.

We write the delayed source term as

$$\begin{array}{l}\begin{array}{ll}S(r,\theta ,t)=\hfill & f(r,\theta ){\displaystyle \sum _{{\tau }_i}\{{\left[1+{\left(B_{\varphi }(r_0,\theta ,t-{\tau }_i)/B_{\mathrm{quench}}\right)}^2\right]}^{-1}}\hfill \\ \hfill & B_{\varphi }(r_0,\theta ,t-{\tau }_i)\},\hfill \end{array}\\ \mathrm{for}\mathrm{\hspace{1em}}{B}_{\varphi }(r_0,\theta ,t-{\tau }_i)>B_{\mathrm{threshold}}, S=0 \mathrm{otherwise}.\end{array}$$(3)

Here the sum is computed over all magnetic flux tubes that reach the surface at a given time t.

On the one hand, tubes with magnetic flux densities higher than B_quench are weakly affected by the Coriolis effect and emerge as untilted active regions (D’Silva & Choudhuri 1993; Jouve et al. 2013; Fournier et al. 2017). Such untilted regions do not provide poloidal magnetic flux to the Babcock-Leighton effect. On the other hand, weak flux tubes may be strongly affected by the turbulent convective motions such that they do not reach the surface. We consider a threshold, B_threshold, below which the toroidal field does not participate in the dynamo mechanism. Such a lower limit prevents the dynamo from growing from an arbitrarily low seed field.

Since the place where the source S operates remains unknown, we used the same arbitrary profile as Jouve et al. (2008), where S is nonzero at r ≳ 0.9 and is maximum at 45° latitude at the surface,

$$f(r,\theta )=\frac{1}{2} \left[1+\mathrm{erf}\left(\frac{r-0.9}{0.02}\right)\right]\cos \theta \sin \theta .$$(4)

2.1. Introducing the 3D results into the model

Global simulations have enabled us to constrain the correlation time, τ_delay, the source-term quenching, B_quench, and the magnetic field threshold, B_threshold. Here τ_delay represents the rise time of magnetic flux tubes, that is, coherent magnetic structures that may form from the destabilization of a previously amplified magnetic layer (Rempel & Schüssler 2001; Hotta et al. 2012) by the buoyancy instability (Parker 1955; Matthews et al. 1995; Wissink et al. 2000; Fan 2001; Kersalé et al. 2007; Favier et al. 2012). The amplification only depends on the stratification of the solar convection zone and has been found to be F_amp = 10 (Rempel & Schüssler 2001; Hotta et al. 2012). As a result, amplified flux tubes can reach up to 15 × 10⁴ G (i.e., 10 B_eq).

Fournier et al. (2017) demonstrated that the rise time of magnetic flux tubes follows the relation

$${\tau }_{\mathrm{delay}}\propto P_{\mathrm{rot}}{\left({\Gamma }_{{\alpha }_1}^{{\alpha }_2}\right)}^{{\alpha }_3} ,$$(5)

where P_rot is the rotation period of the star, and ${\Gamma }_{{\alpha }_1}^{{\alpha }_2}$ is the ratio between the buoyant force and the Coriolis effect, modified by the magnetic tension. The exponent α₃ is a function of the azimuthal mode number with which the magnetic flux tube rises. In the case of the Sun, taking P_rot constant, this relation reduces to

$${\tau }_{\mathrm{delay}}\propto {\tau }_0{\left|\frac{B_{\varphi }}{B_{\mathrm{eq}}}\right|}^{\alpha } ,$$(6)

with α varying between −0.91 and −2.0 depending on the azimuthal mode. The parameter τ₀ is the rise time for a field in equipartition with the convective velocity, B_ϕ = B_eq , and is a free parameter. It depends on the rotation period of the star, the depth of its convection zone, and the details of the destabilization process forming the flux tubes as well as on the profiles of the turbulent thermal conductivity, viscosity, and diffusivity.

In Eq. (6) there is no latitudinal dependence. Recalling that this equation describes the reduction of the rise velocity by the tension force, however, we can model it following the latitudinal dependence of the tension force of 1/sin θ,

$${\tau }_{\mathrm{delay}}=\frac{{\tau }_0}{\sin \theta }{\left|\frac{B_{\varphi }}{B_{\mathrm{eq}}}\right|}^{\alpha }.$$(7)

We further define an effective delay τ_eff, which is a very useful parameter. It corresponds to the shortest delay in any given moment, that is,

$${\tau }_{\mathrm{eff}}={\tau }_0{\left|\frac{B_{\varphi }^{\max }}{B_{\mathrm{eq}}}\right|}^{\alpha },$$(8)

where $B_{\varphi }^{\max }$ is the strongest field strength generated by the Ω-effect at the bottom of the convection zone. The maximum is taken over time and latitude (after saturation).

Rempel (2006) showed that Babcock-Leighton dynamos provide B_ϕ up to 3B_eq. Since the resulting amplified flux tubes are in the buoyancy-dominated regime and therefore weakly participate in the dynamo because the Coriolis effect is not strong enough for significantly tilted active regions (see the appendix), we set B_quench = 3B_eq. Fan et al. (2003) have shown that amplified flux tubes weaker than 3B_eq will not reach the surface as coherent structures. The lower threshold on B_ϕ is B_threshold = 3B_eq/F_amp = 0.3B_eq. Since magnetic flux tubes of threshold flux density may have rise times on the order of a cycle, we implemented a source term buffer that lasted until t + 5 τ_diff for any moment t in the simulation, in order to prevent any issues with very long delays.

Fig. 2. Zero-dimensional toy model illustrating the effect of the delay. In each panel the green curve represents the evolution of the delay with time for a given Bϕ, indicated in the lower panels as the blue curve. The dashed line represents time. The histogram on the right-hand side represents the resulting source term S normalized in arbitrary units. Panels a and b illustrate the case of a weakly nonlinear delay (α = −1), and panels c and d show a nonlinear delay (α = −2). Panels a and c show a short effective delay (e.g., $B_{\varphi }^{\max }=B_{\mathrm{eq}}$ and τ0 = 3 · 10−5 Pcyc) and panels b and d show a long effective delay (e.g., $B_{\varphi }^{\max }=0.5B_{\mathrm{eq}}$ and τ0 = 10−3 Pcyc).

Fig. 3. Color-coding of the toroidal magnetic field strengths of an artificial oscillatory dipole (top panel) and an artificial oscillatory quadrupole (bottom panel). The contours represent the source term (S) at the surface for α = −2 and τ0 = 10−2Pcyc.

2.2. Behavior of the delay

The dependence of the delay on B_ϕ is controlled by α. When α is set to zero, the delay is constant. A constantly delayed source term has the same time-profile as a non-delayed source term, but is shifted in time. However, as soon as α becomes nonzero, the delay becomes time dependent. Weaker flux tubes rise longer than stronger ones. Weak flux tubes may therefore reach the surface at the same time as stronger flux tubes that formed at a later time. The time-dependence of the delay results in an accumulation of the source at the surface at certain times.

We illustrate this behavior in Fig. 2 for four different idealized cases. We consider a single point with a sinusoidally varying toroidal magnetic field,

$$B_{\varphi }=B_{\varphi }^{\max }\sin \left(\omega t\right),$$(9)

with amplitude $B_{\varphi }^{\max }$ and frequency ω. The field is shown as blue curves in Fig. 2, while the red curves are the delay computed from B_ϕ. The histograms show the resulting source of poloidal field accumulated from B_ϕ in the past. The two upper panels show the case of a weakly nonlinear short- and a long-τ_eff, respectively (with α = −1). The two lower panels illustrate a strongly nonlinear short- and long-τ_eff (with α = −2). In the weakly nonlinear short-τ_eff case, the accumulation seen in the diagram is almost negligible and the time-profile of the source looks almost identical to the profile of B_ϕ. Only for a longer τ_eff does the accumulation become visible. A similar accumulation is found for the nonlinear short-τ_eff, demonstrating that the nonlinearity support the accumulation. In the case of the long-τ_eff , the deformation of the profile leads to a front.

Clearly, two aspects of the delay lead to accumulation: the length of the effective delay, τ_eff, controlled by τ₀ and $B_{\varphi }^{\max }$, and its nonlinearity, controlled by α.

Because of the latitudinal dependence of the toroidal magnetic field, the delay naturally varies with latitude. Like its time dependence, the latitudinal dependence of the delay is a function of α and τ₀, but it additionally depends on the sign of the latitudinal gradient of the toroidal field, ∂B_ϕ/∂θ. Figure 3 shows that if the toroidal field decreases toward the equator, the resulting delayed field, illustrated by the contour plots, migrates equatorward, whereas when the gradient is directed poleward, the dynamo wave propagates poleward because weaker fields rise on a longer timescale and the accumulation is delayed, which sets the direction of propagation with the latitudinal gradient of B_ϕ.

Fig. 4. Radial magnetic field at the surface as a function of time from the solutions in the advection-dominated regime (ADV, top panel) and the diffusion-dominated regime (DIFF, bottom panel). Contours also represent Br, but show the 15, 1.5, and 0.15% levels of $B_r^{\max }$ for the top panel (ADV) and the 50, 25, and 10% levels of $B_r^{\max }$ for the bottom panel (DIFF). The parameters are summarized in Table 1.

3. Nondelayed model

We solve Eq. (1) with the pseudo-spectral spherical code by Hollerbach (2000), in which the diffusion term is solved in spectral space, but the induction term is solved in real space. We do not employ the momentum and temperature equations of the code here. The induction equation solution took part in the benchmark by Jouve et al. (2008), and after implementing the delay term for the present study, we confirmed their results for τ_delay = 0, as shown in Fig. 4. Our constant turbulent diffusivity of η_t = 10¹¹ cm² s⁻¹ as compared to the radius-dependent diffusivity in the benchmark, and our theory-based differential rotation versus the closed approximation in the benchmark do not modify the overall picture of the solutions. The slight difference in the amplitude of the meridional flow, due to the theoretical profile, explains the longer magnetic cycle of 40 yr compared to 30 yr of Jouve et al. (2010b). The solution is antisymmetric and oscillatory, with concentrated strong polar regions and rather weak magnetic fields at low latitudes, which are about 10⁻³ of the high-latitude fields, migrating equatorward. The model is labeled ADV in Table 1 and falls short of producing strong enough low-latitude toroidal fields for a realistic solar butterfly diagram. Another issue is the cycle, which remains too long.

This solution is located in the advection-dominated regime, where the magnetic cycle mostly depends on the meridional circulation. The meridional circulation amplitude is defined as a byproduct of the Λ-effect that reproduces the solar differential rotation and only requires the mixing length parameter α_MLT. This means that the only possibility to decrease the magnetic cycle in this setup is to increase the turbulent magnetic diffusion toward the diffusion-dominated regime.

It is well known that dynamo solutions in the advection-dominated regime differ from those in the diffusion-dominated regime. In the lower panel of Fig. 4, we show a model in the diffusion-dominated regime, with η_t = 10¹² cm² s⁻¹. The relevant parameters are listed in Table 1 under the label DIFF. The solution remains antisymmetric and oscillatory, with an activity-cycle period of 8 yr. The polar regions are less concentrated, more similar to the solar characteristics, and the low latitudes exhibit stronger fields (≈10% of the high-latitude fields), but the dynamo wave propagates purely radially, as predicted by the Parker–Yoshimura rule.

Even if an η_t exists for which the nondelayed model gives a solution with an 11-year cycle period, the low-latitude radial fields remain weaker than the polar regions, and the migration becomes radial while moving to the diffusion-dominated regime.

4. Recovering the solar characteristics with the delay model

We modeled the rise time as a temporal nonlocality, the delay, because Fournier et al. (2017) have shown in global simulations of rising magnetic flux tubes that the rise time is independent of the turbulent magnetic diffusion. In Sect. 2.2 we have shown that a time-dependent delay leads to temporal peaks in the source term S of the induction equation near the surface, migrating latitudinally in the direction of decreasing toroidal field. We show below that the solar characteristics can only be recovered in the diffusive regime (model D-DIFF). We present one series of simulations in the advection-dominated regime (D-ADV) and another one in the diffusion-dominated regime (D-DIFF). In both series, C_S and τ₀ are varied. Since we found that the accumulation is largest in the nonlinear regime, we present two series with α = −2. The parameters Re and C_Ω are both fixed by the differential rotation profile, whose shape is determined by the Λ-effect and a standard solar model. B_quench and B_threshold are both constrained by the results of Fournier et al. (2017) and F_amp by Hotta et al. (2012). The two remaining parameters are C_S and τ₀. They are not yet constrained by any global simulation results. We have seen that τ₀ may determine the variability of S, and C_S may determine its amplitude. The parameters are summarized in Table 1.

Fig. 5. D-ADV case. Saturation field strength (left panel), effective delay (middle panel), and cycle period (right panel) as functions of the source term amplitude CS for four different delays. Each line represents a series of runs for a given delay. The vertical solid line marks the critical $C_{S0}^{\mathrm{crit}}$ for the nondelayed case. The circles identify the solutions shown in Fig. 6.

4.1. Advection-dominated case

In this section we consider the ADV setup in the advection-dominated regime, with η_t = 10¹¹ cm s⁻¹. In the left panel of Fig. 5, we illustrate the maximum of the toroidal magnetic field at the bottom of the convection zone, $B_{\varphi }^{\mathrm{cz}}$, against C_S. The nondelayed series with τ₀ = 0 shows a critical source term amplitude of $C_{S0}^{\mathrm{crit}}=3.534$ below which the model produces decaying solutions. It is remarkable that delayed dynamos deliver nondecaying solutions for weaker C_S than $C_{S0}^{\mathrm{crit}}$. The delay has the surprising effect of reducing the criticality of the dynamo, opening a window to unknown solutions.

These solutions have the particular property of self-saturating at low amplitudes. The nonlinearity of the delay acts as a quenching mechanism, which is appealing since the solutions become independent of the model used for the quenching.

The middle panel of Fig. 5 illustrates the dependence of the relative effective delay (τ_eff/P_cyc) on τ₀ and C_S. In this figure we identify two different regimes: the short-τ_eff and the long-τ_eff regime.

In the short τ_eff regime, the effective delay is several orders of magnitude shorter than the cycle period and therefore produces solutions that are almost identical to the nondelayed case. Some accumulation is visible at mid-cycle (compare the top panels of Figs. 4 and 6) for strongly nonlinear delays (α = −2), but disappears for weaker nonlinearity with α = −1. The cycle period is independent of C_S , as also shown in Fig. 5.

The second regime we identify is the long-τ_eff regime, which is characterized by a τ_eff between 10⁻³ P_cyc and 10⁻¹ P_cyc. This domain corresponds to the delay expected in the solar case (between a few days and a few months) and is therefore the regime of interest. The solutions self-saturate, which leads to a strong dependence of the saturation field and of the cycle period on C_S. The cycle period is also shown to be an increasing function of τ_eff.

It appears that the transient peaks in S at the surface provide a sufficient source to maintain a dynamo that would otherwise decay. Therefore the saturation mechanism is not the quenching of the source term, but the balance between the diffusion and the regular peaks of the accumulated source term S. We still do not understand why an increase of the effective delay (lower C_S) increases the cycle period.

Jouve et al. (2010b) studied the effect of τ₀ for a given $C_S=C_{S0}^{\mathrm{crit}}$. Here we demonstrate that the delay reduces $C_S^{\mathrm{crit}}$ and identify two regimes. We studied the effect of τ₀ for each regime and extended the analysis of Jouve et al. (2010b).

In Fig. 5, it is remarkable that the effective delay and the cycle period are independent of τ₀. The long-τ_eff regime results from the nonlinearity of the delay, controlled by α. Fournier et al. (2017) showed that α is a function of the azimuthal mode with which an unstable flux tube rises. Determining under which condition one or the other mode is preferred will provide solid input to constrain this parameter.

Figure 6 illustrates the radial field at the surface in the long- and short-τ_eff regimes. In both regimes the morphological characteristics of the delayed dynamo resemble the nondelay case: at high latitudes, the strong polar fields are concentrated close to the pole, propagating poleward; at low latitude, the radial fields remain weak, showing an equatorward propagation. Even though the accumulation in S increases the field strength at low latitudes, it remains two orders of magnitude weaker than the polar fields.

Additionally, in the short-τ_eff regime, the cycle period remains exclusively controlled by the meridional circulation, like in the nondelayed solutions, but in the long-τ_eff regime, the cycle period clearly increases, with decreasing C_S.

We conclude that the nonlinear delay does not affect the qualitative characteristics of the dynamo, but by increasing the cycle period, renders it a worse quantitative result. We summarize the dynamo characteristics in Table 2.

Fig. 6. Radial magnetic field at the surface as a function of time from the advection-dominated case with delay (D-ADV) for two source term amplitudes CS. The respective parameters are indicated above the plot. The black solid lines represent the 1.5 and 0.15% levels of $B_r^{\max }$ for the short-τeff case (top panel), and the 25 and 2.5% levels for the long-τeff case (bottom panel).

4.2. Diffusion-dominated case

Estimates of the turbulent magnetic diffusivity from observations and mixing-length-based stellar models suggest a value higher than 10¹² cm² s⁻¹. However, as seen in the lower panel of Fig. 4 for the nondelayed case in the diffusive regime, even though the cycle period of about 8 yr fits the observations better, the low latitudes remain weakly active and the dynamo wave propagates radially. The fact that turbulent magnetic diffusivities lower than 10¹² cm² s⁻¹ are required in order to obtain an equatorward migration has been a long-standing issue in Babcock-Leighton dynamos.

In Sect. 2.2 we have shown that for a prescribed sinusoidal magnetic field, the maximum field propagates in the latitudinal direction of the magnetic gradient. Here we discuss this mechanism for the simulated results.

As in the advection-dominated regime, the introduction of the delay reduces the criticality of the dynamo and opens a window for new types of solutions. We could also identify the same regimes, illustrated in Fig. 7. The short-τ_eff regime remains identical to the nondelayed case, and shows the same differences to the solar cycle (upper panel of Fig. 8). In the long-τ_eff regime the overall behavior is similar: the cycle period and the effective delay are independent of τ₀. The solutions saturate before reaching B_quench (left panel), and the cycle period decreases with C_S (right panel). However, the morphologic characteristics of the dynamo are clearly different. The lower panels of Fig. 8 show that the reduction of criticality allows an accumulation of B_r that reaches ≈25% of the polar region field strength. Furthermore, the latitudinal distribution of the initial magnetic field (decreasing toward the equator) leads to an equatorward propagation in the low latitudes. The last panel of Fig. 8 shows a solution that quantitatively agrees with the solar characteristics (see Table 2).

Fig. 7. D-DIFF case. Saturation field strength (left panel), effective delay (middle panel), and cycle period (right panel) as functions of the source term amplitude CS for four different delays. Each line represents a series of runs for a given delay. The vertical solid line marks the critical $C_{S0}^{\mathrm{crit}}$ for the nondelayed case. The circles identify solutions shown in Fig. 8.

In this particular run, τ_eff = 3 × 10⁻³ P_cyc , which is about a few days. Solutions with effective delays of up to a few months show comparable behaviors, however. This shows that the nonlocality introduces an additional timescale on the order of days, which is sufficient to obtain accumulation of radial-flux generation at low latitudes and equatorward migration.

4.3. Comparison to the solar case

The series we carried out in the diffusion-dominated regime shows a good quantitative agreement with the solar case, but we scanned the constrained parameter space with the remaining free parameters, namely η_t, τ₀ , and C_S, and selected a simulation that we call SOLAR, which reproduces the solar characteristics. By varying C_S and η_t, we were able to find dynamo solutions whose butterfly diagram quantitatively matches several characteritics of the solar observations. We summarize these aspects in Table 2. We compared the activity cycle period, the propagation of the active belt and the high latitude, the extent of the polar regions, and the activity level of the low latitudes. We recall that the averaged strength of the observed active belts is about an order of magnitude weaker than the polar field strength.

The SOLAR solution is antisymmetric and oscillatory, with a cycle of 11 yr, an extent of the polar regions of 36° , the amplitude of the low latitudes is a fourth of the amplitudes of the polar regions, and an equatorward propagation of the active belt as well as a poleward propagation of the high-latitude fields. It is also remarkable that the polar reversal occurs at half-cycle of the low latitude. The turbulent magnetic diffusivity required to obtain this solution is η_t = 6.7 × 10¹¹ cm² s⁻¹, which is in the transitional regime between the diffusive and the advective regime. The resulting effective delay τ_eff is on the order of a month (≈33 days). The toroidal field at the bottom of the convection zone saturates below equipartition with the convective motions, at about 0.15 B_eq.

In Fig. 9 we illustrate various aspects of the dynamo mechanism, showing how the toroidal field profile (a) is distorted into a source term sharply peaked in time (b), building up a front that, added to the diffuse field at the surface, leads to the very characteristics of the butterfly diagram (c). Because of the stiff accumulation, the active latitudes are more strongly localized in latitude and time than the solar butterfly diagram. As Weber et al. (2011) have shown that convective motions introduce stochasticity in the emergence characteristics of active regions, this mechanism could explain the broadening of the real solar activity bands in the butterfly diagram. We illustrate the propagation of the active latitudes in the meridional plane in Fig. 10, where the solid and dashed black contours represent the toroidal field, and the color-coding shows the negative and positive radial field strength, respectively. The slices are taken in the beginning of the activity cycle, at maximum, in the decreasing phase, and at minimum. The toroidal field peaks at high latitudes and propagates almost radially, as predicted by the Parker–Yoshimura rule. The strongest radial field is located at the pole. Close to the surface, where the source term is the strongest, the active latitudes, indicated by the dashed ovals, can be seen to migrate to the equator. These meridional sections demonstrate that the equatorward migration of the active belt does not follow the migration of the dynamo wave, which propagates almost radially, but is caused by the longer delay for weaker toroidal fields closer to the equator.

5. Sensitivity of the results to the free parameters

Since some parameters remain weakly or even not constrained, it is important to discuss the reliability of the results of the previous section. The accumulation of the source term S in short-lived peaks as a result of the delay is the key aspect of the model and provides the solar-like characteristics. For a prescribed oscillating magnetic field, the accumulation is controlled by the two parameters τ₀ and α (see Sect. 2.2). Fournier et al. (2017) were able to constrain α in the range of −0.91 to −2 depending on the unstable azimuthal mode of the rising flux tube. However, one of the limitations of the latter work is the lack of constraint on τ₀. It depends on many details of the formation of magnetic flux tubes and will be clearly challenging to constrain. Fortunately, we were able to demonstrate that the actual solutions in the long-τ_eff regime are clearly independent of τ₀. As for the dependence on α, we have varied α from −2 to −0.1 and found that solar-like solutions could be obtained until α = −1, which remains within the above constrained space. For −1 < α, the accumulation is not sufficient and the dynamo decays.

Fig. 8. Radial magnetic field at the surface as a function of time from the diffusion-dominated case with delay (D-DIFF) for $C_S>C_{S0}^{\mathrm{crit}}$ ({\it top panel}) and $C_S<C_{S0}^{\mathrm{crit}}$ (bottom panel). The respective parameters are indicated above the plot. Contours represent the 25% and 10% levels of $B_r^{\max }$ in the short-τeff case (top panel), and the 25% level in the long-τeff case (bottom panel).

Even though we were able to constrain the value of B_quench in the first term in the sum of Eq. (3) based on global numerical simulations, the quadratic quenching we use is quite arbitrary. Because this model possesses the surprising capability of self-quenching through its nonlinear nonlocality, however, the solutions are independent of the chosen quenching model. The level of saturation only depends on the effective delay, τ_eff, which is an outcome of the model and cannot be chosen arbitrarily.

Because we use a lower threshold on the magnetic field, some mode may not grow. The preferred mode of the dynamo is determined by the choice of the threshold and the initial condition. Although the threshold is relatively well constrained from the stability analysis of the buoyancy instability (Ferriz-Mas et al. 1994), it prevents the dynamo from growing from an arbitrary low seed field.

As is expected for the diffusion-dominated regime, the cycle period depends on η_t, but in contrast to the nondelayed dynamos, the delayed model additionally shows a dependence on C_S. The cycle period may change by a factor of two over the relevant range of C_S. This is remarkable because it might explain how solar-like stars with a comparable rotation period and convective envelope (same η_t) could show different magnetic cycles, which might be due to a difference in metalicity, for instance. This dependence needs to be carefully addressed because the interpretation of C_S as a physical quantity is not trivial. More global simulations are required to be conclusive on this issue. The nonlocal dynamos seem to be good candidates to address this particular issue, however.

Fig. 9. Results of the SOLAR setup, with the color-coded toroidal magnetic field at the bottom of the convection zone (panel a), the resulting delayed source term at the surface (panel b), and the radial magnetic field at the surface (panel c). Contours in panel c represent the 25% level of $B_r^{\max }$.

Fig. 10. Slices of the northern hemisphere. The slices are taken at different phases of the cycle, ordered from left to right and top to bottom panels: the rising phase, the maximum phase, the declining phase, and the minimum of the activity cycle. The colored contours represent the radial magnetic field, with red pointing outward and blue inward. Contours are 5% and 2.5% of the maximum Bϕ.

6. Conclusions

Until now, diffusive Babcock-Leighton dynamos were considered to be unable to qualitatively reproduce the solar dynamo. The Parker–Yoshimura rule implies for the internal differential rotation of the Sun that the dynamo wave propagates poleward. We also find that the cycle period is too short and the low-latitude radial fields are too weak.

We here introduced a delay in the source term of the poloidal field. Like in Jouve et al. (2010b), this delay represents the rise time of magnetic flux tubes through the convection zone. In contrast to previous studies, however, we built this model on the results of global numerical simulations of rising magnetic flux tubes in compressible stellar interiors (Fournier et al. 2017). The model consists of a rise time that depends nonlinearly on the magnetic flux density.

We have shown that the nonlinearity of the delay leads to an accumulation of the Babcock-Leighton source term at certain times. When this accumulation becomes sufficiently important, it may prevent the dynamo from decaying, even though the nondelayed model shows no dynamo action. The reduction of the criticality of the dynamo opens a new window to unknown solutions. These delayed dynamos have the peculiar property of self-quenching.

We found that the nonlinear delay can provide a mechanism for generating a migration of the surface fields in the direction of weaker internal fields. In case of a stationary toroidal internal field at mid-latitudes, for example, the generated poloidal fields at the surface migrate toward the equator at low latitudes and toward the poles at high latitudes. This is independent of the sign of the internal differential rotation.

The requirement of a low turbulent magnetic diffusivity, η_t, for Babcock-Leighton dynamos to qualitatively reproduce the solar cycle has been shown to be unnecessary. We demonstrated that our delayed model, with a turbulent magnetic diffusivity of η_t = 6.7 × 10¹¹ cm² s⁻¹, agrees well with the solar butterfly diagram, even though the diffusivity is relatively high throughout the entire convection zone.

We note that one proposed solution of the low-diffusivity problem is to use different values for η_t for the toroidal and for the poloidal components in the induction equation (Chatterjee et al. 2004). While the diffusivity may well be different in the horizontal and vertical directions, the poloidal field has varying components in both the horizontal and vertical directions, rendering the poloidal-field diffusivity location dependent. We have not tried such a setup, and it may not be necessary given the results we presented.

The model we presented is by design a simplified model. It has allowed us to identify the effect of the delay on the dynamo solutions. However, several ingredients are lacking for a state-of-the-art model (Rempel 2006; Cameron & Schüssler 2015; Pipin 2017). We only solved the induction equation for large-scale fields. We ignored turbulent pumping and back-reaction of the magnetic field on the flow. All these elements will increase the complexity of the model and add more free parameters that need to be constrained.

The large-scale field generation based on the Babcock-Leighton effect has not been derived from first principles. Its validity therefore remains uncertain. In the absence of global simulations addressing the formation of magnetic flux tubes, the current models remain quite arbitrary. Nevertheless, the nonlinearities of our solutions are potentially relevant for dynamos other than those of Babcock-Leighton type.

Finally, the relevance of this work for stellar dynamos will be revealed only if this model is proven to reliably reproduce observed dynamo patterns of other solar-like stars.

Acknowledgments

We would like to thank the participants and organizers of the Natural Dynamos (2016) conference for constructive remarks and input.

References

Appendix A: Constraining B_quench and B_threshold

The quenching field strength is defined such that the buoyant force balances the Coriolis force, resulting in a zero tilt at the surface. Fournier et al. (2017, Sect. 3) showed for the axisymmetric case that

$$\frac{F_{\mathrm{buoy}}}{F_{\mathrm{corio}}}=0.7{\Gamma }_1^1$$(A.1)

with

$${\Gamma }_1^1=\frac{v_{\mathrm{A}}}{\varpi \Omega }\cdot$$(A.2)

Presuming that the ratio of the buoyant force over the Coriolis force is unity for the quenching field strength, we obtain

$$B_{\mathrm{quench}}=\frac{1}{0.7}\varpi \Omega \sqrt{\rho {\mu }_0}=32\times {10}^4\mathrm{G}=32B_{\mathrm{eq}}$$(A.3)

with Ω = 430 nHz, ρ = 0.1 g cm⁻³, ϖ = 0.7 ∼ 2π 7 × 10¹⁰ cm and $B_{\mathrm{eq}}=\sqrt{{\mu }_0\rho u_{\mathrm{rms}}^2}={10}^4\mathrm{G}$ with u_rms = 100 m s⁻¹.

The field strength of a magnetic flux tube that reaches the surface needs to be larger than 3B_eq (Fan et al. 2003). The threshold field strength is therefore 3B_eq.

All Tables

All Figures

	Fig. 1. Solar-like differential rotation and meridional circulation obtained from a mean-field hydrodynamic simulation (Küker & Rüdiger 2011).
In the text

Fig. 2. Zero-dimensional toy model illustrating the effect of the delay. In each panel the green curve represents the evolution of the delay with time for a given Bϕ, indicated in the lower panels as the blue curve. The dashed line represents time. The histogram on the right-hand side represents the resulting source term S normalized in arbitrary units. Panels a and b illustrate the case of a weakly nonlinear delay (α = −1), and panels c and d show a nonlinear delay (α = −2). Panels a and c show a short effective delay (e.g., $B_{\varphi }^{\max }=B_{\mathrm{eq}}$ and τ0 = 3 · 10−5 Pcyc) and panels b and d show a long effective delay (e.g., $B_{\varphi }^{\max }=0.5B_{\mathrm{eq}}$ and τ0 = 10−3 Pcyc).

In the text

	Fig. 3. Color-coding of the toroidal magnetic field strengths of an artificial oscillatory dipole (top panel) and an artificial oscillatory quadrupole (bottom panel). The contours represent the source term (S) at the surface for α = −2 and τ0 = 10−2Pcyc.
In the text

Fig. 4. Radial magnetic field at the surface as a function of time from the solutions in the advection-dominated regime (ADV, top panel) and the diffusion-dominated regime (DIFF, bottom panel). Contours also represent Br, but show the 15, 1.5, and 0.15% levels of $B_r^{\max }$ for the top panel (ADV) and the 50, 25, and 10% levels of $B_r^{\max }$ for the bottom panel (DIFF). The parameters are summarized in Table 1.

In the text

	Fig. 5. D-ADV case. Saturation field strength (left panel), effective delay (middle panel), and cycle period (right panel) as functions of the source term amplitude CS for four different delays. Each line represents a series of runs for a given delay. The vertical solid line marks the critical $C_{S0}^{\mathrm{crit}}$ for the nondelayed case. The circles identify the solutions shown in Fig. 6.
In the text

	Fig. 6. Radial magnetic field at the surface as a function of time from the advection-dominated case with delay (D-ADV) for two source term amplitudes CS. The respective parameters are indicated above the plot. The black solid lines represent the 1.5 and 0.15% levels of $B_r^{\max }$ for the short-τeff case (top panel), and the 25 and 2.5% levels for the long-τeff case (bottom panel).
In the text

	Fig. 7. D-DIFF case. Saturation field strength (left panel), effective delay (middle panel), and cycle period (right panel) as functions of the source term amplitude CS for four different delays. Each line represents a series of runs for a given delay. The vertical solid line marks the critical $C_{S0}^{\mathrm{crit}}$ for the nondelayed case. The circles identify solutions shown in Fig. 8.
In the text

Fig. 8. Radial magnetic field at the surface as a function of time from the diffusion-dominated case with delay (D-DIFF) for $C_S>C_{S0}^{\mathrm{crit}}$ ({\it top panel}) and $C_S<C_{S0}^{\mathrm{crit}}$ (bottom panel). The respective parameters are indicated above the plot. Contours represent the 25% and 10% levels of $B_r^{\max }$ in the short-τeff case (top panel), and the 25% level in the long-τeff case (bottom panel).

In the text

	Fig. 9. Results of the SOLAR setup, with the color-coded toroidal magnetic field at the bottom of the convection zone (panel a), the resulting delayed source term at the surface (panel b), and the radial magnetic field at the surface (panel c). Contours in panel c represent the 25% level of $B_r^{\max }$.
In the text

	Fig. 10. Slices of the northern hemisphere. The slices are taken at different phases of the cycle, ordered from left to right and top to bottom panels: the rising phase, the maximum phase, the declining phase, and the minimum of the activity cycle. The colored contours represent the radial magnetic field, with red pointing outward and blue inward. Contours are 5% and 2.5% of the maximum Bϕ.
In the text