© ESO, 2011

1. Introduction

The observed existence of large-scale astrophysical magnetic fields, for example galactic or solar fields, is usually explained by self-excited dynamo action within electrically conducting fluids or plasmas. However, this mechanism of field amplification continues to be a matter of debate as the existing theory encounters problems when extrapolated to the large magnetic Reynolds numbers of astrophysics. Nonetheless, large-scale astrophysical fields are believed to be predominately generated by so-called αΩ dynamos, in which most of the field amplification occurs through the shearing of field lines by ubiquitous differential rotation, a process known as the Ω effect (Steenbeck & Krause 1969). For example, many models of the solar dynamo invoke the strong shear found in the tachocline at the base of the convection zone (see, e.g., Charbonneau 2010). Shear alone cannot drive dynamo action however, and the α effect, caused by helical motions, provides the necessary twist of the sheared field to complete the magnetic field amplification cycle. In the Sun, an α effect is provided via kinetic helicity due to the interaction of stratified convection and solar rotation.

The α effect can drive a dynamo by itself, which is then of the so-called α² type. These dynamos are of great theoretical interest due to their simplicity, but are expected to be outperformed by αΩ dynamos in the wild. Strictly speaking, αΩ dynamos should be referred to as α²Ω dynamos as the α² process of course continues to occur in reality, even in the presence of the Ω effect. However, in the mean-field approach one sometimes makes the so-called “αΩ” approximation by neglecting the production of toroidal field via the α effect entirely in favor of the Ω effect. In such models, the nonlinear competition between different αΩ modes has been thoroughly studied by monitoring, for example, rapid changes of Lyapunov exponents in the bifurcation structure (Covas et al. 1997).

In the present paper we consider numerical solutions of the compressible MHD equations in three dimensions with turbulent helical flows where, of course, the αΩ approximation is not applicable. Nevertheless we will refer to αΩ and α²Ω regimes when shear is dominant or comparable with amplification by the helical turbulence, respectively.

Very often, a linear stability analysis of a given setup reveals that several different dynamo modes are expected to be excited at the same time. While during the linear stage the relative strength of these modes is determined by the initial conditions, the mode or mixture of modes of the final saturated state is decided by the nonlinear interactions between the modes in their backreaction on the flow. The naive guess that the final state should always be characterized by the mode with the highest growth rate, has turned out not to be valid in general. This was first shown in Brandenburg et al. (1989), who found that, for small enough dynamo numbers, the saturated state of α² and αΩ dynamos is instead determined by the solution with the smallest marginal dynamo number. For larger dynamo numbers, however, axisymmetric dipolar and quadrupolar modes have asymptotically identical growth rates, and in the nonlinear regime there can be several stable solutions, including some with mixed parity, that no longer bifurcate from the trivial one. Moreover, certain axisymmetric solutions turned out to be unstable to non-axisymmetric perturbations and evolved eventually toward another axisymmetric solution (Rädler et al. 1990).

In direct numerical simulations of a geodynamo model with stress-free boundary conditions, it has been observed that again two different dynamo solutions, a dipolar and a “hemispherical” one, can both be stable (Christensen et al. 1999; Grote & Busse 2000). Because of the free fluid surface in that model, this might even be taken as a hint for the possibility of non-unique stable states in stellar setups as well.

Fuchs et al. (1999) have demonstrated an even more extreme case with a dynamo powered by a forced laminar flow. In the course of the magnetic field growth, the Lorentz force arranges the flow into a different pattern, which is hydrodynamically stable, but unable to drive a dynamo. As the dynamo dies out subsequently without a chance to recover, it was named “suicidal”.

Hence, the question of the character of the final, saturated stage of a dynamo cannot reliably be answered on the basis of a linear approach and the study of the nonlinear model might unveil very unexpected results. Here, we will show in a simple setup that, while αΩ dynamos do grow faster than α² dynamos, non-linear effects are capable of driving transitions from αΩ modes to α² modes. As the two competing dynamo modes are excited for the same parameter set, i.e., are solutions of the same eigenvalue problem, we refer to them as fratricidal, in reminiscence of the aforementioned suicidal dynamos.

The two astrophysical dynamos for which we have long time-series, the solar dynamo and that of the Earth, both exhibit large fluctuations. The solar dynamo in particular is known to go through prolonged quiescent phases such as the Maunder minimum (Eddy 1976). A conceivable connection with fratricidal dynamos makes understanding how non-linear effects define large-scale dynamo magnetic field strengths and geometries a matter of more than intellectual curiosity.

