Modulated cycles in an illustrative solar dynamo model with competing αeffects
School of Mathematics and Statistics, Newcastle University, Newcastle Upon Tyne, NE1 7RU, UK
email: l.cole@ncl.ac.uk; paul.bushby@ncl.ac.uk
Received: 19 December 2013
Accepted: 17 February 2014
Context. The largescale magnetic field in the Sun varies with a period of approximately 22 years, although the amplitude of the cycle is subject to longterm modulation with recurrent phases of significantly reduced magnetic activity. It is believed that a hydromagnetic dynamo is responsible for producing this largescale field, although this dynamo process is not well understood.
Aims. Within the framework of meanfield dynamo theory, our aim is to investigate how competing mechanisms for poloidal field regeneration (namely a timedelayed BabcockLeighton surface αeffect and an interfacetype αeffect), can lead to the modulation of magnetic activity in a deepseated solar dynamo model.
Methods. We solve the standard αΩ dynamo equations in one spatial dimension, including source terms corresponding to both of the competing αeffects in the evolution equation for the poloidal field. This system is solved using two different methods. In addition to solving the onedimensional partial differential equations directly, using numerical techniques, we also use a local approximation to reduce the governing equations to a set of coupled ordinary differential equations (ODEs), which are studied using a combination of analytical and numerical methods.
Results. In the ODE model, it is straightforward to find parameters such that a series of bifurcations can be identified as the time delay is increased, with the dynamo transitioning from periodic states to chaotic states via multiply periodic solutions. Similar transitions can be observed in the full model, with the chaotically modulated solutions exhibiting solarlike behaviour.
Conclusions. Competing αeffects could explain the observed modulation in the solar cycle.
Key words: dynamo / magnetohydrodynamics (MHD) / Sun: activity / Sun: interior / Sun: magnetic fields
© ESO, 2014
1. Introduction
At the solar photosphere, bipolar active regions are formed when loops of magnetic flux rise to the surface from the base of the convection zone due to the action of magnetic buoyancy (Parker 1955b). This implies that the properties of sunspotbearing active regions can be used to deduce some of the features of the underlying largescale magnetic field. It is well known (see, for example, Stix 2002; Charbonneau 2005; Jones et al. 2010) that zones of active region emergence follow a cyclic pattern with a period of approximately 11 years. At the beginning of each cycle, sunspots tend to be found at midlatitudes, with zones of emergence drifting towards the equator as the cycle progresses. The underlying largescale (predominantly azimuthal) magnetic field changes sign at the end of each cycle, giving a full magnetic period of approximately 22 years. However, the solar cycle is not strictly periodic. In particular, the peak amplitude (measured, for example, by the sunspot coverage) varies from one cycle to the next. Although this modulation does not usually disrupt the cycle, more extreme episodes of modulation have been recorded. For example, during a period known as the Maunder Minimum, very few sunspots were observed between approximately 1650 and 1720 (Eddy 1976; Ribes & NesmeRibes 1993). However, sunspot records are not the only indicators of modulation. Due to the fact that the Sun’s strong magnetic field protects the Earth from cosmic rays, the abundance of certain isotopes in the Earth’s atmosphere is known to be anticorrelated with the solar cycle. Therefore, by analysing Beryllium10 deposits in ice cores (see, for example, Delaygue & Bard 2011) and Carbon14 levels in tree rings (see, for example, Muscheler et al. 2007) it is possible to deduce the history of the solar cycle. Such studies have indicated that cyclic activity did persist throughout the Maunder Minimum, but at a significantly reduced level (Beer et al. 1998). Furthermore, it is clear that the Maunder Minimum is not exceptional – the solar cycle has often been interrupted by recurrent “Grand Minimum” phases of significantly reduced magnetic activity.
It is believed that the largescale magnetic field in the solar interior is generated and maintained by a hydromagnetic dynamo. From a conceptual point of view, the largescale field can usefully be decomposed into its toroidal (azimuthal) and poloidal (meridional) components – a working dynamo requires mechanisms that allow the poloidal field to be regenerated from toroidal field and vice versa. It is widely accepted that differential rotation (usually referred to as the Ωeffect in dynamo theory) is responsible for the generation of toroidal field from poloidal field. Surface observations indicate that equatorial regions rotate more rapidly than the poles, and helioseismological studies (see, for example, Schou et al. 1998) have shown that this rotation profile persists, approximately independently of radius, throughout most of the convection zone. At the base of the convection zone, a region of strong shear (the tachocline) couples the radiative zone, which rotates almost rigidly, to the differentiallyrotating convective envelope. In most solar dynamo models, it is assumed that a significant fraction of the toroidal field is generated in the vicinity of the tachocline (where the Ωeffect should be very efficient due to the presence of strong differential rotation).
Although the Ωeffect is well understood, the reverse process that generates poloidal field from toroidal field is a topic of some debate. In classical interface dynamo models (see, for example, Parker 1993; Charbonneau & MacGregor 1996) the poloidal field is regenerated at the base of the convection zone by the action of cyclonic convection upon toroidal magnetic field lines (Parker 1955a). This process is usually referred to as the αeffect. Strong toroidal fields will tend to inhibit (or quench) the operation of the αeffect, so interface dynamo models are usually constructed in such a way that the αeffect is restricted to the region just above the base of the convection zone, whilst the Ωeffect operates just below the interface. The two layers are coupled by the effects of magnetic diffusion, as well as magnetic buoyancy and turbulent pumping (see, for example, Tobias et al. 2001). Even with strong αquenching, it has been shown that an interface dynamo of this type can operate efficiently (Charbonneau & MacGregor 1996). In BabcockLeighton dynamo models (Babcock 1961; Leighton 1964), the poloidal field is regenerated at the solar surface through the decay of active regions (which tend to emerge with a systematic tilt with respect to the azimuthal direction). This surface αeffect can only contribute to the dynamo if there is some mechanism that is capable of transporting the resultant poloidal field to the tachocline. This could be achieved by diffusion or by pumping, but meridional flows also could play an important role in this respect. A polewards meridional flow is observed at the solar surface (see, for example, Hathaway & Rightmire 2010) and, by mass conservation arguments, there must be a returning circulatory flow somewhere within the solar interior. A singlecell meridional circulation, with an equatorial flow at the base of the convection zone would couple the surface layers to the tachocline in an effective way, thus completing the dynamo loop.
A complete model of the solar dynamo must be able to explain the observed modulation as well as the 22year magnetic cycle. It has been shown that it is possible to induce modulation by introducing stochastic effects into BabcockLeighton models (Charbonneau & Dikpati 2000; Bushby & Tobias 2007), as well as into models of interface type (Ossendrijver 2000). However, fully deterministic models (with no random elements) can also produce modulated dynamo waves. Weiss et al. (1984) and Jones et al. (1985) considered a simple system in which the dynamo was modelled using a set of coupled ordinary differential equations, which included the nonlinear interactions between the magnetic field and the flow. They found that it was possible to generate quasiperiodic and chaoticallymodulated solutions in addition to standard periodic dynamo waves. More recent studies have shown that the full meanfield equations also exhibit significant modulation when dynamical nonlinearities are included in the governing equations (see, for example, Tobias 1996; Brooke et al. 2002; Bushby 2006). An alternative approach was used by Yoshimura (1978) who demonstrated that modulation can arise if explicit time delays are built into the nonlinear terms in a simple system of model dynamo equations. A more sophisticated model was considered by Jouve et al. (2010) who investigated the effects of magnetic buoyancyinduced time delays in the context of a twodimensional BabcockLeighton dynamo. By introducing time delays into the surface αeffect term, they were able to demonstrate the existence of modulated cycles. They then went on to consider a simpler onedimensional dynamo system in which the surface αeffect term was represented by the inclusion of a timedelayed toroidal field (with a parameterised time delay that was dependent upon the magnetic field strength). They were able to demonstrate the existence of a sequence of bifurcations from periodic to chaotically modulated solutions as the time delay parameter was increased.
The aim of this work is to investigate the competition between a deepseated (interface) αeffect and a surface αeffect. Building on the approach described by Jouve et al. (2010), who did not include a deepseated αeffect, the influence of the surface αeffect will be modelled using a timedelayed toroidal field. The use of a time delay is natural in this context: even if flux tubes rise rapidly to the surface, the time taken for the resultant poloidal field to be transported back to the tachocline will, in general, be nonnegligible compared to the period of oscillation of the dynamo. Previous studies have investigated systems with competing αeffects (see, for example, Dikpati & Gilman 2001; Mason et al. 2002; Mann & Proctor 2009), but we believe that this is the first study to consider the effects of explicit time delays in a model of this type. The paper is structured as follows: In Sect. 2, we describe the full model and an idealised system of equations that can derived from it (based upon a local analysis). This is followed in Sect. 3 by an analysis of the stability of the idealised model and then in Sect. 4 by the corresponding numerical results. In Sect. 5, we describe some numerical calculations which demonstrate the existence of modulated solutions in the full onedimensional model. Finally, in Sect. 6, we present our conclusions and discuss the relevance of our results to the solar dynamo.
2. Model setup
Following a similar approach to that adopted by Jouve et al. (2010), we consider a simple, illustrative model of the solar dynamo. This model is based upon the standard meanfield dynamo equation (see, e.g., Moffatt 1978), (1)where U is the largescale velocity field, α represents the standard meanfield αeffect, η_{T} is the turbulent magnetic diffusivity (which we shall assume to be constant), whilst the mean magnetic field, B, satisfies ∇·B = 0. Instead of solving this equation in spherical geometry, we consider the simpler problem of dynamo action in a flat Cartesian domain, with the axes oriented so that the yaxis would correspond to the azimuthal direction on a spherical surface. We can then look for dynamo solutions that depend only on a single spatial variable x (which can be regarded as being analogous to the colatitude) and time t. The solenoidal constraint upon B can then be satisfied by writing the magnetic field in the following form: (2)where B(x,t) is the toroidal field component, whilst A(x,t) corresponds to the poloidal potential.
Our model is based on the assumption that the solar dynamo is operating primarily in the region around the base of the convection zone. For simplicity, we assume that α, which represents a deepseated αeffect, is constant, i.e. α = α_{0}. Furthermore, we adopt a fixed velocity profile of the form , where v_{0} and Ω_{0} are both assumed to be constant in this illustrative model. This velocity field gives a constant meridional flow and a differential rotation profile that is independent of x. We also make the well known αΩ approximation, which assumes that differential rotation is the dominant mechanism for toroidal field regeneration in this region. Following Jouve et al. (2010), we also introduce a delayed toroidal field, Q(x,t), which lags behind the normal toroidal field with a time delay denoted by τ. However, unlike Jouve et al. (2010), who considered a time delay that was dependent upon the toroidal magnetic field strength, we assume τ to be constant throughout this study. The delayed toroidal field is coupled to the other equations via the inclusion of an additional poloidal source term, SQ(x,t), where S is a constant. This source term can be regarded as being the contribution to the local poloidal field from the nonlocal surface αeffect (which must, therefore, depend upon the strength of the toroidal field at earlier times). Finally, we introduce parameterised quenching nonlinearities into both of the αeffect terms in the poloidal field equation.
Having made these assumptions, we can now write down the three scalar partial differential equations for A(x,t), B(x,t) and Q(x,t): where λ is a constant that determines the strength of the nonlinear quenching. In order to reduce the number of parameters that control the system, the variables can be rescaled as follows: where L is a characteristic lengthscale and B_{0} is a representative value of the magnetic field strength (which may be chosen so that the constant coefficient in the quenching terms equals unity in these scaled variables). On dropping the primes, we obtain the following set of partial differential equations (PDEs): Thus the only parameters to control the system are the Reynolds number, Re = v_{0}L/η_{T}, which measures the strength of the meridional flow, and the dynamo number, , which indicates the strength of the dynamo sources relative to magnetic dissipation.
This system can be further simplified by carrying out a local analysis. Because we have the freedom to choose a convenient characteristic lengthscale, local wavelike solutions can be assumed to have a unit wavenumber without any loss of generality. We therefore seek solutions of the form , and , where , and are complex functions of time only. Dropping the tildes, the governing equations for these quantities become: Following the methods used in Jones et al. (1985) and Jouve et al. (2010) it is possible to reduce the order of this system by using the following representation: where ρ and θ are real quantities and y and z are complex numbers. Upon substituting these expressions into the governing Eqs. (9)–(11), the following set of 5 real ODEs is obtained: where y_{1} and y_{2} represent the real and imaginary parts of y respectively and z_{1} and z_{2} represent the corresponding real and imaginary parts of z. This fifthorder system is a useful alternative representation of the local model.
3. Dynamo transitions
In this section, we focus upon the local model that is described by Eqs. (9)–(11). To further our understanding of this system, we have carried out a series of calculations to determine the critical value of the dynamo number as the parameters S and τ are varied (at fixed Re). It is also possible to identify the value of τ that leads to quasiperiodic solutions analytically by studying the stability of the periodic solution.
3.1. Linear theory
The critical dynamo numbers can be calculated by linearising the governing Eqs. (9)–(11) and writing A, B and Q in the following form: and . The following characteristic equation is generated: (17)Setting the real part of the growth rate to be zero and solving the characteristic equation for the imaginary part of σ will determine the critical value of the dynamo number, D_{c}, at which the trivial (nonmagnetic) solution loses stability to oscillatory dynamo waves. Letting σ = iω, we obtain the following: (18)In the case of S = 0 the system corresponds to a standard αΩ dynamo. It is then straightforward to show that the critical dynamo number is 2, regardless of the magnitudes of τ or Re. For S ≠ 0, this equation can be solved numerically using a NewtonRaphson algorithm.
Fig. 1 Critical dynamo number D_{c} as a function of S, for values of τ less than 1 with Re = 10. Here, the solid line corresponds to τ = 0.1, the dotted line represents τ = 0.3 and the dashdotted line shows τ = 0.5. 

