Kinematic solar dynamo models with a deep meridional flow
Abstract
We develop two different solar dynamo models to verify the hypothesis that a deep meridional flow can restrict the apperance of sunspots below 45 degrees, proposed by Nandy & Choudhuri (2002). In the first one, a single polytropic approximation for the density profile was taken, for both radiative and convective zones. In the second one, two polytropes were used to distinguish between both zones (Pinzon & CalvoMozo, 2001). The magnetic buoyancy mechanism proposed by Dikpati & Charbonneau (1999) was chosen in both models. We, actually, have obtained that a deep meridional flow pushes the maxima of toroidal magnetic field toward the solar equator, but in contrast to Nandy & Choudhuri (2002), a second zone of maximal fields remains at the poles. The second model, although closely resembling the solar standard model of Bahcall, Pinsonneault & Wasserbug (1995); Bahcall, Pinsonneault & Basu (2001), gives solar cyles three times longer than observed.
keywords:
MHD  magnetic field  Sun: interior  Sun: magnetic fields.1 Introduction
Despite its irregular appearance, the solar magnetic cycle evolves in a spatially and temporally well organized manner. Some sunspots associated phenomena, such as the 11 years cycle of alternating polarities (Hale law), allows to conclude that the solar magnetic cycle can be explained by a dynamo process involving the transformation of a poloidal magnetic field into a toroidal one, and the later regeneration of the poloidal field, but of opposed polarity to the initial, and so on. In the BabcockLeighton approach (Babcock, 1961; Leighton, 1969) the dynamo operates in the following way: In a first stage, an initial dipolar field, with lines on meridional planes, is dragged in the eastwest direction by the action of the solar differential rotation to form a toroidal field. This happens close to the base of the solar convection zone (SCZ), in a thin layer called tachocline, where helioseismology has discovered a substancial radial shear in the rotation pattern. In a second stage, tubes of toroidal flux emerge to the surface, since these tubes of intense magnetic field are less dense than their sourrondings (magnetic bouyancy force). While the force lines rises they are twisted by the Coriolis force to form the bipolar magnetic regions (BMR) associated with sunspots in the solar surface. In the last stage, the leader portion of each BMR (i.e, that is farther ahead in the direction of solar rotation) migrates equatorward, while the follower portion migrates poleward, thanks to the meridional circulation. The general effect is that the force lines between the two BMRs lay on the poloidal direction, but oppossite to the original field. On the long run, they cancel out the original dipolar field and give rise to a new one in the opposite direction (Choudhuri, Schüsler & Dikpati, 1995; Durney, 1995, 1996, 1997; Dikpati & Charbonneau, 1999).
In the kinematic regime, where a fixed velocity field is assumed and one investigates just the evolution of the magnetic field due to that velocitiy field (without the backreaction of the magnetic field onto the solar plasma) (Choudhuri, 2000), three main ingredients are used to resemble the mass transport in the BabcockLeighton approach: differential rotation, meridional circulation and magnetic buoyancy. The first one can be suitably determined from helioseismology measurements, and resembles closely the radial shear at the tachocline and the latitudinal angular velocity distribution in the convection zone. The second one is the main flux transport agent. It is responsible to lead the poloidal flux generated in the surface to the deeper layers where toroidal flux is regenerated. Obeservational evidence tell us that there is a poleward meridional flow at the surface with an average speed between m s, but the internal return flow remains unknown (Giles et al., 1997). The last ingredient (magnetic bouyancy) summarizes the results of rising magnetic tube simulations (D’Silva & Choudhuri, 1993; Fan, Fisher & DeLuca, 1993; Caligari, MorenoInsertis & Schüsler, 1995, 1998), where a buoyant magnetic force appears, thanks to the density gradient between the inside and the outside of the tube. We will incorporate a simplified form for magnetic buoyancy which resembles an average number of such events.
Most of the recent kinematic models are able to sucessfully explain the phenomena like the reversal of the solar dipolar field in an 11years cycle or the reversal of polarities of sunspot pairs from one cycle to the next (Hale law) (Dikpati & Choudhuri, 1994; Dikpati & Charbonneau, 1999). Neverthless, they fail to predict the absence of sunspots at latitudes above 30 degrees (Spörer Law). Based on recent numerical simulations which suggest that a subadiabatic stratification, like the one in the solar radiative zone, is able to resist the magnetic buoyancy (Rempel, Schüser & Tóth, 2000), Nandy & Choudhuri (2002) proposed a dynamo model which restrics for the first time the sunspots’ appearance below 45 degrees. In this model, the meridional flow penetrates below the solar convection zone, where the toroidal magnetic field can’t erupt to the surface, and is dragged in equatorward direction by the meridional circulation. All these models distinguish for the magnetic diffusivity and the differential rotation profiles between the radiative and convective zones, but use a meridional flow based on a single polytropic approximation to the density profile for both zones.
This work investigates the hypothesis of Nandy & Choudhuri (2002) of a meridional return flow in the radiative zone in two other different models. The first one is just the same model of Nandy & Choudhuri (2002), but with the magnetic bouyancy mechanism of Dikpati & Charbonneau (1999), and is used mainly as reference considering both confined to convective zone and tachocline and deep meridional flow penetrating in the radiative zone. The second one introduces a bipolytropic approximation of the density profile to construct a meridional flow in the radiative and convective zones. This density profile better resembles the solar standard model of Bahcall, Pinsonneault & Wasserbug (1995); bahhcall01, and it is very interesting to explore how the hypotesis of a deeper meridional flow works there. In section 2 we show the mathematical formulation of the problem, based on the MHD induction equation, and define the meridional flow, differential rotation, magnetic diffusivity profiles and magnetig bouyancy mechanism to be used. Section 3 shows the results obtained with both models. Conclusions and a discussion of our results are exposed in section 4. Some details of the numerical implementation of the simulation are shown in the appendix A.
2 Mathematical Formalism
The MHD induction equation governing the evolution of the magnetic field is (Cowling, 1957)
(1) 
By assuming spherical symmetry, the magnetic and velocity fields can be writen as
(2)  
(3) 
where and correspond to the toroidal and the poloidal components of the magnetic field respectively; is the angular velocity, is the velocity in the meridional plane and is the magnetic diffusivity.
By replacing equations (2) and (3) in the induction equation (1) and by separating the poloidal and toroidal components of the magnetic field, we obtain
(4)  
(5)  
where and . It can be observed that a source term has been added (by hand) to the right side of eq. (4). This term is very necessary in our model for two reasons: it represents the magnetic buoyancy mechanism which transports the toroidal field from the base of the solar convection zone to the surface, and it allows the surface regeneration of the poloidal magnetic field. We will discuss later about the functional form and physical content of this term.
2.1 Differential rotation
As we disscused before, the helioseismology gives a good characterization of the radial shear and the latitudinal distribution of the solar angular velocity. It has been found a thin layer of substancial radial shear called tachocline, located at the base of the solar convection zone, where the toroidal magnetic field is generated ( effect). An analytical expression can be inferred from these mesurements, and we will use the one used by Dikpati & Charbonneau (1999), who were the first to include a solarlike differential rotation profile in a kinematic model,
(6) 
Here, is the latitudinal differential rotation in the surface and is an error function that confines the radial shear to a tachocline of thickness . In this expression, a rigid core rotates uniformly with angular velocity . Other values are , , nHz and .
2.2 Meridional circulation
An analytical expression for the velocity in the meridional plane is more difficult to find that the previous one, because there is not enough observational data. On one hand, the observations suggest a poleward flow with an average velocity of m/s (Giles et al., 1997) in the solar surface, but it is little what we know about the return flow, except that it must exist to satisfy mass conservation. On the other hand, numerical simulations of turbulent convection zones show that a structure of plumes in the base of the solar convection zone is able to push down the flow (and the magnetic field) to the radiative zone, where there is a net movement towards the equator, with a mean velocity of around m s (Brumell, Hurlburt & Toomre, 1998; Miesch et al., 2000). In this sense, we consider a single convection cell for each meridional quadrant, based on the next equation introduced by Dikpati & Choudhuri (1994); Choudhuri, Schüsler & Dikpati (1995) and previosly used by Nandy & Choudhuri (2002):
(7) 
where is the stream function given by
and is the density profile for the sun. Most of the models use a profile for an adiabatic gaseous sphere with a specific heat ratio coefficient , which corresponds to a constant polytropic index . So we have
(9) 
where we chose gr cm as the surface density value (Pinzon & CalvoMozo, 2001). The coefficient was chosen in such a way that the maximal latitudinal velocity at middle latitudes is m s. Other values in the eq. (2.2) are: cm, cm, , and cm. Here is the minimal coordinate value for the integration range, and the free parameter is the maximal depth of return flow (See Dikpati & Choudhuri (1994) for more details).
2.3 Magnetic buoyancy (MB)
We mentioned above that a buoyancy magnetic force upwards is generated in the base (or below) of the solar convection zone, due to of the density difference between the inner and the outer regions of a magnetic flux tube. Simulations in this sense have shown that toroidal magnetic tubes of G emerge to the surface and are twisted by Corilis forces to form tilted bipolar active magnetic regions (tilted BMRs). They have also shown that the tilt and the latitud of those BMRs strongly depend of the initial value of the toroidal magnetic field. Initial values greater than G don’t produce tilt in BMR, and values lower than emerge radially with tilts in disagreement with Joy’s law (D’Silva & Choudhuri, 1993; Fan, Fisher & DeLuca, 1993; Caligari, MorenoInsertis & Schüsler, 1995, 1998). We introduce a simplified form of including these results, due to Dikpati & Charbonneau (1999), with an initial value of G. Its analytical form is
It can be seen that this term acts nonlocally in (Durney, 1995, 1996, 1997): values of toroidal magnetic field at the tachocline produce proportional poloidal magnetic fields at a thin layer close to the surface, defined by the error (erf) functions. The last term, , anihillates the subsequent increase of at the surface, beyond some maximal level . It is, in other words, a saturation term, and it is the only source of nonlinearity in the system. The other parameters are , , , so the alpha effect is confined beneath the surface in a layer of thickness .
2.4 Magnetic diffusivity
Magnetic diffusivity is different for the radiative and the convective zones. The turbulent regime present at the convective zone makes two orders of magnitude higher than the radiative one (see fig 1 of Dikpati & Charbonneau (1999)),
(11) 
with and cm s and .
3 Results
We solved the model described above by means of the ADI method on a twodimensional mesh of spatial divitions, with and , and by imposing the following boundary conditions:
at  (12)  
at  (13)  
at  (14)  
at  (15) 
At the upper boundary it must be ensured, in addition, that satisfies the free space condition,
(16) 
The initial condition is not determinant for the evolution of the model, as long as the coefficient is large enough to avoid the magnetic field to diffuse. In this case the system always relaxes to peridic solutions. In all cases, we took the initial condition
(17) 
The code was tested by using all profiles and parameters of Dikpati & Charbonneau (1999), and by reproducing the results therein. The numerical details on discretization and integration are shown in appendix A.
3.1 Meridional flow confined to the convective zone
In this first model, the meridional flow penetrates until , that is, it is confined to the solar convection zone and the tachocline. The free parameters used in this model are shown on the left side of Table LABEL:table1. As the butterfly diagram shows (Fig 1A), the maximal toroidal magnetic field in the base of SCZ is located at latitudes between and degrees, with weak branches migrating in the equatorward direction. Under the primordial assumption that the toroidal field lines emerge to the surface like magnetic flux tubes, this model generates sunspots close to the poles. The rigth side of Table LABEL:table1 summarizes the most important results of this model. We obtain a period slightly larger than the one observed and the maximal intensities of the magnetic fields are in the correct order of magnitude for the toroidal one but one order the magnitude over the one expected for the radial field.
Parameter  Value  Result  Value 