In Sect. 2 we sketch the mean-field theory of α² and α²Ω dynamos. In Sect. 3 we describe our numerical set-up and briefly discuss the test-field method, a technique to extract the turbulent transport coefficients of mean-field theory from direct numerical simulations. In Sects. 5 and 6 we describe different transition types, and we conclude in Sect. 7.

2. Mean field modeling

In the magneto-hydrodynamic approximation, the evolution of magnetic fields is controlled by the induction equation (1)where B is the magnetic field, J = ∇ × B is the current density in units where the vacuum permeability is unity, and η is the microphysical resistivity. A common approach to (1)is mean-field theory, under which physical quantities (upper case) are decomposed into mean (overbars) and fluctuating (lower case) constituents: (2)The mean can be any which obeys the Reynolds averaging rules, but is frequently assumed to be a spatial one filtering out large length-scales (a two-scale approach). Here we will however use planar averaging, either over the xy plane so that or over the yz plane, that is, , where  ⟨ · ⟩ _ξ denotes averaging over all values of the variable ξ (or volume, if not specified). A mean defined by averaging over y only will also be used (⟨ · ⟩ _y). We humbly ask the reader to consider these definitions carefully, given that the two planar averages will in places be used simultaneously. It is important to remember that the superscript refers to the direction of variation, rather than the direction of averaging.

Under Reynolds averaging Eq. (1)becomes where is the mean electromotive force (EMF) associated with correlations of the fluctuating fields.

Symmetry considerations allow one to write the as a function of the mean-fields in the system. In the case of a planar averaging scheme, this relation becomes (5)where α_ij and η_ij are turbulent transport coefficients, and averaged quantities depend on one spatial coordinate only. The traditional α effect is described by the symmetric part of the tensor α_ij, and requires helicity in the flow. The symmetric part of η_ij describes turbulent dissipation, and, in the isotropic case, appears equivalently to the microphysical resistivity η. It is therefore termed the turbulent resistivity, η_t. When assuming that can be completely represented by the mean magnetic field and its first spatial derivatives, the Taylor series in (5)can be truncated after the term . A more complete formula would include higher spatial as well as temporal derivatives.

2.1. Mean-field dynamo action

Let us assume a large-scale shearing flow of the simple form (6)and velocity fluctuations which are isotropic, homogeneous, and statistically stationary. Consequently, if α_ij and η_ij are assumed to be independent of (the kinematic limit), then they reduce to constant scalars α and η_t¹.

If this system were to contain a y-dependent mean field, the shear would induce field constituents which are proportional to x. We restrict ourselves here to periodic spatial dependencies and hence exclude such unbounded fields. When identifying the y direction with the azimuthal one in a spherical system, this is equivalent to restricting to axisymmetric fields. The evolution of harmonic mean magnetic fields is given by the solution of the eigenvalue problem (7)where , η_T = η_t + η, and . The resulting dispersion relation reads (8)with eigenvalues (apart from the always decaying modes with B_y = 0) (9)It can easily be seen that there are two “pure” modes with particularly simple geometries: the α² mode with k_z = 0 does not depend on S and has the form (10)where the growth rate is and is an amplitude factor. The upper (lower) sign corresponds to positive (negative) αk_x.

In contrast, the α²Ω mode with k_x = 0 does depend on S and has, for S ≫ αk_z (the αΩ approximation) the form In the above, is again an amplitude factor, φ represents, for S > 0 (S < 0), the  ± π/4 (± 3π/4) phase shift between the x and y components of the mean field, and upper (lower) signs apply for positive (negative) values of αk_z; see Table 3 of Brandenburg & Subramanian (2005). The corresponding growth rate is (13)For equal |k|, the αΩ mode grows faster than the α² mode².

A key characteristic of α²Ω solutions is that the growth rate λ has a non-vanishing imaginary part k_zc which results in traveling waves with phase speed c. The wave nature of α²Ω solutions is a significant draw in explaining the solar magnetic cycle. For a characteristic α²Ω dynamo found in direct numerical simulations with a setup described below, we show in Fig. 1 the time-series of rms values of B alongside the traveling wave in the z − t plane (“butterfly diagram”). This solution is similar to those considered recently by Käpylä & Brandenburg (2009). There are other sources for such oscillations however. Admittance of a spatially varying α enables oscillatory and hence traveling wave solutions in pure α² dynamos, see Baryshnikova & Shukurov (1987), Rädler & Bräuer (1987), Stefani & Gerbeth (2003), Mitra et al. (2010).

