HELIOSEISMIC DATA INCLUSION IN SOLAR DYNAMO MODELS

The dynamic nature of solar activity can often be traced back to the presence and evolution of magnetic fields in the Sun. The more intense magnetic fields on the order of 1000 Gauss (G) are observed to be concentrated within regions known as sunspots, which often appear in pairs of opposite magnetic polarities (Hale 1908). Sunspots have been observed regularly now for about four centuries starting with the telescopic observations of Galileo Galilei in the early 1600s. These observations point out that the number of sunspots on the solar surface varies in a cyclic fashion with an average periodicity of 11 years (Schwabe 1844), although there are variations both in the amplitude and period of this cycle. At the beginning of a cycle sunspots appear at about midlatitudes in both hemispheres (with opposite polarity orientation across the hemispheres) and then progressively appear at lower and lower latitudes as the cycle progresses (Carrington 1858) until no sunspots are seen (i.e., solar minimum). In the next cycle, the same pattern repeats, but the new cycle spots have their bipolar magnetic orientation reversed relative to the previous cycle (in both hemispheres). So considering sign as well as amplitude, the solar cycle has a period of 22 years.

There is a weaker, more diffuse component of the magnetic field outside of sunspots which is seen to have a somewhat different evolution. This field—whose radial component has been observable at the solar surface since the advent of the magnetograph—was believed to be on the order of 10 G. It is found that this field is concentrated in unipolar patches, which originate at sunspot latitudes at the time of sunspot maxima, and then moves poleward with the progress of the sunspot cycle (Babcock 1959; Bumba & Howard 1965; Howard & LaBonte 1981). The sign of this radial field of any given cycle is opposite to the old cycle polar field, which it cancels and reverses upon reaching the poles. The amplitude of this radial field achieves a maximum (at the poles) at the time of sunspot minimum (i.e., with a 90° phase difference relative to the sunspots). However, the periodicity of the cycle of this field matches the sunspot cycle period, underscoring that they are related. Recent observations by Hinode indicate that this radial field gets concentrated within unipolar flux tubes with field strength on the order of 10³ G (Tsuneta et al. 2008).

Explanations of this observations of the solar magnetic cycle rely on the field of magnetohydrodynamic dynamo theory, which seeks to address the generation and evolution of magnetic fields as a complex nonlinear process involving interactions between the magnetic field and plasma flows within the solar interior. In particular, it is now believed that the solar cycle involves the generation and recycling (feeding on the energy available in plasma motions) of two components of the magnetic field—the toroidal component and the poloidal component. In an axisymmetric spherical coordinate system, the magnetic and velocity fields can be expressed as

$$\begin{eqnarray} {\bf B} &=& B_{\phi } {\bf \hat{e}}_{\phi } + \nabla \times (A {\bf \hat{e}}_{\phi }), \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} {\bf v} &=& r\sin (\theta)\Omega {\bf \hat{e}}_{\phi } + {\bf v}_p, \end{eqnarray} \tag{ 2 }$$

where the first term on the right-hand side of Equations (1) and (2) is the toroidal component (in the case of the velocity field this corresponds to the differential rotation) and the second term is the poloidal component of the field (in the case of the velocity field this corresponds to the meridional circulation). The toroidal component of the magnetic field is thought to be produced by stretching of an initially poloidal field by the differential rotation of the Sun (the dynamo Ω-effect); subsequently, strong toroidal flux loops rise up due to magnetic buoyancy emerging through the solar surface as sunspots (Parker 1955a). To complete the dynamo cycle, the poloidal field (whose radial component is manifested as the observed vertical field on the solar surface) has to be regenerated from this toroidal field in a process that is traditionally called the dynamo α-effect. The first explanation of this α-effect was due to Parker (1955b) who suggested that helical turbulent convection in the solar convection zone (SCZ) would twist rising toroidal flux tubes into the poloidal plane, recreating the poloidal component of the magnetic field. Much has changed since this pioneering description of the first solar dynamo model by Parker, although the basic notion of the recycling of the toroidal and poloidal components remains the same.