cm s  max  G  
cm s  max  G  
cm s  T  years  
3.2 Deep meridional flow
A different situation is observed if the meridional flow penetrates into the radiative zone, as deep as . A better distribution of the toroidal lines is obtained, however, a polar branch still remains, in contrast to Nandy & Choudhuri (2002). There are two zones of maximal intensity: one below degrees and another above degrees (at the poles). This result can be interpreted in the following way: if most of the magnetic field is dragged below the tachocline, only a small portion emerges to the surface by the action of the magnetic bouyancy force. The rest migrates to the equator, dragged by the solar plasma, and will undergo the alpha effect only at middle and low latitudes.
It is important to notice that this model gives a period of years (close to the observed one) and maximal values of both toroidal and radial fields within the right order of magnitudes. The parameters we used and the results are summarized in Table 2. The large scale behaviour of the model can be observed in the butterfly diagrams of Fig 2. The equatorward migration of the branches of toroidal magnetic field is evident.
Parameter  Value  Result  Value 

cm s  max  G  
cm s  max  G  
cm s  T  years  
3.3 A more realistic density profile
As the model of Nandy and Choudhouri does (Nandy & Choudhuri, 2002), the model above employs a single polytropic density profile (Eq. 9), with (). This is just true for the convective zone (), but not for the radiative one. We can introduce a more realistic density profile by using a bipolytropic aproximation due to Pinzon & CalvoMozo (2001), with for the convective zone and for both radiative zone and core (Fig 3). We fit the numerical results of Pinzon & CalvoMozo (2001) by the analytical expression
(18) 
with the same surface density value as (9).
Parameter  Value  Result  Value 