The mean fields of α² modes are force free, while α²Ω modes cause a potential force which has minimal effect as long as the peak Alfvén speed is subsonic so that the magnetic force cannot generate significant density perturbations. Within kinematics, the induction equation allows for superimposed α² and α²Ω modes. Such a superposition can extend approximately into the non-kinematic regime, and in Sect. 5 we will discuss the implications of the Lorentz forces generated when this occurs.

3. Model and methods

3.1. Numerical setup

We have run simulations of helically forced sheared turbulence in homogeneous isothermal triply (shear) periodic cubic domains with sides of length 2π. The box wavenumber, which is also the wavenumber of the observed mean fields, is therefore k₁ = 1. Unless otherwise specified, our simulations have 64³ grid points. For the shear flow we have taken the one defined by (6). We solve the non-dimensionalized system where U is the fluid velocity excluding the shear flow, is the advective derivative, c_s = 1 is the isothermal sound speed, ρ the density, the viscous force, is the rate of strain tensor, ν is the kinematic viscosity and f the forcing term. We use the Pencil Code³, which employs sixth-order explicit finite differences in space and a third order accurate time stepping method. While our code allows full compressibility, the simulations are only weakly compressible (small Mach number). As in earlier work (Brandenburg 2001), in each time step the forcing function is a snapshot of a circularly polarized plane wave. All these waves have the same handedness, but their direction and phase change randomly from one time step to the next. This forcing provides kinetic helicity. The nondimensional forcing amplitude, f₀ (see Haugen et al. 2004, for details), is arranged such that the Mach number associated with the turbulent rms velocity is of the order of 0.1. In practice, this means that f₀ is in the range 0.02 to 0.05. The wavevectors are taken from the set of vectors that satisfy periodicity and whose moduli are adequately close to the target forcing wavenumber k_f.

The magnetic vector potential is initialized with a weak Gaussian random field, the initial velocity is given by and the initial density is uniform. In Table 1 we have collected the control parameters and some key derived quantities of the model. Two parameters of note are the magnetic Reynolds and Prandtl numbers, (17)To characterize the turbulence, we provide values of α and η_t which represent the corresponding tensors as described in Sect. 2. These were determined using the test-field method (see Sect. 3.2) with test-field wavevector or .

For our purposes, we require the helical turbulence to be strong enough that the α² dynamo can safely be excited. Accordingly, we guaranteed that in all our simulations, Re_M is above the critical value (of the order of unity) for α² dynamos in the corresponding shearless setup (Brandenburg 2009). Further, some of the transitions we will study require long simulation times due to their rarity, which constrains us to modest numerical resolutions. This in turn prevents our (explicit) numerical resistivity from being small, so the turbulent velocities must be reasonably large for the stated super-critical values of Re_M. Choosing furthermore subsonic shear speeds, we are restricted to a modest region of parameter space. In light of these limitations we operate mostly in a Pr_M > 1 regime.

3.2. Test-field method

A fundamental difficulty in extracting the tensors α_ij and η_ij from a numerical simulation of (14)–(16)is that (5)is under-determined. However, turbulent transport depends only on the velocity field, so “daughter” simulations of the induction equation, whose velocity fields are continuously copied from the main run, share the same tensors α_ij and η_ij. It is therefore possible to lift the degeneracy by running an adequate number of daughter simulations with suitably chosen “test” mean fields. We employ this test-field method (TFM); for an in-depth overview see Schrinner et al. (2005, 2007) and Brandenburg et al. (2008a,b). Recently the original method has been extended to systems with rapidly evolving mean-fields, requiring a more complicated ansatz than Eq. (5)(Hubbard & Brandenburg 2009) and to the situation with magnetic background turbulence (Rheinhardt & Brandenburg 2010).

In addition to calculating turbulent tensors for planar averages as described in the references above, we will be interested in tensors that depend both on x and z (that is, are y-averages). For this, we generalize (5)to (18)There are 27 tensor components (as ), so nine test-fields are required, which we choose to be of the form (19)where f_q(x,z) is defined, according to the choice of q, to be one of the following functions:

and B^T is, as standard for test-field methods, an arbitrary amplitude factor. Although the wavenumber of the test fields is usually treated as a varying parameter, we need here to consider only the single value k₁ because the fastest growing and also the saturated dynamos in the simulations are dominated by this wavenumber, the smallest possible in our periodic setup. As is often the case in applications of the test-field method, we will occasionally be faced with unstable solutions of the test problems. We treat that difficulty by periodically resetting the test solutions (see Hubbard et al. 2009). Since it takes a finite time for the stable parts of the test solutions to reach their stationary values, and as this time is frequently close to the required reset time, only limited windows in the time series of the data are valid.

4. Kinematic regime