First of all, simulations of the buoyant rise of toroidal flux tubes point out that to match the observed properties of sunspots at the solar surface, the strength of these flux tubes at the base of the SCZ has to be much more than the equipartition field strength of 10⁴ G (D'Silva & Choudhuri 1993; Fan et al. 1993). The classical dynamo α-effect due to helical turbulence is expected to be quenched for superequipartition field strengths and therefore other physical processes have to be invoked as a regeneration mechanism for the poloidal field. One of the alternatives is an idea originally due to Babcock (1961) and Leighton (1969). The Babcock and Leighton (hereby BL) model proposes that the decay and dispersal of tilted bipolar sunspot pairs at the near-surface layers, mediated by diffusion, differential rotation, and meridional circulation, can regenerate the poloidal field. This mechanism is actually observed and is simulated as a surface flux transport process that can reproduce the solar polar field reversals (Wang et al. 1989). Therefore, a synthesis of Parker's original description along with the BL mechanism for poloidal field generation is now widely accepted as a leading contender for explaining the solar dynamo mechanism (Choudhuri et al. 1995; Durney 1997; Dikpati & Charbonneau 1999; Nandy & Choudhuri 2001; Nandy 2003), although there are other alternative suggestions as well. A description of all of those is beyond the scope of this paper and interested readers are referred to the review by Charbonneau (2005).

Second, helioseismology has now mapped the solar internal rotation profile (Schou et al. 1998; Charbonneau et al. 1999), which is observed to be vary mainly in the latitudinal direction in the main body of the SCZ. Helioseismology has also discovered the tachocline—a region of strong radial and latitudinal shear beneath the base of the SCZ which is expected to play an important role in the generation and storage of strong toroidal flux tubes.

Third, more is now known about the meridional circulation, which is observed to be poleward at the surface (Hathaway 1996). To conserve mass, this circulation should turn equatorward in the solar interior. This circulation is deemed to be important for the dynamics of the solar cycle (see Hathaway et al. 2003, and the review by Nandy 2004) but the profile of this in the solar interior remains poorly constrained. Nandy & Choudhuri (2002) proposed a deep equatorward counterflow in the circulation (penetrating into the radiative interior beneath the SCZ) to better reproduce in dynamo simulations the latitudinal distribution of sunspots (because equatorward advective transport and storage of the deep-seated toroidal field is more efficient at these depths where turbulence is greatly reduced). However, Gilman & Miesch (2004), based on a laminar analysis, argued that the penetration of the circulation would be limited to a shallow Ekman layer close to the base of SCZ. A recent and more detailed analysis of the problem by Garaud & Brummell (2008) suggests that the circulation can in fact penetrate deeper down into the radiative interior. At this point there is no consensus on the profile and nature of the meridional circulation in the solar interior. Helioseismic data do provide some information about the depth dependence of this circulation at near-surface layers (Braun & Fan 1998; Giles 2000; Chou & Ladenkov 2005; González-Hernández et al. 2006), which, in conjunction with reasonable theoretical arguments, can be used to construct some plausible interior profiles of this flow.

Numerous kinematic dynamo models have been constructed in recent years (see Charbonneau 2005 for a review) incorporating these large-scale flows (differential rotation and meridional circulation) as drivers of the magnetic evolution. More recently such dynamo models (based on the BL idea of poloidal field generation) have also been utilized to make predictions for the upcoming cycle (Dikpati et al. 2006; Choudhuri et al. 2007). At present, all these kinematic dynamo models incorporate the information on large-scale flows as analytic fits to the differential rotation profile and a theoretically constructed meridional circulation profile that is subject to mass conservation but matches the flow speed only at the solar surface (i.e., without incorporating the depth-dependent information that is available). However, these large-scale flows are crucial to the generation and transport of magnetic fields; the differential rotation is the primary source of the toroidal field that creates solar active regions, and the meridional flow is thought to play a crucial role in coupling the two source regions for the poloidal and toroidal field through advective flux transport. Given this, it is obvious that the next step in constructing more sophisticated dynamo models of the solar cycle is to move toward a more rigorous use of helioseismic data to constrain these models in a way such that they conform more closely to the best available observational constraints; that is the goal of this study.

In Section 2, we describe the basic features of the kinematic dynamo model based on the BL idea that we use for our study; in this model, we use fairly standard parameterization (commonly used in the community) of various processes such as the diffusivity, dynamo α-effect, and magnetic buoyancy. In Sections 3.1 and 3.2, we present the methodologies for using the helioseismically observed solar differential rotation and constraining the meridional circulation profiles within this dynamo model and describe how they improve upon the commonly used analytic profiles. In Section 4, we present results from dynamo simulations using these improved helioseismic constraints and conclude in Section 5 with a summary of our main results and their contextual relevance.

We substitute Equations (1) and (2) into the magnetic induction equation

$$\begin{equation} \frac{\partial \textbf {B}}{\partial t} = \nabla \times \left(\textbf {v}\times \textbf {B} - \eta \nabla \times \textbf {B}\right), \end{equation} \tag{ 3 }$$

and add the phenomenological BL poloidal field source α to obtain the axisymmetric dynamo equations:

$$\begin{eqnarray} \frac{\partial A}{\partial t} + \frac{1}{s}[ \textbf {v}_p \cdot \nabla (sA)] &=& \eta \left(\nabla ^2 - \frac{1}{s^2} \right)A\nonumber\\ && + \alpha (r,\theta)F(B_{\rm av})B_{\rm av}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} && \nonumber \\ &&\frac{\partial B}{\partial t} + s\!\left[\! \textbf {v}_p \cdot \nabla \left(\frac{B}{s} \right) \!\right] + (\nabla \cdot \textbf {v}_p)B\nonumber\\ &&= \eta \!\left(\nabla ^2 - \frac{1}{s^2}\! \right) B + s\left(\left[ \nabla \times (A\bf \hat{e}_\phi) \right]\!\cdot\! \nabla \Omega \right) +\frac{1}{s}\frac{\partial (sB)}{\partial r}\frac{\partial \eta }{\partial r},\nonumber\\ \end{eqnarray} \tag{ 5 }$$

where s = rsin(θ). The terms on the left-hand side of both equations with the poloidal velocity (v_p) correspond to the advection and deformation of the magnetic field by the meridional flow. The first term on the right-hand side of both equations corresponds to the diffusion of the magnetic field. The second term on the right-hand side of both equations is the source of that type of magnetic field (BL mechanism for A and rotational shear for B). Finally, the third term on the right-hand side of Equation (5) corresponds to the advection of toroidal field due to a turbulent diffusivity gradient.

As mentioned before, the generation of poloidal field near the surface due to the decay of active regions is modeled through the inclusion of a source term, which acts as a source for the vector potential A. It is localized both in radius and latitude matching observations of active region emergence patterns and depends on the average toroidal field at the bottom of the convection zone (Dikpati & Charbonneau 1999). The radial and latitudinal dependence of the source is the following:

$$\begin{eqnarray} \alpha (r,\theta) &=& \frac{\alpha _0}{16}\cos (\theta) \left(1 + {\rm erf}\left(\frac{\theta - (90^\circ - \beta)}{\gamma } \right) \right)\nonumber\\ &&\times\left(1 - {\rm erf}\left(\frac{\theta - (90^\circ + \beta)}{\gamma } \right) \right) \nonumber\\ && \times \left(1 + {\rm erf}\left(\frac{r - r_{\rm al}}{d_{\rm al}}\right) \right)\left(1 - {\rm erf}\left(\frac{r - r_{\rm ah}}{d_{\rm ah}} \right) \right),\nonumber\\ \end{eqnarray} \tag{ 6 }$$