cm s  max  G  
cm s  max  G  
cm s  T  years  
Despite that this profile closely resembles the internal structure of the solar standard model of Bahcall, Pinsonneault & Wasserbug (1995); Bahcall, Pinsonneault & Basu (2001), our results are not so good as before.
On one hand, we obtain again two maximal intensity zones: the first one very close to the poles, but the second one at middle latitudes, between and degrees, in disagreement with the observations. On the other hand, since the bipolytropic density profile is around three times larger that the usual density profile at the return point, , (see Fig 3), the counter flow is around three times lower that before in order to hold mass conservation, . Thus, our model gives a new period much larger than observed (actually, three times larger, but this is not a linear relationship). If we allow the flow to go deeper, the second zone of maximum toroidal field approaches the equator, but at the cost of increasing the period very quickly.
4 Conclusions
In this work three different models of solar dynamo are shown. We developed these models by using the velocity field of Nandy & Choudhuri (2002), the magnetic buoyancy mechanism introduced by Dikpati & Charbonneau (1999) and different values of the diffusivity coefficient for the convective and the radiative zones, as in Dikpati & Charbonneau (1999); Nandy & Choudhuri (2002). In other words we have recovered the Nandy & Choudhuri (2002) results even with a different source formulation. This suggests that the good performace of the kinematic models in the papers above is general in nature. Moreover, our results suggests that, at least for this kind of models, the effect of a deep meridional flow to solve the apperance of sunspots at high latitudes is general, too.
The most important results are:

