A BABCOCK–LEIGHTON SOLAR DYNAMO MODEL WITH MULTI-CELLULAR MERIDIONAL CIRCULATION IN ADVECTION- AND DIFFUSION-DOMINATED REGIMES

We use the spherical polar coordinates $({\rm r},\theta ,\phi )$, and assume the axisymmetry. The magnetic field, as the sum of a toroidal component (${B}_{\phi }$) and a poloidal component (${\boldsymbol{B}}_{\boldsymbol{p}}$), can be written as:

$$\begin{eqnarray}&&\boldsymbol{B}(r,\theta ,t)={B}_{\phi }(r,\theta ,t){\hat{\boldsymbol{e}}}_{\phi }+{\bf \nabla }\times \left[A(r,\theta ,t){\hat{\boldsymbol{e}}}_{\phi }\right]\end{eqnarray} \tag{ 1 }$$

where the toroidal component of the magnetic field is ${B}_{\phi }(r,\theta ,t){\hat{\boldsymbol{e}}}_{\phi }$ and the vector potential is $A(r,\theta ,t){\hat{\boldsymbol{e}}}_{\phi }$. Both components can be influenced by a flow. The large-scale flow field $\boldsymbol{U}(r,\theta )$ is expressed as the sum of differential rotation $({\rm \Omega }(r,\theta ))$ and the meridional circulation $(\boldsymbol{u}(r,\theta )={u}_{r}(r,\theta ){\hat{\boldsymbol{e}}}_{r}+{u}_{\theta }(r,\theta ){\hat{\boldsymbol{e}}}_{\theta })$:

$$\begin{eqnarray}&&\boldsymbol{U}(r,\theta )=\boldsymbol{u}(r,\theta )+r\;\mathrm{sin}\;\theta {\rm \Omega }(r,\theta ){\hat{\boldsymbol{e}}}_{\phi }\end{eqnarray} \tag{ 2 }$$

as toroidal and poloidal parts of the total flow field.

The evolution of the large-scale magnetic field B according to:

$$\begin{eqnarray}&&\frac{{\bf \partial }\boldsymbol{B}}{\partial {\rm t}}={\bf \nabla }\times (\boldsymbol{U}\times \boldsymbol{B}-\eta {\bf \nabla }\times \boldsymbol{B}).\end{eqnarray} \tag{ 3 }$$

The toroidal component becomes

$$\begin{eqnarray}&&\begin{array}{c}{\displaystyle \frac{\partial {B}_{\phi }}{\partial t}}+{\displaystyle \frac{1}{r}}\left[{\displaystyle \frac{\partial }{\partial r}}({{ru}}_{r}{B}_{\phi })+{\displaystyle \frac{\partial }{\partial \theta }}({u}_{\theta }{B}_{\phi })\right]\\ \quad =r\;\mathrm{sin}\;\theta ({\boldsymbol{B}}_{\boldsymbol{p}}\cdot {\bf \nabla }){\rm \Omega }-{\bf \nabla }\eta \times {\bf \nabla }\times {B}_{\phi }{\hat{\boldsymbol{e}}}_{\phi }\\ \quad +\eta \left({{\rm \nabla }}^{2}-{\displaystyle \frac{1}{{r}^{2}{\mathrm{sin}}^{2}\theta }}\right){B}_{\phi }\end{array}\end{eqnarray} \tag{ 4 }$$

where $\eta (r)$ is the magnetic diffusivity.

$$\begin{eqnarray}&&\begin{array}{c}{\displaystyle \frac{\partial A}{\partial t}}+{\displaystyle \frac{1}{r\;\mathrm{sin}\;\theta }}(\boldsymbol{u}\cdot {\bf \nabla })(r\;\mathrm{sin}\;\theta A)\\ \quad =\eta \left({{\rm \nabla }}^{2}-{\displaystyle \frac{1}{{r}^{2}{\mathrm{sin}}^{2}\theta }}\right)A+S(r,\theta ;{B}_{\phi })\end{array}\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{B}}_{\boldsymbol{p}}={\bf \nabla }\times ({\rm A}{\hat{\boldsymbol{e}}}_{\phi })$.

There is considerable uncertainty about what is the most realistic profile of turbulent magnetic diffusivity with radius. Direct measurements of magnetic diffusivity as a function of depth are not possible yet. The mixing-length theory gives us a rough estimate of the supergranular diffusivity in the supergranulation layer near the surface with a range ${\eta }_{\mathrm{super}}={10}^{12}-{10}^{14}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$. The magnetic diffusivity in the convective envelope of the Sun is dominated by its turbulent contribution, but below the convection zone there is much less turbulence and the core is stabler, so the diffusivity should be determined essentially from the molecular contribution in the stably stratified deep radiative interior (Dikpati et al. 2006). We assume that the turbulence governs the diffusivity in the convection zone and is significantly reduced in the subadiabatically stratified radiative zone below (Dikpati et al. 2004). The diffusivity profile can be written as:

$$\begin{eqnarray}\begin{array}{ccc}\eta (r) & = & {\eta }_{\mathrm{core}}+{\displaystyle \frac{{\eta }_{T}}{2}}\left[1+\mathrm{erf}\left({\displaystyle \frac{r-{r}_{8}}{{d}_{8}}}\right)\right]\\ & & +{\displaystyle \frac{{\eta }_{\mathrm{super}}}{2}}\left[1+\mathrm{erf}\left({\displaystyle \frac{r-{r}_{9}}{{d}_{9}}}\right)\right].\end{array}\end{eqnarray} \tag{ 6 }$$