Open with DEXTER 
As a specific example, Fig. 1 shows how the critical dynamo number changes as S is varied between −50 and 50 with values of τ less than 1 and for fixed Re = 10. For these values of τ, D_{c} < 2 for most values of S. In these regions of parameter space, it is easier to excite a dynamo than it would be in the corresponding αΩ system, so we can see that the nonlocal αeffect is enhancing the dynamo for these parameter values. However, for larger values of τ (where there is a significant timelag between B and Q) there is a finite range of values of S in which D_{c} > 2, with a local maximum in D_{c} occurring somewhere in this range. In this case, the two competing αeffects appear to be impeding each other, thus making it more difficult to excite a dynamo. So it is clear that, even in linear theory, the interaction between competing αeffects is nontrivial. Additional calculations have been made to determine the critical dynamo number for different values of Re and a broader range of values for τ, but this case (which is of greatest relevance to the present study) is fairly representative in terms of the behaviour that is exhibited as S is varied.
3.2. The transition from periodic to quasiperiodic behaviour
Restoring the nonlinear terms to the governing equations, finite amplitude oscillations can be found when the dynamo number exceeds D_{c}. The stability of the periodic solution can be analysed by expressing the magnetic fields in the following form: where A_{0}, B_{0} and Q_{0} are the complex wave amplitudes, ω is the frequency (with the ω subscript denoting the periodic state). The substitution of these expressions into the governing Eqs. (9)–(11) produces the following simultaneous equations: which can be solved using standard methods. Once the amplitude and frequency of the periodic solution have been determined, it is possible to perturb this solution to study its stability. Following the general method described by Jouve et al. (2010), this can be achieved by setting: where α_{1},α_{2},β_{1},β_{2},γ_{1} and γ_{2} are the coefficients of the perturbed fields, p is the complex growth rate of the perturbation and the symbol ⋆ represents the complex conjugate. Substituting these expressions into the governing Eqs. (9)–(11) results in a system of 6 coupled equations that relates the coefficients of the perturbed fields to the growth rate p for a given set of parameters. After solving this system to find the growth rate of the perturbation, it is then possible to determine the stability of the periodic solutions.
Critical values of τ for D = 1000.
Critical values of τ for D = −1000.
Values of  B_{0} ,  Q_{0}  and ω from both analytical calculations and the numerical simulations.
Tables 1 and 2 illustrate some of the results from this stability analysis. These tables show the parametric dependence of the critical value of τ for the transition from periodic to quasiperiodic solutions. The results in Table 1 correspond to D = 1000, with 0 ≤ Re ≤ 50 and 0 ≤ S ≤ 70. In Table 2, we have used the same values of Re, but D = −1000, whilst 0 ≥ S ≥ − 70. An entry of (...) in either table indicates that no transition exists. Unsurprisingly, no modulation is found for S = 0. In this case the delayed toroidal field Q decouples from the system and we have a standard αΩ dynamo model. More unexpectedly, these results suggest that Re ≠ 0 is a necessary condition for modulation in this system. So the meridional flow seems to play a crucial role in driving the modulation, perhaps by introducing an additional (advective) timescale into the problem. If it is simply the presence of an additional timescale that is the key ingredient here, then it may still be possible to drive modulation in the absence of a flow if some other physical process (such as turbulent pumping) was included in the model. However, it is beyond the scope of this paper to investigate whether or not this is indeed the case. In the case of positive D, no modulation was found for negative S, whilst the same is true for positive values of S in the negative D case. Given the idealised nature of this local model, we should probably not read too much into this result, but (if nothing else) this again illustrates that competing αeffects interact in a rather nontrivial way in this system. Where modulation does occur, some trends can be identified. For example, for fixed Re in the D = 1000 case, the critical value of τ decreases with increasing S (whereas it increases with increasing  S  in the D = −1000 case). At fixed S, the critical value of τ tends to decrease with increasing values of Re, although this trend appears to reverse at low S and high Re in the D = 1000 case. We have no definitive physical explanation for this behaviour but can speculate that this is somehow related to the nonmonotonicity that was observed in the D_{c} calculations in the previous subsection.
4. Numerical simulations of the local model
In this section, we apply a numerical approach to the local model that is described by Eqs. (9)–(11). Decomposing the system into its real and imaginary parts, we use a fourthorder RungeKutta scheme in Fortran to timestep the governing equations.
4.1. Validation of numerical calculations
To validate the code, it is possible to check that the results agree with the critical dynamo numbers that can be obtained from linear theory (see Sect. 3.1). Fixing the values of Re = 10, S = 20 and τ = 10, gives a prediction of D_{c} = 2.499. This is consistent with the numerics: we find decaying oscillations for D = 2.4, whilst D = 2.6 gives a stable periodic solution. We can also compare the amplitude and frequency of the periodic solutions with the corresponding analytical predictions, using a Fourier transform to determine the frequency of oscillation in the numerical case. Table 3 shows the results of such a comparison, for variable Re, using D = 1000, S = 10 and τ = 0.1 (which includes the “reference case” below). All results are accurate to within 1% which clearly validates both the numerical scheme and the analytical calculations. Finally, fixing Re = 10, D = 1000 and S = 10, we find that the solution exhibits a transition from periodic to quasiperiodic dynamo waves at τ = 0.347, which is compares very favourably to the analytic value of τ = 0.348 (see Table 1).
Fig. 2 Reference case for τ = 0.2: this shows the toroidal field B (solid line) and the delayed toroidal field Q (dashed line) as a function of time (which is expressed in dimensionless units). 