where α₀ sets the strength of the source term and we set it to the value in which our system starts having oscillating solutions. The parameters β = 40° and γ = 10° characterize the latitudes in which sunspots appear. On the other hand, r_al = 0.94 R_☉, d_al = 0.04 R_☉, r_ah = R_☉, and d_ah = 0.01 R_☉ characterize the radial extent of the region in which the poloidal field is deposited (see Figure 1(a) and (b)). Besides radial and latitudinal dependence, we also introduce lower and upper operating thresholds on the poloidal source that is dependent on toroidal field amplitude—which is a more realistic representation of the physics involving magnetic buoyancy (as argued in Nandy & Choudhuri 2001; Nandy 2002; see also Charbonneau et al. 2005). The presence of a lower threshold is in response to the fact that the plasma density inside weak flux tubes is not low enough (compared with the density of the surrounding plasma) to make them unstable to buoyancy (Caligari et al. 1995) and those that manage to rise have very long rising times (Fan et al. 1993). On the other hand, if flux tubes are too strong they are not tilted enough when they reach the surface to contribute to poloidal field generation (D'Silva & Choudhuri 1993; Fan et al. 1993). The dependence of the poloidal source on magnetic field is

$$\begin{equation} F(B_{\rm av}) = \frac{K_{\rm ae}}{1 + \left(B_{\rm av}/B_h \right)^2}\left(1 - \frac{1}{1 + \left(B_{\rm av}/B_l \right)^2} \right), \end{equation} \tag{ 7 }$$

where K_ae = 1/max(F(B_av)) is a normalization constant and B_h = 1.5 × 10⁵ G and B_l = 4 × 10⁴ G are the operating thresholds (see Figure 1(c)).

Zoom In Zoom Out Reset image size

Another ingredient of this model is a radially dependent magnetic diffusivity; in this work we use a double-step profile (see Figure 2) given by

$$\begin{eqnarray} \eta (r) &=& \eta _{\rm bcd} + \frac{\eta _{\rm cz} - \eta _{\rm bcd}}{2}\!\left(\!1 + {\rm erf}\left(\frac{r - r_{\rm cz}}{d_{\rm cz}}\! \right)\! \right)\nonumber\\ &&+ \frac{\eta _{\rm sg} - \eta _{\rm cz} - \eta _{\rm bcd}}{2} \left(1 + {\rm erf}\left(\frac{r - r_{\rm sg}}{d_{\rm sg}}\right) \right), \end{eqnarray} \tag{ 8 }$$

where η_bcd = 10⁸ cm² s⁻¹ corresponds to the diffusivity at the bottom of the computational domain, η_cz = 10¹¹ cm² s⁻¹ corresponds to the diffusivity in the convection zone, η_sg = 10¹³ cm² s⁻¹ corresponds to the supergranular diffusivity, and r_cz = 0.73 R_☉, d_cz = 0.03 R_☉, r_sg = 0.95 R_☉, and d_sg = 0.05 R_☉ characterize the transitions from one value of diffusivity to the other. Although it is common these days to use these sort of profile, a point we note in passing is that the exact depth dependence of turbulent diffusivity within the SCZ is poorly, if at all, constrained. Besides the value of the supergranular diffusivity η_sg, which can be estimated by observations (Hathaway & Choudhary 2008), the others parameters are not necessarily well constrained.

Zoom In Zoom Out Reset image size

Once all ingredients are defined (see Section 3 for the flow field), we solve the dynamo equations using a recently developed and novel numerical technique called exponential propagation (see the Appendix). Our computational domain is defined in a 250 × 250 grid comprising only one hemisphere. Since we run our simulations only in one hemisphere our latitudinal boundary conditions at the equator (θ = π/2) are ∂A/∂θ = 0 and B = 0. Furthermore, since the equations we are solving are axisymmetric, both the potential vector and the toroidal field need to be zero (A = 0 and B = 0) at the pole (θ = 0), to avoid singularities in spherical coordinates. For the lower boundary condition (r = 0.55 R_☉), we assume a perfectly conducting core, such that both the poloidal field and the toroidal field vanish there (i.e., A, B = 0 at the lower boundary). For the upper boundary condition, we assume that the magnetic field has only a radial component (B = 0 and ∂(rA)/∂r = 0); this condition has been found necessary for stress balance between subsurface and coronal magnetic fields (for more details refer to van Ballegooijen & Mackay 2007). As initial conditions we set A = 0 throughput our computational domain and B ∝ sin(2θ)*sin(π*((r − 0.55 R_☉)/(R_☉ − 0.55 R_☉))). After a few cycles, all transients related to the initial conditions typically disappear and the dynamo settles into regular oscillatory solutions whose properties are determined by the parameters in the dynamo equations.

The last two ingredients of the dynamo model are the velocity fields (differential rotation and meridional circulation), which are the focus of this work. The differential rotation is probably the best constrained of all dynamo ingredients but the actual helioseismology data have never before been used directly in dynamo model, only an analytical fit to it. We discuss below how the actual rotation data can be used directly within dynamo models through the use of a weighting function to filter out the observational data in the region where it cannot be trusted. On the other hand, the meridional circulation is one of the most loosely constrained ingredients of the dynamo. Traditionally, only the peak surface flow speed is used to constrain the analytical functions that are used to parameterize it, in conjunction with mass conservation. In this work, we take advantage of the properties of such functions and make a fit to the helioseismic data on the meridional flow that constrains the location and extent of the polar downflow and equatorial upflow, as well as the radial dependence of the meridional flow near the surface—thereby taking steps toward better constrained flow profiles.

3.1. Differential Rotation

As opposed to meridional circulation, there are helioseismic measurements of the differential rotation for most of the convective envelope which can be used directly in our simulations. Here we use data from the Global Oscillation Network Group (GONG; courtesy Dr. Rachel Howe) obtained using the Regularized Least Squares (RLS) inversion mapped onto a 51 × 51 grid (see Figure 3(a)). However, these observations cannot be trusted fully in the region within 0.3 R_☉ of the rotation axis (specifically at high latitudes), because the inversion kernels have very little amplitude there. Below, we outline a method to deal with this suspect data by creating a composite rotation profile that replaces these data at high latitudes with plausible synthetic data, which smoothly matches to the observations in the region of trust.

Zoom In Zoom Out Reset image size

3.1.1. Adaptation of the Data to the Model

In the first step, we use a splines interpolation in order to map the data to the resolution of our simulation (a grid of 250 × 250, see Figure 3(a)). The next step is to make a composite with the data and the analytical form of Charbonneau et al. (1999; see Figure 3(b)). The analytical form is defined as

$$\begin{eqnarray} \Omega _A(r,\theta) &=& 2\pi \bigg[\Omega _{c} + \frac{1}{2}\bigg(1 - \texttt {erf}\bigg(\frac{r - r_{\rm tc}}{w_{\rm tc}} \bigg)\bigg)\nonumber\\ &&\times (\Omega _{e} - \Omega _{c} + (\Omega _{p} - \Omega _{e})\Omega _S(\theta))\bigg]\nonumber\\ &&\times \Omega _S(\theta) = a\cos ^2(\theta) + (1-a)\cos ^4(\theta), \end{eqnarray} \tag{ 9 }$$

where Ω_c = 432 nHz is the rotation frequency of the core, Ω_e = 470 nHz is the rotation frequency of the equator, Ω_p = 330 nHz is the rotation frequency of the pole, a = 0.483 is the strength of the cos²(θ) with respect to the cos⁴(θ) term, r_tc = 0.716 the location of the tachocline, and w_tc = 0.03. We use the parameters defining the tachocline's location and thickness as reported by Charbonneau et al. (1999) for a latitude of 60°. This is because they match the data better at high latitudes (which is the place where the data merge with the analytical profile) than those reported for the equator. This composite replaces the suspect data at high latitudes within 0.3 R_☉ of the rotation axis with that of the analytic profile. However, it is important to note that at low latitudes, within the convection zone, the actual helioseismic data are utilized.

In order to make the composite, we create a weighting function m(r,  θ) with values between 0 and 1 for each gridpoint defining how much information will come from the RLS data and how much from the analytical form (see Figure 3(d)). We define the weighting function in the following way:

$$\begin{equation} m(r,\theta) = 1 - \frac{1}{2}\left(1 + \texttt {erf}\left(\frac{r^2\cos (2\theta) - r_{m}^2}{d_{m}^2}\right)\right), \end{equation} \tag{ 10 }$$

where r_m = 0.5 R_☉ is a parameter that controls the center of the transition and d_m = 0.6 R_☉ controls the thickness. The resultant differential rotation profile, which can be seen in Figure 3(c), is then calculated using the following expression:

$$\begin{equation} \Omega (r,\theta) = m(r,\theta)\Omega _{\rm RLS}(r,\theta) + [1 - m(r,\theta)]\Omega _A(r,\theta). \end{equation} \tag{ 11 }$$

3.1.2. Differences Between the Analytical Profile and the Composite Data

It is instructive to compare the analytical and composite profiles with the actual helioseismology data. In Figure 4(a), we present the residual error of subtracting the analytical profile of Charbonneau et al. (1999), with r_tc = 0.7 and w_tc = 0.025, from the RLS data. In Figure 4(b), we present the residual error of subtracting our composite from the RLS data. As is expected, there is no difference between the composite and raw data at low latitude, but the residual increases as we approach the rotation axis—where the RLS data cannot be trusted. For the analytical profile, the residuals errors are more significant, even at low latitudes. This demonstrates the ability of our methodology to usefully integrate the helioseismic data for differential rotation.

Zoom In Zoom Out Reset image size

3.2. Meridional Circulation

The meridional flow profile remains rather poorly constrained in the solar interior even though the available helioseismic data can be used to constrain the analytic flow profiles that are in use currently. Here, we present ways to betters constrain this profile with helioseismic data. In order to do that, we use data from GONG (courtesy Dr. Irene González-Hernández) obtained using the ring-diagrams technique. These data, which we can see in Figure 6(a), correspond to a time average of the meridional flow between 2001 and 2006; and it comprises 19 values of r from R_☉ down to a depth of 0.97 R_☉, and 15 different latitudes between −525 and 525. It is important to note that our work relies heavily in the assumption that the meridional flow is adequately described by a stream function with separable variables. This is consistent with the assumption present implicitly in all work on axisymmetric solar dynamo models up to this date. Below, we use this property of our stream function, along with weighted latitudinal and radial averages of the data, to completely constrain its latitudinal dependence, as well as the topmost 10% of its radial dependence. As these data currently do not constrain the depth of penetration of the flow in the deep solar interior, we explore two different plausible penetration depths of the circulation. For reasons described later, we choose to perform simulations with two different peak meridional flow speeds, therefore exploring four plausible meridional flow profiles altogether.

3.2.1. Constraining the Latitudinal Dependence of the Meridional Flow

Meridional circulation has been typically implemented in these type of dynamo models by using a stream function in combination with mass conservation, i.e.,

$$\begin{equation} \overrightarrow{\textbf {v}}_p\left(r,\theta \right) = \frac{1}{\rho (r)}\vec{\nabla }\times (\Psi (r,\theta)\widehat{\textbf {e}}_{\phi }). \end{equation} \tag{ 12 }$$

The two stream functions that are commonly used were proposed by van Ballegooijen & Choudhuri (1988) and Dikpati & Choudhuri (1995). They have in common the separability of variables and thus can be written in the following way:

$$\begin{equation} \Psi (r,\theta)=v_0F(r)G(\theta), \end{equation} \tag{ 13 }$$

where v₀ is a constant which controls the amplitude of the meridional flow.

Using such a stream function the components of the meridional flow become

$$\begin{equation} v_r(r,\theta) = v_0\frac{F(r)}{r\rho (r)}\frac{1}{\sin (\theta)}\frac{\partial }{\partial \theta }\left(\sin (\theta)G(\theta)\right), \end{equation} \tag{ 14 }$$

$$\begin{equation} v_{\theta }(r,\theta) = -v_0\frac{1}{r\rho (r)}\frac{\partial }{\partial r}(r F(r))G(\theta), \end{equation} \tag{ 15 }$$

which can themselves be separated into the multiplication of exclusively radially and latitudinally dependent functions. This property allows us to constrain the entire latitudinal dependence of this family of functions by using the available helioseismology data for v_θ at the surface. This can be done because the latitudinal dependence of v_θ is exactly the same as that of the stream function and only the amplitude of this functional form changes with depth, see for example Figure 5 for the latitudinal velocity at different depths used by van Ballegooijen & Choudhuri (1988). In this work, we assume a latitudinal dependence like the one they used, i.e.,

$$\begin{equation} G(\theta) = \sin ^{(q+1)}(\theta)\cos (\theta). \end{equation} \tag{ 16 }$$

Zoom In Zoom Out Reset image size

In order to estimate the parameter q, we first take a density-weighted average of the helioseismic data using the values for solar density from the solar model S (Christensen-Dalsgaard et al. 1996), such that

$$\begin{equation} \bar{v}_\theta (\theta _j) = \frac{\sum _i v_{\theta }(r_i,\theta _j)\rho (r_i)}{\sum _i \rho (r_i)}. \end{equation} \tag{ 17 }$$

In Figure 6(b), we plot the meridional flow at different depths weighted by density, which we add for each latitude in order to find the average. We then use this average to make a least-squares fit to the analytical expression (Equation (16)), which we can see in Figure 6(c). We find that a value of q = 1 fits the data best. This therefore constrains the latitudinal (θ) dependence of the flow profile.

Zoom In Zoom Out Reset image size

3.2.2. Constraining the Radial Dependence of the Meridional Flow

As opposed to the latitudinal dependence, the radial dependence of the meridional flow is less constrained since there is no data below 0.97 R_☉. However, at least some of the parameters can be constrained: we start with the solar density, for which we perform a least-squares fit to the solar model S using the following expression:

$$\begin{equation} \rho (r) \sim \left(\frac{R_\odot }{r} - \gamma \right)^{m}, \end{equation} \tag{ 18 }$$

we find that values of γ = 0.9665 and m = 1.911 fit the model best (see Figure 7(c)).

Zoom In Zoom Out Reset image size

In the second step, we constrain the radial dependence of the stream function. We begin with the function

$$\begin{equation} F(r) = \frac{1}{r}(r - R_p)(r - R_\odot)\sin \left(\pi \frac{r - R_p}{R_1-R_p}\right)^a, \end{equation} \tag{ 19 }$$

where R_☉ corresponds to the solar radius, R_p to the maximum penetration depth of the meridional flow, and a and R₁ control the location in radius of the poleward peak and the value of the meridional flow at the surface. In order to constrain them we use the helioseismic data again, but this time we use the radial dependence (see Figure 7(a)). We first remove the latitudinal dependence, which we do by dividing the data of each latitude by G(θ) using the value of q = 1 found in the previous section. In Figure 7(b), we plot the flow data after removing the latitudinal dependence, note that there is no longer any sign difference between the two hemispheres. The next step is to generate the latitudinal average which we can see as black dots in Figure 7(c). It is evident from looking at the radial dependence of the meridional flow that the velocity increases with depth for most latitudes, and that the point of maximum velocity is not within the depth up to which the data extends. Since the exact radial dependence of the data is too complex for our functional form to grasp, the features we concentrate on reproducing are the presence of a maximum inside the convection zone, as well as the amplitude of the flow at the near-surface layers. The logic here is to use the fewest possible parameters and a simple, physically transparent profile that does a reasonable job of matching the data. In view of the lack of better constraints, we assume here that the peak of the return flow is at 0.97 R_☉ (which is the depth at which the current helioseismic data have its peak).

Following the procedures and steps above, we construct profiles to fit the depth dependence pointed out in the available helioseismic data; however, this does not constrain how far the meridional flow can penetrate and therefore we try two different penetration depths, one shallow 0.71 R_☉ (i.e., barely beneath the base of the SCZ) and one deep R_p = 0.64 R_☉ (into the radiative interior)—both of which match the latitudinal and radial constraints as deduced from the near-surface helioseismic data. Note that observed light-element abundance ratios limit the depth of penetration of the circulation to about R_p = 0.62 R_☉ (Charbonneau 2007). Now we know that the meridional flow speed is highly variable, with fluctuations that can be quite significant and the measured flow speed can change depending on the phase of the solar cycle (Hathaway 1996; Gizon & Rempel 2008). The magnetic fields are also expected to feedback on the flow (Rempel 2006). Taken together these considerations point out that the effective meridional flow speed to be used in dynamo simulations could be less than that implied from the González-Hernández et al. data, a possibility that is borne out by the work of Braun & Fan (1998) and Gizon & Rempel (2008), who find much lower peak flow speeds in the range 12–15 m s⁻¹. Keeping this in mind, we use the same latitudinal and radial constraint as deduced earlier, but consider in our simulations two additional profiles with peak flow speeds of 12 m s⁻¹ with deep and shallow penetrations. Therefore, in total we explore four plausible meridional flow profiles in our dynamo simulations (see Figure 7(c) and Table 1 for an overview), the results of which are presented in the next section.

Table 1. Sets of Parameters Characterizing the Different Meridional Flow Profiles used in Our Dynamo Simulations

	vo (m s−1)	Rp (R☉)	a	R1 (R☉)
Set 1	12	0.64	1.795	1.027
Set 2	22
Set 3	12	0.71	2.03	1.03
Set 4	22

Note. v_o corresponds to the meridional flow peak speed, R_p the maximum penetration of the flow, and a and R₁ are the parameters that control the location of the poleward flow as well as the surface speed.

Download table as:  ASCIITypeset image

4.1. Analytic versus Helioseismic Composite Differential Rotation

We first compare dynamo solutions found using the helioseismology composite of the differential rotation (Section 3.1.1) and the analytical profile of Charbonneau et al. (1999). In doing so, we perform simulations with a model as described in Section 2 (see also numerical methods in the Appendix) and generate field evolution maps for the toroidal and poloidal fields. From the simulated butterfly diagrams for the evolution of the toroidal field at the base of the SCZ and the surface radial field evolution (Figure 8), we find that large-scale features of the simulated solar cycle are generally similar across the two different rotation profiles (even with different meridional flow penetration depths and speeds), specially for the shallowest penetration (sets 3 and 4 with R_p = 0.71 R_☉).

Zoom In Zoom Out Reset image size

In order to understand this similarity, it is useful to look at Figures 10 and 11 corresponding to simulations using the composite differential rotation, and the meridional flow sets 1 (R_p = 0.64 R_☉, v_o = 12 m s⁻¹) and 4 (R_p = 0.71 R_☉, v_o = 22 m s⁻¹). The first two columns, from left to right of both figures show the evolution of the shear sources (B_r · ∇_rΩ) and (B_θ · ∇_θΩ)—which contribute to toroidal field generation by stretching of the poloidal field. It is evident that the location and strength of these sources is different—the radial shear source is mainly present near the surface whereas the latitudinal shear is spread throughout the convection zone. It can also be seen that the radial shear source is roughly five times stronger than the latitudinal. However, if attention is paid to the evolution of the toroidal field (third column from the right in both figures), it is clear that this radial shear term has no significant impact on the structure and magnitude of the toroidal field (a similar result was found by Dikpati et al. 2002). The reason is that the upper boundary condition (B = 0 at r = R_☉), in combination with the high turbulent diffusivity (and thus short diffusive timescale) there, imposes itself very quickly on the toroidal field generated by the radial shear—washing it out. This greatly reduces the relative role of the surface shear as a source of toroidal magnetic field, effectively making the surface dynamics very similar across simulations using the analytical rotation profile (without any surface shear) or the composite helioseismic profile.

Once the surface layers are ruled out as important sources of toroidal field generation we are faced with the fact that the strongest source of toroidal field is the latitudinal shear inside the convection zone (and not the tachocline radial shear), as is evident in Figures 10 and 11. This goes against the commonly held perception that the tachocline is where most of the toroidal field is produced. However, the importance of the latitudinal shear term in the SCZ is clearly demonstrated here where we have plotted the shear source terms, which is not normally done by dynamo modelers.

The establishment of the SCZ as an importance source region of toroidal field is of relevance when regarding the similarity of the solutions obtained using the composite data and the analytical profile. This is because the region where the shear of the analytical profile and the helioseismic composite data differ most (both for radial and latitudinal shear) is in the tachocline. This is evident in Figure 9 where we plot the residual of subtracting the radial and latitudinal shear of the analytical profile from the shear of the composite data. The similarity between solutions is specially important for shallow meridional flow profiles with low penetration—which does not transport any poloidal field into the deeper tachocline, thereby further diminishing the role of the tachocline shear.

Zoom In Zoom Out Reset image size

4.2. Shallow versus Deep Penetration of the Meridional Flow

In the second part of our work, we compare dynamo solutions obtained for each of the four different meridional flow profiles with two different penetration depth and with two different peak flow speeds. First, from Figure 8 it is evident that the shape of the solutions changes with varying penetration depth; this is caused by the increasing role of the tachocline shear in generation and storage of the toroidal field as the penetration depth of the meridional flow increases. This is apparent when one compares Figure 10 (deep penetration) to Figure 11 (shallow penetration). If we look at the inductive shear sources it is evident that for the shallowest penetration no field is being generated inside the tachocline. In the poloidal field plots (right column) we see that no poloidal field is advected into the tachocline region for the shallowest penetration, but some is advected for the deepest.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Second, we compare the periods of our solutions in Table 2: we find that most solutions have a sunspot cycle (i.e., half-dynamo cycle) period that is comparable to that of the Sun, with the exception of the fast flow with deep penetration and the slow flow with shallow penetration, which have, respectively, a comparatively smaller and larger period. As the meridional flow is buried deeper, one expects the length of the advective circuit to increase, thereby resulting in larger dynamo periods. However, it is evident that as we increase the penetration, the period decreases—even if the length of the flow loop that supposedly transports magnetic flux increases; this is counterintuitive but has a simple explanation. Our simulations and exhaustive analysis points out that it is not how deep the flow penetrates that governs the cycle period, but it is the magnitude of the meridional counterflow right at the base of the SCZ (R_p = 0.713 R_☉) that is most relevant. This is because most of the poloidal field creation at near-surface layers is coupled to buoyant eruptions of toroidal field from this layer of equatorward migrating toroidal field belt at the base of the SCZ (see third column in Figures 10 and 11) and it is the flow speed at this region that governs the dynamo period. In Figure 7, it is clear that the speed of the counterflow in this convection zone–radiative interior interface increases as the flow becomes more penetrating (a consequence of the constraints set by mass conservation and the fits to the near-surface helioseismic data), thereby reducing the dynamo period.

Table 2. Simulated Sunspot Cycle Period for the Different Sets of Meridional Flow Parameters

	Rp (R☉)	vo (m s−1)	τ—HS Data (yr)	τ—Analytical (yr)
Set 1	0.64	12	9.67	10.00
Set 2		22	5.63	5.67
Set 3	0.71	12	14.67	14.02
Set 4		22	11.85	12.85

Notes. R_p corresponds to the maximum penetration depth of the meridional flow, v_o to the peak speed in the poleward flow, and τ is the period of the solutions in units of years. For the rest of the parameters in each set please refer to Table 1.

Download table as:  ASCIITypeset image

Overall, an evaluation of the butterfly diagram (Figure 8), points out that the toroidal field belt extends to lower latitudes (where sunspots are observed) for deeper penetrating meridional flow, although there is a polar branch as well. For the shallow flow, we find that the toroidal field belt is concentrated around midlatitudes with almost symmetrical polar and equatorial branches—a signature of the convection zone latitudinal shear producing most of the toroidal field (as in the interface dynamo models, see, e.g., Parker 1993; Charbonneau & MacGregor 1997).

4.3. Dependence of the Solutions on Changes in the Turbulent Diffusivity Profile

Although a detailed exploration of the turbulent diffusivity parameter space is outside the scope of this work, we study two special cases in which we vary a single parameter while leaving the rest fixed.

For the first case, we lower the diffusivity in the convection zone, η_cz, from 10¹¹ cm² s⁻¹ to 10¹⁰ cm² s⁻¹. As can be seen in Figure 12, this introduces two important changes in the dynamo solutions: the first one is an overall increase in magnetic field magnitude due to the reduction in diffusive decay while keeping the strength of the field sources constant. The second is a drastic increase in the dynamo period (which can be seen tabulated in Table 3) for the solutions that use a meridional flow with a penetration of R_p = 0.71 R_☉. The reason behind such a change resides in the nature of the transport processes at the bottom of the convection zone, which are a combination of both advection and diffusion. In the case of flow profiles with deep penetration the velocity at the bottom is high enough for downward advection to transport flux into the tachocline. On the other hand, in the case of low penetration, the last bit of downward transport into the tachocline is done by diffusive transport and thus dominated by diffusive timescales. Because of this, by decreasing turbulent diffusivity by an order of magnitude, we drastically increase the period of the solutions.

Zoom In Zoom Out Reset image size

Table 3. Simulated Sunspot Cycle Period for the Different Sets of Meridional Flow Parameters when using a Low Diffusivity in the Convection Zone (η_cz = 10¹⁰ cm² s⁻¹)

	Rp (R☉)	vo (m s−1)	τ—HS Data (yr)	τ—Analytical (yr)
Set 1	0.64	12	12.46	12.86
Set 2		22	6.47	6.55
Set 3	0.71	12	87.78	90.92
Set 4		22	70.48	74.41

Note. R_p corresponds to the maximum penetration depth of the meridional flow, v_o to the peak speed in the poleward flow, and τ is the period of the solutions in units of years.

Download table as:  ASCIITypeset image

In the second parameter space experiment, we lower the supergranular diffusivity η_sg from 10¹³ cm² s⁻¹ to 10¹¹ cm² s⁻¹. This was done to study the impact of the surface radial shear under low diffusivity conditions. However, as can be seen in Figure 13, there is very little difference between the two solutions. This means that even after reducing the supergranular diffusivity by 2 orders of magnitude, the radial shear has very little impact on the solutions and the upper boundary conditions still play an important role in limiting the relative contribution from the near-surface shear layer.

Zoom In Zoom Out Reset image size

In summary, we have presented here methods which can be used to better integrate helioseismic data into kinematic dynamo models. In particular, we have demonstrated that using a composite between helioseismic data and an analytical profile for the differential rotation, we can directly use the helioseismic rotation data in the region of trust and substitute the suspect data by smoothly matching it to the analytical profile where the data are noisy. This paves the way for including the helioseismically inferred rotation profile directly in dynamo simulations. We have also shown how mathematical properties of the commonly used analytic stream functions describing the meridional flow can be fitted to the available near-surface helioseismic data to entirely constrain the latitudinal dependence of the meridional flow, as well as weakly constrain the radial (depth) dependence.

In our simulations, comparing the helioseismic data for the differential rotation with the analytical profile of Charbonneau et al. (1999), with four plausible meridional flow profiles, we find that there is little difference between the solutions using the helioseismic composite and the analytical differential rotation profile—specially for shallow penetrations of the meridional flow and even at reduced supergranular diffusivity. This is because the impact of the surface radial shear, which is present in the helioseismic composite but not the analytic profile, is greatly reduced by the proximity of the upper boundary conditions. Also, for the shallow circulation, the toroidal field generation occurs in a region located above the tachocline with mainly latitudinal shear, where the difference between the composite data and the analytical profile is not significant.

The main result from this comparative analysis is that the latitudinal shear in the rotation is the most dominant source of toroidal field generation in these type of models that are characterized by high diffusivity at near-surface layers, but lower diffusivity within the bulk of the SCZ—specially near the base where most of the toroidal field is being created. Since this latitudinal shear exists throughout the convection zone, an interesting question is whether toroidal fields can be stored there long enough to be amplified to high values by the shear in the rotation, without being removed by magnetic buoyancy. If this were to be the case, i.e., the latitudinal shear is indeed confirmed to be the dominant source of toroidal field induction, we anticipate then that downward flux pumping (Tobias et al. 2001; see also Guerrero & de Gouveia Dal Pino 2008)—which tends to act against buoyant removal of flux, may have an important role to play in this context. This could also call into question the widely held view that the solar tachocline is where most of the toroidal field is created and stored (see Brandenburg 2005 for arguments favoring a more distributed dynamo action throughout the SCZ).

Our attempts to integrate helioseismic meridional flow data into dynamo models and related simulations have uncovered points that are both encouraging and discouraging.

On the discouraging side, we find that the currently available observational data are inadequate to constrain the nature and exact profile of the deep meridional flow, especially the return flow. Neither do the simulation results and their comparison with observed features of the solar cycle clearly support or rule out any possibility. A recent analysis on light-element depletion due to transport by meridional circulation indicates that solar light-element abundance observations restrict the penetration to 0.62 R_☉ (Charbonneau 2007); however, this analysis does not necessarily suggest that the flow does penetrate that deep. Also vexing is the fact that different inversions, involving different helioseismic techniques such as ring-diagram or time–distance analysis recovers different profiles and widely varying peak meridional flow speeds (Giles et al. 1997; Braun & Fan 1998; González-Hernández et al. 2006; Gizon & Rempel 2008). In our analysis, we chose to use the González-Hernández et al. data because at present, this provides the (relatively) deepest full inversion of the flow within the SCZ. Chou & Ladenkov (2005) reported time–distance diagrams reaching a depth of 0.79 R_☉ but have not yet reported a full inversion that could be used on our simulations.

We point out that there is an important consequence of the presence of the flow speed maximum inside the convection zone—which is related to mass conservation: if the maximum poleward flow speed is found to be deeper inside the convection zone this would result in a stronger mass flux poleward, which needs to be balanced by a deeper counterflow subject to mass conservation; the density of the plasma increases rapidly as one goes deeper, e.g., the density at 0.97 R_☉ is 10,000 times larger than at the surface. Although that is not achieved currently, our extensive efforts to fit the data point out that stronger constraints on the return flow may be achieved even with data those do not necessarily go down to where this return flow is located, a fact that may be usefully utilized when better depth-dependent helioseismic data on meridional circulation become available.

Although the depth of penetration of the circulation is an important constraint on the flow itself, our results indicate that the period of the dynamo cycle does not in fact depend on this depth. Rather, our simulations point out that the period of the dynamo cycle is more sensitive to changes on the speed of the counterflow than changes anywhere else in the transport circuit, as this is where the dynamo loop originates. An accurate determination of the average meridional flow speed over this loop closing at the SCZ base is very important in the context of the field transport timescales. As shown by the analysis of Yeates et al. (2008), the relative timescales of circulation and turbulent diffusion determine whether the dynamo operates in the advection or diffusion dominated regime—two regimes which have profoundly different flux transport dynamics and cycle memory (the latter may lead to predictability of future cycle amplitudes). Getting a firm handle on the average meridional flow speed is therefore very important and that is not currently achieved from the diverging helioseismic inversion results on the meridional flow.

This suggests that a concerted effort using different helioseismic techniques on data for the meridional flow over at least a complete solar cycle (over the same period of time) may be necessary to generate a more coherent picture of the observational constraint on this flow profile. It is important to note that even though we used time-averaged data, nothing prevents one from using the same methods to assimilate time-dependent helioseismic data at different phases of the solar circle, allowing us to study the impact of time varying velocity flows on solar cycle properties and their predictability.

On the encouraging side, our dynamo simulations show that it is relatively straightforward to use the available helioseismic data on the differential rotation (on which there is more consensus and agreement across various groups) within dynamo models. Also encouraging is the fact that the type of solar dynamo model presented here is able to handle the real helioseismic differential rotation profile and generate solarlike solutions. Moreover, as evident from our simulations, this dynamo model also generates plausible solarlike solutions over a wide range of meridional flow profiles, both deep and shallow, and with fast and slow peak flow speeds. This certainly bodes well for assimilating helioseismic data to construct better constrained solar dynamo models—building upon the techniques outlined here.

We are grateful to Irene González-Hernández and Rachel Howe at the National Solar Observatory for providing us with helioseismic data and useful counsel regarding its use. We also thank Jørgen Christensen-Dalsgaard of the Danish AsteroSeismology Center for data and discussions related to the solar model S. Finally, we want to thank Paul Charbonneau and an anonymous referee for useful comments and recommendations. The simulations presented here were performed at the computing facilities of the Harvard-Smithsonian Center for Astrophysics and this research was funded by NASA Living With a Star Grant NNX08AW53G to the Smithsonian Astrophysical Observatory. D.N. acknowledges support from the Department of Science and Technology of the Government of India through the Ramanujan Fellowship.

In order to use exponential propagation, we transform our system of partial differential equations (PDEs) into a system of coupled ordinary differential equations (ODEs) by discretizing the spatial operators using the following finite difference operators. For advective terms $\left(\frac{\partial A}{\partial t} = -v \frac{\partial A}{\partial x} + \chi (x) \right)$, where v is the velocity, we use a third-order upwind scheme:

$$\begin{flequation} v\frac{\partial A}{\partial x} = \left\lbrace \begin{array}{@{}ll@{}} \dsty\frac{v}{6\triangle x} (-2 A_{i-1} - 3 A_i \\ \ms \quad +\, 6 A_{i+1} - A_{i+2}) & \mbox{if $v < 0$}\\ \ms \dsty\frac{v}{6\triangle x} (2 A_{i+1} + 3 A_i \\ \ms \quad -\, 6 A_{i-1} + A_{i-2}) & \mbox{if $v \ge 0$}\\ \end{array} \right.\!\!\quad + \mathcal {O}(\triangle x^3). \end{flequation} \tag{ A1 }$$

For diffusive terms $\left(\frac{\partial A}{\partial t} = \eta \frac{\partial ^2 A}{\partial x^2} + \chi (x) \right)$, where η is the diffusion coefficient, we use a second-order space-centered scheme:

$$\begin{equation} \eta \frac{\partial ^2 A}{\partial x^2} = \frac{\eta }{(\triangle x)^2} (A_{i-1} - 2 A_i + A_{i+1}) + \mathcal {O}(\triangle x^2). \end{equation} \tag{ A2 }$$

For other first derivative terms $\left(\frac{\partial A}{\partial t} = \frac{\partial B}{\partial x} + \chi (x) \right)$ we use a second-order space-centered scheme:

$$\begin{equation} \frac{\partial B}{\partial x} = \frac{1}{2\triangle x} (A_{i-1} - A_{i+1}) + \mathcal {O}(\triangle x^2). \end{equation} \tag{ A3 }$$

Here, χ(x) corresponds to all the additional terms a PDE might have on the right-hand side besides the term under discussion and A_i = A(x₀ + i▵x), i = 1,  2,  ...,  N_x is our variable evaluated in a uniform grid of N_x elements separated by a distance ▵x.

A.1. Exponential Propagation

After discretization and inclusion of the boundary conditions we are left with an initial value problem of ODEs:

$$\begin{eqnarray} \frac{\partial \textbf {U}(t)}{\partial t} &=& F(\textbf {U}(t)), \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} \textbf {U}(t_0) &=& \textbf {U}_0,\quad\;\; \end{eqnarray} \tag{ A5 }$$