If the meridional flow is confined to the convective zone (), the emergence latitude of sunspots is just near the poles, as expected. This model also gives wrong maximal values of the radial field intensities.

When the flow penetrates until , the model gives a better solution. However, two regions of maximal toroidal field appear (in contrast with the single region obtained by Nandy & Choudhuri (2002)): one of them is located in the right range of latitudes and the other remains very near to the poles. Except for this last region, the model approximates well the observed butterfly diagram.

In the third model we attempt to improve it by incorporating a more realistic density profile, but the results are not satisfactory. There are two zones of maximal toroidal field as before, but the second one appears at latitudes to high to agree with observations. Moreover, the period of cycle is three times larger than observed.
Our simulations show again the good effects of a deep meridional flow to solve the appearance of sunspots at high latitudes, proposed by Nandy & Choudhuri (2002), but this hypothesis deserves a more extended discussion. First, recent results on magnetic tube flux simulations (Rempel, Schüser & Tóth, 2000) suggests that a subadiabatic stratification may suppress the magnetic buoyancy force and address the equilibrium of a toroidal band. Despite that Rempel, Schüser & Tóth (2000) do not speak of a deep flow, such kind of stratification is the one to be found in the radiative zone and the tachocline, supporting the process suggested by Nandy & Choudhuri (2002) to drag the toridal field towards lower latitudes. However, they do not use a subadiabatic stratification, but an adiabatic background () everywhere. We have introduced a subadiabatic background in a simple way, by using a more realistic solar structure with a bipolytropic density profile, but our results are not in agreement with the observations. Nevertheless, this may suggest that a more realistic meridional circulation profile is needed, not that the idea of a magnetic buoyancy suppressed by a subadiabatic stratification should be rejected.
Second, it should be noticed that the discussion around a deep meridional flow hipothesis is not closed. On one hand, turbulent convection models obtain flow penetration (in a plumes dominated structure) in the stable layer (Miesch et al., 2000; Tobias et al., 2001), and some analysis on helioseismological observations give a poleward meridional flow across the entire convective zone (Giles et al., 1997), supporting the idea of a deep flow. On the other hand, the meridional flow is such a weak flow that it is very unlikely that it can penetrate the strongly subadiabatically stratified radiative core of the Sun. In addition, other problems would have to be reconsidered if a deep flow is assumed, like a larger angular momentum transfer to the radiative zone (Durney, 2000) and the changes in the relative abundance of Lythium and other elements in both radiative and convective zones (Zahn, 2001; Brun, TurckChièze & Zahn, 1999), but these are still open problems. Despite these uncertainties, our aim with this paper is to verify what would happen if this scenario were possible.
Acknowledgments
We thank G. Pinzon for his colaboration with the bipolytropic density data, P.D. Mininni and D. Nandy for help and useful discussions, and an annonymous referee for his invaluable and accurate suggestions. This work was supported by the Bogota Research Division of the National University of Colombia, project DIB803755.
References
 Ames (1977) Ames W. F., 1977, Numerical Methods for partial differential equations
 Babcock (1961) Babcock H.W., 1961, ApJ, 133, 572
 Bahcall, Pinsonneault & Wasserbug (1995) Bahcall, M., Pinsonneault M., Wasserbug G. J., 1995, Rev. Mod. Phys., 67, 781
 Bahcall, Pinsonneault & Basu (2001) Bahcall, M., Pinsonneault M., Basu S., 2001, ApJ, 555, 990
 Brumell, Hurlburt & Toomre (1998) Brumell N. H., Hurlburt N. E., Toomre J., 1998, ApJ, 493, 955
 Brun, TurckChièze & Zahn (1999) Brun A. S., TurckChièze S., Zahn J. P., 1999, ApJ, 525, 1032
 Caligari, MorenoInsertis & Schüsler (1995) Caligari P., MorenoInsertis F., Schüsler M., 1995, ApJ, 441, 886
 Caligari, MorenoInsertis & Schüsler (1998) Caligari P., MorenoInsertis F., Schüsler M., 1998, ApJ, 502, 481
 Choudhuri (2000) Choudhuri A. R., 2000, JA&A, 21, 373
 Choudhuri, Schüsler & Dikpati (1995) Choudhuri A. R., Schüsler M., Dikpati M., 1995, A&A, 303, L29
 Cowling (1957) Cowling T. G., 1957, Magnetohidrodynamics
 Dikpati & Charbonneau (1999) Dikpati M., Charbonneau P., 1999, ApJ, 518, 508
 Dikpati & Choudhuri (1994) Dikpati M., Choudhuri R. A., 1994, A&A, 291, 975
 Dikpati & Choudhuri (1995) Dikpati M., Choudhuri R. A., 1995, Solar Phys, 161, 9
 D’Silva & Choudhuri (1993) D’Silva S., Choudhuri R. A. 1993, A&A, 272, 621
 Durney (1995) Durney B. R., 1995, Solar Phys., 160, 213
 Durney (1996) Durney B. R., 1996, Solar Phys., 166, 231
 Durney (1997) Durney B. R., 1997, ApJ, 486, 1065
 Durney (2000) Durney B. R., 2000, ApJ, 528, 486
 Fan, Fisher & DeLuca (1993) Fan Y., Fisher G. H., DeLuca E.E., 1993, ApJ, 405, 390
 Giles et al. (1997) Giles P. M., Duvall Jr. T. L., Scherrer P. H., Bogart R. S., 1997, Nat, 390,52
 Leighton (1969) Leighton R., 1969, ApJ, 156, 1
 Miesch et al. (2000) Miesch M. S., Elliott J. R., Toomre J., Clune T. L., Glatzmaier G. A., Gilman P. A., 2000,ApJ, 532, 539
 Nandy & Choudhuri (2002) Nandy D., Choudhuri A. R., 2002, Sci, 296, 1671
 Pinzon & CalvoMozo (2001) Pinzon G. CalvoMozo B., 2001, astroph/0102429
 Press et al. (1992) Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical Recipes
 Rempel, Schüser & Tóth (2000) Rempel M., Schüsler M., Tóth G., 2000, A&A, 363, 789
 Tobias et al. (2001) Tobias S. M., Brumell N. H., Clune T. L., Toomre J., 2001, ApJ, 549, 1183
 Zahn (2001) Zahn J. P.,1992, A&A, 265, 115