The parameters in the diffusivity profile: ${\eta }_{T}=3\cdot {10}^{11}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$ is the turbulent diffusivity and ${\eta }_{\mathrm{core}}={10}^{9}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$ is the core diffusivity, ${\eta }_{\mathrm{super}}=3\cdot {10}^{12}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$ is the supergranular diffusivity; ${r}_{8}=0.7\;{R}_{\odot }$, ${d}_{8}=0.0125\;{R}_{\odot }$, ${r}_{9}=0.9562\;{R}_{\odot }$, ${d}_{9}=0.025\;{R}_{\odot }$.

A two-step profile, which leaves out the supergranule layer, will have impact on the surface poloidal fields that are used to create the poloidal part of the butterfly diagram.

Our flux-transport dynamo can be driven by both a tachocline α-effect and Babcock–Leighton-type surface poloidal source. The decay of tilted bipolar active regions produce poloidal fields near the surface, but perhaps that cannot be the sole driver for the large-scale solar dynamo. A Babcock–Leighton dynamo is not a self-excited dynamo, and hence cannot come back to a normally cycling dynamo if it goes to a grand minimum state. Flux-transport solar dynamo models with both Babcock–Leighton-type surface poloidal source and convection zone α-effect have been used (Hazra et al. 2014; Passos et al. 2014), in the context of reviving a Babcock–Leighton dynamo from a grand minima. Furthermore, it is well known that a flux-transport dynamo driven solely by a Babcock–Leighton poloidal source cannot sustain its antisymmetric magnetic field about the equator, as inferred from Hale's polarity rule (Hale et al. 1919), but slowly drifts to symmetric magnetic field solutions. This is true even in the case of a 3D Babcock–Leighton dynamo (see, e.g., Miesch & Dikpati 2014). Inclusion of a tachocline instability-driven α-effect is demonstrated to be one of the solutions to stop such a drift to a nonsolar-like, quadrupolar parity in a pure Babcock–Leighton dynamo models (Dikpati & Gilman 2001; Bonanno et al. 2002). In the present paper, since one of our goals is to investigate the parity produced by a Babcock–Leighton dynamo model operating with different multi-cellular meridional circulation patterns, we assumed zero tachocline α-effect in all cases. We use the following expressions respectively for Babcock–Leighton surface source:

$$\begin{eqnarray}\begin{array}{ccc}S(r,\theta ,{B}_{\phi }) & = & \left\{{S}_{\mathrm{BL}}(r,\theta )+{S}_{\mathrm{tachocline}}(r,\theta )\right\}\\ & & \times {\left[1+{\left({\displaystyle \frac{{B}_{\phi }(r,\theta ,t)}{{B}_{0}}}\right)}^{2}\right]}^{-1}\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{S}_{\mathrm{tachocline}}(r,\theta )=0\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{ccc}{S}_{\mathrm{BL}}(r,\theta ) & = & {\displaystyle \frac{{s}_{1}}{2}}\left[1+\mathrm{erf}\left({\displaystyle \frac{r-{r}_{4}}{{d}_{4}}}\right)\right]\left[1-\mathrm{erf}\left({\displaystyle \frac{r-{r}_{5}}{{d}_{5}}}\right)\right]\\ & & \times \mathrm{sin}\theta \mathrm{cos}\theta \left[{\displaystyle \frac{1}{1+{e}^{{\gamma }_{1}(\pi /4-\theta )}}}\right]\end{array}\end{eqnarray} \tag{ 9 }$$

for $\theta \lt \pi /2$ and for $\theta \gt \pi /2$

$$\begin{eqnarray}\begin{array}{ccc}{S}_{\mathrm{BL}}(r,\theta ) & = & {\displaystyle \frac{{s}_{1}}{2}}\left[1+\mathrm{erf}\left({\displaystyle \frac{r-{r}_{4}}{{d}_{4}}}\right)\right]\left[1-\mathrm{erf}\left({\displaystyle \frac{r-{r}_{5}}{{d}_{5}}}\right)\right]\\ & & \times \mathrm{sin}\theta \mathrm{cos}\theta \left[{\displaystyle \frac{1}{1+{e}^{{\gamma }_{1}(\theta -3\pi /4)}}}\right].\end{array}\end{eqnarray} \tag{ 10 }$$

The parameters used in Equations (9) and (10) are: ${s}_{1}=2.5\;{\mathrm{ms}}^{-1}$, ${r}_{4}=0.95\;{r}_{\odot }$, ${r}_{5}=0.9875\;{R}_{\odot }$, ${d}_{4}={d}_{5}=0.0125\;{R}_{\odot }$, ${\gamma }_{1}=30.0$, $\beta =70.0$, and ${B}_{0}=10\;\mathrm{kG}$.

The solar internal rotation profile includes primarily latitudinal shear in the convection zone, as found by two-dimensional helioseismic inversions. A solar-like internal differential rotation profile (Dikpati & Charbonneau 1999) is given by

$$\begin{eqnarray}&&{\rm \Omega }(r,\theta )={{\rm \Omega }}_{c}+\frac{1}{2}\left[1+\mathrm{erf}\left(2\frac{r-{r}_{c}}{{d}_{1}}\right)\right]\{{{\rm \Omega }}_{s}(\theta )-{{\rm \Omega }}_{c}\}\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{{\rm \Omega }}_{s}={{\rm \Omega }}_{\mathrm{Eq}}+{a}_{2}{\mathrm{cos}}^{2}\theta +{a}_{4}{\mathrm{cos}}^{4}\theta \end{eqnarray} \tag{ 12 }$$