where U is the solution vector in $\mathbb {R}^N$. Provided that the Jacobian ∂F(U(t)) exists and is continuous in the interval [t₀, t₀ + Δt], we can linearize F(U(t₀ + Δt)) around the initial state to obtain

$$\begin{eqnarray} &&\frac{\partial \textbf {U}(t)}{\partial t} = F(\textbf {U}_0) + \partial F(\textbf {U}_0)(\textbf {U}(t_0+\Delta t) - \textbf {U}_0) + R(\textbf {U}(t_0+\Delta t)),\nonumber\\ \end{eqnarray} \tag{ A6 }$$

where R(U(t₀ + Δt) are the residual high-order terms. The solution to this equation can be written as

$$\begin{equation} \textbf {U}(t_0+\Delta t) = \textbf {U}_0 + \frac{e^{A\Delta t}-I}{A} + \mathcal {O}(\Delta t^2), \end{equation} \tag{ A7 }$$

where A = ∂F(U₀). Neglecting higher order terms leaves us with a scheme that is second-order accurate in time and is an exact solution of the linear case. However, there is a way of increasing the time accuracy of this method by following a generalization of Runge–Kutta methods for nonlinear time-advancement operators proposed by Rosenbrock (1963). The combination of exponential propagation with Runge–Kutta methods was first proposed by Hochbruck & Lubich (1997) and then generalized by Hochbruck et al. (1998). In this work, we use a fourth-order algorithm which goes to the following intermediary steps to advance the solution vector between time steps (Uⁿ) → U^{n + 1})):

