A&A 398, 825-844 (2003)
DOI: 10.1051/0004-6361:20021699
B. von Rekowski^{} - A. Brandenburg^{} - W. Dobler^{} - A. Shukurov
Department of Mathematics, University of Newcastle, Newcastle upon Tyne NE1 7RU, UK
Received 25 May 2000 / Accepted 14 November 2002
Abstract
We present an axisymmetric numerical model of a dynamo active accretion disc.
If the dynamo-generated magnetic field in the disc is sufficiently strong
(close to equipartition with thermal energy),
a fast magneto-centrifugally driven outflow develops within a conical
shell near the
rotation axis, together with
a slower pressure driven outflow from the outer
parts of the disc as well as around the axis. Our results show that a dynamo
active accretion disc can
contribute to driving an outflow even without any external magnetic field.
The fast outflow in the conical shell is confined by the
azimuthal field
which is produced by the dynamo in the disc and advected to the disc
corona. This part of the outflow has high angular momentum and
is cooler and less dense than its surroundings.
The conical shell's half-opening angle is typically about
near the disc and decreases slightly with height.
The slow outflow is hotter and denser.
Key words: ISM: jets and outflows - accretion, accretion disks - magnetic fields - MHD
Extensive numerical studies of collimated disc winds have been performed using several types of model. Uchida & Shibata (1985), Matsumoto et al. (1996) and Kudoh et al. (1998) consider an ideal MHD model of an accretion disc embedded in a non-rotating corona, permeated by an external magnetic field initially aligned with the disc rotation axis. Intense accretion develops in the disc (Stone & Norman 1994), accompanied by a centrifugally driven outflow. The wind is eventually collimated by toroidal magnetic field produced in the corona by winding up the poloidal field.
Bell (1994), Bell & Lucek (1995) and Lucek & Bell (1996, 1997) use two- and three-dimensional ideal MHD models, with a polytropic equation of state, to study the formation and stability of pressure driven jets collimated by poloidal magnetic field, which has a minimum inside the jet. The general structure of thermally driven disc winds was studied analytically by, e.g., Fukue (1989) and Takahara et al. (1989).
Another type of model was developed by Ustyugova et al. (1995), Romanova et al. (1997, 1998), Ouyed et al. (1997) and Ouyed & Pudritz (1997a, 1997b, 1999) who consider ideal MHD in a polytropic corona permeated by an external poloidal magnetic field and subsume the physics of the accretion disc into boundary conditions at the base of the corona. The disc is assumed to be in Keplerian rotation (with any accretion neglected). The corona is non-rotating initially, and the system is driven by the injection of matter through the boundary that represents the disc surface. These models develop a steady (or at least statistically steady) state consistent with the analytical models by Blandford & Payne (1982), showing a magneto-centrifugal wind collimated by toroidal magnetic field, which again is produced in the corona by the vertical shear in the angular velocity. Three-dimensional simulations suggest that the resulting collimated outflow does not break due to the kink instability (Ouyed et al. 2003; Thomsen & Nordlund 2003).
It is not quite clear how strong the external magnetic field in accretion discs can be. Dragging of an external field from large radii in a viscous disc requires that magnetic diffusivity is much smaller than the kinetic one (Lubow et al. 1994), i.e. that the magnetic Prandtl number is significantly larger than unity which would be difficult to explain in a turbulent disc (Heyvaerts et al. 1996). This argument neglects, however, the effect of magnetic torques which could produce significant field line dragging even when the magnetic Prandtl number is of order unity (Shalybkov & Rüdiger 2000). On the other hand, the efficiency of trapping an external magnetic field at initial stages of the disc formation is questionable because only a small fraction of the external flux can be retained if the density contrast between the disc and the surrounding medium is large. In that case the magnetic field will be strongly bent and reconnection will remove most of the flux from the disc. Furthermore, the ionization fraction is probably too small in the inner disc to ensure sufficient coupling between the gas and the magnetic field (Fromang et al. 2002). For these reasons it seems appropriate to explore whether or not magnetic fields generated by a dynamo within the disc can produce an outflow with realistic properties.
The models of Ustyugova et al. (1995) and Ouyed & Pudritz (1997a, 1997b, 1999) have a
distributed mass source and imposed poloidal velocity at the (unresolved) disc surface. As a result,
their system develops a persistent outflow which is not just a transient. On the other hand, the models of
Matsumoto et al. (1996) resolve the disc, but have no mass replenishment to compensate losses
via the outflow and accretion, and so the disc disappears with time. Our model is an attempt to
combine advantages of both of these approaches and also to add dynamo action in the disc.
Instead of a rigidly prescribed mass injection, we allow for
self-regulatory replenishment of matter within a resolved disc. Instead of prescribing
poloidal velocity at the disc surface, we resolve the disc and prescribe an
entropy contrast between the disc and
the corona, leaving more freedom for the system. Such an increase of entropy with
height is only natural to expect for a disc surrounded by a hot corona,
and we
parameterize the coronal heating by a (fixed) entropy contrast. We further add self-sustained, intrinsic
magnetic field to our system, as opposed to an external field used in the other models.
Since our model goes beyond ideal MHD, the magnetic field must be maintained against
decay. A simple form of mean-field dynamo action is included for this
purpose.
Figure 1: General structure of the outflows typically obtained in the magnetic runs with mass sink presented in this paper. The cool, dense disc emits (a) a hardly collimated, thermally driven wind (slow, hot, dense, magnetized, rotating), (b) a magneto-centrifugally driven outflow in a conical shell (faster, cooler, less dense, magnetized, quickly rotating), and (c) a thermally driven outflow near the axis (slow, hot, dense, weakly magnetized and weakly rotating). | |
Open with DEXTER |
Like in the model of Ouyed & Pudritz (1997a, 1997b, 1999), the hot, pressure supported corona does not rotate initially. The disc is cool and is therefore centrifugally supported, so its rotation is nearly Keplerian. The corresponding steady-state solution is used as our initial condition. This solution is, however, unstable because of the vertical shear in the angular velocity between the disc and the corona (Urpin & Brandenburg 1998). Angular momentum transfer by viscous and magnetic stresses also leads to a departure from the initial state. As a result, a meridional flow develops, which exchanges matter between the corona and the disc surface layers. Mass losses through the disc surface and to the accreting central object are then replenished in the disc where we allow local mass production whenever and wherever the density decreases below a reference value. Matter is heated as it moves into the corona; this leads to an increase in pressure which drives the wind. Another efficient driver of the outflow in our model is magneto-centrifugal acceleration. A strong toroidal magnetic field produced by the dynamo in the disc is advected into the corona and contributes to confining the wind.
Altogether, our model contains many features included in a range of earlier models. For example, the disc is essentially resistive to admit dynamo action, whereas magnetic diffusion turns out to be relatively unimportant in the corona, similarly to the models of, e.g., Wardle & Königl (1993), Ferreira & Pelletier (1995) and Casse & Ferreira (2000a, 2000b). The synthetic nature of the outflow, driven by both pressure forces and magneto-centrifugally, is characteristic of the self-similar solutions of Ferreira & Pelletier (1995) and Casse & Ferreira (2000b). The latter authors also stress the rôle of the vertical entropy gradient in enhancing the mass flux in the outflow and prescribe a vertical profile of the entropy production rate. Structured outflows outwardly similar to our results have been discussed by Krasnopolsky et al. (1999) and Goodson et al. (1997). However, in the former paper the structure in their solutions is due to the boundary conditions imposed on the disc surface, and in the latter model the structure is due to the interaction between the stellar magnetic field and the disc under the presence of differential rotation between star and disc. On the contrary, the structure in our model results from the dominance of different driving mechanisms in the inner and outer parts of the disc.
While we had initially expected to find collimated outflows similar to those reported by Ouyed & Pudritz (1997a, 1997b), the outflow patterns we obtained turned out to be of quite different nature. However, a significant difference in our model is the resolved disc, which is necessary to model dynamo action in the disc, and the relatively small extension into the corona, where collimation is not yet expected. Figure 1 shows a sketch of the overall, multi-component structure of the outflows in the presence of magnetic fields and mass sink, which needs to be compared to the individual figures in the rest of the paper.
The plan of the paper is as follows. We introduce our model in Sect. 2, and then consider a range of parameters in Sect. 3 to clarify and illustrate the main physical effects. In Sect. 4 we explore the parameter space, and our conclusions are presented in Sect. 5.
Our intention is to make our model as simple as possible and to avoid
detailed modelling of mass supply to the accretion disc, which occurs
differently in different accreting systems. For example,
matter enters the accretion disc
at large radii in a restricted range of azimuthal angles in binary systems
with Roche lobe overflow. On the other hand, matter supply
can be more uniform in both azimuth and radius in active galactic nuclei.
Being interested in other aspects of accretion physics, we prefer to avoid
detailed modelling of these processes. Instead,
similarly to the model of Ouyed & Pudritz
(1997a, 1997b), we inject matter into the system, but an important difference is that
we introduce a self-regulating mass source in the disc.
We also allow for a mass sink at the centre to model
accretion onto the central object (star). With these two effects included,
the continuity equation becomes
(3) |
In Eq. (2), is a response time which is chosen to be significantly shorter than the time scale of the depletion processes, which is equivalent to the time scale of mass replenishment in the disc, (cf. Sect. 3.4), to avoid unphysical influences of the mass source. We do not fix the distribution and magnitude of beforehand, but the system adjusts itself such as to prevent the disc from disappearing.
The self-regulating mass sink at the position of the central star
is modelled in a similar manner,
Apart from the continuity Eq. (1), the mass source also appears
in the Navier-Stokes equation, unless matter is always injected with the
ambient velocity of the gas. In that case, however, a runaway instability
can occur: if matter is already slower than Keplerian, it falls inward, and so
does the newly injected matter. This enhances the need for even more mass
injection. A similar argument applies also if matter is rotating faster
than Keplerian. This is why we inject matter at the Keplerian speed.
This leads to an extra term in the Navier-Stokes equation,
,
which would only
be absent if the gas were rotating at the Keplerian speed. Thus,
the Navier-Stokes equation takes the form
(8) |
(9) |
(12) |
We assume that the magnetic field in the disc is generated by
a standard
-dynamo (e.g., Krause & Rädler 1980),
which implies an extra electromotive force,
,
in the
induction equation for the mean magnetic field, .
To ensure that
is solenoidal, we solve the
induction equation in terms of the vector potential ,
Since
has to be
antisymmetric about the midplane and vanishing outside the disc,
we adopt the form
Depending on the sign of and the vertical distribution of , the dynamo can generate magnetic fields of either dipolar or quadrupolar symmetry. We shall discuss both types of geometry.
Protostellar systems are known to be strong X-ray sources (see, e.g., Glassgold et al. 2000; Feigelson & Montmerle 1999; Grosso et al. 2000). The X-ray emission is generally attributed to coronae of disc-star systems, plausibly heated by small scale magnetic reconnection events (Galeev et al. 1979), for example in the form of nanoflares that are caused by slow footpoint motions (Parker 1994). Heating of disc coronae by fluctuating magnetic fields is indeed quite natural if one accepts that the disc turbulence is caused by the magneto-rotational instability. Estimates for the coronal temperatures of YSOs range from to (see, e.g., Feigelson & Montmerle 1999). For the base of the disc corona, temperatures of at least are to be expected in order to explain the observed mass loss rates (Kwan & Tademaru 1995; Kwan 1997). The discs, on the other hand, have typical temperatures of a few (e.g., Papaloizou & Terquem 1999).
A simple way to implement a dense, relatively cool disc embedded in a rarefied, hot corona without modelling the detailed physics of coronal heating is to prescribe the distribution of specific entropy, , such that s is smaller within the disc and larger in the corona. For a perfect gas this implies (in a dimensionless form), where is a function of position (here p and are gas pressure and density, , and and are the specific heats at constant pressure and constant volume, respectively).
We prescribe the polytrope parameter K to be unity in the
corona and smaller in the disc, so we put
In the present case it is advantageous to use the potential enthalpy,
(17) |
(18) |
The first law of thermodynamics allows us to
express the pressure gradient in terms of
h and s,
We use a softened, spherically symmetric gravitational potential of the form
Our initial state is the hydrostatic equilibrium obtained by
solving, for h, the vertical balance equation obtained from Eq. (22),
The
initial rotation velocity,
,
follows from the radial
balance equation,
As a rough estimate, the value of h in the midplane of the
disc is
(27) |
Our model is scale invariant and can therefore be applied to various
astrophysical objects.
We consider here the range of parameters typical of protostellar discs, for which a typical
surface density is
.
A typical coronal sound speed is
,
which corresponds to a temperature of
.
This allows us to fix relevant units as follows:
Since [h]=[u]^{2},
the dimensionless value h=1 corresponds
to
.
With a mean specific weight ,
the universal gas constant
and
,
we have
(29) |
We choose between 0.1 and 0.005, corresponding to a typical disc temperature (in the model) of to ; see Eq. (26); real protostellar discs have typical temperatures of about a few thousand Kelvin (see Sect. 2.2).
In our models, we use for the disc outer radius, z_{0}=0.15 for the disc semi-thickness and r_{0}=0.05 for the softening (stellar) radius. The disc aspect ratio is . Note that r_{0}=0.05 corresponds to , i.e. one solar radius. Therefore, we shall not reduce it much below this physically meaningful value. Note, however, that smaller values of r_{0} would result in faster outflows (Ouyed & Pudritz 1999).
Furthermore, and . We vary the value of between 0.003 and 0.007.
The mean-field dynamo is characterized by the parameters
,
,
and
.
The total magnetic diffusivity in the disc is therefore
.
In terms of the usual Shakura-Sunyaev viscosity parameter,
this corresponds to
In terms of the usual nondimensional dynamo parameters we have
(31) |
(32) |
We take , , and consider two values of . For , the mass sink at the central object is suppressed, whereas implies instantaneous accretion of any extra matter (relative to the hydrostatic equilibrium) by the central object. A realistic lower limit for can be estimated as , where is the free fall velocity (given by in dimensional quantities). In most cases we used , but we also tried an even smaller value of and obtained very similar results. A finite value of implies that matter is not instantaneously absorbed by the sink. Therefore, some matter can leave the sink if it moves so fast that its residence time in the sink is shorter than . As can be seen below, a small (negligible) fraction of mass does indeed escape from the sink.
Computations have been carried out in domains ranging from to , but the results are hardly different in the overlapping parts of the domains. In our standard computational box, and in the case of the larger computational domain .
We use third order Runge-Kutta time-stepping and a sixth order finite-difference scheme in space. Details and test calculations are discussed by Brandenburg (2003).
On the outer boundaries, the induction equation is evolved using
one-sided derivatives (open boundary conditions).
The normal velocity component has
zero derivative normal to the boundary, but the velocity
is required to be always directed outwards.
The tangential velocity
components and potential enthalpy on the boundaries are similarly obtained
from the next two interior points away from the boundary.
Tests show that the presence of the boundaries does not affect the flow
inside the computational domain. Regularity
conditions are adopted on the axis where
the
and
components of vectorial quantities vanish,
whereas scalar variables and the z components of vectorial quantities
have vanishing radial derivative.
Figure 2: A nonmagnetic model without mass sink at the centre: velocity vectors and logarithmically spaced density contours (from 0.008 to 6000, separation factor of 1.6) superimposed on a grey scale representation of at time t=223. Specific enthalpy h is directly proportional to temperature T, and corresponds to . Note hot gas near the central object and in the near corona, and cooler gas in the disc and the far corona. The dashed line shows the sonic surface where the poloidal velocity equals . The disc boundary is shown with a thin black line; . | |
Open with DEXTER |
Figure 3: As in Fig. 2, but with a mass sink at the centre, , at time t=312. The outflow speed near the axis is strongly reduced in comparison to that in the model without mass sink shown in Fig. 2. | |
Open with DEXTER |
We illustrate in Fig. 2 results obtained for a model without any magnetic fields and without mass sink. A strong outflow develops even in this case, which is driven mostly by the vertical pressure gradient in the transition layer between the disc and the corona, in particular by the term in Eq. (21). The gain in velocity is controlled by the total specific entropy difference between the disc and the corona, but not by the thickness d of the transition layer in the disc profile (4).
The flow is fastest along the rotation axis and within a cone of polar angle of about , where the terminal velocity is reached. The conical shape of the outflow is partly due to obstruction from the dense disc, making it easier for matter to leave the disc near the axis. Both temperature and density in nonmagnetic runs without mass sink are reduced very close to the axis where the flow speed is highest.
The general flow pattern is sensitive to whether or not matter can accrete onto the central object. We show in Fig. 3 results for the same model as in Fig. 2, but with a mass sink given by Eq. (6) with . This can be compared with earlier work on thermally driven winds by Fukue (1989) who also considered polytropic outflows, but the disc was treated as a boundary condition. In Fukue (1989), outflows are driven when the injected energy is above a critical value. The origin of this energy injection may be a hot corona. The critical surface in his model (see the lower panel of his Fig. 2) is quite similar to that found in our simulation (our Fig. 3), although our opening angle was found to be larger than in Fukue's (1989) model. (Below we show, albeit with magnetic fields, that smaller values of do result in smaller opening angles, see Fig. 9, which would then be compatible with the result of Fukue (1989).)
As could be expected, the mass sink hampers the outflow in the cone (but not at ). The flow remains very similar to that of Fig. 3 when is reduced to 0.005. Thus, the nonmagnetic outflows are very sensitive to the presence of the central sink. As we show now, magnetized outflows are different in this respect.
In this section we discuss results obtained with magnetic fields,
first without mass sink at the centre
and then including a sink. We show that the effects of the sink are
significantly weaker than in the nonmagnetic case.
Figure 4: A magnetic model without mass sink at the centre: velocity vectors and poloidal magnetic field lines (white) superimposed on a grey scale representation of temperature (in terms of h) for a run with (resulting in roughly dipolar magnetic symmetry) at time t=269. In contrast to the nonmagnetic model of Fig. 2, a conical shell has developed that is cooler and less dense than its exterior. The conical shell intersects at around . The dashed line shows the surface where the poloidal velocity equals , with the Alfvén speed from the poloidal magnetic field (fast magnetosonic surface with respect to the poloidal field). The disc boundary is shown with a thin black line; . | |
Open with DEXTER |
Figure 5: Profiles of vertical velocity and sound speed along the axis, , for the nonmagnetic model of Fig. 2 with smoothing radius r_{0}=0.05(dotted) and two runs with magnetic field, with r_{0}=0.05 as in Fig. 4 (solid) and r_{0}=0.02 (dashed, t=162). The shaded area marks the location of the disc. Terminal wind speeds are reached after approximately three disc heights. The presence of a magnetic field and a deeper potential well (smaller r_{0}) both result in faster outflows. The models shown here have no mass sink at the centre and . | |
Open with DEXTER |
Figure 6: Results for a larger domain, , with the same parameters as in Fig. 4, averaged over times , . Left panel: velocity vectors, poloidal magnetic field lines and grey scale representation of hin the inner part of the domain. Right panel: velocity vectors, poloidal magnetic field lines and normalized specific enthalpy in the full domain. | |
Open with DEXTER |
A magnetized outflow without the central mass sink, shown in Fig. 4, is similar to that in Fig. 2, but is denser and hotter near the axis, and the high speed cone has a somewhat larger opening angle. In addition, the outflow becomes more structured, with a well pronounced conical shell where temperature and density are smaller than elsewhere (the conical shell reaches at in Fig. 4). Here and in some of the following figures we also show the fast magnetosonic surface with respect to the poloidal field. In Sect. 3.5 we show that this surface is close to the fast magnetosonic surface. As shown in Fig. 5, the outflow becomes faster inside the cone ( on the axis). As expected, we find that deeper potential wells, i.e. smaller values of r_{0} in Eq. (5) and (23), result in even faster flows and in larger opening angles.
Our results are insensitive to the size and symmetry of the computational domain: we
illustrate this in Fig. 6 with a larger
domain of the size
.
The disc midplane in this run is located asymmetrically in zin order to verify that the (approximate) symmetry of the solutions
is not imposed by the symmetry of the computational domain. Figure 6 confirms that our results
are not affected by what happens near the computational boundaries.
Figure 7: As in Fig. 4, but with a mass sink at the centre, . The outflow speed near the axis and the opening angle of the conical shell (now reaching at ) are reduced in comparison to that in the model without mass sink shown in Fig. 4, but the general structure of the outflow is little affected. Averaged over times , . (Reference model.) | |
Open with DEXTER |
Unlike the nonmagnetized system, the magnetized outflow changes only comparatively little when the mass sink is introduced at the centre. We show in Fig. 7 the results with a sink ( ) and otherwise the same parameters as in Fig. 4. As could be expected, the sink leads to a reduction in the outflow speed near the axis; the flow in the high speed cone becomes slower. But apart from that, the most important effects of the sink are the enhancement of the conical structure of the outflow and the smaller opening angle of the conical shell. A decrease in by a factor of 2 to 0.005 has very little effect, as illustrated in Figs. 8 and 10.
Increasing the entropy contrast (while keeping the specific entropy unchanged in the corona)
reduces the opening angle of the conical shell.
Pressure driving is obviously more important in this case, as compared
to magneto-centrifugal driving (see Sect. 3.5).
A model with
(corresponding to a density and inverse temperature
contrast of about 50:1 between the disc and the corona) is shown in
Fig. 9.
At
,
the radial velocity
is slightly
enhanced relative to the case
(contrast 10:1);
see Fig. 10.
At
,
on the other hand, the flow with the larger entropy contrast
reaches the Alfvén point close to the disc (at
as opposed to in the other case), which leads to a smaller terminal velocity.
Figure 8: As in Fig. 7, but with a stronger mass sink at the centre, . The flow pattern is very similar to that of Fig. 7 where the time scale of the sink is twice as large. Averaged over times , . | |
Open with DEXTER |
Figure 9: As in Fig. 7, but with a larger entropy contrast, . The opening angle of the conical shell is reduced in comparison to that of Figs. 7 and 8 where the entropy contrast is smaller (the shell crosses at ). Averaged over times t=140...240, . | |
Open with DEXTER |
We conclude that the general structure of the
magnetized flow and its typical parameters
remain largely unaffected by the sink, provided its efficiency
does not exceed a certain threshold.
It is plausibly the build-up of magnetic pressure at the centre that shields the central
object to make the central accretion inefficient. This shielding would be
even stronger if we included a magnetosphere of the central object. We discuss in
Sect. 4.1 the dependence of our solution on the geometrical size of the sink and
show that the general structure of the outflow persists as long as the size of the sink does
not exceed the disk thickness.
Figure 10: Spherical radial velocity component, (dashed), and poloidal Alfvén speed, (solid), as functions of spherical radius at polar angles (thick lines) and (thin lines) for the models of Figs. 7 (top panel), 8 (middle panel) and 9 (bottom panel). | |
Open with DEXTER |
The dynamo in most of our models has , consistent with results from simulations of disc turbulence driven by the magneto-rotational instability (Brandenburg et al. 1995; Ziegler & Rüdiger 2000). The resulting field symmetry is roughly dipolar, which seems to be typical of disc dynamos with in a conducting corona (e.g., Brandenburg et al. 1990). We note that the dominant parity of the magnetic field is sensitive to the magnetic diffusivity in the corona: a quadrupolar oscillatory magnetic field dominates for if the disc is surrounded by vacuum (Stepinski & Levy 1988).
For
,
the critical value of
for dynamo action is
about 0.2, which is a factor of about 50 larger than without
outflows.
Our dynamo is then only less than twice supercritical. A
survey of the dynamo regimes for similar models is given by
Bardou et al. (2001).
Figure 11: Radial profiles of magnetic pressure from the toroidal (solid) and poloidal (dotted) magnetic fields and thermal pressure (dashed) for the model of Fig. 6. Shown are the averages over the disc volume (lower panel) and over a region of the same size around z=8 in the corona (upper panel). . | |
Open with DEXTER |
The initial magnetic field (poloidal, mixed parity) is weak ( ), cf. Fig. 11 for comparison with the gas pressure], but the dynamo soon amplifies the field in the disc to , and then supplies it to the corona. As a result, the corona is filled with a predominantly azimuthal field with at larger radii; see Fig. 11. We note, however, that the flow in the corona varies significantly in both space and time^{}. Magnetic pressure due to the toroidal field , , exceeds gas pressure in the corona outside the inner cone and confines the outflow in the conical shell. The main mechanisms producing in the corona are advection by the wind and magnetic buoyancy (cf. Moss et al. 1999). Magnetic diffusion and stretching of the poloidal field by vertical shear play a relatively unimportant rôle.
The field in the inner parts of the disc is dominated by the toroidal component; at ; this ratio is larger in the corona at all . However, as shown in Fig. 11, this ratio is closer to unity at larger radii in the disc.
As expected, results in mostly quadrupolar fields (e.g., Ruzmaikin et al. 1988). As shown in Fig. 12, the magnetic field in the corona is now mainly restricted to a narrow conical shell that crosses at . Comparing this figure with the results obtained with dipolar magnetic fields (Fig. 4), one sees that the quadrupolar field has a weaker effect on the outflow than the dipolar field; the conical shell is less pronounced. However, the structures within the inner cone are qualitatively similar to each other.
The magnitude and distribution of in Eq. (14) only weakly affect magnetic field properties as far as the dynamo is saturated. For a saturated dynamo, the field distribution in the dynamo region ( , |z|<0.15) roughly follows from the equipartition field, with . In other words, nonlinear states of disc dynamos are almost insensitive to the detailed properties of (e.g., Beck et al. 1996; Ruzmaikin et al. 1988).
A discussion of disc dynamos with outflows, motivated by the present model, can
be found in Bardou et al. (2001). It is shown that the value of magnetic
diffusivity in the corona does not affect the dynamo solutions strongly.
Moreover, the outflow is fast enough to have the magnetic Reynolds number
in the corona larger than unity, which implies that ideal integrals
of motion are very nearly constant along field lines;
see Sect. 3.6. The most important
property is the sign of
as it controls the global symmetry of the
magnetic field.
Figure 12: As in Fig. 4, but with , at time t=132. The magnetic field geometry is now mostly quadrupolar because . | |
Open with DEXTER |
The mass injection and loss rates
due to the source, sink and wind are defined as
Figure 13: Azimuthally integrated mass flux density, represented as a vector , in the simulation of Fig. 7 with a dipolar magnetic field and a mass sink at the centre. Shades of grey show the distribution of the mass source in the disc, . The disc boundary is shown with a white line. | |
Open with DEXTER |
The mass loss rate in the wind fluctuates on a time scale of 5 time units, but remains constant on average at about , corresponding to , in the models of Figs. 7 and 8. The mass in the disc, , remains roughly constant.
The rate at which mass needs to be replenished in the disc, , is about 0.4. This rate is not controlled by the imposed response rate of the mass source, , which is 25 times larger. So, the mass source adjusts itself to the disc evolution and does not directly control the outflow. We show in Fig. 13 trajectories that start in and around the mass injection region. The spatial distribution of the mass replenishment rate shown in Fig. 13 indicates that the mass is mainly injected close to the mass sink, and remains moderate in the outer parts of the disc. (Note that the reduced effect of the mass sink in the magnetized flow is due to magnetic shielding rather than to mass replenishment near the sink - see Sect. 3.2.)
The angular structure of the outflow can be characterized by the following
quantities calculated for a particular spherical
radius, r=8,
for the model of Fig. 6: the azimuthally integrated normalized
radial mass flux density,
,
where
(38) |
(39) |
(40) |
(41) |
Figure 14: Dependence, on polar angle , of azimuthally integrated radial mass flux density through a sphere r=8 (solid, normalized by the disc mass ), azimuthally integrated radial angular momentum flux density (dashed-dotted, normalized by the disc angular momentum ), azimuthally integrated radial Poynting flux density (divided by 5, dashed, normalized by the magnetic energy in the disc ), and azimuthally integrated radial kinetic energy flux density (divided by 10, dash-3dots, normalized by the kinetic energy in the disc ), for the model of Fig. 6. The unit of all the quantities is . | |
Open with DEXTER |
Figure 15: Angular momentum (normalized by the maximum angular momentum in the disc) for the model with and , shown in Fig. 7. The maximum value is . The dashed line shows the fast magnetosonic surface with respect to the poloidal field (cf. Fig. 4), the solid line the Alfvén surface where the poloidal velocity equals the poloidal Alfvén speed, and the dotted line the sonic surface. Averaged over times t=130...140. | |
Open with DEXTER |
Figure 16: The ratio of the poloidal magneto-centrifugal and pressure forces, as defined in Eq. (42), is shown with shades of grey, with larger values corresponding to lighter shades. The maximum value is . Superimposed are the poloidal magnetic field lines. The dashed line shows the fast magnetosonic surface with respect to the poloidal field (cf. Fig. 4). The parameters are as in the model of Fig. 7. Averaged over times t=130...140. | |
Open with DEXTER |
Figure 17: The ratio in the corona for the model with and , shown in Fig. 7. For a force-free magnetic field, , this ratio is . Averaged over times t=130...140. | |
Open with DEXTER |
Figure 18: As in Fig. 16, but for the model of Fig. 9, i.e. with larger entropy contrast, . Averaged over times t=140...240. | |
Open with DEXTER |
The radial kinetic and magnetic energy flux densities, integrated over the whole sphere, are and , respectively, where is the fast magnetosonic speed (with respect to the poloidal field) at the critical surface (where u_{z}=c_{*}) on the axis. Thus, can be taken as a good indicator of the kinetic energy loss, and the magnetic energy loss into the exterior is about 3% of this value. These surface-integrated flux densities (or luminosities) are, as expected, roughly independent of the distance from the central object.
The magnetized outflows in our models with central mass sink have a well-pronounced structure, with a fast, cool and low-density flow in a conical shell, and a slower, hotter and denser flow near the axis and in the outer parts of the domain. Without central mass sink, there is a high speed, hot and dense cone around the axis.
The magnetic field geometry (e.g., Fig. 7) is such that for , the angle between the disc surface and the field lines is less than , reaching at -1.5, which is favourable for magneto-centrifugal acceleration (Blandford & Payne 1982; Campbell 1999, 2000, 2001). However, the Alfvén surface is so close to the disc surface in the outer parts of the disc that acceleration there is mainly due to pressure gradient. The situation is, however, different in the conical shell containing the fast wind. As can be seen from Fig. 15, the Alfvén surface is far away from the disc in that region and, on a given field line, the Alfvén radius is at least a few times larger than the radius of the footpoint in the disc. This is also seen in simulations of the magneto-centrifugally driven jets of Krasnopolsky et al. (1999); see their Fig. 4. The lever arm of about 3 is sufficient for magneto-centrifugal driving to dominate. As can be seen also from Fig. 10, the flow at the polar angle is mainly accelerated by pressure gradient near the disc surface (where the Alfvén surface is close to the disc surface). However, acceleration remains efficient out to at least r=1within the conical shell at . This can be seen in the upper and middle panels of Fig. 10 (note that the conical shell is thinner and at a smaller in the model with larger entropy contrast, and so it cannot be seen in this figure, cf. Fig. 9). These facts strongly indicate that magneto-centrifugal acceleration dominates within the conical shell.
Another indicator of magneto-centrifugal acceleration in the conical shell is the distribution
of angular momentum (see Fig. 15), which is significantly larger in the outer parts of
the conical shell than in the disc, which suggests that the magnetic field plays
an important rôle in the flow dynamics. We
show in Fig. 16 the ratio of the "magneto-centrifugal force''
to pressure gradient,
,
where
the subscript "pol'' denotes the poloidal components.
Here, the "magneto-centrifugal force'' includes all terms in the
poloidal equation of motion, except for the pressure gradient (but we ignore
the viscous term and the mass production term, the latter being
restricted to the disc only),
As further evidence of a significant contribution from
magneto-centrifugal driving in the conical shell, we show in Fig. 17
that the magnetic field is close to a force-free configuration in regions where
angular momentum is enhanced, i.e. in the conical shell and in the corona surrounding
the outer parts of the disc. These are the regions where the Lorentz force contributes
significantly to the flow dynamics, so that the magnetic field performs work and therefore relaxes
to a force-free configuration.
The radial variation in the sign of the current helicity
is due to
a variation in the sign of the azimuthal magnetic field and of the current
density.
Such changes originate from the disc where they
imprint a corresponding variation in the sign of the
angular momentum constant, see Eq. (46).
These variations are then carried along magnetic lines into the corona.
The locations where the azimuthal field reverses are still relatively
close to the axis, and there the azimuthal field relative to the poloidal field
is weak compared to
regions further away from the axis.
Figure 19: Flow parameters as functions of height z, along the field line with its footpoint at , for the model of Fig. 7. Upper left: poloidal Alfvén Mach number, (dashed), and a similar quantity that includes the poloidal Alfvén speed as well as the sound speed (fast magnetosonic Mach number with respect to the poloidal field), (dotted) - the two curves are practically identical; upper right: ratio of toroidal ( ) and poloidal ( ) magnetic fields; lower left: toroidal velocity (solid) and poloidal velocity (dash-3dots), in units of the toroidal velocity at the footpoint, ; lower right: density in units of the density at the footpoint, . The position of the Alfvén (and fast magnetosonic) point on the field line is indicated by the vertical line. This figure is useful to compare with Fig. 3 of Ouyed et al. (1997). | |
Open with DEXTER |
Pressure driving is more important if the entropy contrast between the disc and the corona is larger (i.e. is smaller): the white conical shell in Fig. 16 indicative of stronger magneto-centrifugal driving shifts to larger heights for , as shown in Fig. 18. We note, however, that there are periods when magneto-centrifugal driving is dominant even in this model with higher entropy contrast, but pressure driving dominates in the time averaged picture (at least within our computational domain).
We show in Fig. 19 the variation of several quantities along
a magnetic field line that has its footpoint at the disc
surface at
and lies around the conical shell.
Since this is where magneto-centrifugal driving is still dominant, it is useful
to compare Fig. 19 with Fig. 3 of Ouyed et al. (1997), where a
well-collimated magneto-centrifugal jet is studied.
Since our outflow is collimated only weakly within our computational domain,
the quantities are plotted against height z, rather than
as in
Ouyed et al. (1997) (
is nearly constant along a field line for weakly collimated
flows, whereas approximately
along a magnetic line for well-collimated flows).
The results are qualitatively similar, with the main difference that
the fast magnetosonic surface in our model almost coincides with the Alfvén surface
in the region around the conical shell where the outflow is highly supersonic.
Since we include finite diffusivity, the curves in Fig. 19 are smoother
than in Ouyed et al. (1997), who consider ideal MHD.
A peculiar feature of the conical shell is that the flow at
is sub-Alfvénic
but strongly supersonic.
The fast magnetosonic surface is where the poloidal
velocity
equals the fast magnetosonic speed for the
direction parallel to the field lines,
(43) |
Axisymmetric ideal magnetized outflows are governed by five
Lagrangian invariants,
the flux ratio,
(44) |
(45) |
(47) |
Figure 20: The four Lagrangian invariants k(a), , , and e(a)as a function of the flux function afor the model of Fig. 6. All points in the domain , are shown (provided ). The dots deviating from well-defined curves mostly originate in . Note that specific entropy as fifth Lagrangian invariant is trivially conserved in our model, because entropy is spatially constant throughout the corona and temporally fixed. | |
Open with DEXTER |
We show in Fig. 20 scatter plots of k(a), , , and e(a)for the model of Fig. 6. Points from the region collapse into a single line, confirming that the flow in the corona is nearly ideal^{}. This is not surprising since the magnetic Reynolds number is much greater than unity in the corona for the parameters adopted here. For , there are departures from perfect MHD; in particular, the angular velocity of magnetic field lines, , is somewhat decreased in the upper parts of the domain (indicated by the vertical scatter in the data points). This is plausibly due to the finite magnetic diffusivity which still allows matter to slightly lag behind the magnetic field. As this lag accumulates along a stream line, the departures increase with height z. Since this is a "secular'' effect only, and accumulates with height, we locally still have little variation of k and , which explains why magneto-centrifugal acceleration can operate quite efficiently.
The corona in our model has (turbulent) magnetic diffusivity comparable to that in the disc. This is consistent with, e.g., Ouyed & Pudritz (1999) who argue that turbulence should be significant in coronae of accretion discs. Nevertheless, it turns out that ideal MHD is a reasonable approximation for the corona (see Fig. 20), but not for the disc. Therefore, magnetic diffusivity is physically significant in the disc and insignificant in the corona (due to different velocity and length scales involved), as in most models of disc outflows (see Ferreira & Pelletier 1995 for a discussion). Thus, our model confirms this widely adopted idealization.
The models presented so far have some clear deficiencies when compared with characteristic features relevant to protostellar discs. The first deficiency concerns the relative amount of matter accreted onto the central object compared with what goes into the wind. A typical figure from the models presented above was that as much as 30% of all the matter released in the disc goes into the wind and only about 70% is accreted by the central object. Earlier estimates (e.g., Pelletier & Pudritz 1992) indicate that only about 10% of the matter joins the wind. Another possible deficiency is the fact that in the models presented so far the low-temperature region extends all the way to the stellar surface whilst in reality the cool disc breaks up close to the star because of the stellar magnetic field. Finally, the overall temperature of the disc is generally too high compared with real protostellar discs which are known to have typical temperatures of about a few thousand Kelvin.
The aim of this section is to assess the significance of various improvements in the model related to the above mentioned characteristics. We consider the effect of each of them separately by improving the model step by step.
As discussed in Sects. 3.1 and 3.2, properties of the outflow,
especially of the nonmagnetic ones, are sensitive to the parameters of the
central sink. It is clear that a very efficient sink would completely inhibit
the outflow near the axis. On the other hand, a magnetosphere of the central
object can affect the sink efficiency by channelling the flow along the stellar
magnetic field (Shu et al. 1994; see also Fendt & Elstner 1999, 2000).
Simulating a magnetosphere turned out
to be a difficult task and some preliminary attempts proved unsatisfactory.
Figure 21: As in Fig. 7, but with a mass sink at the centre which is five times larger, r_{0} = 0.25, and . Note the absence of outflows along the axis and a larger opening angle of the conical shell, which crosses at . Averaged over times , . | |
Open with DEXTER |
Instead, we have considered a model with a geometrically larger sink, and illustrate it in Fig. 21. Here, r_{0}=0.25 in Eqs. (5) and (23), which is five times larger than the sink used in our reference model of Fig. 7. The relaxation time of the sink in Eq. (6), , was rescaled in proportion to the free-fall time at r_{0} as , which yields . The resulting mass loss rate into the wind is , which corresponds to in dimensional units. Although this is indeed smaller than the value for the reference model ( ), the overall mass released from the disc is also smaller, resulting in a larger fraction of about 40% of matter that goes into the wind; only about 60% is accreted by the central object. As we show below, however, larger accretion fractions can more easily be achieved by making the disc cooler. We conclude that even a sink as large as almost twice the disc half-thickness does not destroy the outflow outside the inner cone. However, the outflow along the axis is nearly completely suppressed.
As can be seen from Fig. 22, the gap in the disc translates directly into a corresponding gap in the resulting outflow pattern. At the same time, however, only about 20% of the released disc material is accreted by the sink and the rest is ejected into the wind. A small fraction of the mass that accretes toward the sink reaches the region near the axis at some distance away from the origin and can then still be ejected as in the models without the gap. Note that our figures (e.g., Fig. 22) show the velocity rather than the azimuthally integrated mass flux density (cf. Fig. 13); the relative magnitude of the latter is much smaller, and so the mass flux from the sink region is not significant.
We now turn to the discussion of the case with a cooler disc. Numerical
constraints prevent us from making the entropy gradient between disc and
corona too steep. Nevertheless, we were able to reduce the value of down to 0.005, which is 20 times smaller than the value used for the
reference model in Fig. 7. The
value for the sink was
reduced to 0.02; smaller values proved numerically difficult.
Figure 22: As in Fig. 7, but with a gap between the sink and the disc with the disc inner radius at . The radius of the stellar mass sink is three times larger, r_{0} = 0.15, and . Outflow is absent along field lines passing through the gap. Averaged over times , in both the disc and the sink. | |
Open with DEXTER |
Figure 23: As in Fig. 22, but with in the disc and in the sink. Averaged over times when the overall magnetic activity in the disc is relatively low. | |
Open with DEXTER |
Figure 24: As in Fig. 23, but averaged over a period of enhanced magnetic activity in the disc, . A conical shell develops which crosses at . | |
Open with DEXTER |
The value of results in a disc temperature of K in the outer parts and K in the inner parts. As in Sect. 4.2, the disc terminates at an inner radius of . At , the density in dimensional units is about 10^{-9} g cm^{-3}, which is also the order of magnitude found for protostellar discs. The resulting magnetic field and outflow geometries are shown in Figs. 23 and 24.
A characteristic feature of models with a cooler disc is a more vigorous temporal behaviour with prolonged episodes of reduced overall magnetic activity in the disc during which the Alfvénic surface is closer to the disc surface, and phases of enhanced magnetic activity where the Alfvénic surface has moved further away. Figures 23 and 24 are representative of these two states. It is notable that the structured outflow of the type seen in the reference model occurs in states with strong magnetic field and disappears during periods with weak magnetic field.
Another interesting property of the cooler discs is that now a smaller fraction of matter goes into the wind (10-20%), and 80-90% is accreted by the central object, in a better agreement with the estimates of Pelletier & Pudritz (1992).
We note in passing that channel flow solutions typical of two-dimensional simulations of the magneto-rotational instability (Hawley & Balbus 1991) are generally absent in the simulations presented here. This is because the magnetic field saturates at a level close to equipartition between magnetic and thermal energies. The vertical wavelength of the instability can then exceed the half-thickness of the disc. In some of our simulations, indications of channel flow behaviour still can be seen. An example is Fig. 23 where the magnetic energy is weak enough so that the magneto-rotational instability is not yet suppressed.
According to the Shakura-Sunyaev prescription, turbulent viscosity and magnetic diffusivity are reduced in a cooler disc because of the smaller sound speed (cf. Eq. (10)). Since in the corona the dominant contributions to the artificial advection viscosity come from the poloidal velocity and poloidal Alfvén speed (and not from the sound speed), has to be reduced. Here we choose and reduce by a factor of 10 to . Since we do not explicitly parameterize the turbulent magnetic diffusivity with the sound speed (cf. Eq. (15)), also needs to be decreased, together with the background diffusivity . We choose here values that are 25 times smaller compared to the previous runs, i.e. and , which corresponds to ranging between 0.001 and 0.007. With this choice of coefficients, the total viscosity and magnetic diffusivity in the disc have comparable values of a few times 10^{-5}.
The effect of reduced viscosity and magnetic diffusivity is shown in
Fig. 25. Features characteristic of the
channel flow solution are now present,
because the vertical wavelength of the magneto-rotational instability is here
less than the half-thickness of the disc.
Figure 25: As in Figs. 23 and 24, but with , , , and . Averaged over times (where the magnetic energy is enhanced), . | |
Open with DEXTER |
If the disc dynamo is sufficiently strong, our model develops a clearly structured outflow which is fast, cool and rarefied within a conical shell near the rotation axis where most of the angular momentum and magnetic energy is carried, and is slower, hotter and denser in the region around the axis as well as in the outer parts of the domain. The slower outflow is driven mostly by the entropy contrast between the disc and the corona, but the faster wind within the conical shell is mostly driven magneto-centrifugally. Without a central mass sink, the flow near the axis is faster, but otherwise the flow structure is similar to that with the sink.
The half-opening angle of the cone with hot, dense gas around the axis is about -; this quantity somewhat changes with model parameters but remains close to that range. The outflow in our models does not show any signs of collimation. It should be noted, however, that not all outflows from protostellar discs are actually collimated, especially not at such small distance from the source. An example is the Becklin-Neugebauer/Kleinmann-Low (BN/KL) region in the Orion Nebula (Greenhill et al. 1998), which has a conical outflow with a half-opening angle of out to a distance of 25-60 AU from its origin. Therefore, collimation within a few AU (the size of our computational domain) is expected to be only weak.
The region around the fast, cool and rarefied conical shell seen in Fig. 7 is similar to the flow structure reported by Krasnopolsky et al. (1999); see their Fig. 1. In their model, however, the thin axial jet was caused by an explicit injection of matter from the inner parts of the disc which was treated as a boundary. In our reference model the fast outflow is sub-Alfvénic because of the presence of a relatively strong poloidal field, whereas in Krasnopolsky et al. (1999) the outflow becomes super-Alfvénic at smaller heights. Outside the conical shell the outflow is mainly pressure driven, even though the criterion of Blandford & Payne (1982) is fulfilled. However, as Casse & Ferreira (2000b) pointed out, pressure driven outflows might dominate over centrifugally driven outflows if thermal effects are strong enough.
In our model, matter is replenished in the resolved disc in a self-regulatory manner where and when needed. We believe that this is an improvement in comparison to the models of Ouyed & Pudritz (1997a, 1997b, 1999) and Ustyugova et al. (1995), where mass inflow is prescribed as a boundary condition at the base of the corona. If we put in Eqs. (1) and (7), the disc mass soon drops to low values and the outflow ceases. This is qualitatively the same behaviour as in the models of, e.g., Kudoh et al. (1998).
We should stress the importance of finite magnetic diffusivity in the disc: although poloidal velocity and poloidal magnetic field are well aligned in most of the corona, dynamo action in the disc is only possible in the presence of finite magnetic diffusivity, and the flow can enter the corona only by crossing magnetic field lines in the disc.
An outflow occurs in the presence of both dipolar and quadrupolar type magnetic fields, even though fields with dipolar symmetry seem to be more efficient in magneto-centrifugal driving (cf. von Rekowski et al. 2000). The effects of the magnetic parity on the outflow structure deserves further analysis.
The dynamo active accretion disc drives a significant outward Poynting flux in our model. Assuming that this applies equally to protostellar and AGN discs, this result could be important for understanding the origin of seed magnetic fields in galaxies and galaxy clusters; see Jafelice & Opher (1992) and Brandenburg (2000) for a discussion. We note, however, that the pressure of the intracluster gas may prevent the magnetized plasma from active galactic nuclei to spread over a significant volume (Goldshmidt & Rephaeli 1994).
Our model can be improved in several respects. In many systems, both dynamo-generated and external magnetic fields may be present, so a more comprehensive model should include both. We used an dynamo to parameterize magnetic field generation in the disc because we have restricted ourselves to axisymmetric models. As argued by Brandenburg (1998), dynamo action of turbulence which is driven by the magneto-rotational instability can be roughly described as an dynamo. But this parameterization can be relaxed in three-dimensional simulations where one may expect that turbulence will be generated to drive dynamo action. Such simulations will be discussed elsewhere (von Rekowski et al. 2002).
Since our model includes angular momentum transport by both viscous and magnetic stresses, it is natural that the accreted matter is eventually diverted into an outflow near the axis; this is further facilitated by our prescribed entropy gradient at the disc surface. We believe that this picture is physically well motivated (Bell & Lucek 1995), with the only reservation that we do not incorporate the (more complicated) physics of coronal heating and disc cooling, but rather parameterize it with a fixed entropy contrast. We include a mass sink at the centre which could have prevented the outflow, and indeed the sink strongly affects nonmagnetized outflows. We have shown, however, that the magnetic field can efficiently shield the sink and thereby support a vigorous disc wind.
The assumption of a prescribed entropy distribution is a useful tool to control the size of the disc and to parameterize the heating of the disc corona. However, it should be relaxed as soon as the disc physics can be described more fully. The energy equation, possibly with radiation transfer, should be included. This would lead to a self-consistent entropy distribution and would admit the deposition of viscous and Ohmic heat in the outflow. In the simulations by von Rekowski et al. (2002), entropy is evolved.
We believe that a mass source is a necessary feature of any model of this kind if one wishes to obtain a steady state. In the present paper the mass source is distributed throughout the whole disc to represent replenishment of matter from the midplane of the disc. Alternatively, a mass source could be located near to or on the domain boundary.
Acknowledgements
We are grateful to C. G. Campbell, R. Henriksen, J. Heyvaerts, R. Ouyed and R. E. Pudritz for fruitful discussions. We acknowledge many useful comments of the anonymous referee. This work was partially supported by PPARC (Grants PPA/G/S/1997/00284 and 2000/00528) and the Leverhulme Trust (Grant F/125/AL). Use of the PPARC supported supercomputers in St Andrews and Leicester is acknowledged.