is the surface latitudinal differential rotation. Values were chosen to closely resemble the best fit to the helioseismic solution of Charbonneau et al. (1998). The angular velocity of the rigidly rotating core is ${{\rm \Omega }}_{c}/2\pi =432.8\;\mathrm{nHz}$. ${{\rm \Omega }}_{\mathrm{Eq}}/2\pi =460.7\;{nHz}$ is the rotation rate at the equator. The other parameters are set to be ${a}_{2}/2\pi =62.69\;\mathrm{nHz}$ and ${a}_{4}/2\pi =67.13\;\mathrm{nHz}$. ${r}_{c}=0.7\;{R}_{\odot }$ indicates the central radius of the tachocline thickness ${d}_{1}=0.025\;{R}_{\odot }$.

The differential rotation profile, prescribed by three terms as shown in Equation (12), is well-formulated up to $60^\circ$, but to adequately fit the measured rotation rate including more poleward latitudes requires additional terms (Schou et al. 1998; Dikpati et al. 2002). We will discuss how the implementation of such differential rotation profile affects our results.

The circulation is represented in the spherical shell by the stream function of Forgács-Dajka & Petrovay (2002) for both hemispheres. To study the effects of the meridional circulation, we use a simple and easily adjustable spatial structure, with realistic amplitude. The components of the meridional circulation can be written as

$$\begin{eqnarray}&&{u}_{r}(r,\theta )=\frac{1}{\rho (r)\cdot {r}^{2}\mathrm{sin}\theta }\cdot \frac{\partial }{\partial \theta }{\bf \Psi }(r,\theta )\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{u}_{\theta }(r,\theta )=\frac{-1}{\rho (r)r\;\mathrm{sin}\;\theta }\cdot \frac{\partial }{\partial r}{\bf \Psi }(r,\theta )\end{eqnarray} \tag{ 14 }$$

Note that the spherical polar geometrical factor ($r\;\mathrm{sin}\theta$) is absorbed in ${\bf \Psi }(r,\theta )$.

We reproduce the stream function of Forgács-Dajka & Petrovay (2002), which has the following form:

$$\begin{eqnarray}&&{\bf \Psi }(r,\theta )=\boldsymbol{\psi }(r){\mathrm{sin}}^{2}{n}_{1}\theta \mathrm{cos}{n}_{2}\theta \end{eqnarray} \tag{ 15 }$$

where the form of the given function $\psi (r)$ specifies the flow.

$$\begin{eqnarray}&&\psi (r)={\psi }_{0}\mathrm{sin}\left[\frac{k\pi \left(r-{r}_{\mathrm{mc}}\right)}{{R}_{\odot }-{r}_{\mathrm{mc}}}\right]\mathrm{exp}\left[\frac{{(r-{r}_{0})}^{2}}{{{\rm \Gamma }}^{2}}\right]\end{eqnarray} \tag{ 16 }$$

where ${\psi }_{0}$ sets the amplitude of the meridional circulation. Two parameters play the primary role in setting the meridional circulation for the simulations. We fix the number of cells in latitude by changing the parameter n and the number of cells in radius by changing the parameter k. ${r}_{0}$ and ${\rm \Gamma }$ are geometric parameters, ${r}_{0}=({R}_{\odot }-{r}_{\mathrm{mc}})/30\;\mathrm{cm}$, and ${\rm \Gamma }=6.4\cdot {10}^{10}\;\mathrm{cm}$. r_mc is the radius to which the meridional flow penetrates from the base of the convective zone. Given the observed depth of the tachocline, we set ${r}_{\mathrm{mc}}=4.79\cdot {10}^{10}\;\mathrm{cm}$. We use the same density profile ($\rho (r)$) as used by Dikpati & Charbonneau (1999).

With all parameters specified, we can compute the stream functions for the five circulation patterns we use in the dynamo. Figure 1 displays the circulation patterns that we will use for our analysis of the differences in dynamo behavior that arise from different circulation patterns. These include the reference case (a) a single cell in each hemisphere, (b) a long, primary cell associated with a weak reverse cell in high latitudes, (c) two cells with depth, (d) two cells in latitude and two in depth, and (e) four cells in latitude. There is some observational support for each of these patterns, as reviewed in the introduction. The parameters chosen in the five cases are given respectively as: (a) k = 1, ${n}_{1}={n}_{2}=1$, and ${{\rm \Psi }}_{0}=-6.1771\times {10}^{21}$; (b) k = 1, ${n}_{1}=1$, ${n}_{2}=3$, and ${{\rm \Psi }}_{0}=3.519\times {10}^{21}$; (c) k = 2, ${n}_{1}=1$, ${n}_{2}=1$, and ${{\rm \Psi }}_{0}=3.0937\times {10}^{21}$; (d) k = 2, ${n}_{1}=1$, ${n}_{2}=3$, and ${{\rm \Psi }}_{0}=-1.7625\times {10}^{21}$; and (e) k = 1, ${n}_{1}=4$, ${n}_{2}=4$, and ${{\rm \Psi }}_{0}={{{\rm \Psi }}_{0}}_{a}+{{{\rm \Psi }}_{0}}_{b}+{{{\rm \Psi }}_{0}}_{c}+{{{\rm \Psi }}_{0}}_{d}$, in which we assign values for ${{{\rm \Psi }}_{0}}_{a},{{{\rm \Psi }}_{0}}_{b},{{{\rm \Psi }}_{0}}_{c},{{{\rm \Psi }}_{0}}_{d}$ for a certain range in θ and zero elsewhere, namely ${{{\rm \Psi }}_{0}}_{a}=1.7043\times {10}^{21}$ for 0 ≤ θ ≤ 225, ${{{\rm \Psi }}_{0}}_{b}=4.1340\times {10}^{21}$ for 225 ≤ θ ≤ 45°, ${{{\rm \Psi }}_{0}}_{c}\;=-6.8972\times {10}^{21}$ for 45° ≤ θ ≤ 675, and ${{{\rm \Psi }}_{0}}_{d}=-7.8002\times {10}^{21}$ for 675 ≤ θ ≤ 90°. The units of ${{\rm \Psi }}_{0}$ are in c.g.s, i.e., ${\mathrm{cm}}^{3}\;{{\rm s}}^{-1}\;\mathrm{unit}\;\mathrm{of}\;{\rho }^{-1}$.