The fruits of the linear analysis of Sect. 2.1 are not always clear in the kinematic regime of direct numerical simulations. Consider first a low value of Re_M ≈ 5, slightly above its supposed marginal values for both the α² and the α²Ω modes. When starting with a Beltrami field as in (10)with k_x = k₁, the initialized dynamo mode indeed starts to grow exponentially after a short adjustment stage. Its corresponding sibling mode, that is, the α²Ω one, is fed via fluctuations and after a certain delay starts to grow exponentially in turn. Even without including memory effects, the estimated growth rates are in good agreement with those obtained from the transport coefficients measured by the test-field method. Further, the expected geometries of both modes can be identified satisfactorily. If instead one uses the measured transport coefficients to generate an appropriate initial α²Ω field whose geometry is described by (12)with k_z = k₁, it grows again at the right rate. Its α² sibling, however, never appears. Instead a field follows the α²Ω mode, albeit with an rms value one order of magnitude lower: we interpret this as “enslavement” of the field by the growing α²Ω mode: i.e., much like in Sect. 6.3, the magnetic field of the α²Ω mode forces via turbulent motions that can be represented by transport coefficients depending on both x and z. Only the geometry of the α²Ω mode can be safely identified, while the field shows a strong imprint of the forcing wavenumber, furthering the enslavement hypothesis.

Even for higher Re_M the kinematic picture should still be clear in the shearless case, S = 0: when starting with random initial conditions one expects all three (Beltrami) α² modes possible in a box, that is, with |k| = k₁, to grow at similar rates. In Fig. 2 we present time series, a butterfly diagram and spatial cuts for Run A of Table 2 (Re_M ≈ 37) along with its shearless counterpart, where the spatial cuts show the mean fields during the kinematic regime. For S = 0 all three fields , and indeed appear with similar amplitudes and grow at a rate extracted with the test-field method. Their identification as Beltrami fields, however, is all but clear, see Fig. 2, in lower panel. At least a strong |k| = k₁ harmonic is detectable in each of these fields.

The sheared case is even more obscure. During the kinematic stage there are some indications of a dynamo wave, albeit strongly entangled by fluctuating fields, see Fig. 2, middle panel. More worrying however, the growth rates and amplitudes of the and fields are very similar and no pronounced features of a Beltrami field can be found in . The spatial structure of the field, expected to be characterized by k_z = k₁, is even less clear in this respect than that of the field in the shearless case. However, an identifiable dynamo wave is markedly visible after entering the saturation stage where the growth rate of falls below that of .

For an explanation we refer to the occurrence of a small-scale dynamo at the given values of Re_M and S, checked by integrating Eq. (4)with . It dominates the linear stage and we hypothesize that it “enslaves” the and fields, and is responsible for the ragged character of their geometries. The fact that with shear there is slower growth than without (see Fig. 2, upper panel), contradictory to the linear analysis, can perhaps be traced down to the reduction of α by the growing magnetic turbulence, cf. the magnetokinetic contribution to α of Rheinhardt & Brandenburg (2010). After the saturation of the small-scale dynamo, the α²Ω mode can establish itself on the now settled MHD turbulence, whereas the α² mode had perhaps not enough time to take shape. In the shearless case, in contrast, the small-scale dynamo, if it exists at all, has at least a much smaller growth rate and is hence less capable of influencing the mean field growth. The appearance of a small-scale dynamo in this parameter range is plausible in view of the results in Schekochihin et al. (2005) and Käpylä et al. (2008), although those works use different setups.

5. Deterministic interactions of α² and α²Ω modes

5.1. Numerical results

Here we report on the results of our simulations, a first set of which is characterized in Table 2. In Fig. 3 we show time series for Run A, which saw a transition from a z varying α²Ω dynamo () to an x varying α²  dynamo (). As is made clear in the bottom panel, there was a prolonged period where the two modes were coexisting while their relative strengths were changing monotonically. However, note that is stronger than , that is, the α² field is distorted during the transition. Run A was repeated 16 times with the same parameters, but different random seeds, and all these runs exhibited similar behavior. Likewise we performed runs where both the value of η and the numerical resolution (cf. Runs B–D, I, J) were varied. As these additional runs also showed the same transition pattern, we conclude that it is deterministic for this level of shear and forcing. More, we conclude that for these cases the α²Ω mode is unstable to the growth of an α² mode due to non-linear effects. Runs with the dynamical parameters (S, u_rms) of Table 2 inevitably generate α² fields from α²Ω fields after modest times, runs with significantly different shear will usually (for most of the random seeds) exit the kinematic regime into an α²Ω mode, and stay in that mode for a prolonged time with no sign of an α² field (but with random transitions into the α² mode, see Sect. 6).

