2 The model

2.1 Basic equations

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

$${\partial\varrho\over\partial t}+\vec{\nabla}\cdot(\varrho{\v... ..._\varrho^{\rm disc}+q_\varrho^{\rm star} \equiv \dot{\varrho},$$ (1)

where ${\vec{u}}$ is the velocity field and $\varrho$ is the gas density. The mass source, $q_\varrho ^{\rm disc}$, is localized in the disc (excluding the star) and turns on once the local density drops below a reference value $\varrho_0$ (chosen to be equal to the initial density, see Sect. 2.4), i.e.

$$q_\varrho^{\rm disc}=\frac{\xi_{\rm disc}(\vec{r})-\xi_{\rm star}(\vec{r})} {\tau_{\rm disc}}~(\varrho_0-\varrho)_+,$$ (2)

where the subscript "+'' indicates that only positive values are used, i.e.

$$x_+= \left\{ \begin{array}{ll} x\quad&\mbox{if $x>0$ ,}\\ 0&\mbox{otherwise}. \end{array}\right.$$ (3)

This means that matter is injected only if $\varrho<\varrho_0$, and that the strength of the mass source is proportional to the gas density deficit. Throughout this work we use cylindrical polar coordinates $(\varpi,\varphi,z)$, assuming axisymmetry. The shapes of the disc and the central object are specified with the profiles $\xi_{\rm disc}$ and $\xi_{\rm star}$ via

  

$$\xi_{\rm disc}(\varpi,z) \Theta\left({\varpi_0{-}\varpi},{d}\right)~\Theta\left({z_0{-}\vert z\vert},{d}\right),$$ = (4)

$$\xi_{\rm star}(\varpi,z) \Theta\left({r_0{-}r},{d}\right).$$ = (5)

Here, $\Theta(x,d)$ is a smoothed Heaviside step function with a smoothing half-width d set to 8 grid zones in $\xi_{\rm disc}$ and d=r₀ in $\xi_{\rm star}$; $r=\sqrt{\varpi^2+z^2}$ is the spherical radius and r₀ is the stellar radius introduced in Eq. (23); $\varpi_0$ and z₀ are the disc outer radius and disc half-thickness, respectively. So, $\xi_{\rm disc}$ is equal to unity at the disc midplane and vanishes in the corona, whereas  $\xi_{\rm star}$ is unity at the centre of the central object and vanishes outside the star. Note that $\xi_{\rm disc}-\xi_{\rm star}\ge0$ everywhere, since $\xi_{\rm disc}>0$) all the way to the origin, $z=\varpi=0$.

In Eq. (2), $\tau_{\rm disc}$ 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, $M_{\rm disc}/\dot M_{\rm source}$(cf. Sect. 3.4), to avoid unphysical influences of the mass source. We do not fix the distribution and magnitude of $q^{\rm disc}_\varrho$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,

$$q_\varrho^{\rm star} =-\frac{\xi_{\rm star}(\vec{r})}{\tau_{\rm star}}~(\varrho-\varrho_0)_+,$$ (6)

where $\xi_{\rm star}$ is defined in Eq. (5) and $\tau_{\rm star}$is a central accretion time scale that controls the efficiency of the sink. We discuss physically meaningful values of $\tau_{\rm star}$ in Sect. 2.5.

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, $(\vec{{\vec{u}}}_{\rm K}-\vec{{\vec{u}}}) q_\varrho^{\rm disc}$, which would only be absent if the gas were rotating at the Keplerian speed. Thus, the Navier-Stokes equation takes the form

$$\frac{{\rm D} {}{\vec{u}}}{{\rm D} {}t} = - {1\over\varrho}... ...vec{u}}}_{\rm K}-\vec{{\vec{u}}}) q_\varrho^{\rm disc}\right],$$ (7)

where ${\rm D}/{\rm D}t=\partial/\partial t+{\vec{u}}\cdot\vec{\nabla}$ is the advective derivative, p is the gas pressure, $\Phi$ is the gravitational potential, ${\vec{F}}={\vec{J}}\times{\vec{B}}+\vec{\nabla}\cdot{\vec{\tau}}$ is the sum of the Lorentz and viscous forces, ${\vec{J}}=\vec{\nabla}\times{\vec{B}}/\mu_0$ is the current density due to the (mean) magnetic field ${\vec{B}}$, $\mu_0$ is the magnetic permeability, and ${\vec{\tau}}$is the viscous stress tensor,