Zoom In Zoom Out Reset image size

Equations (4) and (5) are solved in a full spherical shell ($0\geqslant \theta \geqslant \pi$), extending radially from below the bottom of convection zone ($r/{R}_{\odot }=0.6$) to the surface ($r/{R}_{\odot }=1$). We use a 4th order Runge–Kutta method with the central finite difference scheme for r and θ derivatives to solve the equations. A validity test of the scheme was performed by reproducing the results the benchmark cases of a Babcock–Leighton flux-transport dynamo (Jouve et al. 2008). We use the boundary conditions as used previously by Dikpati & Charbonneau (1999). We set ${B}_{\phi }=0$ and $A=0$ at the bottom boundary ($r/{R}_{\odot }=0.6$). Note that the bottom boundary condition for the toroidal field is modified from the condition for a perfectly conducting core (see, e.g., Dikpati & Gilman 2001 for physical explanation). Both ${B}_{\phi }=0$ and $A=0$ are also set to zero along the polar axis ($\theta =0{\rm }\;{\rm or}\;{\rm }\pi$) to ensure physicality. At the surface, the toroidal fields are set to zero, and we demand that the poloidal field lines match smoothly with the potential field solutions of the free space above the surface (see Dikpati & Choudhuri 1994). The vector potential $A$ satisfies the following equation above the surface:

$$\begin{eqnarray}&&\left({{\rm \nabla }}^{2}-\frac{1}{{r}^{2}{\mathrm{sin}}^{2}\theta }\right)A=0.\end{eqnarray} \tag{ 17 }$$

A general solution for this equation can be written as:

$$\begin{eqnarray}&&A(r\geqslant {R}_{\odot },\theta ,t)=\underset{n}{\sum }\frac{{a}_{n}(t)}{{r}^{n+1}}{P}_{n}^{l}(\mathrm{cos}\theta )\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{ccc}{a}_{n}(t) & = & {\displaystyle \frac{(2n+1){R}_{\odot }^{n+1}}{n(n+1)}}\\ & & \times {{\displaystyle \int }}_{0}^{\pi /2}A(r={R}_{\odot },\theta ,t){P}_{n}^{l}(\mathrm{cos}\theta )\mathrm{sin}\theta d\theta \end{array}\end{eqnarray} \tag{ 19 }$$

where ${P}_{n}^{l}(\mathrm{cos}\theta )$ is the associated Legendre polynomial. The derivative of $A$ at the solar surface is

$$\begin{eqnarray}&&{\left.\frac{\partial A}{\partial r}\right|}_{(r={R}_{\odot })}=-\underset{n}{{\displaystyle \sum }}\frac{(n+1){a}_{n}(t)}{{R}_{\odot }^{n+2}}{P}_{n}^{l}(\mathrm{cos}\theta ).\end{eqnarray} \tag{ 20 }$$

First we establish a reference case dynamo solution, to which other individual solutions can be compared, to see what changes are created by changing the meridional circulation pattern. We choose for reference the frequently used single cell meridional circulation that has poleward flow near the outer boundary, and return flow at the base of the convection zone (see Figures 1(a) and 2(a)). This flow penetrates slightly below the base of the convection zone, given by $r=0.7{R}_{\;\odot }$. The thin black dashed semicircular arc represents the location of the center of the tachocline. Two frames, (b) and (c) on the right panel of Figure 2, show the time–latitude diagrams respectively of toroidal field at the bottom of convection zone (${B}_{\phi }{| }_{r=0.7\;{R}_{\odot }}$) and surface radial field (${B}_{r}{| }_{r={R}_{\odot }}$). To get these results we took a maximum surface flow speed of $15\;{\mathrm{ms}}^{-1}$, poloidal source amplitude of ${s}_{1}=3.0\;{\mathrm{ms}}^{-1}$, and turbulent diffusivity ${\eta }_{T}=3\cdot {10}^{11}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$.

Zoom In Zoom Out Reset image size

Figure 3 shows how the poloidal and toroidal fields evolve through a sunspot cycle for the entire meridional cross-section. When the toroidal field is strong, it is confined to the lower layers of the domain, and its peak is clearly migrating toward the equator, along with the poloidal field lines that are sheared by differential rotation to produce it. We expect sunspots to emerge from at or near the latitude of maximum in the toroidal field. The reversal of polar fields is seen here in the interval between 5 and 6.25 years, at which time the toroidal is a maximum with its peak near $20^\circ$.

Zoom In Zoom Out Reset image size

Babcock–Leighton flux-transport dynamo models with one-celled meridional circulation in each hemisphere can reproduce many features of the solar cycle. These include (i) the equatorward migration of toroidal flux at lower latitudes, (ii) the 11 year sunspot cycle, (iii) the observed phase shift between poloidal and toroidal components, (iv) the short rise of toroidal field to maximum followed by the long decline to minimum—in the reference case, the ascending phase is $16.77\%$ of the whole cycle, (v) peak tachocline toroidal fields of $37\;\mathrm{kG}$, (vi) and peak surface radial fields of $\pm 204\;{\rm G}$, similar to values obtained by many previous authors (Durney 1995; Dikpati & Charbonneau 1999; Dikpati & Gilman 2001; De Rosa 2005; Rempel et al. 2005; Dikpati & Gilman 2006; Jouve & Brun 2007; Dikpati et al. 2010; Belucz & Dikpati 2013). This solution also reproduces the observed phase shift between the surface poloidal field and the toroidal field at the tachocline; the poloidal field polarity changes from positive to negative when the toroidal field is near maximum and its polarity is negative.