Nonetheless all our simulation parameters allow runs to occasionally fail to fully enter the α²Ω regime, instead exiting the kinematic regime via a short transient α²Ω phase into the α² mode, as shown in Fig. 4. One might speculate that this “weakness” of the α²Ω dynamo is due to strong fluctuations driven by powerful shear (compare the fluctuations of the α²Ω and α² modes in Fig. 7) but the direct exit into the α² mode has also been seen with lower shear as well. We emphasize that this “direct exit” seems to belong neither to the deterministic transitions discussed here nor to the random transitions discussed later.

5.2. Analysis: mean-field approach

Clearly, the transition from the α²Ω mode to the α² one must be a consequence of the back-reaction of onto the flow. Within the mean-field picture, there are two channels available for it: (i) the back-reaction onto the fluctuating flow, usually described as a dependence of α_ij (more seldom η_ij) on the mean field and (ii) the back-reaction onto to the mean flow by the mean Lorentz force, which might again be decomposed into a part resulting from the fluctuating field, , and one resulting from the mean field, . Here, we will deal with a flow generated by the latter force that straddles the distinction of means and fluctuations: it survives under y-averaging, but vanishes under the xy and yz averaging that reveals the α²Ω and α² dynamos respectively. For simplicity we consider magnetic field configurations that would result from a superposition of linear modes of the αΩ and α² dynamos, given in Eqs. (11) and (10) respectively. Such a situation will inevitably occur during the kinematic growth phase if both dynamos are supercritical, but is only relevant for analyzing the back-reaction onto the flow if it at least to some extent continues into the non-linear regime. Our analysis is linear in nature, so while it provides a qualitative framework for understanding the transition process, it is surely not quantitatively accurate.

In order to be able to consider both and as mean fields under one and the same averaging, we have now to resort to y averaging. Moreover, for the sake of clarity we will occasionally subject the resulting x and z dependent mean fields further to spectral filtering with respect to these coordinates. That is, we will consider only their first harmonics  ~e^{ik₁(x + z)} as mean fields.

Let us represent the mean field  ⟨ B ⟩ _y as superposition of a resembling the (x varying) α² mode (Eq. (10)) and a resembling the (z varying) αΩ mode (Eq. (11)): (20)with z′ ≡ z − ct, recalling that c is the speed of the dynamo wave (Eq. (12)). In the above, π/4 ≤ φ ≤ 3π/4 and G,H,k₁ > 0 are appropriate for α > 0. The parameters G and H capture the difference in the strengths of the y and z components of or the x or y components of , respectively. We expect G > 1 as shear amplifies the y component of an α²Ω mode well above its x component. The inclusion of the parameter H, which is unity for pure α² modes will be justified below, see the different strengths of and in Fig. 3, lower panel.

The mean Lorentz force for the superimposed fields can be written as (21)As the Mach numbers were found to be small throughout, we assume incompressibility and hence drop the potential component ∇Φ. Further, we assume that and the mean velocity driven by it are simply linked by a coefficient , where the total viscosity ν_T is the sum of the molecular ν, and the turbulent viscosity ν_t. Thus we can approximate the mean velocity due to the interaction of the superimposed mean fields as (22)where . Clearly, this flow, having merely a y component, shows quadrupolar geometry in the x − z plane as can be rewritten in the form with a new amplitude and phase, and φ′.

The simulations show indeed a dominant part of that shape in the Lorentz-force generated mean flow as can be seen from Fig. 5. There the y averaged U_y is shown together with its Fourier constituent  ~e^{ik₁(x + z)}. The latter contains approximately 10% of the energy in this component, or , indicating that the assumptions made in deriving (22)are reasonably well justified in a non-linear system.

Upon interaction with a or a of the form (20), the mean flow in (22)generates an and , respectively.

5.2.1. Dominating α²Ω mode

If , then can be treated as a perturbation, and we can drop higher order terms in . Accordingly, the z-averaged EMF due to the flow is (23)The curl of this EMF is (24)If Gsinφ > H, then I > 0 and for H > 0 this EMF reinforces . Thus we see that the inclusion of the parameter H in the ansatz for , Eq. (20), was justified as receives enhanced forcing in comparison to .

5.2.2. Dominating α² mode

If then we can in turn treat as a perturbation. Further, as the system is dominated by the α² mode, we will have H ~ 1. In this case we find (25)and (26)We can write (27)If Hsinφ − G < 0, as expected since H ≈ 1, G > 1, term  [] ₁ in (27)will act to damp , that is, the perturbative α²Ω wave. Further, term  [] ₂ is opposite in sign to the time-derivative of such a wave, so it slows or reverses the direction of wave-propagation.

