A DoubleRing Algorithm for Modeling Solar Active Regions: Unifying Kinematic Dynamo Models and Surface FluxTransport Simulations
Abstract
The emergence of tilted bipolar active regions (ARs) and the dispersal of their flux, mediated via processes such as diffusion, differential rotation and meridional circulation is believed to be responsible for the reversal of the Sun’s polar field. This process (commonly known as the BabcockLeighton mechanism) is usually modeled as a nearsurface, spatially distributed effect in kinematic meanfield dynamo models. However, this formulation leads to a relationship between polar field strength and meridional flow speed which is opposite to that suggested by physical insight and predicted by surface fluxtransport simulations. With this in mind, we present an improved doublering algorithm for modeling the BabcockLeighton mechanism based on AR eruption, within the framework of an axisymmetric dynamo model. Using surface fluxtransport simulations we first show that an axisymmetric formulation – which is usually invoked in kinematic dynamo models – can reasonably approximate the surface flux dynamics. Finally, we demonstrate that our treatment of the BabcockLeighton mechanism through doublering eruption leads to an inverse relationship between polar field strength and meridional flow speed as expected, reconciling the discrepancy between surface fluxtransport simulations and kinematic dynamo models.
1 Introduction
Currently, some of the best tools for understanding the solar magnetic cycle are axisymmetric kinematic dynamo models and surface fluxtransport simulations. On the one hand kinematic dynamo models (which are usually based on an axisymmetric formulation), attempt to model the magnetic cycle selfconsistently by using a prescribed meridional flow, differential rotation, turbulent diffusivity and poloidal source (see Sec. 5). They have been successful in reproducing several of the characteristics of the solar cycle (see for example: Choudhuri, Schüssler & Dikpati 1995; Durney 1997; Dikpati & Charbonneau 1999; Covas et al. 2000; Nandy & Choudhuri 2001; Rempel 2006; Guerrero & de Gouveia Dal Pino 2007; Jouve & Brun 2007; MuñozJaramillo, Nandy & Martens 2009, MNM09 from here on; for more information about kinematic dynamo models see review by Charbonneau 2005). On the other hand, surface fluxtransport simulations study the evolution of the photospheric magnetic field by integrating the induction equation using a prescribed meridional flow, differential rotation and turbulent diffusivity. There are two main differences between surface fluxtransport simulations and kinematic dynamo models: in the former the computational domain is restricted to the surface (without imposing axisymmetry) and they are not selfexcited, but driven by the deposition of active region (AR) bipolar pairs. This type of models has proved a successful tool for understanding surface dynamics on long timescales (see, for example, Mackay, Priest & Lockwood 2002; Wang, Lean & Sheeley 2002; Schrijver, De Rosa & Title 2002) and the evolution of coronal and interplanetary magnetic field (see for example Lean, Wang & Sheeley 2002; Yeates, Mackay & van Ballegooijen 2008). However, a discrepancy exists between kinematic dynamo models and surface fluxtransport simulations regarding the relationship between meridional flow amplitude and the strength of the polar field (Schrijver & Liu 2008; Hathaway & Rightmire 2010; Jiang et al. 2010). On the one hand kinematic dynamo models find that a stronger meridional flow results in stronger polar field (Dikpati, de Toma & Gilman 2008), on the other hand surface fluxtransport simulations find an inverse relationship (Wang, Sheeley & Lean 2002; Jiang et al. 2010). In this work we improve upon the idea proposed by Durney (1997) and further elucidated by Nandy & Choudhuri (2001) of using axisymmetric ring doublets to model individual ARs. We show that this captures the surface dynamics better than the effect formulation and resolves the discrepancy between dynamo models and surface fluxtransport simulations regarding the relationship between meridional flow speed and polar field strength.
2 Evolution of the Axisymmetric Component of the Magnetic Field on Timescales Comparable to the Solar Cycle
As mentioned before, kinematic dynamo models are usually based on an axisymmetric formulation and our model is not an exception. Given that here we introduce an improved axisymmetric doublering algorithm for modeling AR eruptions (see below), but AR emergence is strictly a nonaxisymmetric process, it is important to study the amount of information lost by averaging over the longitudinal dimension. We do this by performing surface transport simulations driven by a synthetic set of AR cycles based on Kitt Peak data using the model of Yeates, Mackay & van Ballegooijen (2007). We perform a regular surface fluxtransport simulation in which the bipolar ARs are distributed all across the surface of the Sun (Case 1) and another in which the same set of ARs is deposited at the same Carrington longitude while leaving other properties (time, tilt, latitude of emergence and flux) intact (Case 2). The difference between both simulations is clear from the top row of Fig. 1, where we show a snapshot of the surface magnetic field at the peak of the cycle for Case 1 (Fig. 1a) and Case 2 (Fig. 1b). Obviously these cases have entirely different magnetic configurations at the time of deposition. However, when the magnetic field is averaged in longitude and stacked in time to create a magnetic synoptic map (also know as butterfly diagram; Figs. 1c & 1d), a careful examination shows that the results are essentially the same within a margin of 1% (Figs. 1e & f). The reason the simulations have identical outcomes is that the differential rotation and the meridional flow are both independent of longitude in the simulations. Note that nonaxisymmetry is essential for the evolution of the corona and interplanetary magnetic field. This result simply indicates that an axisymmetric representation of surface dynamics is a reasonable approximation if we are only concerned with the general properties of the magnetic field at the surface over solar cycle timescales in the context of dynamo models.
3 Modeling Individual Active Regions as Axisymmetric DoubleRings
The initial implementation of the doublering algorithm by Durney (1997) and Nandy & Choudhuri (2001) consisted in searching the bottom of the convection zone (CZ) for places in which the toroidal field exceeds a buoyant threshold and placing two axisymmetric rings of constant radial flux directly above them. This implementation had two important deficiencies: strong sensitivity to changes in grid resolution and the introduction of sharp discontinuities in the component of the vector potential. The first necessary step to address these problems is a careful mathematical definition of the vector potential associated with each ring doublet, which ensures a continuous first derivative in the computational domain. We do so by building a separable function:
(1) 
where is a constant we introduce to ensure supercritical solutions and defines the strength of the ring doublet. is defined as
(2) 
where m corresponds to the radius of the Sun and represents the penetration depth of the AR. This depth is motivated from results indicating that the disconnection of an AR fluxtube happens deep down in the CZ (Longcope & Choudhuri 2002). , on the other hand, is easier to define in integral form:
(3) 
where () defines the positive (negative) ring:
(4) 
Here is the colatitude of emergence, is the diameter of each polarity of the doublet , for which we use a fixed value of (heliocentric degrees) and is the latitudinal distance between the centers, which in turn depends on the angular distance between polarity centers and the AR tilt angle ; is calculated using the spherical law of sines (see Fig. 2a for a diagram illustrating these quantities). Fig. 2b shows the axisymmetric signature of one of such axisymmetric ARs.
4 Recreating the Poloidal Field
Given that the accumulated effect of all ARs is what regenerates the poloidal field, we need to specify an algorithm for AR eruption and decay in the context of the solar cycle. On each solar day of our simulation we randomly chose one of the latitudes with fields higher than a buoyancy threshold of Gauss at the bottom of the CZ (), and calculate the amount of magnetic flux present within its associated toroidal ring. The probability distribution we use is not uniform, but is restricted to observed active latitudes. We do this by making the probability function drop steadily to zero between 30 (30) and 40 (40) in the northern (southern) hemisphere:
(5) 
We then calculate the corresponding AR tilt, using the local field strength , the calculated flux and the latitude of emergence . For this we use the expression found by Fan, Fisher & McClymont (1994):
(6) 
reducing the magnetic field of the toroidal ring from which the AR originates. In order to do this, we first estimate how much magnetic energy is present on a partial toroidal ring (after removing a chunk with the same angular size as the emerging AR). Given that this energy is smaller than the one calculated with a full ring, we set the value of the toroidal field such that the energy of a full toroidal ring filled with the new magnetic field strength is the same as the one calculated with the old magnitude for a partial ring. Finally, we deposit a doublering (as defined in Section 3) with these calculated properties, at the chosen latitude.
5 The Kinematic MeanField Dynamo Model
We perform dynamo simulations to explore how the doublering formulation compares to the near surface effect formulation. In particular we focus on the relationship between meridional flow speed and polar field strength. Our model is based one the axisymmetric dynamo equations:
(7) 
(8) 
where A is the component of the vector potential (from which and can be obtained), B is the toroidal field (), is the meridional flow, the differential rotation, the turbulent magnetic diffusivity and . The second term on the righthand side of Equation 7 corresponds to the poloidal source in the meanfield formulation. In this formulation is a constant that sets the strength of the source term and is usually used to ensure supercritical solutions; captures the spatial properties of the BL mechanism: confinement to the surface, observed active latitudes and latitudinal dependence of tilt, while adds nonlinearity to the dynamo by quenching the source term for values of the toroidal field at the bottom of the CZ that are too strong or too weak. More information about this source can be found in MNM09. Note that for simulations using the doublering algorithm this term is not present in the equations ().
In order to integrate these equations, we need to prescribe four ingredients: meridional flow, differential rotation, the poloidal field regeneration mechanism, and turbulent magnetic diffusivity. For the differential rotation, we use the analytical form of Charbonneau et al. (1999), with a tachocline centered at whose thickness is and we use the meridional flow profile defined in MNM09. This meridional flow better captures the features present in helioseismic data, specially the latitudinal dependence. We use an amplitude of m/s for the results shown in Figure 3 and a variable amplitude for the results shown in Figure 5 (see below). We use a double stepped diffusivity profile as described in MNM09. It starts with a diffusivity value cm/s at the bottom of the CZ, jumps to a value of cm/s in the CZ, and then to a value of cm/s in the nearsurface layers. The first step is centered at and has a halfwidth of and the second step is centered at and has a halfwidth of .
For the poloidal field regeneration mechanism we use the improved ringdoublet algorithm described above, using a value of , in order to insure supercriticality (for a meridional flow of m/s). For those simulations which use an effect formulation, we use the nonlocal poloidal source described above (more information in MNM09) using a value of , in order to insure supercriticality (for a meridional flow of m/s).
6 Addressing the Discrepancy Between Kinematic Dynamo Models and Surface FluxTransport Simulations
In order to have a net accumulation of unipolar field at the poles, it is necessary to have an equal amount of flux cancellation across the equator. Since the meridional flow is poleward in the top part of the convection zone, it essentially acts as a barrier against flux cancellation by sweeping both positive and negative AR polarities towards the poles resulting in weak polar fields. This leads to an inverse correlation between flow speed and polar field strength which is accurately captured in surface flux transport simulations. Contrarily, dynamo simulations in typically used parameter regimes obtain an opposite relationship not consistent with the above physics. This is because if there is already a strong separation of flux, a fast meridional flow will lead to an enhancement of the polar field due to flux concentration. This unrealistically strong separation is typical of kinematic dynamo models that use a nonlocal effect BL source (see Fig. 3c). The reason is that by increasing the vector potential proportionally to the toroidal field at the bottom of the CZ (Eq. 7), one creates strong gradients in the vector potential above the edges of the toroidal field belt; this ends up immediately producing poloidal field which is as large in length scale as the toroidal field itself, circumventing the whole process of flux transport by circulation and diffusion. Figure 3 illustrates this fundamental difference: The top row shows the evolution of the surface magnetic field for a dynamo model using the doublering algorithm (Fig. 3a) versus one using the effect formulation (Fig. 3b). The different way in which each formulation handles the surface dynamics is evident. The doublering simulation clearly shows a mixture of polarities and smallscale features which migrate to the poles (very much like the observed evolution of the surface magnetic field). On the other hand, the mean field formulation only shows two large scale polarities whose centroids drift apart as the cycle progresses. The bottom row depicts a snapshot of the poloidal field for the doublering algorithm (Fig. 3c) and the effect formulation (Fig. 3d) – both snapshots taken at solar max. Here the presence of smallscale features and a mixture of polarities is evident for the double ring, whereas the effect formulation only shows a largescale magnetic field with two polarities. It is clear that although the large scale internal field is similar for both, the doublering algorithm does a much better job of capturing the surface dynamics.
6.1 Polar Field Strength vs. Meridional Flow Speed
In order to study the relationship between meridional flow and polar field strength, we perform simulations in which we randomly change the meridional flow amplitude from one sunspot cycle to another (between m/s). This is illustrated in Fig. 4 where a series of sunspot cycles is plotted along with their associated meridional flow. We then evaluate the correlation between the amplitude of the meridional flow of a given cycle and the polar field strength at the end of it. Since we want to evaluate the relative performance of the doublering algorithm as opposed to the nonlocal BL source, we perform the same simulation for both types of sources. Aside from the varying meridional flow amplitude and the poloidal source, the rest of the ingredients are the same. It is important to note that partly due to difficulties in tracking the exact occurrence of solar minimum, the two hemispheres eventually drift out of phase in long simulations – sometimes this phase difference leads to quadrupolar solutions which often go back to the observed dipolar solution. This parity issue only appears when the meridional flow is changed at solar minimum: if there are no variations, or if the variation takes place at solar maximum, the cycle is always locked in phase with dipolar parity. Nevertheless, to compare our simulations with surface fluxtransport models, we change the flow speed only at solar minimum. To be consistent, we accumulate statistics only from cycles in which the two hemispheres are in dipolar phase. The statistics performed for both types of source contain about 200 sunspot cycles.
The values of polar field we find using the kinematic dynamo simulations are of the order of 10 kG which is a common feature of dynamo models, which are successful in simulating the strong toroidal field necessary to produce sunspots and sustain the solar cycle (Dikpati & Charbonneau 1999; Chatterjee, Nandy & Choudhuri 2004; Jiang & Wang 2007; Jouve at al. 2008). Recent high resolution observations of the polar region have now confirmed the existence of such strong kiloGauss unipolar flux tubes (Tsuneta et al. 2008). Fig. 5 shows the results of both simulations. We find a weak positive correlation between meridional flow and polar field strength for the simulations using the nonlocal effect formulation (Fig. 5top), which is in general agreement with the results of Dikpati, de Toma & Gilman (2008). On the other hand, the simulations using the doublering formulation distinctively show a negative correlation (Fig. 5bottom), in agreement with surface fluxtransport simulations (Wang, Sheeley & Lean 2002). This clearly establishes that the discrepancy between the models is resolved by introducing the doublering algorithm and that the doublering formalism does a better job at capturing the observed surface magnetic field dynamics than the nonlocal effect formalism.
7 Concluding Remarks
In the first half of this work, we perform surface fluxtransport simulations to test the validity of the axisymmetric formulation of the kinematic dynamo problem. Our results suggest that this axisymmetric formulation captures well the surface flux dynamics over spatial and temporal scales that are relevant for the solar cycle. Building upon this we introduce an improved version of the doublering algorithm to model the BabcockLeighton mechanism for poloidal field regeneration in axisymmetric, kinematic dynamo models. We show that this new doublering formulation generates surface field evolution and polar field reversal which is in close agreement with observations. Additionally, we find that this improved treatment of the BabcockLeighton process generates an inverse relationship between meridional flow speed and polar field strength – which is suggested by simple physical arguments and also predicted by surfaceflux transport simulations. This resolves the discrepancy between kinematic dynamo models and surface fluxtransport simulations regarding the dynamics of the surface magnetic field. Since the latter drives the evolution of the corona and the heliosphere, our work opens up the possibility of coupling dynamo models of the solar cycle with coronal and heliospheric field evolution models.
8 Acknowledgements
We want to thank Aad van Ballegooijen for useful discussions that were crucial for the development of the algorithm mentioned in Section 4. The computations required for this work were performed using the resources of Montana State University and the HarvardSmithsonian Center for Astrophysics. We thank Keiji Yoshimura at MSU, and Alisdair Davey and Henry (Trae) Winter at the CfA for much appreciated technical support. This research was funded by NASA Living With a Star grant NNG05GE47G and has made extensive use of NASA’s Astrophysics Data System. D.N. acknowledges support from the Government of India through the Ramanujan Fellowship. A.R.Y. thanks the UK STFC for financial support.
Case 1  Case 2 
(a)  (b) 
(c)  (d) 
(e)  (f) 
(a)  (b) 
Doublering Algorithm  effect Formulation 
(a)  (b) 
(c)  (d) 
effect BL formulation  