In the next simulations we study how the characteristic features of butterfly diagram change when the meridional circulation cell contains a second, high-latitude, reverse cell. The form of this meridional circulation is shown in Figures 1(b) and 4(a). The peak flow speed of the primary cell is still $15\;{\mathrm{ms}}^{-1}$, poleward at the surface, but the latitude of this peak is slightly lower, at 253. The peak flow speed of the secondary cell is $3\;{\mathrm{ms}}^{-1}$, equatorward at the surface. The boundary between cells is near $61^\circ$ latitude.

Zoom In Zoom Out Reset image size

The right panel of Figure 4 shows the time–latitude diagrams of ${B}_{\phi }{| }_{r=0.7{R}_{\odot }}$ in panel (b) and ${B}_{r}{| }_{r={R}_{\odot }}$ in panel (c). Not surprisingly, due to the effect of the second cell, a more pronounced poleward branch can be seen in the butterfly diagram of toroidal field in panel (b) compared with that in Figure 2. The sunspot cycle length (i.e., half magnetic cycle) is just $8.3$ years, due to the shorter primary conveyor belt. The strength of toroidal field is similar to that of reference case, $37\;\mathrm{kG}$. Comparing the time–latitude diagrams of Figures 2(c) and 4(c), we find that the polar field peaks around $50^\circ$ latitude instead of peaking near the pole as in the case of a single cell in each hemisphere. This is due to flow convergence at $61^\circ$ latitude instead of at the pole. The second cell also causes a delay in the polarity change by advecting polar fields away from the pole.

The rise of the cycle from minimum to maximum in this case is slightly longer compared with that in the single cell case. This is probably because some of the poloidal flux advected to the bottom in between the primary and secondary cells is advected toward the poles, retarding the early production of the equatorward migrating branch of toroidal field there.

For the next simulation, we add a second, reversed meridional cell below the primary cell (see Figure 1(c)). The two cells are similar in amplitude and radial extent. The results for this meridional flow are shown in Figure 5.

Zoom In Zoom Out Reset image size

This dynamo solution is very different from the reference case. At all latitudes, the tachocline toroidal field is migrating with time toward the poles rather than the equator. The rate of this migration increases with latitude. This is due to the poleward flow at the bottom of the second cell. Thus, based on tachocline toroidal fields, this meridional circulation pattern produces an "antisolar" butterfly diagram. The toroidal field at the bottom of the domain is weaker than in the reference case, because polar fields are not brought down from the top, but instead are advected from lower latitudes near the bottom.

If we plot instead the toroidal field near the middle depth, for example, $0.83125\;{\rm R}$, shown in Figure 6, we get a more solar-like butterfly, with both poleward and equatorward branches. This is because at these depths both circulation cells have equatorward flow, so they advect toroidal field toward the equator in lower latitudes. The relatively high speed total flow there also makes the dynamo period slightly shorter. In addition, since in this case the poleward flow near the outer boundary reaches to a shallower depth, less poloidal flux is advected toward to the pole. These weaker polar fields lead to weaker toroidal fields at all latitudes.

Zoom In Zoom Out Reset image size

The next dynamo simulation is for the case of two meridional cells in both latitude and depth in each hemisphere. The amplitudes of upper and lower cells are about the same. These results are shown in Figure 7. Panel (b) again shows the toroidal fields near the bottom of the convective zone. Panel (c) depicts the surface radial field. As in the case of two cells in radius, we again get an "antisolar" butterfly diagram when we plot the toroidal field at the bottom. Again this is because the flow is toward the pole at the bottom of the stack of cells at all latitudes except the highest. If we plot the toroidal field contours at mid-depth (not shown), we will again get a solar type butterfly, as seen in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the last simulation, we also have four cells, but these cells are located side by side, as seen in the Figure 1(e). In this case, the fields are confined to lower latitudes because the multiple cells in latitude prevent poloidal field transport all the way to the poles as in the reference case. We get a solar-like time–latitude diagram up to about $22^\circ$; the toroidal field migration is equatorward. We do not see dynamo activity beyond about $50^\circ$. The cycle length is very short, 3.125 years, due to the very short conveyor belts represented by the two circulation cells closest to the equator. The strength of the toroidal field is just half of the reference case, because the dynamo is confined to the lower latitudes where the differential rotation is smaller so the production of toroidal field from a given poloidal field is smaller.

How dependent are the results described above on the particular choice of turbulent magnetic diffusivity? We address that question by repeating the dynamo simulations for the five meridional circulations using a lower value, $7\times {10}^{10}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$. We selected this value of diffusivity (${\eta }_{T}=7\times {10}^{10}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$) to represent the low-diffusivity cases, and in Section 3.1 a diffusivity value of $3\times {10}^{11}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$ to denote the representative value for high diffusivities. However, these values are not unique; neighboring values can be considered also.