Open with DEXTER 
4.2. Results
Initially, the parameters are chosen such that D = 1000, Re = 10, S = 10 and τ is varied (we will refer to this as the “reference case”). Figure 2 shows that a periodic solution can be found provided that τ is sufficiently small. Both B(t) and Q(t) oscillate with constant amplitude although Q(t) has a smaller amplitude of oscillation and, as expected, lags behind B(t). The effects of increasing the value of τ are shown in Fig. 3. As τ is increased through the threshold value of τ = 0.347, the lag between B(t) and Q(t) increases to such an extent that we see a transition to a quasiperiodic state. Further increases in τ lead to further transitions, from multiply periodic to chaotically modulated states. Figure 4 shows the time evolution of the toroidal field energy B^{2} for τ = 0.86, at which point the solution is chaotically modulated. It is clear that there are several phases of significantly reduced magnetic activity, and it is tempting to compare these to grand minima. The extent to which this behaviour is “solarlike” is a matter of some debate – this is, after all, a highly idealised model. Nevertheless, it is encouraging that this simple model, with competing αeffects, is capable of producing highly modulated dynamo waves when the time delay is large.
Fig. 3 Reference case for τ = 0.35 (top), τ = 0.61 (middle) and τ = 0.86 (bottom). The plots on the left show the timedependence of the toroidal field, whilst the plots on the right show the phase portraits of the amplitudes of B(t) against Q(t) (as derived from the 5thorder system). 

