A&A 423, 1101-1107 (2004)
DOI: 10.1051/0004-6361:20040435
S. B. F. Dorch
The Niels Bohr Institute for Astronomy, Physics and Geophysics, Juliane Maries Vej 30, 2100 Copenhagen Ø, Denmark
Received 13 March 2004 / Accepted 10 May 2004
Abstract
Evidence is presented from numerical
magneto-hydrodynamical simulations for the existence of magnetic
activity in late-type giant stars. A red supergiant with stellar
parameters similar to that of Betelgeuse ( Orionis) is
modeled as a "star-in-a-box'' with the high-order "Pencil Code''.
Both linear kinematic and non-linear saturated dynamo action are
possible: the non-linear magnetic field saturates at a
super-equipartition value corresponding to surface magnetic field
of field strengths up to 500 Gauss. In the linear regime two
different modes of dynamo action are found with exponential growth
rates of 4 and 25 years, respectively. It is speculated that
magnetic activity of late-type giants may influence dust and wind
formation and possibly lead to the heating of the outer
atmospheres of these stars.
Key words: stars: AGB and post-AGB - stars: late-type - stars: activity - stars: individual: Betelgeuse - stars: magnetic fields - magnetohydrodynamic: (MHD)
There are indications from both dynamo theory and observations that some late-type giant stars such as red supergiants and asymptotic-giant-branch stars (AGB stars) may harbor magnetic fields. On the theoretical side, it has been suggested that non-spherically symmetric planetary nebulae (PNe) formed during late stages of AGB star evolution may be a result of the collimating effect of a strong magnetic field: Blackman et al. (2001) studied interface dynamo models similar to mean field theory's solar -dynamo and found that the generated magnetic surface fields typically could be 400 Gauss, strong enough to shape bipolar outflows, producing bipolar PNe, while also braking the stellar core thereby explaining the slow rotation of many white dwarf stars. Also using mean field dynamo theory Soker & Zoabi (2002) propose instead an dynamo due to the slow rotation of AGB stars rendering the -effect ineffective. They find that the magnetic field may reach strengths of 100 Gauss, significantly less than that found by Blackman et al. (2001). On the one hand, they believe that the large-scale field is strong enough for the formation of magnetic cool spots (see also Soker & Kastner 2003, on AGB star flaring). These spots in turn may regulate dust formation, and hence the mass-loss rate, but the authors argue that they cannot explain the formation of non-spherical PNe (see also Soker 2002): on the other hand, the locally strong magnetic tension could enforce a coherent flow that may favor a maser process.
On the observational side of things, maser polarization is known to exist in circumstellar envelopes of AGB stars (e.g., Gray et al. 1999; Vlemmings et al. 2003, and recently Sivagnanam 2004) and X-ray emission has been observed from some cool giant stars (e.g., Hünsch et al. 1998; Ayres et al. 2003). These observations are generally taken as evidence for the existence of magnetic activity in late-type giant stars (cf. Soker & Kastner 2003).
The cool star Betelgeuse (a.k.a. Orionis) is an example of an abundantly observed late-type supergiant that displays irregular brightness variations interpreted as large-scale surface structures (e.g., Lim et al. 1998; Gray 2000). It is one of the stars with the largest apparent sizes on the sky - corresponding to a radius in the interval 600-800 . Freytag et al. (2002) performed detailed numerical 3-d radiation-hydrodynamic (RHD) simulations of the convective envelope of the star under realistic physical assumptions, while trying to determine if the star's known brightness fluctuations may be understood as convective motions within the star's atmosphere: the resulting models were largely successful in explaining the observations as a consequence of giant-cell convection on the stellar surface, very dissimilar to solar convection. Dorch & Freytag (2002) performed a kinematic dynamo analysis of the convective motions in the above model (i.e. not including the back-reaction of the Lorentz force on the flow) and found that a weak seed magnetic field could indeed be exponentially amplified by the giant-cell convection on a time-scale of about 25 years.
This paper reports on full non-linear magneto-hydrodynamical (MHD) numerical simulations of dynamo action in a late-type supergiant star with fundamental stellar parameters set equal to that of Betelgeuse. The paper is organized as follows: Sect. 2 contains a description of the numerical model, code and setup, Sect. 3 describes and discusses the results, the convective flows and dynamo action, and Sect. 4 contains a short summary and conclusion.
The 3-d MHD equations are solved for a fully convective star. This
is an example of a "star-in-a-box'' simulation, where the entire
star is contained within the computational box. The computer code
used is the "Pencil Code'' by Brandenburg & Dobler^{}. The code has
been employed in several astrophysical contexts including e.g. hydromagnetic turbulence (see Brandenburg & Dobler
2002; Dobler et al. 2003; Haugen et al. 2003). The numerical method
performs well on many different computer architectures especially
on MPI machines, and uses 6th-order spatial derivatives and a
2N-type 3rd-order Runge-Kutta scheme. The code has a
"convective star'' module that allows the solution of the
non-linear MHD equations by the numerical pencil scheme in a star
with a fixed radius R and mass M. The code solves the following
general form of the compressible MHD equations:
(1) | |||
(2) | |||
(3) | |||
(4) |
Variables are computed in terms of R and M so that e.g. the unit of the star's luminosity L becomes . In the present case these fundamental parameters are set to and yielding a luminosity of , consistent with current estimates of Betelgeuse's size, mass and luminosity. The model employs a fixed gravitational potential , an inner tiny heating core entering into Eq. (4) through , and an outer thin isothermally cooling spherical surface at r = R (corresponding to in Eq. (4)), with a Newtonian cooling time scale set to year corresponding to the typical convective turn-over time in the model of Freytag et al. (2002).
Due to constraints on computer time and resolution the true thermodynamic range of the star can not be represented and the surface is cooled at a temperature that is 4.5 times higher than Betelgeuse's effective temperature which is about K (various estimates exist in the literature, see e.g. Freytag et al. 2002). The fixed gravity is approximately given by a 1/r-potential (as in Freytag et al. 2002) and initially the star expands resulting in a slight decrease of 1.4% of the mean mass density, subsequently it re-contracts by roughly 0.1%. Betelgeuse is only slowly rotating and a rotational frequency was chosen corresponding to a surface rotational velocity of 5 km s^{-1}, yielding a large Rossby number. Although a strong differential rotation deep within the star can not be ruled out, this type of internal rotation is not included in the present model.
Models with different numerical resolutions have been run for testing: in this paper the results come from a model with 128^{3}uniformly distributed grid points yielding a spatial resolution of the physical size of the box being R^{3}. The models were computed at the Danish Center for Scientific Computing Horseshoe 512 cpu Linux cluster typically allocating between 16 and 32 cpus.
Boundary conditions on the computational box are anti-symmetric for components of the vector fields and in the direction of the component (yielding a vanishing value at the boundary) and symmetric across the boundary in the direction perpendicular to the component (yielding a vanishing gradient across the boundary). Boundary conditions for the density are anti-symmetric with respect to an arbitrary value across all boundaries (yielding a vanishing second derivative), and the boundary condition for the entropy corresponds to a constant temperature at the boundary.
Dynamo action by flows are often studied in the limit of
increasingly large magnetic Reynolds numbers Re
,
where
and U are characteristic length and
velocity scales. Most astrophysical systems are highly conducting
yielding small magnetic diffusivities ,
and their dimensions
are huge resulting in huge values of Re.
Betelgeuse is
not an exception and most parts of the star are better
conducting than the solar photosphere that has a magnetic
diffusivity of the order of
m^{2}/s:
Figure 1: Diffusivity: average radial magnetic diffusivity (m^{2}/s) in the model of Betelgeuse as a function of radial distance in solar radii . Adopted from Dorch & Freytag (2002). | |
Open with DEXTER |
Figure 1 shows the average Spitzer's resistivity as a function of radius in the model of Betelgeuse by Freytag et al. (2002). Spitzer's formula (e.g., Schrijver & Zwaan 2000) assumes complete ionization and hence the precise values of are uncertain in the outer parts of the star, where the atmosphere borders on neutral. There is some uncertainty connected also with defining the most important length scale of the system, but preliminarily taking to be 10% of the radial distance R from the center (a typical scale of the giant cells), and along the radial direction yields Re -10^{12} in the interior part of the star where .
In the present case we cannot use that large values of Re (partly due to the fact that we are employing a uniform fixed diffusivity ), but rely on the results from generic dynamo simulations indicating that results converge already at Reynolds numbers of a few hundred (e.g., Archontis et al. 2003a,b). Furthermore, Dorch & Freytag (2002) obtained kinematic dynamo action in their model of a magnetic Betelgeuse at Re . Here we use an that is about 10^{8} times too large compared to the estimated surface value in Betelgeuse, leading to a magnetic Reynolds number of Re .
There is some disagreement as to what one should require for a
system to be a "true'' astrophysical dynamo. Several ingredients
seem to be necessary: the flows must stretch, twist and fold the
magnetic field lines (e.g., Childress & Gilbert 1995);
reconnection at
must take place to render
the processes irreversible; weak magnetic field must be circulated
to the locations where flow can do work upon it (cf. Dorch
2000); and finally, the total volume magnetic energy
must
increase (the linear regime) or remain at a constant saturation
amplitude on a long time scale (the non-linear regime). These
points are based largely on experience from idealized kinematic
and non-linear MHD dynamo models; e.g. Archontis et al.
(2003a,b). This paper deals mainly
with the question of the exponential growth and saturation of
.
Figure 2: Simulated surface intensity snapshots at four different instants, time = 256, 347, 457 and 494 years (from upper left to lower right). | |
Open with DEXTER |
Although not the main topic of this paper, it is appropriate to discuss here also the properties of the convective flows in the model, since these ultimately supply the kinetic energy forming the basic energy reservoir for any dynamo action that might be present. It is not expected that the flows match exactly what is found in more realistic RHD simulations, but at least a qualitative agreement should be inferred since the fundamental parameters of this MHD model and the RHD model of Freytag et al. (2002) are the same.
The velocity is initialized with a random flow with a small
amplitude. Rapid large scale convection cells develop
throughout the star: the giant cell convection is evident in both
the thermodynamic variables, such as temperature and gas pressure,
as well as in the flow field. The observational equivalent
however, is the surface intensity. Since the model does not
incorporate realistic radiative transfer (as opposed to the model
of Freytag et al. 2002), only a simulated
intensity can be derived. Using a frequency independent LTE source
function a simulated intensity
can be defined as:
(5) | |||
(6) |
Figure 3: Powerspectrum of u^{2}(k) at time = 457 years, and a line corresponding to Kolmogorov scaling (dashed line). The vertical lines, from left to right, denote the wavenumbers corresponding to the computational box size k_{0}, the stellar radius , and the Nyquist frequency (the numerical resolution). | |
Open with DEXTER |
More quantitatively, the kinetic power-spectrum (Fig. 3) illustrates that there is much more power on large scales than on the small scales of the velocity field: below a wavenumber of (based on the box size in units of the star's radius) power is decreasing fast, but at larger scales the power is proportional to k^{-2/3} corresponding to normal Kolmogorov scaling, the inertial range spans however only roughly one order of magnitude. In conclusion the large-scale convective patterns are then typically larger than 15-30% of the radius, and are actually often on the order of the radius in size. The corresponding radial velocities range between 1-10 km s^{-1} in both up and down flowing regions.
There are at least three different evolutionary phases of convection in the simulations, depending on the level of the total kinetic energy of the convection motions: initially there is a transient of about 30 years after which the rms velocity field reaches a level where it fluctuates around a value of about 800 m/s (this corresponds to the kinematic phase of the dynamo, where the flow is unaffected by the presence of the still weak magnetic field, see below). During the rest of the simulation after about 290 years, the rms speed measured in the entire box decreases to 500 m/s (when the energy in the magnetic field becomes comparable to the kinetic energy density). During the stretch of the simulation however, the maximum speed in the computational box fluctuates around a constant value of about 90 km s^{-1}. The flows are not particularly helical and the mean kinetic helicity is on the order of 10^{-6} m/s^{2}. Mean field -type solar dynamos do not produce large-scale fields if the kinetic helicity is less than a certain value (cf. Maron & Blackman 2002) and hence we cannot expect a large-scale toroidal field in the solar sense to be generated.
In an earlier kinematic study of Betelgeuse using a completely
different numerical approach (Freytag et al. 2002;
Dorch & Freytag 2002), dynamo action was
obtained when the specified value of Re
was larger than
approximately 500 and at lower values of Re
the total
magnetic energy decayed. In the present case Re
is of
the same order of magnitude and we find an initial clear
exponential growth over several turn-over times, and many orders
of magnitude in energy. Figure 4 shows the evolution of
as a function of time, for the first 225 years (in
Betelgeusian time):
once the giant cell convection has properly begun the magnetic
field is amplified and we enter a linear regime of exponential
growth. There are two modes of amplification in the linear regime;
the initial mode with a growth rate of about 4 years, which in the
end gives way to a mode with a smaller growth rate corresponding
to a time-scale of 25 years. This is a slightly unusual situation,
since normally modes with smaller growth rates are overtaken by
modes with larger growth rates (cf. Dorch 2000); the explanation
is that while both modes are growing modes, only the one with the
largest growth rate is a purely kinematic mode - while the
exponential growth of the second mode is linear it is not
kinematic - the presence of the magnetic field is felt by the
fluid through the back-reaction of the Lorentz force in Eq. (2)
becoming important. This quenches the growth slightly and
henceforth one can refer to the second mode as a "pseudo-linear''
mode.
Figure 4: Linear regime: energy as a function of Betelgeusian time in years. The upper almost horizontal line is total thermal energy (dashed-dotted), middle curve with wiggles is total kinetic energy (dotted), and lower full curve is . The three thin dashed lines corresponding to exponential growth with characteristic time-scales of 3.8, 25 and years. | |
Open with DEXTER |
No exponential growth can go on forever and eventually the
magnetic energy amplification must come to a halt: the question is
then whether the magnetic field retains a more or less constant
saturation value, or if it dissipates. The latter is only possible
if the non-linear mode is a decaying mode that corresponds to a
negative growth rate. In case of saturation the typical field
strength is expected to be on the order of the equipartition value
corresponding to equal magnetic and kinetic energy densities.
Figure 5 shows the rms magnetic field strength
within the entire model star as a function of time for 800 Betelgeusian years: the pseudo-linear mode as well as the
mode in the non-linear regime are shown. The rms magnetic field
saturates at a value slightly above the rms equipartition field
strength
90-100 Gauss, corresponding to a value of about
120-130 Gauss. In terms of total energy this means that the
magnetic energy
is above equipartition with the
kinetic energy
by approximately a factor of two.
Hence the field cannot be said to be extremely strong, but it is
not particularly weak in most parts of the star either. What may
be interesting from an observational point of view is the strength
of the field at the surface. In the non-linear regime the field
strength at the sphere with radius r = R can be up to
500 Gauss, while in the interior going downwards in the
star the field strength rises and can be as high as a few kG: the
strongest intermittent magnetic structures almost completely
quenches the velocity field in these regions that are small-scale
compared to the scale of the convection; i.e. the local field can
be far above equipartition. The downward increase of the field
strength is analogue to the flux pumping effect that has been
found in the solar context (cf. Dorch & Nordlund
2001).
Figure 5: Transition to the non-linear regime: rms magnetic field B in the whole computational box in Gauss as a function of Betelgeusian time in years. The upper curve is the equipartition field strength corresponding to the average kinetic energy density of the fluid motions (dotted curve) and the lower full curve is the actual rms field strength B (full curve). The dashed thin curve correspond to growth times of 25 years and a horizontal reference line at a field strength of 100 Gauss. | |
Open with DEXTER |
Figure 6: Distribution of magnetic field strength: PDF for the magnetic field |B| (in kG) at time = 732 years within the non-linear saturated regime. | |
Open with DEXTER |
It is interesting to examine the geometry of the magnetic field
that the saturating non-linear dynamo generates since this could
be relevant for the influence of the field on e.g. asymmetric
dust and wind formation. Qualitatively speaking the field becomes
concentrated into elongated structures much thinner than the scale
of the giant convection cells, but perhaps due to the very
irregular nature of the convective flows, no "intergranular
network'' is formed in the solar sense. On the one hand, at times magnetic
structures coincide with downflows, but not as a general rule. On the
other hand, strong fields are seldomly located within the general
upflow regions.
Figure 7: Energy spectra of the magnetic energy density (solid curve) and kinetic energy density (dashed curve) at time = 512 years, and a line corresponding to a power-law with an exponent of -2/3 (dashed-dotted line). The vertical lines, from left to right, denote the wavenumbers corresponding to the computational box size k_{0}, the stellar radius , and the Nyquist frequency (the numerical resolution). | |
Open with DEXTER |
Figure 8: Magnetic filling factor as a function of radius (solid curve): shown is the relative volume that is occupied by a field stronger than 1% of the maximum magnetic field strength at time = 695 years, which is 1.4 kG. Also shown is the average field scaled to fit the maximum of the filling factor (dashed curve). | |
Open with DEXTER |
Figure 6 shows the PDF of the magnetic field: the distribution is a typical signature of highly intermittent structures, i.e. only a very small fraction of the volume carries the strongest structures and the probability of finding a vanishing field strength at a random point in space, is far greater than finding strong fields.
An energy spectrum reveals that the magnetic structures are well resolved with little power at the Nyquist wavenumber , and that the power at the largest wavenumbers is two-orders of magnitude smaller than that at the largest scales (see Fig. 7). Maximum power is obtained on the largest scales corresponding to wavenumbers of a few, while there is a dip at corresponding to the scale of the radius, where the power is minimum. The power on scales 10-20 is flat leaning towards being proportional to k^{-2/3}corresponding to Kolmogorov scaling, and at small-scales power steeply drops: magnetic structures in the non-linear regime are then large by solar standards, but smaller than the giant convection cells that show increasing power towards large scales.
Blackman (1996) argues that in turbulent dynamos the magnetic filling factor would be of the order ( being the ratio of gas pressure to magnetic pressure) so that for relatively weak fields the filling factor would be small. This is partly the case in the present model where the average field is weak field above the stellar surface where it has a low filling factor (see Fig. 8): however, the filling factor scales very well with the average field strength , not with .
Figure 9 is a map of the spherical surface at r = R in
terms of magnetic energy in the non-linear regime: there are both
dark patches of strong magnetic field, e.g. a 500 Gauss region at
longitude around
and latitude
north of
the equator, and large areas with a vanishing field (e.g. at
longitude
on the equator). On this map the areal
magnetic filling factor is 55% for |B|>50 Gauss,
while it is only 0.6% for |B|>500 Gauss. However,
this mapping does not represent any physical surface of the star.
Due to the fact that the actual upper boundary consists of a few
large cells the surface cannot be captured by a simple sphere with
radius R; this is illustrated by a volume rendering (Fig. 10) of the isosurface at the cooling temperature. E.g. in the upper left corner of Fig. 10 there is a "hill''
in the temperature isosurface and further large "slopes'' can be
seen across the star in this illustration. In the latter figure,
the magnetic field lines illustrate that there is a slight trend
inside the star, looking through the partly transparent surface,
towards a radial orientation of the magnetic field, while the
strong fields near the surface of the star are predominantly
horizontally aligned. This was also observed by Dorch & Freytag
(2002) and may be a generic trade of
giant-cell fully convective slowly rotating stars.
Figure 9: An illustration of the unsigned magnetic field strength |B|at the spherical surface r = R of the model star using an orthographic map projection. The darkest patches correspond to a maximum field strength of 500 Gauss (black on the continuous scale bar). From a snapshot at time = 695 years. The views are centered on longitudes of ( left) and ( right). The grid indicated has a longitudinal spacing of and a latitudinal spacing of . The numerical resolution of the map is 180^{2} grid points. | |
Open with DEXTER |
In summary three different modes of dynamo action are recognized:
Based on the results presented here, it is not possible to state conclusively if Betelgeuse actually has a magnetic field, since such a field is unobserved. However, one may conclude that it seems possible that late-type giant stars such as Betelgeuse can indeed have presently undetected magnetic fields. These magnetic fields are likely to be close to or stronger than equipartition yielding surface strengths on the order of 500 Gauss at maximum; this may be difficult to detect directly, due to the relatively small filling factors of the strong fields, but even the moderately strong fields may have influence on their immediate surroundings through altered dust, wind and mass-loss properties. The formation of dust in the presence of a magnetic field will be the subject of a subsequent paper along the lines presented here: the "Pencil Code'' has recently been augmented with modules for radiation and dust modeling.
The dynamo of the late-type giant studied here may be characterized as belonging to the class called "local small-scale dynamos'' another example of which is the proposed dynamo action in the solar photosphere that is sometimes claimed to be responsible for the formation of small-scale flux tubes (cf. Cattaneo 1999). However, in the case of Betelgeuse this designation is less meaningful since the generated magnetic field is both global and large-scale, but because of the slow and non-differential rotation, no large-scale solar-like toroidal field is formed although the situation might be different in more rapidly rotating AGB stars.
It is interesting to note that very recently Lobel et al. (2004) published spatially resolved spectra of the
upper chromosphere and dust envelope of Betelgeuse. Based on
various emission lines they provide evidence for the presence of
warm chromospheric plasma away from the star at around 40 R. The
spectra reveal that Betelgeuse's upper chromosphere extends far
beyond the circumstellar envelope. They compute that temperatures
of the warm chromospheric gas exceed 2600 K. The presence of a hot
chromosphere lead this author to speculate on the possible connection
to coronal heating in the Sun, which is likely to be magnetic in
origin and caused by flux braiding motions in the solar
photosphere (cf. Gudiksen & Nordlund
2002): it remains to be proven whether a
similar process could be operating in late-type giant stars.
Figure 10: A 3-d volume rendering of magnetic field lines (white) and flux ropes (dark structures). Also show is a isosurface (transparent) at the surface temperature value. | |
Open with DEXTER |
Acknowledgements
This work was supported by the Danish Natural Science Research Council. Access to computational resources granted by the Danish Center for Scientific Computing in particular the Horseshoe cluster at Odense University.