The results for low diffusivity are displayed in Figures 9–13. Figure 9 displays the reference case with single circulation cell for the lower magnetic diffusivity. As in earlier figures, panel (b) gives the toroidal field amplitude near the bottom and panel (c) shows the surface poloidal field. We use the same peak surface flow speed ($15\;{\mathrm{ms}}^{-1}$) and poloidal source amplitude (${{\rm s}}_{1}=3.0\;{\mathrm{ms}}^{-1}$) as for the solution seen in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this lower diffusivity case, the model still reproduces the observed phase shift between the surface poloidal field and the tachocline toroidal field. Due to the lower turbulent diffusivity, much stronger fields are produced. The highest tachocline toroidal field is $52\;\mathrm{kG}$ and the maximum value of the radial fields is $\pm 207\;{\rm G}$ near the poles. The period is also somewhat longer than in Figure 2, namely $12.7$ years, as measured by the time between adjacent peaks in toroidal and poloidal field at the same latitude. But there is also significant overlap between adjacent cycles, so the time between the high latitude peak of a new sign of toroidal field to its disappearance at low latitudes is more than $20$ years. How much overlap there is in the Sun itself is unclear, since from observations we know where the tachocline toroidal field is only from the latitude of sunspots seen. The slight poleward migration of tachocline toroidal fields seen at high latitudes comes from the strong negative radial shear there overcoming the relatively weak equatorward advection of toroidal field by the meridional flow there.

Figure 10 shows the time–latitude diagrams of toroidal and radial fields in the two cell case (Figure 10(a)) with $7\cdot {10}^{10}\;{\mathrm{cm}}^{2}\;{{\rm s}}^{-1}$ turbulent diffusivity. For this case we see that the high and low latitude branches of the butterfly diagram for tachocline toroidal field are about the same, despite the unequal amplitudes of the two circulation cells. Both the poleward meridional circulation at the bottom and the radial gradient of rotation in high latitudes at tachocline depths are contributing to this pattern, which is more pronounced than seen in Figure 9, for which the bottom meridional circulation is toward the equator at high latitudes. The sunspot cycle is shorter, about $10.7$ years, due to the shorter conveyor belt that is the primary circulation cell. Here too, the induced toroidal field is stronger than in the high diffusivity solution with the same meridional circulation. The maximum value of the toroidal field is $78\;\mathrm{kG}$; the maximum radial field is $77\;{\rm G}$. In this case we get substantially less overlap of adjacent cycles than in the low-diffusivity reference case (Figure 9).

In the next simulation with low diffusivity, we turn on a second meridional cell below the primary cell, as for the case previously displayed in Figure 5. The results are shown in Figure 11. The most prominent feature of the butterfly diagram is the longer sunspot cycle, about $50$ years. What little migration of toroidal field there is, is toward the poles, leading to a slightly antisolar butterfly. Here again, as in the higher diffusivity case shown in Figure 6, toroidal field contours (not shown) at mid-depth of the convection zone, yields a butterfly diagram with both poleward and equatorward branches.

The results of the next simulation is shown in Figure 12. Here, with two strong cells in low and mid-latitudes, and two reversed, relatively weak cells in polar latitudes, we again get a predominantly antisolar butterfly but with a small domain of equatorward migration of surface poloidal fields in low latitudes. This must be caused by the equatorward flow at mid-depth coupled with upward magnetic diffusion to the surface. Finally, here again we get a solar-like butterfly for low and middle latitudes when we plot toroidal field contours at mid-depth in the dynamo domain.

In the last simulation, we show results for the four cell case for low diffusivity in Figure 13. Due to the short conveyor belts, decreasing the turbulent diffusivity does not significantly change the cycle length. But the lower turbulent diffusivity has other effects. First, the fields are stronger, as we should expect. Second, we can see dynamo activity at the higher latitudes than in the high diffusivity case, though it is still low compared with that in low latitudes. The orientation of the wings of the time–latitude diagram at the different latitudes is determined by the latitudinal direction of the flow near the bottom; equatorward flow leads locally to migration toward the equator, and poleward flow to migration toward the poles, as seen in Figure 13(b).

From all of the time–latitude plots shown above, we can see that the most solar-like diagrams are produced if there is a single primary circulation cell in each hemisphere, with possibly a weaker secondary, reversed cell in polar latitudes. Circulation with two cells in depth, or two cells in both latitude and depth, give solar-like butterflies only from toroidal fields at mid-depth, not the bottom. For these to be correct for the Sun, the toroidal fields at the bottom must not come to the surface because of their magnetic buoyancy or for any other reason, and a mechanism must exist that keeps mid-depth toroidal fields from rising buoyantly too fast to be amplified to produce sunspots. Neither requirement is easily satisfied using known MHD processes.

We have shown how solutions from a solar flux-transport dynamo model differ for different forms of meridional circulation. The solutions we obtained are all for the same differential rotation, since from helioseismic measurements that is relatively well known for the Sun. These were found for fixed meridional circulation amplitude, but that amplitude is not well constrained from observations, and it is bound to have time variations. In addition, there are uncertainties in the amplitude and form of the surface poloidal flux source as well as the turbulent magnetic diffusivity. It is of interest to know how basic characteristics of a simulated sunspot cycle differ for different values of these uncertain parameters. Three prominent features of simulated cycles to focus on are its period, amplitude, and shape—the times spent in ascending and descending phases. We focus here on the first two of these.

We choose the range of meridional circulation amplitudes to cover the range observed, and somewhat beyond. The choice of surface poloidal source amplitude range is guided by estimates of relation to photospheric magnetograms as discussed in Dikpati & Gilman (2006). While the choice of ranges of values for the meridional circulation amplitude and surface poloidal source is guided by observations, and hence adequately maps the plausible parameter space for the Sun, the range of the turbulent diffusivity amplitudes is selected based on the values that produce a sustaining dynamo. Recent measurements lead to a range of magnetic diffusivity from 1 × 10¹⁰ cm² s⁻¹ to 2 × 10¹¹ cm² s⁻¹ (Chae et al. 2008) while much order results show larger values from 2.5 × 10¹² cm² s⁻¹ (Schrijver & Martin 1990) to 6.0 × 10¹² cm² s⁻¹ (Wang et al. 1989). Given the turbulent diffusivity is one of the biggest unknown ingredients in the solar interior, we performed extensive numerical experiments to choose the range of ${\eta }_{T}$ so that the dynamo does not die due to too much diffusive decay, or does not produce unusually large magnetic field amplitudes due to too low diffusivity.