5.2.3. Mean-field evolution

Here we assume again domination of the α²Ω mode, that is, . With Eqs. (7)and (24)the eigenvalue problem for the modified α² field is then (adopting k_x = k₁, k_z = 0) (28)with eigenvalues (29)Making the approximation I ≫ αk₁, similar to the αΩ approximation S ≫ αk₁, we find (30)The above should be compared with the growth rate of the αΩ dynamo, λ^αΩ from (13)which is not touched by the occurrence of I. The αΩ dynamo saturates when α has been quenched such that the product αS settles at the marginal value . If the parameter I becomes comparable with the shear, i.e., I ~ S, then might grow even when the α²Ω field is saturated, i.e. . In other terms, the saturated α²Ω mode is unstable to the growth of a fratricidal α² field, so the transition will take a well defined time from the onset of the non-linear stage, determined by λ^X.

We test this theory for Run A at the time of Fig. 5, t = 1.45t_res, extracting G and H from the relative strengths of the x and y or y and z components of the mean fields or , respectively, after a projection onto the first harmonics; see Eq. (20). The parameter I is calculated from the magnetic and velocity fields using (31)with , where is the amplitude of the quadrupolar constituent of the velocity field seen in Fig. 5. We find , H ≃ 2.9, G ≃ 4.9, I ≃ 0.09, and confirm that φ ≃ π/4. As I > S = 0.05, the growth of the x varying mode even when the α²Ω mode is saturated is not surprising. Repeating this run16 times with different random seeds (keeping the control parameters fixed) changed the occurrence time of the transition by only one resistive time, suggesting that the transition is an essentially deterministic process.

We have never seen a reverse transition from the α² state back to the α²Ω state. This may be understood in terms of interacting modes, with the α²Ω mode being suppressed once the α² mode is dominating; see Sect. 5.2.2.

6. Random transitions

6.1. Non-interacting modes

Not all transitions fit into the above deterministic picture of interacting α²Ω and α² modes. The behavior above and below the shear value S = −0.05 that applies to the runs in the previous section is qualitatively different. In Fig. 6 we present a set of time series of the rms values of and , all related to Run F of Table 3. Secondary runs were performed by branching off from the original simulation either at t = 5t_res, when the α²Ω mode is well established and stationary, or at t = 8.25t_res, immediately before the transition to the α² mode in the primary run is launched. The only difference between all these runs is in the random seed, which is used by the forcing algorithm. In all, the time until the transition starts varies by  ≈ 2.5t_res, and many more turbulent turnover times ( turbulent turnover times per resistive time). The time elapsed during a transition is always of the order of t_res/2, unlike 3t_res for the process seen in Fig. 3. Thus it is suggestive to assume that there might be a very slow, still essentially deterministic process, preparing the transition, which is likely resistive in nature as that is the longest obvious “native” timescale of the system. Slow resistive effects are known to exist in dynamos, for example the slow resistive growth of α² dynamos in periodic systems. However, transitions can indeed occur at very different times including the extreme case in which a run never develops a quasi-stationary α²Ω mode, but instead enters the α² state almost immediately after the end of the kinematic phase, see Fig. 7 (run G of Table 3). We believe therefore that under certain circumstances the transition process is not a deterministic one, in that it is impossible to predict or at least estimate the time until the transition. Figure 8 is a synopsis of simulations that belong to that type, hence do not show the instability discussed in Sect. 5. Note that, while corresponding setups without shear are known to enable α² modes for the entire Re_M range studied (Brandenburg 2009), we cannot rule out that the α² mode is sub-critical for Re_M = 10,S = −0.1 (i.e., it is possible that the run denoted by the bottom left diamond in the upper panel did not transition because it simply could not).