$$\begin{eqnarray} k_1 &=& \Phi \left(\dsty\frac{1}{2}\Delta t A\right)F(\textbf {U}^n),\nonumber \\ k_2 &=& \Phi \left(\Delta t A\right)F(\textbf {U}^n),\nonumber \\ w_3 &=& \dsty\frac{3}{8}(k_1+k_2),\nonumber \\ u_3 &=& \textbf {U}^n + \Delta t w_3,\nonumber \\ d_3 &=& F(u_3) - F(\textbf {U}^n) - \Delta t A w_3,\nonumber \\ k_3 &=& \Phi \left(\dsty\frac{1}{2}\Delta t A\right)d_3,\nonumber \\ \textbf {U}^{n+1} &=& \textbf {U}^n + \Delta t \left(k_2 + \dsty\frac{16}{27}k_3 \right),\nonumber \\ A &=& F(\textbf {U}^n),\qquad \Phi \left(\Delta t A\right) = \dsty\frac{e^{A\Delta t}-I}{A}.\nonumber \end{eqnarray}$$

A.2. Krylov Approximation to the Exponential Operator

Without any further approximation, this method is very expensive computationally due to the need of continuously evaluating the matrix exponential. However, it is possible to make a good approximation by projecting the operator into a finite dimensional Krylov subspace