Appendix A Numerical methods
The equations (4) and (5) are two coupled partial differential equations of advectiondiffusion type, with a nonlinear term (). There are two important aspects to have into account: numerical dicretization and numerical integration. In the first one, the advectivediffusive character of the equations must be considered. In order to ensure precision and stability for all times the advective terms are discretized by using the LaxWendroff mechanisms, which is a secondorder discretization in both space and time, and the diffusive terms are discretized by a simple forwardtime centeredspace (FTCS) mechanism (Ames, 1977; Press et al., 1992). In the second one, the two dimmensional character and the nonlinearity are considered. Alternating DirectionImplicit (ADI) method is a good numerical technique for these characteristics (Ames, 1977; Press et al., 1992).
A first step is to write the equations in operator notation, as follows:
(19)  
(20) 
where
and and are the crossed terms. is given by equation (2.3) and is given by
(25) 
In the ADI method, the time step is divided into two steps of size . In each half step, one spatial dimmension is treated implicitly and the other is treated explicitly. We treated the terms in implicit form and the terms in explicit form in the first half step, as follows:
(26)  
(27) 
where are the spatial divisions and is the time step. The next half step is treated in the reverse manner,
(28)  
(29) 
The source terms are always treated explicitly, to guarantee the linearity of equations in and . These equations can be organized in such a way that they can be solved by a standard tridiagonal algorithm. The numerical treatment of the boundary conditions is developed in detail in Dikpati & Choudhuri (1995).