$$\tau_{ij}=\varrho\nu(\partial u_i/\partial x_j+\partial u_j/\partial x_i).$$ (8)

Here, $\nu$ is the kinematic viscosity, which has been subdivided into three contributions,

$$\nu=\nu_{\rm t}+\nu_{\rm adv}+\nu_{\rm shock}.$$ (9)

The first term is a turbulent (Shakura-Sunyaev) viscosity in the disc,

$$\nu_{\rm t}=\alpha_{\rm SS}c_{\rm s}z_0\xi_{\rm disc}(\vec{r}),$$ (10)

where $c_{\rm s}=(\gamma p/\varrho)^{1/2}$ is the sound speed and $\gamma$ is the ratio of specific heats. The second term is an artificial advection viscosity,

$$\nu_{\rm adv}=c_\nu^{\rm adv}\delta x~({\vec{u}}_{\rm pol}^2+c_{\rm s}^2+v_{\rm A,pol}^2)^{1/2},$$ (11)

which is required to stabilize rapidly moving patterns. Here, $\delta x = \min(\delta\varpi,\delta z)$ is the mesh size, and $c_\nu^{\rm adv}$ is a constant specified in Sect. 2.5. In Eq. (11) we have only used the poloidal velocity, ${\vec{u}}_{\rm pol}$, and the Alfvén speed due to the poloidal magnetic field, $v_{\rm A,pol}=(\vec{B}_{\rm pol}^2/\varrho\mu_0)^{1/2}$, because in an axisymmetric calculation advection in the $\varphi$-direction is unimportant. Finally,

$$\nu_{\rm shock}=c_\nu^{\rm shock}\delta x^2(-\vec{\nabla}\cdot{\vec{u}})_+$$ (12)

is an artificial shock viscosity, with $c_\nu^{\rm shock}$ a constant specified in Sect. 2.5. Note that $\nu_{\rm adv}$ and $\nu_{\rm shock}$ are needed for numerical reasons; they tend to zero for increasing resolution, $\delta x \to 0$.

We assume that the magnetic field in the disc is generated by a standard $\alpha^2\Omega$-dynamo (e.g., Krause & Rädler 1980), which implies an extra electromotive force, $\alpha{\vec{B}}$, in the induction equation for the mean magnetic field, ${\vec{B}}$. To ensure that ${\vec{B}}$ is solenoidal, we solve the induction equation in terms of the vector potential ${\vec{A}}$,

$${\partial{\vec{A}}\over\partial t}={\vec{u}}\times{\vec{B}}+\alpha{\vec{B}}-\eta\mu_0{\vec{J}},$$ (13)

where ${\vec{B}}=\vec{\nabla}\times{\vec{A}}$, and $\eta$ is the magnetic diffusivity.

Since $\alpha(\vec{r})$ has to be antisymmetric about the midplane and vanishing outside the disc, we adopt the form

$$\alpha=\alpha_0~{z\over z_0}~{\xi_{\rm disc}(\vec{r})\over1+v_{\rm A}^2/v_0^2},$$ (14)

where $v_{\rm A}$ is the Alfvén speed based on the total magnetic field, and $\alpha_0$ and v₀ are parameters that control the intensity of dynamo action and the field strength in the disc, respectively. The $\alpha$-effect has been truncated near the axis, so that $\alpha=0$ for $\varpi \le 0.2$. For the magnetic diffusivity we assume

$$\eta=\eta_0+\eta_{\rm t},$$ (15)

where $\eta_0$ is a uniform background diffusivity and the turbulent part, $\eta_{\rm t}= \eta_{{\rm t}0}\xi_{\rm disc}(\vec{r})$, vanishes outside the disc. Thus, magnetic diffusivity in the corona is smaller than in the disc.