Figure 14 displays the variation of cycle period with circulation amplitude (panel (a)), turbulent diffusivity (panel (b)), and poloidal source amplitude (panel (c)). We see from panel (a) that as the circulation amplitude is increased, in almost all cases the period declines. This is to be expected, because in all cases, unless diffusion dominates, the period is set by the speed of the conveyor belt. The primary exception we see is, that for low speeds, decreasing the circulation in the single cell case (and to much lesser degree, the case with a second weak cell at high latitudes) leads to a decrease in period. This happens because, while the circulation is decreased, turbulent diffusion starts to "short circuit" the conveyor belt, since some poloidal flux is diffused toward the bottom from the top before it reaches polar latitudes.

Zoom In Zoom Out Reset image size

This short circuiting effect is even more evident in panel (b), where we have plotted cycle length versus turbulent diffusivity. For the same circulation amplitude, the solutions become more affected by diffusion to the right in the figure. The periods decline, in some cases by factors of five or more. This result shows that to have a flux-transport dynamo calibrated to the observed sunspot cycle period requires careful choice of the turbulent diffusivity, no matter what circulation pattern is assumed. By contrast, panel (c) shows that the cycle period is almost independent of the amplitude of the surface poloidal source. This is also expected, because the dynamo is nearly linear. Changing the poloidal source amplitude should change primarily the peak amplitude of the cycle, as we shall see below.