The picture of random transitions is different from the interacting mode one of Sect. 5 in several interesting ways. Crucially, the α²Ω mode is here at least meta-stable against growth of the α² mode, as evinced by its prolonged life-time (hundreds of turbulent times) and the small magnitude of , which further is not dominated by an α² mode. A reasonable working hypothesis for the cases of Sect. 5 is then that, there, the α² mode is the only stable solution and, as soon as the nonlinear stage has been entered, it starts to devour the α²Ω one, settling after a time which is related to basic parameters of the system and hence not random. In contrast, for the cases considered here, we conclude that both the α²and the α²Ω solutions are indeed stable (not only metastable) and the latter has a well extended basin of entrainment. Due to its higher growth rate the system settles first in the α²Ω mode and suppresses the α² mode efficiently. A transition to the latter can only occur if a random fluctuation in the forcing is strong enough to push the system over the separatrix into the basin of entrainment of the α² mode. This may happen after only a long time or even never during a simulation’s run time, cf. diamond and square symbols, respectively, in Fig. 8, see also the absence of a transition in Fig. 1 where S = −0.1, Re_M = 20. From these data one can also derive that the tendency to leave the α²Ω state by a random transition increases both with growing shear and Re_M (note the corresponding reduction in the blue squares’ size  =  time until transition in Fig. 8). Although one would expect the basin of entrainment of the α²Ω mode to enlarge with growing S, it is conceivable that at the same time the level of velocity fluctuations is raised, making transitions eventually more likely. Growing Re_M can readily be thought of increasing the fluctuation level and hence the transition probability likewise. As already mentioned in Sect. 5, we see in addition “direct exits” from the kinematic stage into a saturated α² dynamo; see Fig. 8 (cross symbols) and Fig. 9.

Given that the examples for the first scenario (Table 2) differ from those for the second (Table 3) mainly in their lower rate of shear, our conclusion seems reasonable as stronger shear should result in a clearer preference of the α²Ω mode because the α² mode does not feel the shear. Or, in other terms, from a certain shear rate S on, the α²Ω mode should acquire a basin of entrainment with a finite “volume” that grows with S. If this picture is true, transitions in the two scenarios should have clearly different characteristics, and indeed, the transition in Fig. 6 is markedly faster than that seen in Fig. 3.

Once the shear rate is dropped markedly below S = 0.05, the typical value for deterministic transitions, the system again seems to settle in a domain similar to the high-shear () one, namely one with both the α² and the α²Ω modes being stable. However, entering into the α² state immediately after the kinematic stage (the “direct exit”) appears now to be more likely, while at the same time random transitions away from a settled α²Ω state become extremely rare: Over a run time of 3 × 10³   t_res accumulated over several runs with S = 0.01,0.02 as well as target values of Re_M around 30, 60, and 90, only one transition was observed (for Re_M = 80,Pr_M = 5, S = 0.02, after 8 × 10²   t_res); see also Fig. 8, lower panel. Plausibility arguments for the occurrence of this “low shear” domain, in particular for the apparent regaining of a finite basin of entrainment by the α²Ω mode which has to be implied, are not in sight.

As in the transitions discussed in Sect. 5, we have not here seen the α² mode transition back into the α²Ω mode. Some attempts were made to provoke this reverse transition by perturbing the α² state with a (sufficiently strong) α²Ω mode. While in some runs it indeed took over, velocities were attained for which the numerics are unreliable, and often proved numerically unstable, making the results inconclusive. However, such a behavior is not entirely surprising as the α²Ω saturation process can anyway be somewhat wild, cf. Fig. 4.

The absence of spontaneous reverse transitions appears plausible insofar the time variability of the α² mode is much smaller than that of the α²Ω mode, which can clearly be seen in Fig. 10 for Run H. That is, events capable of pushing the system over the separatrix are simply much rarer. Significantly longer integration times are likely needed for their eventual detection, but it is also conceivable that the triggering event never shows up.

6.2. Large scale patterns

Run H will be examined here in more detail. Curiously,  ⟨ U_y ⟩ _y taken just during the transition as shown in Fig. 11 does not show the quadrupolar pattern of Fig. 5. It is therefore not surprising that the butterfly diagrams in Fig. 12 do not show a direct transition from the α²Ω to the α² dynamo, as develops significantly later than . This is clearly visible in Fig. 10, lower panel. As consideration of the mean flow due to the Lorentz force of the mean field alone is obviously not fruitful in explaining this transition, we recall that the back-reaction of the mean field onto the turbulence opens another channel of nonlinear interaction.

According to elementary mean-field dynamo theory, the α effect is caused by the helicity in the flow: α ~  ⟨ w·u ⟩ , where w ≡ ∇ × u is the fluctuating vorticity. Further, the back-reaction of the mean field on the turbulence, which saturates the dynamo, is assumed to be captured by the current helicity  ⟨ j·b ⟩ . It is often related to the magnetic helicity  ⟨ a·b ⟩  and thought to reduce the original α by producing a magnetic contribution of opposite sign. In Fig. 13 we present time-series of the power spectra of these helicity correlators across the transition. However, we see no clear signal around the transition event.

6.3. Mean-field analysis with y averaging