Depending on the sign of $\alpha$ and the vertical distribution of $\eta$, the dynamo can generate magnetic fields of either dipolar or quadrupolar symmetry. We shall discuss both types of geometry.

   
2.2 Implementation of a cool disc embedded in a hot corona

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 $10^6~{\rm K}$to $10^8~{\rm K}$ (see, e.g., Feigelson & Montmerle 1999). For the base of the disc corona, temperatures of at least $8\times10^3$ 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 $10^3~{\rm K}$ (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, $s(\vec{r})$, such that s is smaller within the disc and larger in the corona. For a perfect gas this implies $p=K\varrho^\gamma$ (in a dimensionless form), where $K = {\rm e}^{s/c_{\rm v}}$ is a function of position (here p and $\varrho$ are gas pressure and density, $\gamma=c_{\rm p}/c_{\rm v}$, and $c_{\rm p}$ and $c_{\rm v}$ 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

$$[K(\vec{r})]^{1/\gamma} = {\rm e}^{s/c_{\rm p}} = 1-(1{-}\beta)\xi_{\rm disc}(\vec{r}),$$ (16)

where $0<\beta<1$ is a free parameter with $\beta={\rm e}^{-\Delta s/c_{\rm p}}$, that controls the entropy contrast, $\Delta s>0$, between corona and disc. We consider values of $\beta$ between 0.1 and 0.005, which yields an entropy contrast, $\Delta s/c_{\rm p}$, between 2.3 and 5.3. The temperature ratio between disc and corona is roughly $\beta$; see Eq. (26). Assuming pressure equilibrium between disc and corona, and $p\propto\rho T$ for a perfect gas, the corresponding density ratio is $\beta^{-1}$.

2.3 Formulation in terms of potential enthalpy

In the present case it is advantageous to use the potential enthalpy,

$$H=h+\Phi,$$ (17)

as a variable. Here, h is the specific enthalpy, $h = c_{\rm v} T + p/\varrho = c_{\rm p} T$for a perfect gas (with constant specific heats), and T is temperature. Therefore, specific enthalpy h is related to p and $\varrho$ via $h = \gamma(\gamma{-}1)^{-1} p/\varrho$. Specific entropy s is related to p and $\varrho$ (up to an additive constant) through $s=c_{\rm v}\ln p-c_{\rm p}\ln\varrho$ for a perfect gas. We have therefore

$$\frac{{\rm D} {}\ln h}{{\rm D} {}t} = \frac{{\rm D} {}\ln p}... ... \frac{\gamma}{c_{\rm p}}\frac{{\rm D} {}s}{{\rm D} {}t} \cdot$$ (18)

Since s is independent of time, ${\rm D} {}s/{\rm D} {}t={\vec{u}}\cdot\vec{\nabla} s$. Together with Eq. (1), this yields an evolution equation for h, which can be written in terms of H as

$$\frac{{\rm D} {}H}{{\rm D} {}t}={\vec{u}}\cdot\vec{\nabla}\Ph... ...right) + {\gamma h\over c_{\rm p}}{\vec{u}}\cdot\vec{\nabla}s.$$ (19)

In the following we solve Eq. (19) instead of Eq. (1). In terms of h, density and sound speed are given by

$$\varrho^{\gamma-1}=\frac{\gamma-1}{\gamma}~\frac{h}{K},\quad c_{\rm s}^2=(\gamma-1)h= \gamma\frac{p}{\varrho}\cdot$$ (20)

It proved necessary to include an artificial diffusion term in Eq. (19) with a diffusion coefficient proportional to $\nu$.

The first law of thermodynamics allows us to express the pressure gradient in terms of h and s,

$$-{1\over\varrho}\vec{\nabla}p=-\vec{\nabla}h+T\vec{\nabla} s,$$ (21)

and with $T=h/c_{\rm p}$ we obtain

$$\frac{{\rm D} {}{\vec{u}}}{{\rm D} {}t} = - \vec{\nabla}H +... ...vec{u}}}_{\rm K}-\vec{{\vec{u}}}) q_\varrho^{\rm disc}\right],$$ (22)

which now replaces Eq. (7).

We use a softened, spherically symmetric gravitational potential of the form

$$\Phi=-GM_* \left(r^n+r_0^n\right)^{-1/n},$$ (23)

where G is the gravitational constant, M_* is the mass of the central object, r is the spherical radius, r₀ is the softening radius, and n=5; tentatively, r₀ can be identified with the stellar radius.

   
2.4 The initial state

Our initial state is the hydrostatic equilibrium obtained by solving, for h, the vertical balance equation obtained from Eq. (22),

$$-{\partial \over\partial z}(h+\Phi) + h{\partial s/c_{\rm p}\over\partial z} = 0,$$ (24)

from large z (where $h+\Phi=0$) down to z=0. The initial density distribution $\varrho_0$ is then obtained using Eq. (20); it decreases monotonically with both z and $\varpi$ in the equilibrium state.

The initial rotation velocity, $u_{\varphi0}$, follows from the radial balance equation,

