A&A 376, 713-726 (2001)
DOI: 10.1051/0004-6361:20011041

## Magnetoconvection and dynamo coefficients:

### Dependence of the effect on rotation and magnetic field

M. Ossendrijver1 - M. Stix1 - A. Brandenburg2,3

1 - Kiepenheuer-Institut für Sonnenphysik, Schöneckstraße 6, 79104 Freiburg, Germany
2 - Department of Mathematics, University of Newcastle upon Tyne, NE1 7RU, UK
3 - Nordita, Blegdamsvej 17, 2100 Copenhagen Ø, Denmark

Received 24 October 2000 / Accepted 12 July 2001

Abstract
We present numerical simulations of three-dimensional compressible magnetoconvection in a rotating rectangular box that represents a section of the solar convection zone. The box contains a convectively unstable layer, surrounded by stably stratified layers with overshooting convection. The magnetic Reynolds number, Rm, is chosen subcritical, thus excluding spontaneous growth of the magnetic field through dynamo action, and the magnetic energy is maintained by introducing a constant magnetic field into the box, once convection has attained a statistically stationary state. Under the influence of the Coriolis force, the advection of the magnetic field results in a non-vanishing contribution to the mean electric field, given by . From this electric field, we calculate the -effect, separately for the stably and the unstably stratified layers, by averaging over time and over suitably defined volumes. From the variation of we derive an error estimate, and the dependence of on rotation and magnetic field strength is studied. Evidence is found for rotational quenching of the vertical effect, and for a monotonic increase of the horizontal effect with increasing rotation. For , our results for both vertical and horizontal effect are consistent with magnetic quenching by a factor . The signs of the small-scale current helicity and of the vertical component of are found to be opposite to those for isotropic turbulence.

Key words: convection - magnetic fields

### 1 Introduction

Magnetic fields are observed on a wide variety of cosmical bodies, among which are planets, stars and galaxies. With the exception of a few types of objects whose magnetic fields are thought to be frozen-in relic fields, cosmical magnetic fields are attributed to dynamo action. Dynamo theory concerns the generation of magnetic fields in electrically conducting fluids. In stars and planets, dynamo action is the result of an interplay between convection and rotation. The present simulations focus on dynamo action in late-type stars, which are characterized by an outer convective zone on top of a stably stratified radiative interior. Many late-type stars are magnetically active, and some exhibit magnetic cycles, as does the Sun.

A successful solar model was first proposed by Parker (1955) who recognized that shear and helicity cooperate in such a way that the mean magnetic field oscillates in a migratory manner; later Yoshimura (1975) has shown that the field generally propagates along the surfaces of constant angular velocity. Since then, magnetic cycles and butterfly diagrams have been produced by mean-field models using spherical geometry (Steenbeck & Krause 1969), with depth-dependent magnetic diffusivity (Roberts & Stix 1972), and with the character of an interface wave at the base of the convection zone (Parker 1993; Charbonneau & MacGregor 1997). As an alternative to the helical effect of convection (the  effect), the systematic tilt of flux tubes that become unstable and erupt to the solar surface has been employed as a regenerative agent for the mean poloidal field (Leighton 1969; Durney 1995, 1997; Schmitt et al. 1996; Dikpati & Charbonneau 1999), with equal success regarding the existence of Sun-like cycles and butterfly diagrams. In all cases, however, the regenerative terms in the mean-field equations were based either on approximations such as first-order smoothing, or on the limited available knowledge about the behavior of magnetic flux tubes in the convection zone; the review of Stix (2001) summarizes some of the current problems.

In principle, the problem can be attacked by direct numerical integration of the equations of magnetohydrodynamics (MHD) for the star's convection zone (for a review, see e.g. Nordlund et al. 1994). But severe problems are encountered. Firstly, a stellar convection zone is a highly turbulent plasma ( ), so that the kinetic energy spectrum encompasses a very wide range of length scales. If these were all to be resolved, a prohibitive number of mesh points would be required. Secondly, since the magnetic Reynolds number is also large ( ), a similar argument holds for the length scales of the magnetic field. For these reasons three-dimensional MHD simulations often have been restricted to a rectangular box that represents only a small section of a stellar convection zone. In the present work we do this, too. Thirdly, the Prandtl number, the ratio of the kinematic viscosity to the radiative diffusivity, is very small in the solar convection zone ( -10-7), so that the flow can vary on much smaller scales than the temperature. This is known to have consequences for the topology and temporal behavior of convective turbulence (Cattaneo et al. 1991; Brandenburg et al. 1996; Brummell et al. 1996), but in the present simulations, we ignore such effects, and set . Fourthly, the Mach number, , decreases with depth in the solar convection zone from about 0.1 in the photosphere to less than 0.001 near the base. The corresponding high sound speed in the lower convection zone necessitates the use of very small time steps in order to fulfill the CFL condition for compressible hydrodynamics ( ), while it is often desirable to continue simulations for at least several convective turnover times. One possibility to circumvent this problem is to adopt the anelastic approximation whereby sound waves are excluded, so that larger time steps can be used (Ogura & Charney 1962; Lantz & Fan 1999). Here we maintain full compressibility, and attain a mean Mach number of the order 0.1; in a subsequent paper we shall study effects of a smaller Mach number.