To examine the problem more closely, we recall Eq. (18) for a mean defined by the y average: (32)It is clear that the Fourier constituents of α_ij and η_ijk with wavenumber k₁ in both x and z (the quadrupolar constituents) can create an EMF out of a field , both with the same wavenumber: (33)where the superscripts indicate the coefficients to be the Fourier constituents  ~e^{ik₁(x + y)} and is assumed to vanish. Note that each of them is actually given by four values, e.g., the two amplitudes and phases in: (34)The coefficients relevant for the generation of (from alone) are and . We have used the test-field method (see Sect. 3.2) to find them and present the results in Fig. 14. They turn out to be surprisingly large, when compared to the rms velocity (e.g., ) and some may show a trend across the transition from the α²Ω to the α² field. This overall trend is hypothesized to be due to the increase in u_rms that accompanies the transition from a stronger α²Ω field to a weaker α² field with less potential to inhibit the flow. It is interesting that the large transport coefficients (ignoring and which are without effect) are all those which generate an out of , i.e., out of the component that feels the effect of shear. We speculate that these coefficients feel the shear quite strongly.

7. Discussion and conclusions

We have demonstrated that, while α²Ω modes are kinematically preferred to α² modes in homogeneous systems that support both, the α² mode acts in a fratricidal manner against the former after the nonlinear stage has been reached. This transition can occur in at least two different fashions. Further, we have not observed the reverse process. One of the two transition processes, based on superposed α²Ω and α² modes, operates in a basically deterministic fashion through a large-scale velocity pattern generated by the interaction of the modes. In contrast, we interpret the mechanism of the second process, which may start only many resistive times past the saturation of the α²Ω dynamo, by assuming that both the α²Ω and the α² modes are stable solutions of the nonlinear system. Transitions occur if due to the random forcing a sufficiently strong perturbation builds up which tosses the system out of the basin of entrainment of the α²Ω mode into that of the α² mode. This hypothesis is bolstered by both the random timing of these transitions and by the large time-variability seen in the amplitude of the α²Ω field. A return seems to be much less likely as the level of fluctuations of the α² mode is, by contrast, greatly reduced. Random transitions are influenced by the shear rate. Above the shear level of S = 0.05 associated with deterministic transitions, eventual random transitions appear to be inevitable. Markedly below that level they are, at best, extremely rare. This suggests strongly that the vulnerability of an α²Ω mode to transitions increases with its dynamo number.

These results fit with earlier work studying dynamos whose non-linear behavior is fundamentally different from their linear one (e.g., Fuchs et al. 1999). While our simulations are limited to Cartesian, cubic, shearing-periodic domains, they are particularly exciting given that the only dynamo which has been observed over a long baseline and which could be either α²Ω or α²  , the solar dynamo, indeed shows differing modes of operation (regular cycles vs. deep minima). The results are also disturbing in that we have evidence for non-deterministic, rare (as they occur on scales of multiple resistive times or hundreds of turbulent turnovers) mode changes that show no evidence for a return. While a bifurcation between different stable modes has long been an acknowledged possibility for dynamos (Brandenburg et al. 1989; Jennings 1991; Covas & Tavakol 1997), a rare, stochastic, possibly uni-directional transition is perhaps the most troublesome consequence of such bifurcations except for the ultimate self-extinction.

Given that the mean-field transport coefficients have been systematically derived for many of the runs, it may sound plausible to study the nonlinear stability of different dynamo regimes by means of mean field models. This would allow to construct a stability map for the system to test the hypothesis made in the paper. However, there are two difficulties: firstly, the derived transport coefficients apply only for the specific mean field for which they have been calculated, i.e. they are not the functional derivatives of the mean EMF with respect to the mean field or current density and to our knowledge the technology to do better does not yet exist. This deficiency applies only to the non-linear regime, of course. Secondly, to overcome this restriction one could assume that the transport coefficients are smooth functionals of the mean field and so can be obtained from a finite set of data points by interpolation. This would require only to perturb the system into neighboring states, from which, however, it would tend to return to the original state in a self-regulated way. As the transport coefficients can only be calculated for sufficiently slowly changing flows, the set of obtainable data points would be strongly limited.

The αΩ dynamo is believed to be common and important for systems like the Sun or accretion disks, which all have long life-times compared to turbulent turnover times. It is then a daunting possibility that we could be forced to stretch our simulations over very long temporal base-lines to find the actual long-lasting field configuration. More positively, our result, while in a different geometry, increases the importance of recent work on non-oscillatory αΩ and oscillatory α² modes in spherical shells for the solar dynamo (Mitra et al. 2010; Schrinner et al. 2011).

---

Acknowledgments

We thank the anonymous referee for having made detailed suggestions that have helped improving the clarity of the paper. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the National Supercomputer Centers in Linköping. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952 and the Swedish Research Council Grant No. 621-2007-4064.

References

All Tables

All Figures