$$\begin{equation} S_{\rm Kr} = {\rm span}\lbrace \textbf {U}_0,A\textbf {U}_0,A^2\textbf {U}_0,...,A^{m-1}\textbf {U}_0\rbrace. \end{equation} \tag{ A8 }$$

In order to do this, we first compute an orthonormal basis for this subspace using the Arnoldi algorithm (Arnoldi 1951).

• 1.  
  v₁ = U₀/|U₀|².
• 2.  
  For j = 1, ..., m do

  • a)  
    for i = 1, ..., i compute h_i,j = v^T_i*Av_j,
  • b)  
    calculate $w=Av_j-\displaystyle \sum _{i=0}^j h_{i,j}v_i$,
  • c)  
    evaluate h(j + 1, j) = |w|²,
  • d)  
    if h_{j + 1,j} < stop, else compute the next basis vector v_{j + 1} = w/h_{j + 1,j},

where is the parameter that sets the error tolerance in this approximation. Once we have finished computing the algorithm, we have the relationship

$$\begin{equation} AV_m \approx V_mH, \end{equation} \tag{ A9 }$$

where V_m = [v₁, ..., v_m] and H is a matrix whose elements are h_i,j = v^T_i*Av_j. The validity of this approximation depends on the dimension of the Krylov subspace but numerical experiments have found that 15–30 Krylov vectors are usually enough (see Hochbruck & Lubich 1997; Tokman 2001). Since {v₁, ..., v_m} is an orthonormal basis of the Krylov subspace S_Kr then V^T_mV_M is a m × m identity matrix and V_mV^T_m is a projector from $\mathbb {R}^N$ onto S_Kr. Using this projector, we can find an approximation to any matrix vector multiplication by projecting them onto the Krylov subspace by calculating Ab ≈ V_mV^T_mAV_mV^T_mb. After using Equation (A9) this becomes