Brandenburg et al. (1990) extended a 3D hydrodynamical code to the case of magnetoconvection, including the effect of rotation. Subsequently, a spontaneous dynamo instability was observed in the simulations (Nordlund et al. 1992; Brandenburg et al. 1996). Recent simulations of isotropically forced helical turbulence (Brandenburg 2001) have verified the existence of large scale dynamo action, and it was possible to identify this as the result of an effect (in the sense of a non-local inverse cascade). However, the large scale field generated possesses magnetic helicity and, since for closed or periodic boundaries the magnetic helicity can change only resistively, the growth of the large scale field is slowed down as the magnetic Reynolds number increases. This translates inevitably to a magnetic-Reynolds-number dependent effect and turbulent magnetic diffusivity, as suggested by Vainshtein & Cattaneo (1992). The argument of Vainshtein & Cattaneo is however only phenomenological, not based upon the fundamental concept of magnetic helicity conservation, and hence their conclusion that strong large-scale fields are impossible is not borne out by the simulations of Brandenburg (2001). Furthermore, in that paper it was shown that the obtained by imposing a magnetic field is indeed a reasonable approximation to the that results naturally even if no field is imposed.

In the present simulations we do not calculate box dynamos, but rather concentrate on the dynamo coefficients that occur in mean-field theory. Therefore, all simulations are performed with an imposed magnetic field, and the magnetic Reynolds number is chosen subcritical, so that the magnetic field is not self-sustained. We think that progress in stellar dynamo theory can be made through a combination of exact MHD simulation and mean-field dynamo theory. Furthermore, the existence of systematic, large-scale magnetic fields and stellar cycles with periods far in excess of convective time scales suggests that a mean-field description is possible in some form. Also, we argue that mean-field theory (in the wide sense that includes the above-mentioned models such as Leighton's) is the only current model that reproduces large-scale magnetic fields and cycles, including such an outstanding feature as the solar butterfly diagram, even though the diverse approximations made for calculating the dynamo coefficients may not be valid under stellar conditions. The case of self-excited dynamo action would of course be very interesting too, especially in view of the question of whether -quenching depends on the magnetic Reynolds number, as suggested by Vainshtein & Cattaneo (1992). However, we postpone this until a later study, partly because the measurement of the dynamo coefficients in the presence of a large-scale field that is different from the imposed one is less straightforward.

The present investigation is limited to the regime where the coefficients are determined directly by the flow, and we do not address the question of whether the large-scale magnetic field of late-type stars is in fact generated by such an ordinary effect in the convection zone, or by a magnetic instability in an underlying stably stratified layer (e.g., Brandenburg & Schmitt 1998). We do include a stably stratified layer with overshooting convection, but its main purpose here is to provide realistic conditions for the flow in the unstable region. In most runs, the strength of the imposed field is set to a value which amounts to typically of the equipartition field with respect to the kinetic energy density. During the subsequent evolution the field strength also remains small compared to the equipartition value. Furthermore we shall explore the influence of the imposed magnetic field by increasing its strength up to values somewhat in excess of the equipartition value. Therefore the results may have relevance for the solar convection zone, where the magnetic field is no stronger than the equipartition value. The effect produced in the galactic gas by supernova explosions has been calculated by Ziegler et al. (1996), in a similar spirit as in the present work.

The relevance of the transport coefficients for mean-field dynamo theory becomes clear from the equation for the mean magnetic field,

 (1)

where is the mean magnetic field, is the mean flow, is the fluctuating component of the flow, and is a contribution to the mean electric field, sometimes referred to as the electromotive force. The dynamo coefficients appear in an expansion of this mean electric field in terms of spatial derivatives of the mean magnetic field. In general they can be represented as kernels of an integral equation (Brandenburg & Sokoloff 2001). In the simplest case, the coefficients are treated as local tensors, which leads to

 (2)

The first term in the expansion is the effect. The alpha tensor is a pseudo tensor that exists only in non mirror-symmetric flows, as they occur in stellar convection. For the solar dynamo, the main significance of the effect lies in generating the poloidal mean magnetic field, which is achieved predominantly by . Toroidal fields are generated by the strong differential rotation that exists in the tachocline near the base of the solar convection zone. In other solar-type stars and fully convective stars though, differential rotation could be weak (Küker & Stix 2001), so that the dynamo becomes of the -type. In that case, other components of the alpha tensor, such as , are important for maintaining the toroidal magnetic field.

The remainder of the paper is structured as follows. After introducing the model and the equations, we focus on the dependence of on rotation. In the following section, effects of the strength and orientation of the imposed field are studied. In the final section, a discussion of the results is presented.

### 2 The model

The model used in the simulations is the same as that of Brandenburg et al. (1996). The simulation domain consists of a rectangular box which is defined on a Cartesian grid, with x and y denoting the two horizontal coordinates (corresponding to latitude and longitude in spherical coordinates), and z denoting depth, respectively (Fig. 1). In most cases, convectively stable layers are included below and above the unstable layer. The upper layer (region 1) is stabilized by a cooling term in the energy equation, which leads to an almost isothermal, highly stable stratification. The lower stable layer (region 3) represents an overshoot zone, whose thickness is chosen such that overshooting plumes do not reach the lower impenetrable boundary. The box rotates about an axis , that makes an angle with the z-axis. Hence corresponds to a situation where the box is located at the south pole of the Sun. In the present paper we treat only this special case.

 Figure 1: Geometry of the simulation domain. The layer boundaries are at z1=-0.15, z2=0, z3=1, and z4=2.85, except for run 10 that has no stable regions, so that z1=z2=0 and z3=z4=1. The horizontal extent is 4. Open with DEXTER

The governing equations are those describing magnetic induction, mass continuity, and the balance of momentum and energy:

 (3) (4) (5) (6)

where is the current density, is the magnetic diffusivity, is the kinematic viscosity, is the adiabatic index, and is the radiative conductivity. In all simulations, and are taken constant, and is a prescribed function of depth. The stress tensor is given by

 (7)

and stands for . The equation of state is that of an ideal gas, i.e.,

 (8)

The term Q, given by

 (9)

represents cooling or heating, depending on whether the internal energy density exceeds e1 or falls below it, respectively; represents the cooling/heating rate. The depth-dependent function f ensures that Q is non-vanishing only in region 1, and f(z2)=0. The inclusion of this term leads to the presence of a highly stable, thin overshoot layer, thereby providing a realistic upper boundary condition for the convection zone. In the horizontal directions, periodic boundary conditions are imposed. The upper and lower boundaries of the domain are impenetrable and stress-free, and the horizontal components of the magnetic field variation are set to zero. At the lower boundary, the energy flux is prescribed, and at the top boundary the internal energy (which is proportional to the temperature) is fixed.

 (10)

In the actual runs we set B0y=0; B0x is either 0 or the imposed horizontal field.

All quantities are made dimensionless by setting

 (11)

where is the initial density at a depth z=z3, the bottom of the unstable layer. Thus length is expressed in terms of d=z3-z2, the thickness of the convectively unstable zone (region 2). It follows that time is measured in terms of , velocity in terms of , the magnetic field strength in terms of , and entropy in terms of .

On several occasions we shall consider the energy balance, which is governed by a conservation law,

 (12)

where is the total specific energy, and the fluxes, , are given by
 (13) (14) (15) (16) (17)

Apart from the geometry, boundary conditions, and initial conditions, solutions are governed by the following 10 independent dimensionless parameters. The Prandtl number and the magnetic Prandtl number are defined as

 (18)

where denotes a reference value of the thermometric (radiative) diffusivity for region 2. The parameter determines the pressure scale height at the top of the box, :

 (19)

The initial thermal structure of the (maximally) three regions is characterized by the radiative temperature gradients,

 (20)

In addition, one often employs the polytropic index, . The adiabatic temperature gradient is given by . A measure for instability is provided by the superadiabaticity, , which is positive in an unstably stratified medium. The Rayleigh number is defined as

 (21)

where is the pressure scale height in the middle of the unstable layer, as can be shown using the hydrostatic equilibrium (25). The parameters and refer to the unperturbed stratification, i.e. that before the onset of convection. The Rayleigh number is a measure of the strength of convection compared to that of viscous and thermal (radiative) dissipation, as can be seen by writing , where , , and . Alternatively, one may express in terms of the entropy gradient, ds/d . According to the Schwarzschild criterion, a positive value for , i.e. , signifies instability. In reality, Ra must exceed a finite threshold value for convection to set in. It should be noted that the local value of the Rayleigh number, , varies with depth within the unstable layer, because , and, as a result of convection also , are z-dependent. Typically, the local value is several times smaller than , mainly because can be strongly reduced by convection. The Taylor number,

 (22)

measures the importance of rotation relative to viscous dissipation. Finally, represents the rate at which internal energy is lost from the upper stable layer, and is the angle between the rotation vector and the z-axis.

All other parameters are secondary. The Coriolis number, or inverse Rossby number, measures the importance of the Coriolis force and is defined as

 (23)

where is the turnover time. The correlation length, , is taken to be d in the unstable region. The Chandrasekhar number,

 (24)

measures the strength of the imposed magnetic field. In the initial state, the z-component of the radiative energy flux, de/dz, is assumed to be constant throughout the domain. This determines the radiative conductivities in the three regions according to . In fact, is turned into a smooth function of depth by allowing it to change continuously across thin intermediate layers between the three regions. An approximate initial stratification, unperturbed by convection, is calculated on the assumption of hydrostatic equilibrium, and this is done iteratively until the condition is satisfied. It is also assumed that region 1 is cooled efficiently enough to become isothermal. The result is then basically a smoothed version of

 e = (25) = (26)

where stand for the three regions, , , and , respectively. The actual initial stratification in region 2 is calculated numerically using the mixing-length formalism of convection. The advantage of this approach is that it reduces the amount of time required for relaxation to a fully convective state. The internal energy density at the top equals . Using (25), it is easily shown that (de/d . The radiative, kinetic and magnetic diffusivities follow from Eqs. (18) and (21). Reynolds numbers are defined as

 (27)

where is the rms velocity defined in a suitable way (e.g., by averaging over time and over a partial volume of the box).

 run 1 2 3 4 5 6 7 8 9 10 Ta 2 10 30 100 300 2000 T ( ) 206 822 206 841 698 545 940 362 376 554 T ( ) - - - 564 564 456 508 346 290 521 Ma 0.089 0.088 0.088 0.084 0.091 0.089 0.080 0.070 0.056 0.10 Re 34 34 34 34 35 34 30 30 20 38 Co 0.042 0.093 0.16 0.3 0.5 1.3 3.3 5.4 11 2.6

We employ a finite difference scheme, according to which spatial derivatives are calculated with 6th-order accuracy (Lele 1992). Time-stepping is done using a third-order Hyman predictor-corrector method. Table 1 gives a list of the parameters used for a first series of runs in which the influence of rotation is investigated by varying the Taylor number.

 Figure 2: Averaging procedure (run 6). Top: at each time step, is calculated as a function of depth by averaging over the horizontal coordinates (grey curves). The thick drawn curve is the time average, and the surrounding thin drawn curves indicate the error in the mean. Bottom: is calculated separately for two partial volumes, namely for the depth range -0.15

### 3 The effect

Essentially, the effect amounts to a proportionality between the mean electromotive force, , and the mean magnetic field. Magnetoconvection in a rotating, stratified medium provides the necessary anisotropies for a non-zero . In the present case, with the rotation axis chosen parallel to the direction of gravity, the most elementary, non-isotropic form of the expression for the effect is one that distinguishes vertical and horizontal components:

 (28)

The components and correspond to and in spherical coordinates, respectively. We recall that plays the dominant role in the generation of the poloidal magnetic field. The two components of are extracted from the simulations by constructing and . In practice, in these expressions is replaced by , because the difference, , was found to be negligible in every case. Depending on which component of , i.e. or , is to be determined, a magnetic field is imposed in the z-direction or in the x-direction, respectively.

 Figure 3: Row 1: vertical effect. Row 2: horizontal effect. Row 3: kinetic helicity. Row 4: gradients of and , and their sum (rows 3 and 4 are for a vertical imposed field, but since the field is weak the results in these rows would be the same for a horizontal field). Row 5: current helicity for a vertical imposed magnetic field. Left column: run 6 ( ). Middle column: run 7 ( ). Right column: run 10 ( ; no stable layers). In all cases the unstable region is between z=0 and z=1. All calculations are for the solar south pole. Error margins are given for rows 1 and 2 as thin lines enclosing the mean; those in rows 3-5 would be comparable or narrower, but are omitted in order not to overload the drawings. Open with DEXTER

#### 3.1 Statistics

It is well known from other numerical simulations that is an extremely noisy quantity, so that one should do as much averaging as possible in order to optimize the statistics.

For this paper we have chosen two averaging procedures. The first, which is considered in the present section, consists of an average over the horizontal coordinates and over time (Fig. 2, upper panel). The time average is initiated from the moment when the magnetoconvection has attained a statistically stationary state, i.e. when the kinetic and magnetic energy densities are more or less constant. The second is an additional average over two suitably chosen depth ranges (Fig. 2, lower panel): since the effect depends on depth, and is expected to change sign near the bottom of the unstable layer, the averaging over depth is performed separately for the unstably and the stably stratified layers. This second procedure will be used mainly in later sections.

The statistical deviation of from its mean value generally is large, often larger than the difference between that mean and zero, as illustrated in Fig. 3. Therefore, extended simulations are necessary in order to obtain a significant result. In the present paper the error estimate is based on the assumption that a single simulation can be divided into a number of time intervals containing independent realizations of the small-scale flow and magnetic field. The length of those time intervals should then exceed the correlation time, which is 20 (in dimensionless units), as determined from the width of the autocorrelation function. To be on the safe side, we smoothed curves of and for several runs by performing a box average, and then calculated the autocorrelation function. This procedure yields a typical coherence time , which is several times longer than the value derived from the unsmoothed curves.

Our procedure of estimating errors is more primitive than that applied earlier to results of numerical simulations (Pulkkinen et al. 1993). Moreover, the errors of the mean values determined in this manner can only be crude estimates because firstly we cannot be sure that we deal with Gaussian statistics (although inspection of the examples of Fig. 3 suggests that varies in a nearly symmetric manner around its mean), and secondly the number of time intervals is not very large. Nevertheless we think that the procedure indicates whether or not the obtained mean values are significantly different form zero. Test simulations with identical parameters but different initializations confirm this indication. Of course, the resulting error bars do not imply anything about the physical assumptions made. In particular, we cannot expect that models based on different assumptions, such as ours and that of Rüdiger & Kitchatinov (1993), will yield results that lie within error bars obtained from statistics (these bars shrink to zero for ). In contrast, we must see whether we can find qualitatively similar or dissimilar results from different models, especially as none of the models yet meets the real Sun.

#### 3.2 Relations between and turbulence properties

The existence of the effect is attributed to the helical nature of convective flows in a rotating medium. Figure 3 shows , and several quantities that, as is shown below, have a relation to the effect. Three different runs are selected in order to demonstrate the effects of rotation and boundary conditions.

The coefficients are depth-dependent, and undergo a sign change within the unstable layer. For , the vertical alpha coefficient exceeds the horizontal coefficient, and has the opposite sign. If rotation is stronger, the vertical effect is strongly reduced.

For isotropic turbulence, the first-order smoothing approximation (FOSA) yields that can be described by a single scalar

 (29)

where is the correlation time of the convection, and is the kinetic helicity (Krause & Rädler 1980). An alternative representation, valid under the assumption of FOSA, relates to the current helicity (Keinigs 1983; Rädler & Seehafer 1990; note the extra minus sign in Keinigs' definition of ):

 (30)

It should be noted that (30) does not reflect the nonlinear feedback of the Lorentz force on the flow, but is a purely linear result. If one allows for mild anisotropy in the radial direction, alpha assumes the form of a tensor, and a detailed calculation (Steenbeck & Krause 1969) reveals that, for slow rotation ( ), the diagonal term is given by

 (31)

While conditions under which expressions (29)-(31) hold are fulfilled neither in the simulations nor in the Sun (mainly because of the requirements for a short correlation time and only mild anisotropies in the turbulence), it is still instructive to compare them with the numerical results.

The kinetic helicity is positive at the top of the convective layer (Fig. 3, row 3), as is expected for the southern hemisphere based on the direction of the Coriolis force for contracting sinking parcels, or expanding rising parcels. Brummell et al. (1998) find kinetic helicity with the same sign for a box on the northern hemisphere, but their coordinate system turns out to have a left-handed orientation. Near the center of the convective layer, the sign of the kinetic helicity reverses if rotation is sufficiently strong, because below a certain point sinking parcels expand laterally while approaching the level z=1 (i.e., the stably stratified overshoot layer or the impenetrable lower boundary if an overshoot layer is absent), while rising parcels contract while moving away from it. For , no reverse-helicity region seems to exist. Nevertheless, in these cases alpha still changes its sign near the bottom of the convective layer. For comparison we have calculated a case (run 10) with the same Taylor number as in run 7, but where the stably stratified layers are replaced by impenetrable boundaries at z=0 and z=1. In this case the kinetic helicity near the bottom of the unstable layer is more strongly negative, and has a magnitude comparable to that in the upper part of the unstable layer (Fig. 3, row 3, Col. 3). This can be attributed to the fact that an impenetrable boundary forces a sinking parcel to diverge more strongly than does a stably stratified layer. Thus, except for the stable region in some cases, the sign of as predicted by Eq. (29) roughly agrees with that of , while it is opposite to that of . The deviating sign of was explained by Brandenburg et al. (1990). We note that the same unconventional sign of was also found by Ferrière (1992), and Rüdiger & Kitchatinov (1993), although not for the same values of the Coriolis parameter as in our simulations.

The gradient of is positive at the top of the convective layer and changes sign at a point near the bottom of the unstable layer in the runs with an overshoot layer (Fig. 3, row 4, Cols. 1, 2). Hence the sign of as predicted by (31) agrees with that of in the simulations, also if rotation is weak. In this respect, agrees qualitatively more closely than does the kinetic or the current helicity. In the upper region, is positive because the density gradient dominates, while near the bottom of the unstable domain the influence of the convective stability is felt, resulting in a reduction of the turbulent velocity sufficient to produce a negative sign. With regard to , the sign and the dependence on Co differ strongly from (31). Also, in the run without the stably stratified layer, the gradient of the density dominates throughout the box, and no sign change of is observed (Fig. 3, row 4, Col. 3). Both components of do exhibit a sign change, though. This inconsistency may be a result of the impenetrable boundary condition at z=1, which causes a transfer of kinetic energy from the vertical to the horizontal components, an effect which is not accounted for by Eq. (31). A detailed comparison of the rotational dependence of and with Eq. (31), and with the analytical results of Rüdiger & Kitchatinov (1993), will be presented in Sect. 4.2.

 Figure 4: Solid: kinetic helicity, ; dotted: current helicity, , for a vertical imposed magnetic field at the solar south pole. ( a,c) average over ; ( b,d) average over . Open with DEXTER

The current helicity is negative in the upper part of the unstable layer for a vertical imposed magnetic field (Fig. 3, row 5), which confirms earlier results of Brandenburg et al. (1990). A reverse-polarity layer is not observed, except in run 10 with the impenetrable boundary at z=1. Thus, the sign of in the unstable layer as predicted by Eq. (30) is the same as that of . The negative sign of , which corresponds to a positive sign in the bulk of the convection zone on the northern hemisphere, is at odds with forced isotropic turbulence calculations by Brandenburg (2001) as well as with observations of the solar surface (Seehafer 1990; Low 1996) and mean-field calculations (Rädler & Seehafer 1990; Rüdiger et al. 2001). For a horizontal imposed magnetic field, the current helicity is about one order of magnitude smaller than in the case of a vertical imposed field. It is highly fluctuating, and its mean value can have either sign (not shown). The smaller amplitude is consistent with the smaller amplitude of the magnetic fluctuations (Fig. 6).

### 4 The effect as a function of

In this section we investigate the dependence of the effect on the angular velocity . The relevant dimensionless parameters are the Taylor number and the Coriolis number, as defined by (22) and (23). Averages are calculated over time, and separately over the unstable region and the stable region, as explained above.

#### 4.1 Kinetic helicity and current helicity

Figure 4 shows the rotational dependence of the kinetic helicity and the current helicity (for the vertical imposed magnetic field). The kinetic helicity and the current helicity have been multiplied with and , respectively, in order to allow a comparison with , as suggested by Eqs. (29) and (30) if one sets . The kinetic helicity in the unstable layer increases with the rotation rate up to ; its value in the stable layer is much smaller. Hence the general trend and the correct sign of (Fig. 5, bottom) are recovered by (29) in the unstable layer, but not in the stable layer. The current helicity for the vertical imposed magnetic field increases with rotation up to -300, while it decreases for stronger rotation, as do the magnetic fluctuations (Fig. 6). It seems that in the unstable layer roughly traces the current helicity (Figs. 4 and 5). Note however that both the normalized kinetic helicity, , and the normalized current helicity, , increase monotonically with (not shown). This reflects an increasing degree of alignment of vorticity with velocity, and of current with magnetic field, respectively.

The main conclusion to be drawn with regard to the validity of (29)-(31) is that qualitative agreement exists only between , (29), and (31) on the one hand, and between and (30) in the unstable layer on the other hand. For , such qualitative similarities are no longer evident, but in the following sub-section we shall compare the behavior of and with mean-field results that were derived for the case of arbitrary rates of rotation.

 Figure 5: -coefficients, normalized by , as functions of Ta. Here is an average over time and over the unstable layer. Top: vertical effect. Bottom: horizontal effect. The errors are standard deviations divided by the square root of the number of turnover times covered by the simulation. All simulations are for the solar south pole. The dashed curves are analytical results based on Rüdiger & Kichatinov (1993). In each figure, a suitable uniform scaling factor was applied to the analytical curves. Open with DEXTER

 Figure 6: Strength of the velocity and magnetic field perturbations as a function of Ta, in units of , and , respectively. The vertical and horizontal orientation of the imposed magnetic field are used for and , respectively. Open with DEXTER

#### 4.2 Rotational quenching of

Figure 5 shows and , normalized by , as a function of Ta (or Co). For zero rotation the effect vanishes. For weak rotation increases with increasing rotation rate. It reaches a maximum near ( ); for stronger rotation it is quenched. The decrease of occurs in spite of the monotonic increase of the kinetic helicity, and in spite of the decrease of with increasing rotation rate (Fig. 6). The horizontal effect sets in at a higher rotation rate. Its magnitude increases monotonically with the rotation rate, and exceeds that of the vertical effect if . For , Coriolis forces are so strong that they stifle the convection (see below), and hence also alpha. Two effects that are not shown in detail appear to be responsible for the reduction of : an increasing number of sign changes of , resulting in cancelations when the volume average is calculated, and a decreasing amplitude of the magnetic fluctuations (Fig. 6). Both effects are absent in the case of  .

Clearly, the alpha effect is highly anisotropic, and its dependence on rotation is more complicated than suggested by Eqs. (29) and (31). Rüdiger & Kichatinov (1993) present analytical expressions for the rotational dependence of and , based on the first-order smoothing approximation (FOSA), but for arbitrary rotation rates. They employ a linearized equation of motion for the fluctuating quantities which includes a random driving force with prescribed spectral properties to model turbulence, and which takes into account a stratification of the density and of the turbulent velocity. Due to the stratification as well as rotation, the electromotive force becomes anisotropic. In the simplest case of a driving force with zero frequency and containing a single wavelength (the case designated mixing-length approximation by the authors), and for a weak magnetic field, a manageable expression can be derived. For the present geometry, i.e., a local Cartesian grid situated at the south pole, one obtains

 = (32) = (33)

These expressions correspond to those of Rüdiger & Kichatinov (1993) if one sets , , , and (cf. their Eq. (3.11)). The functions are given by
 = (34) = (35)

where . Several features relevant for our discussion should be noted. As is shown in Fig. 7, the functions are almost equal and, for , they are linear in , so that Eq. (33) reduces to Eq. (31), apart from a factor of order unity. Secondly, , and therefore , are subject to rotational quenching. This unexpected behavior, which is not predicted by Eqs. (29)-(31), is roughly confirmed by the simulations, as is shown by Fig. 5. In order to produce the dashed curves, representative values, estimated from the simulations, were inserted for the gradients and , namely 1.5and -1.0 in the unstable layer, and 1.0 and -3.5 in the stable layer (Fig. 3, row 4).
 Figure 7: Functions , which describe the Coriolis number dependence of contributions to and due to density and turbulence stratification, respectively. Open with DEXTER

For a given Taylor number, the same value of was used for both the unstable layer and the stable layer, i.e. the turnover time, , was assumed to be the same. This is motivated by the fact that both the effective thickness, , and for the stable layer are smaller by a factor 3-4 compared to their values in the unstable layer. Secondly, the temporal variations of the flow are very similar in both layers (Fig. 2, bottom). After dividing Eqs. (32) and (33) by , a factor remains in front which can be assumed to be independent of . However, it was found that setting this factor to overestimates the amplitude of . For the dashed curves shown in Fig. 5, it was replaced by a suitable uniform scaling factor, namely 0.19 for and 0.01 for . Thus, our simulations produce much smaller coefficients than predicted by (32) and (33), especially which is smaller by two orders of magnitude. On the other hand, the main point is the agreement with regard to the rotational dependence. Except for the sign of in the unstable layer, the rotational quenching is reproduced, including the position of the maxima. For weak rotation, , so that the contribution of dominates, while for the density gradient dominates; in both cases the mean-field formula yields the wrong negative sign. Clearly, the sign of as predicted by Eq. (32) can have both values in principle, and depends in a delicate way on the relative importance of the two gradients; even a slight change of the coefficients can generate a sign chance. With regard to , we confirm the saturation predicted by Eq. (33) only in the stably stratified layer, but not (yet) in the unstable layer (Fig. 5, bottom). Apart from this, the general trend as well as the correct signs are reproduced by (33) in both layers. Due to the shape of the functions , this feature should be rather robust.

 run 1 2 3 4 5 6 7 8 9 2 10 30 100 300 2000 104 eddy size 2.5 2.5 2.5 2.5 2 1.5 1.3 1.0 0.8 2.8 2.8 2.7 2.4 2.1 1.5 1.1 1.0 0.82

 Figure 8: Snapshot of the vertical velocity (connected lanes: downward, isolated patches: upward). Top: run 6, moderate rotation. Bottom: run 9, strong rotation. The white curve denotes the zero level of the vertical velocity. The top and bottom surfaces correspond to z=0.09 and z=1.81, respectively. Open with DEXTER

In order to shed further light on the influence of rotation, Fig. 8 shows snapshots of the vertical velocity for two simulations with different rotation rates. Clearly, stronger rotation causes the eddies to have a smaller diameter in the horizontal plane (i.e. perpendicular to the rotation axis). Hence the number of eddies within the simulation domain increases with the rotation rate, and this leads to a more accurate result for , as is obvious from the length of the error bars in Fig. 5, taking into account the duration of the runs (Table 1). In the simpler case of Boussinesq convection in a rotating layer, the scale at which the instability sets in also decreases with increasing Taylor number (Chandrasekhar 1961). Although the present simulations are done at supercritical Rayleigh numbers whereas the analysis of Chandrasekhar applies to the marginally stable case, we may for the moment set aside this difference, as well as possible effects of compressibility and different boundary conditions. In fact, the smallest value of the local Rayleigh number in the unstable layer was found to be typically 104, so that the actual supercriticality is less than assumed. The present runs, which are characterized by and , belong to a regime where the onset of instability occurs in the form of non-oscillatory convection. Since the value of the critical Rayleigh number depends on the horizontal wavenumber, k, there is a critical wavenumber, , corresponding to a preferred eddy size , for which the critical Rayleigh number attains a minimal value, . The eddy sizes observed in the simulations are in quite good agreement with the analytical values of for convection between two free boundary surfaces (Table 2). The results suggest that the preferred mode of convection is the same as that pertaining to the onset of convection. The critical Rayleigh number increases with from 661 for up to for . This explains why for a fixed value of the intensity of convection, as measured by , decreases with increasing rotation rate. Since the actual Rayleigh number in the unstable layer has a local minimum of about 104, this also explains why no convection was observed for .

 Figure 9: Energy fluxes. Top: run 6; Bottom: run 9. Only the three most important contributions are shown, the viscous and electromagnetic fluxes being negligible. The total flux does not include the flux associated with the cooling term Q, Eq. (9); hence the steep rise in region 1. Open with DEXTER

As the scale of the individual eddies decreases, the corresponding decrease of the effective Coriolis number, , may be interpreted as a readjustment of the convection in response to increased Coriolis forces, which tend to hamper the convective energy transport. Still the convective flux of the run with the strongest rotation (run 9) is reduced by a factor of about two compared to runs with weak rotation. Apart from a possible effect of the increased Coriolis force, the reduced flux may result from the fact that smaller eddies lead to larger temperature gradients in horizontal planes, and this causes increased horizontal radiative transport between adjacent eddies which reduces the efficiency of convective energy transport in the vertical direction (Fig. 9).

### 5 Alpha and magnetic field strength

If the strength of the imposed magnetic field is increased, effects of the Lorentz force become more noticeable. Moreover, these effects depend on the orientation of the imposed field, i.e. vertical or horizontal. From the previously discussed runs, a representative case (run 7) was selected as the starting point for runs with increasingly strong imposed magnetic fields of both orientations (Table 3). A comprehensive study of the mode structure of magnetoconvection for different imposed magnetic fields is beyond the scope of the present paper, though. Various cases for a vertical imposed field, without rotation, are treated by Matthews et al. (1995). Wissink et al. (2000) consider the breakup of a horizontal magnetic layer into a number of flux tubes, thereby including the effect of rotation.

 Figure 10: Top: rms velocity perturbations as a function of the strength of the imposed magnetic field (normalized to the value of for small B0). Bottom: rms magnetic field and final equipartition field strength. Open with DEXTER

 run 7 7a 7b 7c 7d 7e B0 0.001 0.01 0.02 0.04 0.06 0.08 T ( ) 940 205 203 200 414 96 T ( ) 508 280 278 559 381 268

 Figure 11: Top: and as functions of for run 7, where . Bottom: and as functions of for run 7, where . Here is the equipartition field strength, averaged over the unstable layer, for convection without a magnetic field. Analytical quenching results ( dash-dotted) are given for both signs. All calculations are for the solar south pole. Open with DEXTER

For an imposed magnetic field in the vertical direction, convection is increasingly hampered with increasing field strength, as is indicated by a decrease of the rms turbulent velocity (Fig. 10, top) and of the convective energy flux. Since the magnetic energy is enhanced by advection of the imposed field, while the kinetic energy decreases as a result of the Lorentz force, after a while the rms magnetic field can be much larger than the equipartition field strength (Fig. 10, bottom). Moreover, if the initial magnetic field is too strong, then convection dies out completely. This in fact occured during run 7e ( ). After 100 time units, the volume-averaged energy density of the velocity and magnetic-field perturbations had both fallen by about 2 orders of magnitude. The decrease of is more rapid than that of , because advection of the imposed field by convection is stifled; this explains why alpha decreases also when measured in dynamical units (i.e., when divided by ).

A debate is ongoing in the literature about how and whether the magnetic Reynolds number affects the magnetic quenching of the effect. The dependence of on the magnetic field strength is often schematically represented in the form . While some numerical and analytical results suggest that (Kraichnan 1979a, 1979b; Brandenburg & Donner 1997), others suggest (Cattaneo & Vainshtein 1991; Vainshtein & Cattaneo 1992; Tao et al. 1993; Cattaneo 1994; Cattaneo & Hughes 1996; Vainshtein 1998; Brandenburg 2001). Although the present set of simulations is not well-suited to answer this question because of the relatively small Reynolds numbers ( ), the results for both and are consistent with (Fig. 11). Blackman & Field (2000) argued that the quenching with observed in simulations with periodic boundaries may in fact be a consequence of the boundary conditions rather than a dynamic effect, and that p=0 should be possible if a flux of magnetic helicity through the boundaries is allowed. Our conditions (10) do allow for a flux of magnetic helicity through the top and bottom boundaries. Thus, the fact that quenching sets in when suggests that an Rm-dependence of the quenching is possible in the present case; such quenching could be dynamic according to the reasoning of Blackman and Field. However, this conclusion may change if one allows for shear near the surface. This is the case between the solar surface and the corona, and so a strong helicity flux is possible in that case (Berger & Ruzmaikin 2000).

It appears that for boundaries that allow for a flux of helicity diverse Rm dependencies of the quenching are possible, depending, e.g., on the character of the forcing of the flow. Brandenburg & Dobler (2001) report a case where . We also note that even with the strong -quenching of a large-scale field may build up because of an equally strong quenching of the turbulent magnetic diffusivity (Brandenburg 2001).

The magnetic-field dependence of the coefficients is essentially the same for the cases of the vertical and horizontal imposed fields, even though convection is hampered less in the case of a horizontal imposed field (Fig. 11, top). For field strengths comparable to the equipartition value, the convective pattern becomes very different. If , convection exhibits a pattern of irregular oblique elongated cells that make an angle with the direction of the imposed magnetic field, and a small horizontal effect is still measured. For runs with field strengths and respectively, longitudinal rolls are formed that are increasingly aligned along the x-direction, and the coefficient is practically zero. For still stronger fields ( ), convection assumes the form of highly regular oblique lanes, and a small effect is again measured.

### 6 Discussion

The following discussion of our results in the solar context may be considered as complimentary to that of Brandenburg et al. (1990). First-order smoothing is a justified approximation if either the magnetic Reynolds number is small, or if the flow is slow in the sense . In our dimensionless variables, the latter condition would mean but, with the rms velocities and coherence times obtained in our simulations, we rather have . Our magnetic Reynolds number (20-35) is quite moderate but not small. Hence the simulations clearly go beyond first-order smoothing, although we are far from the conditions met in the solar convection zone. The large fluctuations that occur at large Reynolds number are evident in our calculations. In particular we see that the components of the tensor undergo large fluctuations. Only after averaging over many coherence times and over the horizontal coordinates of the box a significant depth variation is obtained. With this variation we confirm some results obtained in the case of weakly anisotropic turbulence, and extend them beyond the limits of first-order smoothing: in the southern hemisphere the horizontal coefficient (Fig. 3, row 2) is negative (positive) in the upper (lower) part of the box, as is plausible from the rotational effect on convergent and divergent velocities that arise as the flow is forced into the horizontal direction. On the other hand we find the opposite sign for the vertical coefficient (Figs. 2 and 3, row 1). Also of opposite sign is the small-scale current helicity. Indeed, both the sign and the rotational dependence of in the unstable layer roughly agree with that predicted by the simple FOSA-expression which relates to the current helicity (30).

For cases with convective overshooting in the lower part of the box the components of reverse their sign near the transition to the stable stratification. For an dynamo the horizontal component is the essential ingredient; possibly the sign change of has interesting implications for a dynamo that operates as an interface wave at the base of the convection zone. In the layer of convective overshooting the sign of is appropriate for the equatorward migration of the mean field, if the positive shear , as inferred from helioseismology for the depth range around and the lower heliographic latitudes (Kosovichev et al. 1997), is used as the other essential ingredient. These conclusions are drawn from simulations with an imposed magnetic field B0=10-3, and may still be valid for B0=10-2, as the results of run 7a, shown in Fig. 11, suggest. The latter field strength corresponds to 10 T, which is the strength of the field in the stable layer below the convection zone, as deduced from the properties of rising flux tubes (Caligari et al. 1995). For stronger fields we obtain a marked quenching of consistent with a factor , although we have checked this only for . The simulations have also revealed a striking difference between the convective patterns for a vertical and a horizontal orientation of the imposed magnetic field, with rolls being formed for strong horizontal fields. In view of a strong toroidal field at the base of the convection zone a magnetic Rayleigh-Taylor instability (Brandenburg & Schmitt 1998) or an instability of toroidal flux tubes (Ferriz-Mas et al. 1994) should be employed to obtain an effect. Nevertheless, since the sign of follows plausible arguments, we may hope to obtain the same sign for such a dynamically determined (Brandenburg & Schmitt 1998).

As far as the dependence of on rotation is concerned, the simulation may yield more reliable results than for the other parameters. In the solar convection zone the Coriolis number is small near the surface, and increases to values around 1 or somewhat larger near its base. This is the regime shown in Fig. 5. The rotational quenching of and the saturation of (Fig. 5, top) may therefore just begin in the depth region where the dynamo operates, but in solar-type stars which are younger and therefore rotating faster, these effects may be important. We have observed some agreement between the numerical results for and the predictions of Rüdiger & Kichatinov (1993), as far as the dependence on rotation and on depth is concerned. In both the unstable and the stable layer, the approximate Coriolis number dependence and the correct sign for are reproduced. The main discrepancies concern the amplitude of (the scaling factors of the theoretical curves in Fig. 5) and the sign of in the unstable layer, as well as the sign of the small-scale current helicity.

Even within the limited range accessible to numerical simulation, the parameter space is enormous. We have not yet considered the latitude dependence of (although some earlier work of Brandenburg 1994 indicates that is maximum at about of latitude), and we have not yet explored the dependencies on the Reynolds and Rayleigh numbers. This must be done by further simulations.