Open with DEXTER 
Fig. 4 Reference solution for τ = 0.86. A plot of B^{2} against time. 

Open with DEXTER 
As indicated by the results in Table 1, the analysis of the stability of the periodic solution indicates that it is not possible to find a transition to a quasiperiodic solution for negative values of S, when D is positive. This tendency for the periodic state to be stable (for S ≤ 0) regardless of the value of τ has been confirmed numerically. However, for positive values of S it always appears to be possible to find a transition to quasiperiodic solutions, provided that the Reynolds number is nonzero, and these transitions are consistent with those predicted in Table 1. Once quasiperiodic solutions have been found it is usually possible to find chaoticallymodulated states for sufficiently large values of the time delay. For negative values of the dynamo number, the results are again consistent with those predicted analytically. No modulation is found for positive S or for Re = 0. For negative S and positive Re, it is possible to find quasiperiodic and chaotically modulated solutions as the time delay is increased. One such solution is illustrated in Fig. 5. As in Fig. 4, it should again be noted that the chaotically modulated solution that is illustrated in the lower part of Fig. 5 is characterised by phases of significantly reduced magnetic activity.
Fig. 5 Effects of increasing τ for D = −1000, Re = 10, S = −10. The upper plot shows the time evolution of the toroidal field B (solid line) and the delayed toroidal field Q (dashed line) for a periodic solution at τ = 0.2. The lower plot shows the toroidal field energy, B^{2}, as a function of time for τ = 1.08. 