$$\begin{equation} Ab \approx V_mHV_m^Tb. \end{equation} \tag{ A10 }$$

In fact we can approximate the action of any operator Φ(A), that can be expanded on a Taylor series, on the vector U₀ using the Krylov subspace projection:

$$\begin{equation} \Phi (A)\textbf {U}_0 \approx V_m\Phi (H)V_m^T\textbf {U}_0 = |\textbf {U}_0|^2V_m\Phi (H)e_1, \end{equation} \tag{ A11 }$$

where e₁ = [1 0... 0] and we used v₁ = U₀/|U₀|² and V^T_mv₁ = e₁. As we can see in Equation (A11), the use of the Krylov approximation effectively reduces the size of the matrix operator; this makes the use of the exponential operator relatively inexpensive computationally.

These two algorithms, in combination with a robust error control and an adaptative time-step mechanism strategy, form the core of the SD-Exp4 integrator (for more details see Hochbruck & Lubich 1997; Hochbruck et al. 1998; Tokman 2001).

In order to verify the performance standards of the SD-Exp4 code, we ran case C' of the dynamo benchmark by Jouve et al. (2008). This case is similar in nature to the simulations done in this work with the following differences.

• 1.  
  The meridional flow profile has different radial and latitudinal dependence.
• 2.  
  The poloidal source term has no quenching term and a different radial and latitudinal dependence.