$$- {u_{\varphi0}^2\over\varpi} = - {\partial\over\partial\varpi}(h+\Phi) + h{\partial s/c_{\rm p}\over\partial\varpi}\cdot$$ (25)

In the disc, $h=c_{\rm p} T$ is small, so $u_{\varphi0}$ is close to the Keplerian velocity, while the corona does not rotate initially, and so is supported by the pressure gradient.

As a rough estimate, the value of h in the midplane of the disc is

$$h_{\rm disc}\approx\beta\ h_{\rm corona} \equiv\beta GM_*/\varpi,$$ (26)

as can be seen by integrating Eq. (24), ignoring the $\partial\Phi/\partial z$ term. Thus, the initial toroidal velocity in the midplane can be obtained from Eq. (25) using Eq. (26) and recalling that $\partial s/\partial\varpi=0$in the midplane, which gives

$$u_{\varphi0}\approx\sqrt{(1-\beta)GM_*/\varpi}=\sqrt{1-\beta}\vec{{\vec{u}}}_{\rm K}.$$ (27)

For $\beta =0.1$, for example, the toroidal velocity is within 5% of the Keplerian velocity.

   
2.5 Dimensionless variables and choice of parameters

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 $\Sigma_0\approx1~{\rm g}~{\rm cm}^{-2}$. A typical coronal sound speed is $c_{\rm s0} \approx 10^2~{\rm km}~{\rm s}^{-1}$, which corresponds to a temperature of $T \approx 4\times10^5~{\rm K}$. This allows us to fix relevant units as follows:

$$[\Sigma]=\Sigma_0,\ \; [u]=c_{\rm s0},\ \; [x]=GM_* c_{\rm s0}^{-2},\ \; [s]=c_{\rm p}.$$ (28)

With $M_*=1~M_\odot$ we have $[x] \approx 0.1~{\rm AU}$, and $[t]=[x]/[u] \approx 1.5~{\rm day}$. The unit for density is $[\varrho]=[\Sigma]/[x] \approx 7\times10^{-13}~{\rm g}~{\rm cm}^{-3}$, the unit for the mass accretion rate is $[\dot{M}]=[\Sigma] [u] [x] \approx 2\times10^{-7}~M_\odot~{\rm yr}^{-1}$, and the unit for the magnetic field is $[B]=[u] (\mu_0[\varrho])^{1/2} \approx 30~{\rm G}$.

Since [h]=[u]², the dimensionless value h=1 corresponds to $10^{14}~{\rm cm}^2~{\rm s}^{-2}$. With a mean specific weight $\mu=0.6$, the universal gas constant ${\cal R}=8.3\times10^7~{\rm cm}^2~{\rm s}^{-2}~{\rm K}^{-1}$ and $\gamma=5/3$, we have

$$c_{\rm p}={\gamma\over\gamma-1}{{\cal R}\over\mu} \approx3.5\times10^8~{\rm cm}^2~{\rm s}^{-2}~{\rm K}^{-1}.$$ (29)

Therefore, h=1 corresponds to $T=[u]^2/[s]=[h]/c_{\rm p}\approx3\times10^5~{\rm K}$. Using $h=-\Phi$ in the corona, this corresponds to a temperature of $3\times10^5~{\rm K}$ at $r=[x]\equiv0.1~{\rm AU}$.

We choose $\beta$ between 0.1 and 0.005, corresponding to a typical disc temperature (in the model) of $3\times10^4~{\rm K}$ to $1.5\times10^3~{\rm K}$; see Eq. (26); real protostellar discs have typical temperatures of about a few thousand Kelvin (see Sect. 2.2).

In our models, we use $\varpi_0=1.5$ for the disc outer radius, z₀=0.15 for the disc semi-thickness and r₀=0.05 for the softening (stellar) radius. The disc aspect ratio is $z_0/\varpi_0=0.1$. Note that r₀=0.05 corresponds to $7\times10^{10}~{\rm cm}$, i.e. one solar radius. Therefore, we shall not reduce it much below this physically meaningful value. Note, however, that smaller values of r₀ would result in faster outflows (Ouyed & Pudritz 1999).

Furthermore, $c_\nu^{\rm adv}=0.02$ and $c_\nu^{\rm shock}=1.2$. We vary the value of $\alpha_{\rm SS}$ between 0.003 and 0.007.