Open with DEXTER 
5. Solving the PDE system
Although the results from the local model are promising, it is important to verify that they are not crucially dependent upon the simplifying assumptions that have been made when deriving the model. In this section, we return to the original model of partial differential equations, as defined by Eqs. (3)–(5). In dimensionless units, we assume that 0 ≤ x ≤ π/2 (recalling that we interpret x as being analogous to the colatitude on a spherical surface), imposing the boundary conditions that A = B = Q = 0 at x = 0 (the “North pole”) and B = Q = ∂A/∂x = 0 at x = π/2 (the “Equator”). These boundary conditions correspond to the assumption that the global magnetic field has dipolar symmetry. Having neglected the effects of curvature, and having assumed constant α_{0}, v_{0} and Ω_{0}, we should stress again that this should still be regarded as an illustrative model. Nevertheless, it contains the key physical ingredient of two competing αeffects with a surface αeffect contribution that depends upon a timedelayed toroidal field. In order to obtain dynamo waves that propagate towards the Equator, we focus primarily upon the D < 0 parameter regime (which would correspond to a negative deepseated αeffect in the northern hemisphere). We solve the governing equations numerically, approximating derivatives using secondorder finite differences. A 4thorder RungeKutta scheme is again used to timestep the governing equations.
Given that we are investigating the D < 0 regime, the local model suggests that we should be able to find modulation for negative values of S. However, in this region of parameter space there is an overwhelming tendency for steady modes to be preferred at onset (recall that wavelike solutions were assumed when the local model was derived). It is well known that steady and oscillatory modes can bifurcate from the trivial state at similar values of D in global αΩ dynamos (see, for example, Jennings & Weiss 1991), so this behaviour is not entirely unsurprising. However, it is almost certainly rather model specific – experimentation with the inclusion of different nonlinear quenching mechanisms suggests that it is possible to obtain oscillatory solutions in these parameter regimes. Furthermore, oscillatory solutions can be found for positive dynamo numbers and therefore, despite some differences, the results from the local model should not be discarded.
In fact, in the case of this global model, interesting solutions can be found for negative values of D and positive values of S. This is illustrated by Fig. 6 which shows solutions for D = −6000, S = 1 and Re = 10. A periodic solution can be found at τ = 0.01. This is characterised by an oscillatory magnetic field which propagates towards the Equator (note that these contour plots have been plotted as a function of latitude and time, for ease of comparison with observations). Increasing the timedelay leads to a transition to a quasiperiodic solution. Further increases in τ eventually lead to chaotically modulated oscillations (as illustrated in Fig. 7). This solution is rather “solarlike” in many respects, with the dynamo confined to low latitudes, and with strong variations in the amplitudes of successive cycles. Furthermore, the modulation is characterised by periods of reduced magnetic activity. So although the modulation due to these competing αeffects was not in the expected parameter regime, it is clearly a robust feature of this system.
Fig. 6 Dynamo solutions from the full PDE system (D = −6000, S = 1 and Re = 10). Top: contours of toroidal field as a function of latitude and time (a latitude of 90° corresponds to the pole, 0° to the equator) for τ = 0.01. Middle: as above, but for τ = 0.05. Bottom: a plot of the energy in the toroidal field as a function of time for the quasiperiodic solution that is obtained for τ = 0.05. 