• 3.  
  The turbulent diffusivity profile consists of only one step and reaches a peak value of 10¹¹ cm² s⁻¹.
• 4.  
  The analytic differential rotation has no cos⁴(θ) dependence and uses a thinner tachocline.

In order to compare the performance of the different codes in the benchmark study, two quantities are used: C^crit_s = α^crit₀R_☉/η_cz which quantifies the minimum strength of the source relative to dissipation that yields stable oscillations and ω = 2πR²_☉/(Tη_cz) which quantifies the frequency of the magnetic cycle relative to the diffusion timescale. The dependence of this quantities is then plotted versus resolution for all the different codes (see Figure 11 of Jouve et al. 2008). In order to compare our code, we perform simulations with the same parameters and model ingredients as in case C' of Jouve et al. (2008) and plot this two quantities versus resolution. We also plot the location of the mean found by the other codes and a shaded area encompassing one standard deviation around the mean. As can be seen in Figure 14, we find lower values of C^crit_s than the values obtained in Jouve et al. (2008; 1.86 as opposed to an average 2.46) for oscillatory solutions; this may be due to a lower amount of numerical diffusion in our code. Moreover, we find values of ω in very good agreement with those found in Jouve et al. (2008; 540 as opposed to an average of 538.2). We also find that for our code this quantities are less sensitive to resolution when compared with other codes used in the benchmark. This is likely due to the fact that we try to avoid numerical differentiation whenever we can do it analytically.

Zoom In Zoom Out Reset image size

As a final remark, Table 4 of Jouve et al. (2008) specifies the different time steps that are used by the different codes, which range from values of 10⁻⁸ to 10⁻⁵ code time units. Thanks to the use of Krylov approximations, we are able to use time steps as large as 10⁻² in code time units while solving case C'. However, this does not translate to a performance improvement of several orders of magnitude (due to the added computations in the Krylov approximation). In order to quantify the relative performance improvement, we also made comparisons with the Surya code, which has been studied extensively in different contexts (e.g., Nandy & Choudhuri 2002; Chatterjee et al. 2004; Choudhuri et al. 2007), and is made available to the public on request. In comparison to the Surya code, we find that the SD-Exp4 code achieves a performance improvement that reduces runtime from a half to a tenth of the total runtime depending on the particularities of the simulation.