A&A 365, 562-570 (2001)
S. B. F. Dorch1 - Å. Nordlund2
Send offprint request: S. B. F. Dorch
1 - The Royal Swedish Academy of Sciences, Stockholm Observatory, 13336 Saltsjöbaden, Sweden
2 - The Niels Bohr Institute for Astronomy, Physics and Geophysics, Juliane Maries Vej 30, 2100 Copenhagen Ø, Denmark
Received 23 June 2000 / Accepted 26 September 2000
The interaction of magnetic fields and stratified convection was studied in the context of the solar and late type stellar dynamos by using numerical 3D MHD simulations. The topology of stratified asymmetric and over-turning convection enables a pumping mechanism that may render the magnetic flux storage problem obsolete. The inclusion of open boundary conditions leads to a considerable flux loss unless the open boundary is placed close to the physical boundary. Simulations including solar-like latitudinal shear indicates that a toroidal field of several tens of kilo-Gauss may be held down by the pumping mechanism.
Key words: Sun: MHD - magnetic fields - stratified convection - activity-rotation: stars - magnetic fields
Author for correspondance: email@example.com
Magnetic fields play an important rôle for the formation of the spectra of active late type G, K and M dwarf stars (e.g. Schrijver & Rutten 1987; Rutten et al. 1989; Johns-Krull & Valenti 1996). The most studied star in that respect is the Sun. It is generally assumed that solar active regions are manifestations of a strong toroidal flux system, that is generated and anchored deep below the surface of the Sun, possibly in the undershoot layer below the convection zone proper. Toroidal magnetic strands ascend through the convection zone because of buoyancy and the on the average super-adiabatic stratification. The much weaker poloidal field is assumed to be generated in the convection zone from this toroidal field by a cyclonic effect.
In the mean field dynamo context, the mechanism that takes care of communicating the poloidal field back to the region where the toroidal field is generated, is assumed to be a diffusive coupling of the regions where the -effect and the cyclonic -effect operate (e.g. Parker 1993). Choudhuri & Dikpati (1995) have shown that meridional circulation may also couple the two regions, if the time scale of the circulation is shorter than the diffusive time scale, and if the circulation is such that the flow is equator-ward at the bottom of the convection zone. The modeling by diffusive coupling of the regions may be considered to be somewhat unsatisfactory because it relies on a rather ad hoc approach and the meridional circulation approach is equally unsafe since it is not well observed -- at best the amplitudes and directions of the flow are indicated (van Ballegooijen 1998).
The question of how the poloidal field may return to the region where the generation of the toroidal field supposedly takes place is related to the well-known ``buoyancy dilemma'' and the generally assumed solution to this problem: any magnetic field in the convection zone will escape because of its buoyancy but a magnetic field may be stored in the stably stratified region below the convection zone; i.e., in the undershoot layer. In that scenario the convection zone is considered to be a passive one-dimensional medium while the magnetic fields are treated more or less as solid objects that move in it. It is conceivable, however, that some kind of balance may occur between the drag of descending plasma and the buoyancy of the magnetic field embedded in the plasma.
Along this line of thinking Drobyshevski & Yuferev (1974) proposed that a downward ``topological pumping'' of the magnetic field could be occurring, because of the asymmetric nature of the topology of 3D convective flows, i.e. that they consist of networks of descending material embedding regions of ascending material. They investigated the kinematic case at low magnetic Reynolds number by assuming an incompressible simple geometrical velocity flow pattern. Criticism by Parker (1975) made Drobyshevski et al. (1980) redo the experiments with a more suitable upper boundary condition. Arter et al. (1982), Arter (1983) and Galloway & Proctor (1983) extended the work to higher magnetic Reynolds numbers and several different compressible flows. They found that the magnetic energy did indeed increase at the bottom of the domain.
During the last decade several groups have performed more detailed magneto-convective numerical simulations (e.g. Hurlburt & Toomre 1988; Brandenburg et al. 1990; Jennings et al. 1992; Nordlund et al. 1992; Hurlburt et al. 1994; Nordlund et al. 1994; Brandenburg et al. 1996). Moreover, numerical simulations of stratified convection have shown that trace particles initially placed in a horizontal layer of a highly stratified model on the average are transported downwards, as a result of the asymmetric topology of stratified convection (Stein & Nordlund 1989). In the context of magnetic fields such a tendency for downwards transport of magnetic fields has also been seen in convective dynamo simulations (Nordlund et al. 1992; Brandenburg et al. 1996). Most recently Tobias et al. (1998) have investigated this effect using their standard method for studying stellar convection and Mcleod (1998) presented speculations along the same line of thought that shall be followed here.
The objective of the numerical experiments presented here is to study the interaction of magnetic fields and solar-like stratified over-turning convection and differential rotation. The model of the Sun is a ``local box'' model of a convectively unstable layer (henceforth referred to as the ``convection zone'') sandwiched between two stable layers. In order to circumvent problems associated with the very disparate thermal and dynamical time scales, the model has a much higher luminosity than the Sun, and all variables are scaled accordingly. To compare with solar values, the results must be re-scaled as follows. With a flux scale of times the solar flux, the velocity scale factor becomes , the temperature fluctuation scale factor , and the magnetic field strength scale factor . For convenience, the variables in the model are given in units of time s, length Mm, and density gcm-3. A Kramers' opacity scaled inversely with the scale factor for the total luminosity is adopted. This ensures that the boundary between the stable layer and the convection zone is at about the same depth in the model as it is in the Sun.
The full resistive and compressible MHD-equations are solved
using the staggered mesh
method by Galsgaard and others (e.g. Galsgaard & Nordlund 1997;
Stein & Nordlund 1989;
Nordlund et al. 1992, 1994):
The code uses a finite difference staggered mesh with 6th order derivative operators, 5th order centering operators and a 3rd order time-stepping routine (Hyman 1979).
In most stars envelope convection is essentially driven
by surface cooling. The entropy contrast at the surface is far
larger than that at the bottom of the convection zone, if the convection
zone extends over several or many pressure scale heights. To model
this situation, without having to actually include all layers up to the
solar surface, a simple
expression for an isothermal cooling layer at the upper boundary of the model
Both experiments with closed and open upper boundaries were performed. In order to implement a stable open upper boundary a buffer zone was allocated, where a fiducial electric field is gradually turned on. The sense of the electric field is such that it drags the magnetic field out of the buffer zone, and the magnitude is increased from zero to of the order of , where is the maximum velocity in the buffer zone. A layer that is about 10 grid zones from the numerical upper boundary may thus be considered as the physical open upper boundary. This layer is far below the real boundary of the solar convection zone (even the numerical upper boundary of the model is far below the photosphere).
The hydrodynamic part of the initial condition is a snapshot from a well developed stage of a numerical model of the solar convection zone. The initial condition for the magnetic field is given by a unidirectional (``poloidal'') sheet that is placed in the middle of the convection zone. The sheet is initially in a state of isentropic pressure equilibrium with the surroundings.
The model has a high degree of stratification with a density contrast of roughly in the convection zone alone. The thickness of the undershoot layer in the model is approximately equal to 0.8 , where is the pressure scale height at the bottom of the convection zone. This is much larger than the helioseismological upper limits of 0.1-0.2 (see e.g. Christensen-Dalsgaard 1995). Scaling the extension of the undershoot layer naively in proportion to the velocity (Hooke's law -- the retarding force increases approximately linearly with distance), the resulting undershoot thickness, of the order of 0.01 , falls well within the upper limits from observations.
Both experiments with and without differential rotation
in the convection zone were performed, with
the purpose of
illustrating the effects of radial and latitudinal shear and of the
Coriolis force on the magnetic field.
In the rotational cases background differential rotation and the Coriolis force
are included by using the force
in the equation
of motion Eq. (2)
The differential rotation that is implemented in the model may not be
identified as rotation around e.g. one of the horizontal axes.
Rather the mapping of the rotation is such that the rotation axis is
given by a vector in the ``meridional'' (x,y) plane pointing towards the
north pole, i.e. a rotation vector
An azimuthal background
velocity corresponding to the observed rotation is implemented by adding
a Newtonian term in the equation of motion:
The differential rotation is then measured relative to the latitude
In the model that latitude corresponds to the position
at fractions of 0.25, 0.75, 1.25 and 1.75 of the latitudinal extension
computational box counting the total length as 2 Ly(see Fig. 2).
At that latitude there is no radial shear since we set the internal solid
body rotation equal to the in situ azimuthal speed, so that
The strength of the background rotation may also be changed by varying the radius R in Eq. (12). Values from 500 to 1400 Mm were used, the higher values in combination with larger values of , in order to still produce the desired amplitude of the differential rotation.
Several experiments with varying initial magnetic field strengths, numerical resolutions, and upper boundary conditions were performed. First, results from simulations without rotational effects are discussed (see Table 1), and then results from a simulation including latitudinal shear are reported. Although some of the experiments initially had a maximum field strength of the order of the formal equipartition value, they have been denoted ``sub-equipartition'' in Table 1 because the distribution of the kinetic energy density was much broader than the distribution of the magnetic energy density , i.e. the peak magnetic energy density was smaller than the peak kinetic energy density even though the most likely value of was similar to that of .
A moderate number of grid points were used in these experiments (a few experiments were also ran at higher resolutions, up to , to check for resolution effects). The advantage of performing relatively small numerical experiments, is that it is possible to perform a larger number of experiments, with different setups and with a variety of parameter values. Since it is impossible, with the limits of currently available computer power to accurately reproduce solar conditions, it is necessary to experiment with trade-offs between various constraints.
|Experiment||Initial field||Plasma||Grid points||Top boundary|
Figure 1 shows two sets of magnetic field isosurfaces at two
different instants in time for an experiment with a sub-equipartition magnetic
field (experiment [i], in Table 1):
The poloidal sheet of magnetic field initially
placed in the middle of the convection zone quickly starts to interact with
the convection (Fig. 1, left panel).
At a subsequent
time (Fig. 1, right panel) the magnetic field
more or less fills the whole volume of the convection zone and
penetrates into the stable layer below.
|Figure 1: Snapshot of magnetic field isosurfaces: light gray (yellow), low field strength -- dark grey (blue), 10 times higher field strength (colors refer to color versions of the figure, available at http://www.astro.ku.dk/aake). Left: early in the experiment. Right: at a late time|
|Open with DEXTER|
|Figure 2: The background rotational velocity Eq. (12) (in kms-1) as a function of latitude (in units of ); the two equators of the model are at the latitudes 0.5 and 1.5 (corresponding to and , respectively), while the poles are at 0, 2 (north pole) and 1.5 (south pole)|
|Open with DEXTER|
Figure 3 shows the horizontally averaged field at 7 equidistant instants of time for experiment [i]: the poloidal sheet is spread out, and the distribution of poloidal magnetic flux settles to a characteristic distribution. The highest (horizontally averaged) poloidal flux density occurs in the overshoot layer, and in these particular experiments a significant fraction of the total flux also resides there. In the real Sun, this fraction may be expected to be much smaller, because the undershoot layer is much thinner there. Given the shape of the distribution, with no particular enhancement in the undershoot layer, it is likely the Sun has a correspondingly smooth distribution, with the majority of the poloidal flux residing inside the convection zone.
The ``pumping effect'' described in the above takes place because of the topology of the over-turning stratified convection: of all the fluid parcels threaded by magnetic field lines in the initial state about half are initially ascending. However, because of the stratification, most of the ascending fluid parcels have to over-turn and descend and most of these keep descending down to the bottom of the convection zone (Stein & Nordlund 1989). The fluid parcels drag the threading field lines along and hence an appreciable fraction of the field is transported downwards. Fragments of the field that are caught in ascending flows are advected upwards. In the experiments with an open upper boundary some of these fragments escape through the top of the computational domain and flux is systematically lost (see the discussion below).
Even in the cases of super-equipartition fields (e.g. experiments [iv] and [v]) the magnetic field is pushed downwards in the initial phases until it reaches the undershoot layer where the kinetic energy flux decreases. Because the overall magnetic flux decreases through flux loss at the surface the magnetic field eventually enters a state in which it is below equipartition.
In the case of an initial field below equipartition the pumping effect is -- not surprisingly -- even more pronounced than for the super-equipartition fields: Fig. 4 shows the results for three cases where the magnetic field was below equipartition (experiments [i], [ii] and [iii]) and the dynamics thus were dominated by the convective motions. The average magnetic field is distributed over the entire convection zone, with maximum (horizontally averaged) flux density in the undershoot layer.
The radial distribution of the magnetic flux present
in the convection zone
is not very different for the cases
with a closed upper boundary (left panel in Fig. 4) and the cases with an open upper boundary
(right panel in Fig. 4), but the total amount of magnetic flux
is rapidly reduced during the
first few turn-over times, when there is a significant loss of
poloidal flux through the upper boundary.
This does not influence the distribution of flux over depth much, however,
it only influences the amount of flux that is available for distribution.
|Figure 3: The horizontally averaged magnetic field as a function of radius for 7 different times (experiment [i]): t= 0 (initial field - solid curve), 3.3, 6.5, 9.8, 13.0, 16.3 and 19.5 (thin dashed curves), and time-average (solid curve) turn-over times. The approximate bottom and top of the undershoot layer are indicated by two vertical dashed lines|
|Open with DEXTER|
The flux loss in the models with an open upper boundary is strongly exaggerated in comparison to the Sun. The real flux loss may be expected to be significantly smaller than in the models with open upper boundaries discussed here: Much of the weak ascending flux must over-turn rather than reach the solar surface since it is embedded in a fluid of which only a tiny fraction (0.1% estimated from mass flux amplitudes) reaches the solar surface. Note that the magnetic field considered is weak and incoherent and does not have sufficient buoyancy to overcome the drag of the fluid motions.
A supplementary numerical experiment confirms the exaggeration
of the flux loss: moving the open magnetic boundary four mesh
point upwards, corresponding to about 0.6 density scale heights,
reduces the flux loss by about 20%. The new boundary (at
gcm-3) is still over ten density
scale heights (sic) away from the real solar boundary (at
gcm-3). Almost all of the flux that is still lost in the
modified experiment may thus in the case of the real Sun be expected
to turn over before reaching the solar surface.
|Figure 4: Horizontal average of the poloidal By field for two experiments (left panel: solid line, experiment [i] and dashed line, experiment [ii]) that had a closed upper boundary (after 25 and 30 turn-over times respectively), and for experiment [iii] (right panel) that had an open upper boundary (after 18 turn-over times). The amplitudes of the curves have been scaled to have the same maximum as experiment [i] to compare to experiment [i]|
|Open with DEXTER|
Figure 6 shows the average poloidal field strength for the experiments with an open upper boundary (experiments [iii], [iv] and [v]). The relative reduction in the poloidal flux by the escape of magnetic structures through the upper boundary does not depend much on the field strength, for strengths up to several times formal equipartition values.
The rate of flux loss does depend on the field strength but in a rather counterintuitive way: The stronger the initial field the smaller the initial reduction rate. Naively one might expect that if the field strength is large the higher buoyancy would make the field escape faster. The explanation is that the over-turning of the fluid is slowed down in the cases with a stronger initial poloidal flux sheet. The fragmentation of the sheet is also slower in the case with a stronger field.
The magnetic field in the convection zone rapidly becomes very fragmented: while the field in the beginning is uniformly distributed in the sheet, quickly a picture develops where the field in the convection zone becomes very intermittent, while the field that is pumped into the undershoot layer is much more uniform.
The reason for this difference in topology between the magnetic field in the convection zone and in the undershoot layer is that once the magnetic field has been pumped down into the undershoot layer it is less susceptible to fragmentation, since the motions in the stable layer have a much smaller amplitude and are not as systematic as the motions in the convection zone.
The degree of intermittency and fragmentation of the field in the convection zone depends on the strength of the initial poloidal magnetic sheet. The poloidal magnetic sheet is more stable towards the initial fragmentation if its field is stronger since in that case a larger force is needed to over-power the tension in the poloidal field lines.
In what follows, results are presented for a particular simulation that includes shear and the Coriolis force (through the force in the equation of motion Eq. (2)). The physical size of this experiment was Mm and the numerical resolution grid points.
The initial condition of a uniform poloidal sheet placed at a certain depth may be said to be rather arbitrary. However, an arbitrary state is both appropriate and useful in this case, since the system rather quickly enters a state that is independent of the initial condition i.e. there is no internal long-time ``memory'' in the system.
Initially the pumping effect ensures that an on the average poloidal field is injected into the bottom undershoot layer, and after about four turn-over times the poloidal flux density has a maximum below the convection zone of the model.
After an initial transient process where the differential rotation establishes itself and the over-turning convection distorts the poloidal sheet, the system enters a well-developed state, where the magnetic field displays the structure of the differential rotation. At this point, the magnetic energy and the toroidal field strength have increased to a significantly higher level than their initial values (see the two top panels in Fig. 5), and the rms toroidal field is larger than the rms poloidal field. The average poloidal flux has dropped to about half the initial value due to flux loss at the surface (see the bottom two panels in Fig. 5) but this loss of magnetic flux through the upper boundary is somewhat halted after the formation of a large-scale toroidal field.
At first the growth of the toroidal field is fast while the
background rotation increases toward the profile determined by Eq. (12) on the time scale
(cf. Eq. (11)).
background rotation profile is fully attained by the fluid, the
toroidal magnetic field begins to increase linearly, with the rate of
increase given by the latitudinal shear and the original poloidal
The structure of the magnetic field directly
reflects the latitudinal dependence of
the background azimuthal velocity field
as a result of the latitudinal shear.
|Figure 5: The case with latitudinal shear: four panels showing the average magnetic energy , the rms field strengths Bz (toroidal field -- solid line) and By (poloidal field -- dashed line), the ratio of toroidal to poloidal field, and the average poloidal field as functions of turn-over times|
|Open with DEXTER|
In the well-developed state after the main flux loss has taken place, the dominant magnetic field is a strong toroidal flux system located near the bottom of the convection zone, i.e. it is not pumped into the undershoot layer but it does not escape from the convection zone either: the ``center of gravity'' of the magnetic field is above the bottom of the convection zone (see Fig. 7).
Figure 8 shows a view of the toroidal flux system: four toroidal flux streets are formed on each side of the equators, in the regions where there is a maximum shear.
The magnetic field just below the convection zone rotates more slowly than the field just above. This is a mechanism that may twist the magnetic field lines that connect across the equator regions.
As the initial poloidal sheet fragments while being wound up, toroidal flux structures leave the convection zone through the open upper boundary. In a rotating frame of reference vertical motions lead to horizontal motions through the action of the Coriolis force. One special case is the rise of flux loops that are rotated so that they emerge tilted with respect to the latitudinal circle, and another is the excitation of meridional circulation, as a result of the transport of angular momentum: in the simulation both of these effects are found.
The surface of the Sun is an extremely important place: not only may flux be lost there, but the surface constitutes a ``reconnection region'' where vertical field lines are effectively ``cut over'', and where the remaining ``stubs'' may be advected passively by the horizontal surface motions. Field lines of course actually continue into the corona and either connect back into the photosphere or connect out into interplanetary space. However, reconnection is observed to proceed so rapidly in the corona and the mass density is so low there, that from the point of view of the sub-surface dynamics the connections above the surface are of little importance.
Figure 9 shows two panels with a small section of the box, with snapshots of the magnetic field, looking down along the two toroidal flux streets. A magnetic flux structure is sticking out through the surface, and is moving towards the equator between the two toroidal flux streets, but at the same time it is connected through a subsurface structure to one of them. As a result of the drift of the surface-structure, the magnetic field lines that make up the subsurface structure have a tendency to become inclined with respect to the toroidal street to which it connect, and hence a poloidal field component is formed, with opposite polarity relative to the original one.
|Figure 6: The normalized average poloidal field strength as function of time in turn-over times for the experiments with an open upper boundary condition: (experiment [iii] dashed curve, experiment [iv] solid curve and experiment [v] dashed-dotted curve)|
|Open with DEXTER|
|Figure 7: The depth of the location of the maximum magnetic energy as a function of turn-over times|
|Open with DEXTER|
The main conclusions that may be drawn from the study of the interaction of a magnetic sheet with stratified over-turning convection in the absence of rotation are the following:
The toroidal magnetic field in these simulations reach peak field strengths of a several tens of kG (scaled to the Sun). The peak field strengths occur near the bottom, but still inside, the convection zone. It is conceivable that emerging flux regions form when buoyancy finally becomes dominant, and that this occurs at field strengths of the order of 100 kG, as has been deduced from emergence and tilt patterns by several investigations. Because of numerical limitations we were not able to study that process with the current series of experiments--future experiments with higher resolution and even larger density contrasts are needed here.
Besides offering a clue to the operation of the solar dynamo the results
presented here may also contribute to the understanding of magnetic field
generation in other late type stars. For example; the fact that no
undershoot storage is available in the case of
magnetically active fully convective M dwarf stars
(e.g. Chabrier & Baraffe 1997)
is generally considered to be a problem (see e.g. Allard et al. 1997; Küker & Rüdiger 1999) -- this ceases
to be a problem in the scenario presented here.
|Figure 8: A snapshot of strong magnetic field lines (blue) and magnetic field isosurfaces (yellow)|
|Open with DEXTER|
We note that, even though the loss of magnetic flux at the upper boundary is
exaggerated in the models presented here, such a loss is certainly a real
effect that is important to include, since the Sun is known to loose a
considerable amount of toroidal magnetic flux during an activity cycle
(e.g. van Ballegooijen 1998).
|Figure 9: A view of one quarter of the box, centered at one of the equators, where a magnetic structure is emerging through the surface and is moving towards the equator. Field lines in white and grey (purple) isosurfaces of magnetic field strength. The snapshot to the right is taken half a turnover time later than the one on the left|
|Open with DEXTER|
Lastly we find it important to emphasize that the visualization of ``emerging'' flux structures and field lines (Fig. 9) illustrates a mechanism first pointed out by van Ballegooijen (1995, private communication) that may be crucial for the reversal of the poloidal field: flux structures that rise and penetrate the surface effectively results in the field lines being ``cut'' at the surface, with the leading polarity tending to drift towards the equator, and the following polarity tending to drift towards the pole. The result is that subsurface connections between the following polarity of one emerging structure and the leading polarity of another trailing structure may become tilted in the sense opposite to the tilt associated with the normal winding of the field. When such reversed tilts are caught by the differential rotation, they will effectively lead to ``unwinding'' and reversal of the poloidal field component. This again is an important topic for future studies.
This work was supported in part by the Danish Research Foundation, through its establishment of the Theoretical Astrophysics Center. Computing time at the UNIC computing center was provided by the Danish Natural Science Research Council. SBFD acknowledges support through an EC-TMR grant to the European Solar Magnetometry Network.