Open with DEXTER 
Fig. 7 Chaoticallymodulated solution from the full PDE system (D = −6000, S = 1, Re = 10 and τ = 0.09). Top: contours of toroidal field as a function of latitude and time. Bottom: a plot of the energy in the toroidal field as a function of time. 

Open with DEXTER 
6. Conclusions
In this paper, we have investigated the properties of an illustrative meanfield dynamo model which includes two competing αeffects. The first of these is the standard deepseated αeffect, the second is due to a surface αeffect (of BabcockLeighton type). Following the approach described by Jouve et al. (2010), who did not consider competing αeffects, the contribution from the surface αeffect was modelled by assuming that it depends upon a timedelayed toroidal field (with a constant parameterised time delay τ). Two different approaches were applied to this model. Initially, a local approximation was made to reduce the governing equations to a system of coupled ordinary differential equations. A linearised version of these equations was used to determine the dependence of the critical dynamo number upon S (the magnitude of the surface αeffect) and τ. Generally, the larger the magnitude of S, the easier it becomes to excite the dynamo. However, there are some regions of parameter space in which the two competing αeffects appear to impede each other, thus inhibiting the dynamo. Moving beyond linear theory, it was found that there are significant regions of parameter space in which the periodic solution becomes unstable with increasing τ, leading to quasiperiodicity. This was verified numerically, where further increases in the time delay were shown to produce chaotically modulated states with phases of significantly reduced activity. The full PDE model was then investigated. Although modulation was found, this occurs in a different parameter regime to that predicted by the ODE model. This discrepancy could be model specific, although we expected to see some differences between the two models due to the fact that significant simplifications were made when deriving the set of coupled ODEs. Nevertheless, it was possible to find chaotically modulated solutions in the PDE model, and these solutions exhibit certain features that are (at least qualitatively) “solarlike”.
There are many possible areas of future work. In particular, more could be done to explore the robustness of the PDE model to variations in the boundary conditions and the nonlinear quenching mechanisms. As has already been mentioned, preliminary calculations suggest that the adoption of different nonlinearities may make a very significant difference to the behaviour of the model. It may also be possible to improve the existing model by refining the way in which the time delay is implemented – the current approach is simple and effective, but is derived by truncating a Taylor series expansion at lowest order. Retaining higher order terms may make a difference to the behaviour of the model. Moving beyond the onedimensional Cartesian system, it would be natural to explore a twodimensional version of this model in (axisymmetric) spherical geometry. This would open up the possibility of including a more realistic flow geometry (in both the meridional and azimuthal directions) as well as spatially dependent meanfield coefficients. Although this would still be within the framework of meanfield theory, a more realistic model would enable more detailed comparisons to be made between our results and the solar dynamo.
References
 Babcock, H. W. 1961, ApJ, 133, 572 [NASA ADS] [CrossRef] (In the text)
 Beer, J., Tobias, S., & Weiss, N. 1998, Sol. Phys., 181, 237 [NASA ADS] [CrossRef] (In the text)
 Brooke, J., Moss, D., & Phillips, A. 2002, A&A, 395, 1013 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Bushby, P. J. 2006, MNRAS, 371, 772 [NASA ADS] [CrossRef] (In the text)
 Bushby, P. J., & Tobias, S. M. 2007, ApJ, 661, 1289 [NASA ADS] [CrossRef] (In the text)
 Charbonneau, P. 2005, Liv. Rev. Sol. Phys., 2, 2 (In the text)
 Charbonneau, P., & Dikpati, M. 2000, ApJ, 543, 1027 [NASA ADS] [CrossRef] (In the text)
 Charbonneau, P., & MacGregor, K. B. 1996, ApJ, 473, L59 [NASA ADS] [CrossRef] (In the text)
 Delaygue, G., & Bard, E. 2011, Clim. Dyn., 36, 2201 [NASA ADS] [CrossRef] (In the text)
 Dikpati, M., & Gilman, P. A. 2001, ApJ, 559, 428 [NASA ADS] [CrossRef] (In the text)
 Eddy, J. A. 1976, Science, 192, 1189 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Hathaway, D. H., & Rightmire, L. 2010, Science, 327, 1350 [NASA ADS] [CrossRef] (In the text)
 Jennings, R. L., & Weiss, N. O. 1991, MNRAS, 252, 249 [NASA ADS] [CrossRef] (In the text)
 Jones, C. A., Weiss, N. O., & Cattaneo, F. 1985, Physica D Nonlinear Phenomena, 14, 161 [NASA ADS] [CrossRef] (In the text)
 Jones, C. A., Thompson, M. J., & Tobias, S. M. 2010, Space Sci. Rev., 152, 591 [NASA ADS] [CrossRef] (In the text)
 Jouve, L., Proctor, M. R. E., & Lesur, G. 2010, A&A, 519, A68 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Leighton, R. B. 1964, ApJ, 140, 1547 [NASA ADS] [CrossRef] (In the text)
 Mann, P. D., & Proctor, M. R. E. 2009, MNRAS, 399, L99 [NASA ADS] [CrossRef] (In the text)
 Mason, J., Hughes, D. W., & Tobias, S. M. 2002, ApJ, 580, L89 [NASA ADS] [CrossRef] (In the text)
 Moffatt, H. K. 1978, Magnetic field generation in electrically conducting fluids (Cambridge: Cambridge University Press) (In the text)
 Muscheler, R., Joos, F., Beer, J., et al. 2007, Quat. Sci. Rev., 26, 82 [NASA ADS] [CrossRef] (In the text)
 Ossendrijver, M. A. J. H. 2000, A&A, 359, 364 [NASA ADS] (In the text)
 Parker, E. N. 1955a, ApJ, 122, 293 [NASA ADS] [CrossRef] (In the text)
 Parker, E. N. 1955b, ApJ, 121, 491 [NASA ADS] [CrossRef] (In the text)
 Parker, E. N. 1993, ApJ, 408, 707 [NASA ADS] [CrossRef] (In the text)
 Ribes, J. C., & NesmeRibes, E. 1993, A&A, 276, 549 [NASA ADS] (In the text)
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [NASA ADS] [CrossRef] (In the text)
 Stix, M. 2002, The sun: an introduction (Berlin: SpringerVerlag) (In the text)
 Tobias, S. M. 1996, A&A, 307, L21 [NASA ADS] (In the text)
 Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 2001, ApJ, 549, 1183 [NASA ADS] [CrossRef] (In the text)
 Weiss, N. O., Cattaneo, F., & Jones, C. A. 1984, Geophys. Astrophys. Fluid Dynamics, 30, 305 [NASA ADS] [CrossRef] (In the text)
 Yoshimura, H. 1978, ApJ, 226, 706 [NASA ADS] [CrossRef] (In the text)