The mean-field dynamo is characterized by the parameters $\vert\alpha_0\vert=0.3$, $v_0={c_{\rm s}}$, $\eta_{{\rm t}0}=10^{-3}$ and $\eta_0=5\times10^{-4}$. The total magnetic diffusivity in the disc is therefore $\eta_{\rm T0}=\eta_{\rm t0}+\eta_0 = 0.0015$. In terms of the usual Shakura-Sunyaev viscosity parameter, this corresponds to

$$\alpha_{\rm SS}^{(\eta)}\equiv {\eta_{\rm T0}\over c_{\rm s,d... ..._{\rm T0}\over0.0015}\right) \left({z_0\over0.15}\right)^{-1},$$ (30)

where we have used $c_{\rm s,disc}^2\approx\beta c_{\rm s,corona}^2$(cf. Eqs. (20) and (26)) and $c_{\rm s,corona} \approx \varpi^{-1/2}$.

In terms of the usual nondimensional dynamo parameters we have

$$\vert C_\alpha\vert=\vert\alpha_0\vert z_0/\eta_{\rm T0}=30$$ (31)

and, for Keplerian rotation,

$$C_\omega=(\varpi\partial\Omega/\partial\varpi)z_0^2/\eta_{\rm T0} =-22.5\varpi^{-3/2},$$ (32)

so that the dynamo number is given by

$$\vert{\cal D}\vert=\vert C_\alpha C_\omega\vert=675\varpi^{-3/2}.$$ (33)

Note that the value of the dynamo number expected for accretion discs is given by

$$\vert{\cal D}\vert=\frac{ \vert\alpha_0~\varpi\partial\Omega/... ... T0}^2} \la{3\over2}\left(\alpha_{\rm SS}^{(\eta)}\right)^{-2}$$ (34)

for $\eta_{\rm T0}=\alpha_{\rm SS}^{(\eta)}c_{\rm s,disc}z_0$, $\vert\alpha_0\vert\la c_{\rm s,disc}$ and $c_{\rm s,disc}=\Omega z_0$. As follows from Eqs. (30)-(33), for Keplerian rotation,

$${2\over 3} \vert{\cal D}\vert \left(\alpha_{\rm SS}^{(\eta)}\right)^2 = \vert\alpha_0\vert z_0 \beta^{-1} \varpi^{-1/2}.$$ (35)

(Note that this expression is independent of $\eta_{\rm T0}$.) For our choice of parameters of $\vert\alpha_0\vert=0.3$ and z₀=0.15, expression (35) is $\la$1 if $\varpi \ga 0.002\ \beta^{-2} = 0.2$ for $\beta =0.1$, which corresponds to the truncation radius of the $\alpha$-effect. Therefore, our choice of parameters is consistent with the constraint (34).

We take $\gamma=5/3$, $\tau_{\rm disc}=0.1$, and consider two values of $\tau_{\rm star}$. For $\tau_{\rm star}\to\infty$, the mass sink at the central object is suppressed, whereas $\tau_{\rm star}\to0$ implies instantaneous accretion of any extra matter (relative to the hydrostatic equilibrium) by the central object. A realistic lower limit for $\tau_{\rm star}$ can be estimated as $\tau_{\rm star}\ga r_0/v_{\rm ff}\approx0.008$, where $v_{\rm ff}$ is the free fall velocity (given by ${\textstyle{1\over2}}v_{\rm ff}^2=GM_*/r_0$ in dimensional quantities). In most cases we used $\tau _{\rm star}=0.01$, but we also tried an even smaller value of  $\tau _{\rm star}=0.005$ and obtained very similar results. A finite value of $\tau_{\rm star}$ 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 $\tau_{\rm star}$. 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 $(\varpi,z)\in[0,2]\times[-1,1]$ to $[0,8]\times [-2,30]$, but the results are hardly different in the overlapping parts of the domains. In our standard computational box, $\delta \varpi=\delta z=0.01$and in the case of the larger computational domain $\delta \varpi=\delta z=0.02$.

2.6 Numerical method and boundary conditions

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 $\varpi$ and $\varphi$ 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 $\log_{10}h$ at time t=223. Specific enthalpy h is directly proportional to temperature T, and $\log_{10}h=(-1,0,1)$ corresponds to $T\approx (3\times10^{4},3\times10^{5},3\times10^{6})~\mbox{K}$. 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 ${c_{\rm s}}$. The disc boundary is shown with a thin black line; $\beta =0.1$.

Figure 3: As in Fig. 2, but with a mass sink at the centre, $\tau _{\rm star}=0.01$, 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.