Doublering algorithm  

References
 Charbonneau, P. 2005, Living Reviews in Solar Physics, 2, 2
 Charbonneau, P., ChristensenDalsgaard, J., Henning, R., Larsen, R. M., Schou, J., Thompson, M. J., & Tomczyk, S. 1999, Astrophys. J., 527, 445
 Chatterjee, P., Nandy, D., & Choudhuri, A. R. 2004, Astron. Astrophys., 427, 1019
 Choudhuri, A. R., Schussler, M., & Dikpati, M. 1995, Astron. Astrophys., 303, L29+
 Covas, E., Tavakol, R., Moss, D., & Tworkowski, A. 2000, Astron. Astrophys., 360, L21
 Dikpati, M. & Charbonneau, P. 1999, Astrophys. J., 518, 508
 Dikpati, M., de Toma, G., & Gilman, P. A. 2008, Astrophys. J., 675, 920
 Durney, B. R. 1997, Astrophys. J., 486, 1065
 Fan, Y., Fisher, G. H., & McClymont, A. N. 1994, Astrophys. J., 436, 907
 Guerrero, G. & de Gouveia Dal Pino, E. M. 2007, Astron. Astrophys., 464, 341
 Hathaway, D. H. & Rightmire, L. 2010, Science, 327, 1350
 Jiang, J., Işik, E., Cameron, R. H., Schmitt, D., & Schüssler, M. 2010, Astrophys. J., 717, 597
 Jiang, J. & Wang, J. X. 2007, Mon. Not. R. Astron. Soc., 377, 711
 Jouve, L. & Brun, A. S. 2007, Astron. Astrophys., 474, 239
 Jouve, L., Brun, A. S., Arlt, R., Brandenburg, A., Dikpati, M., Bonanno, A., Käpylä, P. J., Moss, D., Rempel, M., Gilman, P., Korpi, M. J., & Kosovichev, A. G. 2008, Astron. Astrophys., 483, 949
 Lean, J. L., Wang, Y., & Sheeley, N. R. 2002, Geophys. Res. Lett., 29, 240000
 Longcope, D. & Choudhuri, A. R. 2002, Sol. Phys., 205, 63
 Mackay, D. H., Priest, E. R., & Lockwood, M. 2002, Sol. Phys., 209, 287
 MuñozJaramillo, A., Nandy, D., & Martens, P. C. H. 2009, Astrophys. J., 698, 461
 Nandy, D. & Choudhuri, A. R. 2001, Astrophys. J., 551, 576
 Rempel, M. 2006, Astrophys. J., 647, 662
 Schrijver, C. J., De Rosa, M. L., & Title, A. M. 2002, Astrophys. J., 577, 1006
 Schrijver, C. J. & Liu, Y. 2008, Sol. Phys., 252, 19
 Tsuneta, S., Ichimoto, K., Katsukawa, Y., Lites, B. W., Matsuzaki, K., Nagata, S., Orozco Suárez, D., Shimizu, T., Shimojo, M., Shine, R. A., Suematsu, Y., Suzuki, T. K., Tarbell, T. D., & Title, A. M. 2008, Astrophys. J., 688, 1374
 Wang, Y., Lean, J., & Sheeley, Jr., N. R. 2002a, Astrophys. J., 577, L53
 Wang, Y., Sheeley, Jr., N. R., & Lean, J. 2002b, Astrophys. J., 580, 1188
 Yeates, A. R., Mackay, D. H., & van Ballegooijen, A. A. 2007, Sol. Phys., 245, 87
 —. 2008, Sol. Phys., 247, 103