All Tables
Values of  B_{0} ,  Q_{0}  and ω from both analytical calculations and the numerical simulations.
All Figures
Fig. 1 Critical dynamo number D_{c} as a function of S, for values of τ less than 1 with Re = 10. Here, the solid line corresponds to τ = 0.1, the dotted line represents τ = 0.3 and the dashdotted line shows τ = 0.5. 

Open with DEXTER  
In the text 
Fig. 2 Reference case for τ = 0.2: this shows the toroidal field B (solid line) and the delayed toroidal field Q (dashed line) as a function of time (which is expressed in dimensionless units). 

Open with DEXTER  
In the text 
Fig. 3 Reference case for τ = 0.35 (top), τ = 0.61 (middle) and τ = 0.86 (bottom). The plots on the left show the timedependence of the toroidal field, whilst the plots on the right show the phase portraits of the amplitudes of B(t) against Q(t) (as derived from the 5thorder system). 

Open with DEXTER  
In the text 
Fig. 4 Reference solution for τ = 0.86. A plot of B^{2} against time. 

Open with DEXTER  
In the text 
Fig. 5 Effects of increasing τ for D = −1000, Re = 10, S = −10. The upper plot shows the time evolution of the toroidal field B (solid line) and the delayed toroidal field Q (dashed line) for a periodic solution at τ = 0.2. The lower plot shows the toroidal field energy, B^{2}, as a function of time for τ = 1.08. 

Open with DEXTER  
In the text 
Fig. 6 Dynamo solutions from the full PDE system (D = −6000, S = 1 and Re = 10). Top: contours of toroidal field as a function of latitude and time (a latitude of 90° corresponds to the pole, 0° to the equator) for τ = 0.01. Middle: as above, but for τ = 0.05. Bottom: a plot of the energy in the toroidal field as a function of time for the quasiperiodic solution that is obtained for τ = 0.05. 

Open with DEXTER  
In the text 
Fig. 7 Chaoticallymodulated solution from the full PDE system (D = −6000, S = 1, Re = 10 and τ = 0.09). Top: contours of toroidal field as a function of latitude and time. Bottom: a plot of the energy in the toroidal field as a function of time. 

Open with DEXTER  
In the text 