Figure 15 shows how the maximum toroidal field varies with meridional circulation (panel (a), turbulent diffusivity (panel (b)), and amplitude of the poloidal source term (panel (c)) for the five circulation patterns. From panel (c) we see that, as we should expect, raising the amplitude of the poloidal source raises the peak toroidal field amplitude. Because of the nonlinear quenching of the source term internal to the model, however, the amplitudes are beginning to approach asymptotic limits. There are also significant differences in the efficiency of different circulation patterns; two cells in depth and four cells in latitude both generate much less toroidal field than the other patterns, which give almost the same amplitudes. In the four cell case, this is because the shearing of poloidal field is largely confined to low latitudes, where the latitudinal rotation gradient is weakest, by the short latitudinal extent of the conveyor belt, so less toroidal field is generated. In the case of two cells in depth, less of the poloidal flux gets to the bottom where the radial shear is strongest, again reducing the dynamo's ability to amplify toroidal field.

Zoom In Zoom Out Reset image size

In panel (b), the result that the peak toroidal field declines with increasing magnetic diffusivity is due to the fact that there is more dissipation in the system that the induction processes in the dynamo have to overcome. Here too we see that, for all diffusivities, configurations with four cells in latitude and two cells in depth generate less toroidal field for the same assumed diffusivity.

From panel (a), we see that the variation in peak toroidal field with changes in circulation amplitude are more complex. With two cells in depth the peak amplitude changes relatively little with circulation amplitude. With two cells in both latitude and depth, we see a similar result, except near circulation amplitude $15\;{\mathrm{ms}}^{-1}$ where the amplitude drops by more than $50\%$ when the magnetic field configuration switches from dipolar to quadrupolar. In both cases, with the relatively fine scale meridional circulation pattern, the solutions are in the diffusivity dominated regime, so the amplitude does not change. Meridional circulation acts mainly as a transporter of toroidal and poloidal flux, rather than as an amplifier. The faster the toroidal field is transported, the weaker is the resulting field amplitude, because the toroidal field peak is displaced before it amplified as much as it could be if it were stationary.

The other three cases each show the largest value of toroidal field amplitude at a meridional circulation amplitude between about $8$ and $15\;{\mathrm{ms}}^{-1}$. At these speeds, there is an optimum balance between amplification of toroidal field by differential rotation shearing, diffusive decay, and meridional transport of toroidal and poloidal flux. In each case, for larger than optimum circulation, the toroidal fields are moved in latitude and/or depth too fast to be as fully amplified as it would have been if moving more slowly; at a less than optimum advection rate, more time is allowed at a given latitude and depth for the toroidal and poloidal fields to decay due to diffusion.

Jouve & Brun (2007) computed a power-law relationship among the various parameters varied for the case of two cells in latitude and two in radius. However, for variation in flow structures in terms of not only number of cells in latitude and depth, but also the node locations in latitude and depth, it will not be possible to reach a unique answer for the power law relationship. Thus we presented here the curves as function of flow speed, diffusivity, and poloidal source amplitudes. Furthermore, some of these flow patterns do not produce solar-like features, and hence power law relationships may not be meaningful.

In all the cases we have presented so far, we have found dipolar parity during the simulation time of interest, i.e., about up to 500 years. We know that the growth rates of quadrupolar parity solutions are slightly higher than the dipolar parity in a Babcock–Leighton dynamo and the dipolar parity slowly drifts to quadrupolar if the dynamo simulations run for more than 2000 years (Dikpati & Gilman 2001; Bonanno et al. 2002; Hotta & Yokoyama 2010; Miesch & Dikpati 2014). However, Jouve & Brun (2007) have shown that this switching from dipolar to quadrupolar parity is very fast in the case of a four-celled meridional circulation that consists of two cells in latitude and two in depth. In our four-celled case consisting of two cells in latitude and two in depth (see Figure 7), the parity change did not occur so quickly. This is because the ratio of poleward surface flow speed of the top cell to that of the bottom cell was too high, ∼50, whereas in Jouve & Brun (2007) that ratio was ∼6. In order to investigate the fast change of parity in the four-celled case, we consider a ratio of poleward surface flow speed to poleward bottom flow speed to be ∼5 and simulate that case and present our results in Figures 16–19.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even though the four-celled case with two cells in latitude and two in depth produce least solar-like solutions, it is worth examining the parity issue, because such a four-celled profile could appear intermittently in one or both hemispheres of the Sun. Fast change in parity in one hemisphere with respect to the other could lead to a large phase shift between north and south cycles, as observed in the recent cycle.

In Figure 16, as in earlier figures' panel (b), shows the toroidal fields near the bottom of convective zone. Panel (c) depicts the surface radial field. We see here a radical change from the earlier examples. The toroidal and poloidal field patterns are now symmetric rather than antisymmetric about the equator. In other words, we have found quadrupolar type rather than dipole type parity. This difference in parity about the equator develops in just a few cycles, so the system in this case has a strong preference for quadrupole parity.

Jouve & Brun (2007) demonstrated, in their simulations of dynamo with four-celled (two in latitude and two in depth) meridional circulation, for what combinations of meridional flow speed and turbulent magnetic diffusivity quadrupolar structure is favored. In the case of a meridional circulation consisting of single cell in each hemisphere, Dikpati & Gilman (2001) reasoned that slow switching to quadrupolar from dipolar parity occurs when bottom poloidal fields become weak enough after a long traversal via the conveyor belt from surface to the bottom, and hence cannot connect with their opposite-hemisphere counterparts about the equator. A global statement would be that the quadrupolar mode is selected when, for a particular meridional circulation, it is dissipated at a substantially lower rate than is the dipolar mode. In other words, the growth rate for a quadrupolar mode is higher than for a dipolar mode. Since a circulation pattern that has two cells each in both latitude and depth is inherently more complicated than a simple, single-celled pattern, the poloidal and toroidal flux of opposite signs and different amplitude are converged together in more places away from the equator. So the dynamo favors quadrupolar symmetry, unless the upper and lower cells are very unequal in amplitude.

We find also in the low diffusivity case that having two cells in both latitude and radius leads to a quadrupolar solution, if the ratio of surface flow to bottom flow is 5:1 instead of 50:1 (we do not produce here the time–latitude diagram for the low-diffusivity case). Therefore, this fast switching from dipolar to quadrupolar parity is a typical phenomena in the four-celled meridional circulation with a much stronger bottom flow of similar order of magnitude as the surface flow, irrespective of turbulent diffusivity value.

Figures 17 and 18 give details of the actual transition from dipolar to quadrupolar symmetry. In Figure 17, panels (a) and (c), we see by eye that the transition appears to occur in about 120 years (16–18 sunspot cycles in this simulation) in the tachocline toroidal field and surface poloidal field respectively. Frames (b) and (d) show finer detail for each in the middle of this transition. The simulation was started from a previous dipole simulation, which inevitably has some slight departures from dipole symmetry at the truncation error level. The quadrupolar symmetry is so strongly preferred for the parameters chosen that even these small differences are enough to start the process of symmetry switching.

What we see in Figure 17 is a very simple process in which switching occurs by the northern hemisphere developing a phase lag relative to the southern hemisphere, which grows until the north lags by a sunspot or half magnetic cycle, with very little change in pattern in each hemisphere. What is actually happening is that the periods of both hemispheres are getting longer, from about 6.5–8.5 years, but the north reaches the longer period faster than the south, so the south gains on it in phase. Presumably in another simulation with different truncation errors, it could be the northern hemisphere that lags, ending up with the same final state.

In Figure 18, we show meridional cross sections of both toroidal field (blue/yellow shading) and poloidal field (solid black and dashed red lines) in 1 year intervals in the middle of the transition. We can see particularly that the peak toroidal fields are moving up in the domain to a mid-depth where the flow in low latitudes is toward the equator. In addition, the polar field reversals are evolving to a state in which they go from positive to negative at both poles at the same time.

We can quantify how long the transition takes by measuring the difference in cycle periods between south (P_S) and north (P_N) as a function of time. The normalized difference in cycle periods between south and north $({P}_{S}-{P}_{N})/[({P}_{S}+{P}_{N})/2]$ is shown as function of time in Figure 19. We see that during this transition, the north adjusts to the longer period faster than the south. The difference in cycle period is positive, and reaches a peak of about $12\%$ midway through the transition.

This switching is possible because with the circulation pattern chosen, there is much less linking of flux between hemispheres at the equator. With dipole symmetry there is strong diffusion across the equator, whereas with quadrupole symmetry there is much less since both sides have the same sign of field there. Furthermore, dipole symmetry is best preserved when the toroidal fields of both signs are strongest, where they are being brought into close proximity near the equator. This is guaranteed when there is single cell in depth in low and mid latitudes, with equatorward flow near the bottom where the turbulent magnetic diffusion is smaller.

The parity switching example shown above is very different from that found in Jouve & Brun (2007), particularly in latitudes equatorward of $45^\circ$. This is not surprising, because there is a significant structure in the flow structures—we have two strong primary cells at the top and bottom, extending from the equator to $60^\circ$ latitudes, associated with two weak secondary cells at polar regions, whereas Jouve & Brun (2007) had four equally strong cells (two in top and two in bottom) of equal latitudinal extents. We found a simple drift of phase in one hemisphere with respect to the other at all latitudes, whereas in Jouve & Brun (2007, their Figure 10) there is a radical rearranging of poloidal flux in low latitudes, with a switch in symmetries by the antisymmetric mode dying out, replaced by the symmetric mode without much change in phase locally. This difference suggests there may be multiple ways for symmetry switching to occur, which should be explored in the future.