A&A 384, 329-342 (2002)
DOI: 10.1051/0004-6361:20011726

Numerical simulations of astrophysical jets from randomly perturbed Keplerian disks

B. F. Vitorino - V. Jatenco-Pereira - R. Opher

Instituto Astronômico e Geofísico, Universidade de São Paulo, Av. Miguel Stéfano 4200,
04301-904 São Paulo, Brazil

Received 29 November 2000 / Accepted 26 November 2001

Abstract
We study the outflows from randomly perturbed Keplerian accretion disks. In the velocity field of the gas at the surface of the disk, we introduce a random perturbation proportional to r^a, where r is the radial distance and a is a fixed parameter in the simulation (with values of -3/2, -1, -1/2, 0, 1/2 or 1). In the simulations, the disk is a fixed boundary from which the gas is ejected with continuous and random velocities added together. The continuous ejection velocity is taken to be a thousandth of that of the local Keplerian disk; the maximum value of the random velocity is 10 times the continuous velocity. We observe the formation of periodic structures along the outflow axis. The distance of the axial separation between the structures is found to be about 11 times the innermost radius of the accretion disk. For a negative radial slope a, the structures are well defined and situated near the jet axis, whereas for a positive slope, decollimation of the jet occurs.

Key words: ISM: jets and outflows - galaxies: jets - accretion, accretion disks - magnetohydrodynamics (MHD)

1 Introduction

Two possible mechanisms have been suggested to explain the production of astrophysical jets: magnetic forces associated with an accretion disk (Königl & Pudritz 1999) and the interaction between an accretion disk and a magnetized star, in the case of protostellar jets (Shu et al. 1999). In the first mechanism, the disk-wind driven model, the gas is centrifugally flung out from the surface of a Keplerian disk when the field lines thread the disk at an angle of $60^\circ$or less with respect to the surface. When the gas reaches the local Alfvén velocity, the ( $\mbox{\boldmath$J$ }_{\rm p} \times \mbox{\boldmath$B$ }_{\phi}$) term becomes dominant, providing the collimation for the outflow. In the second mechanism, the X-wind model, the wind emanates near the region where the disk corotates with the central magnetic star. Some magnetospheric lines, connecting the central star and the disk, are partially open. The wind emanates from the inner edge of the disk and flows out along the open magnetic field lines. Excess angular momentum brought in by the accreting disk is extracted by the magnetocentrifugally driven wind. The first mechanism could be a universal model since it seeks to explain astrophysical jets independently of the detailed nature of the central object, from active galactic nuclei (AGNs) to young stellar objects.

Episodic outbursts and knots have been observed in astrophysical jets. Recent theories and models try to provide a natural explanation for these features, common to all astrophysical jets.

There are two approaches in making the numerical simulations of axisymmetric magnetized flows from accrection disks. One approach attempts to follow the internal dynamics of the disk (Shibata & Uchida 1986; Stone & Norman 1994). The second approach does not follow the internal dynamics of the disk, but treats it as aboundary condition in Keplerian equilibrium at the base of the wind (Bell & Lucek 1995; Ustyugova et al. 1995).

Adopting the latter aproach, Ouyed & Pudritz (1997a,b, 1999) made 2.5-dimensional, high-resolution numerical MHD simulations of the onset and collimation of outflows from the surface of a Keplerian accretion disk around a central object. They obtained steady or episodic jets, depending upon the assumed values of the constant (in time) ejection velocity from the disk to the corona. For the smallest values of the ejection velocity, their results showed the episodic production of knots near the surface of the disk, which propagate along the jet axis. These knots are produced by MHD shocks.

We made 2.5-dimensional time-dependent simulations with the ZEUS-3D code (Stone & Norman 1992a,b) in the axisymmetry option, introducing a random velocity perturbation in the accretion disk. We know that accretion disks are subject to instabilities, such as the Balbus-Hawley instability. However, the origin of instabilities, in general, is unknown. We investigate here the behaviour of the produced jets when the perturbations have a radial dependence in the disk. The perturbation was random in time and its amplitude, proportional to r^a, where r is the radial distance and a defines the radial dependence of the amplitude of the perturbation. For each time interval, the velocity amplitude was chosen by a random generator at all radial positions. The time interval is defined by the time step controled by the code which satisfies the Courant-Friedrichs-Lewy (CFL) stability condition (this condition guaranties causality). In our simulations, the time step was $\sim$ $6.65 \times 10^{-2}$ periods of rotation at the innermost radius of the accretion disk $r_{\rm i}$. The radial position, where the perturbation is introduced, is defined by the radial resolution used. A minimum interval of $0.25~r_{\rm i}$ was used at the radial position $r_{\rm i}$, while a maximum interval of $\sim$1.16 $r_{\rm i}$ was used at $r=20~r_{\rm i}$.

A random component was added to the continuous component with a maximum value at $r_{\rm i}$ chosen to be 10 times the continuous ejection velocity. The simulations were run with different radial dependencies for the perturbation. We simulated the region 1- $80~r_{\rm i}$ of the accretion disk.

Introducing the random perturbation with different radial slopes (i.e., values of a), we observed the fomation of periodic structures along the outflow axis. The distance of the axial separation of the structures was found to be about $11~r_{\rm i}$. For a negative radial slope for the perturbation (i.e., negative values of a), the density structures are well defined near the jet axis. For a positive slope, decollimation of the jet occurs. As the slope from zero to positive values is increased, the region near the jet axis is found to be less and less well-defined.

In Sect. 2 the basic equations that were solved are presented. The model is described in Sect. 3. In Sect. 4, the numerical results are presented and discussed. Finally, in Sect. 5, we present a summary and the final conclusions.

2 Basic equations

Cylindrical coordinates (r, $\phi$, z) were adopted, with the central object at the origin. The z-axis is perpendicular to the disk, such that the surface of the disk lies in the z=0 plane.

The equations that describe the conservation of mass, momentum, energy and magnetic flux are:

$$\frac{\partial \rho }{\partial t}+{\bf\nabla} \cdot \left( \rho \mbox{\boldmath$v$ }\right) =0,$$ (1)

$$\rho \left[ \frac{\partial \mbox{\boldmath$v$ }}{\partial t}+... ...a} \Phi -\mbox{\boldmath$j$ }\times {\mbox{\boldmath$B$ }} =0,$$ (2)

$$\rho \left[ \frac{\partial e}{\partial t}+\left( \mbox{\boldm... ...ht] +p\left( {\bf\nabla} \cdot \mbox{\boldmath$v$ }\right) =0,$$ (3)

$$\frac{\partial \mbox{\boldmath$B$ }}{\partial t}-{\bf\nabla} ... ...t( \mbox{\boldmath$v$ }\times \mbox{\boldmath$B$ } \right) =0,$$ (4)

$${\bf\nabla} \cdot \mbox{\boldmath$B$ }=0,$$ (5)

where $\rho$ is the gas density, $\mbox{\boldmath$v$ }$ the velocity, $\mbox{\boldmath$B$ }$ the magnetic field, p the gas pressure, e the internal energy, $\mbox{\boldmath$j$ }={\bf\nabla} \times \mbox{\boldmath$B$ }/4\pi$the electric current density and $\Phi =-GM/\left( r^2+z^2\right) ^{1/2}$, the gravitational potential.

We assumed a polytropic gas $p=K\rho ^{\gamma}$ with the polytropic index $\gamma=5/3$. Thus, ignoring radiative transfer effects, the energy Eq. (3) need not be solved. It was also assumed that an Alfvénic turbulent pressure $p_{\rm A}$ provides, together with the thermal pressure, the main support against gravity for the cold corona. Here $p_{\rm A}=p/\beta _{\rm T}$, where $\beta _{\rm T}$ is the ratio of the gas pressure to the Alfvénic turbulent pressure, assumed to be constant (OP97a).

All parameters are given in units of their values at the innermost radius of the accretion disk $r_{\rm i}$. Thus, the normalized equation of motion (Eq. (2)) is

 

$$\left[ \frac{\partial \mbox{\boldmath$v$ }^{\prime }}{\partial t}... ...boldmath$J$ }^{\prime }\times \mbox{\boldmath$B$ }^{\prime }}{\rho ^{\prime }},$$ (6)

where the dimensionless variables are: $r^{\prime }=r/r_{\rm i}$, $z^{\prime }=z/r_{\rm i}$, $v^{\prime } =v/v_{{\rm K},{\rm i}}$, $\rho ^{\prime }=\rho /\rho _{\rm i}$, $p^{\prime }=p/p_{\rm i}$, $B^{\prime }=B/B_{\rm i}$, $\Phi ^{\prime }=-1/\left( r^{\prime 2}+z^{^{\prime }2}\right) ^{1/2}$ and $\nabla ^{\prime }=r_{\rm i}\nabla$. The parameter  $\delta _{\rm i}=\rho_{\rm i}v_{{\rm K},{\rm i}}^2/p_{\rm i}$ is the ratio of the Keplerian to the thermal energy density and $\beta _{\rm i}=8\pi p_{\rm i}/B_{\rm i}^2$ is the ratio of the gas to the magnetic pressure in the corona at $r_{\rm i}$. At $r_{\rm i}$, the Keplerian speed is $v_{{\rm K},{\rm i}}=\left( GM/r_{\rm i}\right) ^{1/2}$ and $B_{\rm i}$ is the poloidal field at $r_{\rm i}$. For convenience, we drop the primes and refer only to the normalized variables. Time is given in units of $t_{\rm i}=r_{\rm i}/v_{{\rm K},{\rm i}}$, so that the dimensionless time is

$$\tau =\frac t{t_{\rm i}} \cdot$$ (7)

3 The model

We study a model in which the disk is a fixed boundary condition in pressure equilibrium with the corona. There is a point mass at the center of the disk that represents a star (or a black hole, a neutron star or a white dwarf). An initial potential (zero current) poloidal field threads the disk and the corona.

In the numerical simulations, the computational domain defines the spatial region where the MHD equations are time evolved. ZEUS-3D, a grid code, divides the computational domain into cells, called "active" zones. Finite differencing the evolution near the boundaries of the grid requires that values for the dependent variables be specified beyond the computational domain. Therefore, at each boundary, "ghost" zones are added. Values for the dependent variables in the ghost zones are specified, using boundary conditions appropriate to the geometry and the physics of the problem being solved. Thus, the evolution equations are not solved for the ghost zones.

Inflow and outflow are used to specify the boundary conditions which mean that a physical quantity is flowing inside or outside of the computational domain, respectively. We use the expression injection velocity to have the same meaning as the term inflow.

The values of all the variables in the ghost zones in the inflow boundary condition are held equal to a set of predetermined values (which may not be allowed to vary in time). Outflow is not permited. In the outflow boundary condition, the values of the variables in the ghost zones are equal to those in the corresponding active zones (extrapolating the flow beyond the boundary).

In the reflecting boundary condition, zone-centered variables (scalars), as well as the tangential components of the velocity and the magnetic field in the ghost zones are set equal to the corresponding values of their images in the active zones. The normal component of the velocity and the magnetic field are set equal to zero on the boundary and reflected into the adjacent zone.

We use inflow boundary conditions at the disk surface and open outflow conditions on the remaining boundaries. Reflecting boundaries are used along the axis of symmetry, coinciding with the disk axis. The initial setup and boundaries are shown in Fig. 1.

	Figure 1: The adopted model. The central object is placed at the origin of the coordinate system. The ghost zones (z<0) define the surface of the Keplerian disk, while the active zones (which evolve in time) define the corona. The poloidal magnetic field continues smoothly into the disk surface (Fig. 2 of OP97a).
Open with DEXTER

Initially, the corona has no toroidal magnetic field. In order to assure that an outflow can only arise due to the magnetocentrifugal mechanism, magnetostatic equilibrium is numerically maintained. Resolving Eq. (6) for magnetostatic equilibrium, we obtain the expression for the coronal density distribution,

$$\rho =\left( 1/\sqrt{r^2+z^2}\right) ^{3/2} ,$$ (8)

which is used in our numerical simulations.

Assuming the surface of the accretion disk to be a fixed boundary condition, there is no back reaction by the jet. Time evolution occurs only in the coronal region ($z\geq 0$) where the jet exists. However, we fix the density at the base of the corona (first column of cells), assuming that the disk is a suitable dynamic source of matter. The density distribution of the disk surface has the form

$$\rho _{\rm d}=\eta _{\rm i}r^{-3/2} ,$$ (9)

where $\eta _{\rm i}=\rho _{\rm d}/\rho _0$ is the ratio between the disk density ( $\rho _{\rm d}$) and the density at the base of the corona ($\rho _0$).

A poloidal magnetic field in a potential configuration ( $\mbox{\boldmath$J$ }=0$) threads the corona and the disk. This field is described by the stream function,

$$\Phi_{\rm m} = \sqrt {r^2+\left( z_{\rm d}+z \right)^2}- \left( z_{\rm d}+z \right),$$ (10)

where $z_{\rm d}$ is the dimensionless disk thickness. The field components are given as

$$B_z= \frac {1}{r}\frac {\partial {\Phi_{\rm m}}} {\partial {r... ...rac {1}{r}\frac {\partial {\Phi_{\rm m}}} {\partial {z}} \cdot$$ (11)

Due to the torsion of the poloidal field at the disk surface, a toroidal magnetic field is produced at the base of the corona, which could propagate into the disk leaving it with a nonzero toroidal field. However, our boundary conditions do not allow this to happen. Thus, we antecipate that disks must have nonzero toroidal fields. We chose a current-free configuration for such a field to be

$$B_\phi =\frac{\mu _{\rm i}}r ,$$ (12)

where the input parameter $\mu _{\rm i}=B_{\phi {\rm i}}/B_{\rm i}$ is the ratio of the toroidal to the poloidal magnetic field at $r_{\rm i}$. The disk-corona field was matched by the ZEUS-3D code. There are three important initial timescales defined at the disk's surface: the Keplerian time $\tau _{\rm K}$, the Alfvén time $\tau _{\rm A}$ and the magnetic braking time $\tau _{\rm B}$. Thus, we have (e.g., Eq. (4.64) of OP97a)

$$\left( \tau _{\rm K},\;\tau _{\rm A},\;\tau _{\rm B}\right) =... ...; \eta _{\rm i}\sqrt{\delta _{\rm i}\beta _{\rm i}/2}\right) .$$ (13)

To estimate the turbulent pressure needed to support the corona for a given gas pressure, we used the equation for the ratio of the turbulent to the magnetic pressure $\beta _{\rm T}$,

$$\beta _{\rm T}= \gamma / \left( \left( \gamma -1\right) \delta _{\rm i}-\gamma \right).$$ (14)

As in Eq. (6), $\delta_{\rm i}$ is the ratio of the Keplerian to the thermal energy density.

The gas is continuously injected from the disk into the corona at a velocity $\mbox{\boldmath$v$ }_{\rm p}=v_{\rm inj} v_{\rm K} \mbox{\boldmath$B$ }_{\rm p}/B_{\rm p}$, where $v_{\rm K}=\left( 1/r \right)^{1/2}$ is the local Keplerian velocity at the disk surface, $\mbox{\boldmath$B$ }_{\rm p}$ the local poloidal magnetic field that threads the disk and $v_{\rm inj}$ is the injection velocity, an input parameter. For each time interval, a random perturbation was introduced in the velocity field at each point of the disk surface. The random perturbation was constructed, assuming a normal distribution for each component of the velocity field. This was done using the function GASDEV which appears in Numerical Recipes, 1st edition for FORTRAN, page 203. The function GASDEV gives a normally distributed deviation D with a zero mean and a unit variance. The normal distribution for the z-component of the velocity is given by

$$n\left( v_z \right) = \frac {1}{\sigma_z \sqrt{2\pi}} {\rm e}^{-\left( {v_z}^2 / 2{\sigma_z}^2 \right) },$$ (15)

The z-component of the standard deviation is given by

$$\sigma_z = 10\;v_{\rm inj}\;\frac {{\rm d}r\;{\rm d}\phi} {\sqrt{{\rm d}r^2 + {{\rm d}\phi}^2}} ,$$ (16)

where dr and d$\phi$ are the grid spacings. Similar expressions hold for $\sigma_{\rm r}$ and $\sigma_{\phi}$ for the radial and toroidal components of the velocity, using the pair $({\rm d}\phi, {\rm d}z)$ and $({\rm d}z, {\rm d}r)$, respectively. Introducing the perturbation, the velocity field at the disk surface can be expressed by the poloidal component $v_{\rm p}$ and the Keplerian component $v_{\rm K}$:

$$v_{\rm p} = v_{\rm p0} + \;\sigma\;D\;r^a,$$ (17)

$$v_{\rm K} = v_{\rm k0} + \;\sigma\;D\;r^a,$$ (18)

where, $v_{\rm p0}=v_{\rm inj}v_{\rm K0}$ and $v_{\rm k0}$ are the prescribed injection and Keplerian velocities, respectively, at the disk surface without perturbation. It is to be noted that $v_{\rm inj}$ is dimensionless. The exponent a is an input parameter which defines the radial dependence of the amplitude of the perturbation. We chose $10v_{\rm inj}$ for the maximum value of the perturbation at $r_{\rm i}$. Thus, there are six free parameters in our simulation: $\delta_{\rm i}$, $\beta_{\rm i}$, $\eta_{\rm i}$, $\mu_{\rm i}$, $v_{\rm inj}$ and a.

4 Numerical results and discussion

We used the following parameters to investigate the time evolution of a jet from a randomly perturbed Keplerian disk:

$$\left( \delta _{\rm i}, \beta _{\rm i}, \eta _{\rm i}, \mu_{\... ...{\rm inj}\right) = \left( 100.0,1.0,100.0,1.0,10^{-3}\right) ,$$ (19)

implying the initial timescales

$$\left( \tau _{\rm K},\tau _{\rm A},\tau _{\rm B}\right) = \left(1.0,7.07,707.1\right) .$$ (20)

As discussed above, we did the simulations with $v_{\rm inj}=10^{-3}$ and introduced random perturbations for given values of a in the velocity field of the disk surface. We made seven simulations, A - G. In simulation A, for comparison purposes, there is no perturbation. In simulations B - G we used a=-3/2, -1, -1/2, 0, 1/2, 1, respectively, for the random component.

Figure 2 shows the maximum poloidal velocity $v_{\rm p}$ (perturbed plus continuous velocities) for each simulation, using different values of the parameter a.

Figure 2: Radial profile of the poloidal velocity $v_{\rm p}$, in units of  $v_{{\rm K},{\rm i}}$, for each simulation, using different values of the parameter a. The constant poloidal velocity $0.001~v_{\rm K}=0.001/r^{1/2}$is added to the maximum value of the perturbed poloidal velocity 0.01 ra. The top panel shows these profiles for simulations B - E, in which the parameter a=-3/2, -1, -1/2, 0, respectively. The bottom panel shows these distributions for simulations F and G, in which the parameter a=1/2 and 1, respectively. For comparison, we show in the top panel the $10~v_{\rm p}$ profile, 10 times the unperturbed poloidal velocity.

Open with DEXTER

Similarly, Fig. 3 shows the maximum toroidal velocity (perturbed plus Keplerian velocities).

Figure 3: Radial profile of the toroidal velocity, in units of $v_{{\rm K},{\rm i}}$, for each simulation, using different values of the parameter a. The constant toroidal velocity $v_{\rm K}=1/r^{1/2}$is added to the maximum value of the perturbed toroidal velocity 0.01 ra. Simulations B - G use a=-3/2, -1, -1/2, 0, 1/2, 1, respectively. Simulations B to D show roughly the same toroidal velocity profile as simulation A, where there is no perturbation.

Open with DEXTER

Simulations were made in the domain $\left( z,r\right)=\left(80,20\right)$ with a resolution of $\left( 120,45\right)$ for 500 time units $\left( <\tau _{\rm B}\right)$. Along the z-axis, 30 zones span a distance of $10~r_{\rm i}$ outwards from the disk surface. An additional 90 ratioed zones span an additional distance of $70~r_{\rm i}$, with each subsequent zone increasing in size by a factor of 1.017. Twenty uniform zones span an initial distance of five disk radius in the radial direction. Outside of these uniform grid zones, there are another 25 zones, with each subsequent zone increasing in size by a factor 1.065.

Initial conditions in the atmosphere are shown in Fig. 4.

Figure 4: The initial conditions in the atmosphere. At the top: density distribution in units of $\rho _{\rm i}$ (Eq. (8)); at the bottom: the magnetic potential (Eq. (10)). Forty logarithmically spaced isocontour lines are shown for the density. Forty linearly spaced contours are shown for the magnetic potential. The absolute value for the steps between the isocontours of the logarithm of the density is 0.098. The toroidal magnetic field is initially zero.

Open with DEXTER

At the top of the figure we have the density distribution (Eq. (8)) and at the bottom, the initial magnetic (potential) configuration. Initially, the velocity, and the toroidal magnetic fields are zero.

The figures described below show the results of the simulations. Figures 5-7 show the distributions of the density, the poloidal and toroidal velocities, and the poloidal and toroidal magnetic fields at 500 time units for simulations A, B and G, our extrene cases. Contour lines are isocontours in the r-z plane.

Figure 5: Simulation A (without perturbation). From top to bottom: the distributions of density, poloidal velocity and magnetic field, and toroidal magnetic field and velocity for $\tau =500$ are shown. The density distribution is described by forty logarithmically spaced isocontour lines; the other quantities are described by 40 linearly spaced contours. In the graphs, the toroidal magnetic and velocity contour lines are 100 times their nominal values, in units of $B_{\rm i}$and $v_{{\rm K},{\rm i}}$, respectively. The absolute values for the steps between the isocontours of the logarithm of the density, toroidal magnetic field and toroidal velocity are 0.096, $0.0408~B_{\rm i}$ and $0.0118~v_{{\rm K},{\rm i}}$, respectively. The maximum value of the poloidal velocity is $0.63~v_{{\rm K},{\rm i}}$. We note the formation of a highly collimated steady state jet.

Open with DEXTER

When no perturbation was introduced in the disk (simulation A: Fig. 5), we obtained at $\tau =500$, the formation of a highly collimated steady state jet along the z-axis (see the density distribution at the top of the figure). The strong radial pinch of a dominant toroidal field (fourth panel from top) propagates into the corona and collimates the outflow. The maximum value for the poloidal velocity (second panel from top) is on the order of the Keplerian velocity at $r_{\rm i}$.

In simulations B (Fig.  6) - E , where a is negative or null, we observe the formation of regularly spaced structures along the z-axis with a mean separation distance of about $12~r_{\rm i}$ (see the density distribution at the top panel of the Fig. 6).

Figure 6: Simulation B (a=-3/2). From top to bottom: the distributions of density, poloidal velocity and magnetic field, and toroidal magnetic field and velocity for $\tau =500$ are shown. The density distribution is described by forty logarithmically spaced isocontour lines; the other quantities are described by 40 linearly spaced contours. In the graphs, the toroidal magnetic and velocity contour lines are 100 times their nominal values, in units of $B_{\rm i}$ and $v_{{\rm K},{\rm i}}$, respectively. The absolute values for the steps between the isocontours of the logarithm of the density, toroidal magnetic field and toroidal velocity are 0.096, $0.0312~B_{\rm i}$and $0.0174~v_{{\rm K},{\rm i}}$, respectively. The maximum value of the poloidal velocity is $0.65~v_{{\rm K},{\rm i}}$. We note the formation of regularly spaced structures along the jet axis (see the density distribution at the top).

Open with DEXTER

However, for simulation E, in which a=0, the region near the z-axis becomes less well-defined. The separation between the structures is about $11~r_{\rm i}$.

When we use a positive slope, a=1/2 for simulation F and a=1 for simulation G (Fig. 7), the formation of the regularly spaced structures occurs, but with an ensuing decollimation of the jet.

Figure 7: Simulation G (a=1). From top to bottom: the distributions of density, poloidal velocity and magnetic field, and toroidal magnetic field and velocity for $\tau =500$ are shown. The density distribution is described by forty logarithmically spaced isocontour lines; the other quantities by 40 linearly spaced contours. In the graphs, the toroidal magnetic and velocity contour lines are 100 times their nominal values, in units of $B_{\rm i}$ and $v_{{\rm K},{\rm i}}$, respectively. The absolute values for the steps between the isocontours of density (log $\rho$), toroidal magnetic field and toroidal velocity are 0.106, 0.0345 and 0.0517, respectively. The maximum value of the poloidal velocity is $0.70~ v_{{\rm K},{\rm i}}$. Observe the high decollimation of the jet (see the bottom panel for the toroidal velocity).

Open with DEXTER

The separation between the structures is about $11~r_{\rm i}$. We can see the decollimation clearly when we observe Fig. 7. This effect can be seen by an increase of the poloidal velocity (second panel from the top) in the region away from the jet axis, with a corresponding decrease of the poloidal velocity near the axis. The region near the axis becomes less and less well-defined (see the isocontours of the toroidal velocity at the bottom). We also observe this last effect in simulation E, using a=0. The decollimation is also seen by the bow lines of the poloidal magnetic field in the middle panel (compare with simulation A, where there is no perturbation (Fig. 5)). We observe the increase of the gas ejection from the outer radial regions of the disk surface (see poloidal velocity vectors in the second panel). The decollimation is more pronounced in simulation G (Fig. 7) where a=1.

In Fig. 8, the axial distributions of the density, axial and toroidal velocities and the toroidal magnetic field along a cut at $r=5~r_{\rm i}$ at 500 time units is shown for all the simulations.

	Figure 8: Simulations A - G. The density (solid line), axial velocity (dotted line), toroidal velocity (dashed line) and toroidal magnetic field (dashed-dotted line) along the z-axis at $r=5~r_{\rm i}$ are shown. The curves are normalized by their maximum values.
Open with DEXTER

Figure 9: The density (thicker solid line), pressure (thinner solid line), axial velocity (dotted line), toroidal velocity (dot-dashed line), axial magnetic field (dashed line) and toroidal magnetic field (dot-dot-dashed line) along the r-axis at a z=40.5are shown for simulation A. The following quantities are scaled: $10 \times$ density; $10^4 \times$pressure; $-1.0 \times$ toroidal velocity. The other curves denote respective nominal values.

Open with DEXTER

	Figure 10: The density (thicker solid line), the pressure (thinner solid line), axial velocity (dotted line), toroidal velocity (dot-dashed line), axial magnetic field (dashed line) and toroidal magnetic field (dot-dot-dashed line) along the r-axis at a z=40.5 are shown for simulation B. The scaled quantities are: $10 \times$ density; $10^4 \times$ pressure; $-1.0 \times$ toroidal velocity.
Open with DEXTER

These quantities are normalized to their maximum values along this cut. Observing the axial velocity distribution (dotted line), we see that the flow is rapidly accelerated along the field line from the disk to roughly $z=15~r_{\rm i}$. After this, for simulation A (left top panel), there is a "Hubble flow" constant acceleration, with v_z proportional to z. However, for simulations B - G, this constant acceleration character is modified, reflecting the pulsed nature as well as the decollimation of the jet. The effect of the decollimation of the jet in simulations E - G is seen in the decrease of the axial (dotted line) and toroidal (dashed line) velocities along the z-axis. The less defined features of the distribution along this cut for these simulations also reflect this effect. From the toroidal magnetic field distribution (dashed-dotted line), we observe the effect of torsion of the poloidal field lines due to the rotation of the disk. The structures, in simulations B - G, have their minimum mass concentration (solid line) at the maximum values of the toroidal magnetic field (dashed-dotted line), reflecting the jet constriction by this field. The shift between the toroidal magnetic field and the density in the pinched structures is quite a general result of long standing (e.g., Chan & Henriksen 1980).

Figures 9 to 11 show the radial distributions of some physical quantities in a cut at $z=40.5~r_{\rm i}$ for simulations A, B and G.

	Figure 11: The density (thicker solid line), pressure (thinner solid line), axial velocity (dotted line), toroidal velocity (dot-dashed line), axial magnetic field (dashed line) and toroidal magnetic field (dot-dot-dashed line) along the r-axis at a z=40.5are shown for simulation G. Scaled quantities: $10 \times \mbox{density}$; $10^4 \times$ pressure; $-1.0 \times$ toroidal velocity.
Open with DEXTER

From the density distribution (thicker solid line) we see that the jet is hollow (on radial scales <$r_{\rm i}$) and that no motion (outflow or inflow) is present within this hollow region. The density contrast between the interior of this region and the jet is on the order of 7 for the simulation A (no perturbation) and B (a=1). However, for simulation G, a contrast on the order of 3 shows the production of an outflow less denser.

Figure  12 shows the axial and radial mass fluxes for simulations A to G at $\tau =500$.

Figure 12: Simulations A - G: axial and radial mass fluxes (Eq. (21)). Axial (solid line) and radial (dotted line) mass fluxes at $\tau =500$. The axial mass flux, at each axial distance z, uses the transverse area $A=\pi r^2$ with $r=20~r_{\rm i}$. The radial mass flux, uses the cylindrical surface area $A=2 \pi r z$ with $r=20~r_{\rm i}$, at each axial distance z. The fluxes are given in units of $\dot{m}_{\rm i}=2 \pi r_{\rm i}^2 \rho_{\rm i} v_{{\rm K},{\rm i}}$. The scale on the vertical axis represents a tenth of the nominal value.

Open with DEXTER

The mass flux is given by

$$\dot{m}= \int \rho \, v \, {\rm d}A ,$$ (21)

where A is the area. The axial mass flux refers to the quantity of mass per unit time which passes through a circular area $A=\pi r^2$ and is calculated for each z, using $r=20~r_{\rm i}$; the radial mass flux refers to the quantity of mass per unit time which passes through the area $A=2 \pi r z$ and is shown for each z, using $r=20~r_{\rm i}$. These fluxes, given in units of $\dot{m}_{\rm i}=2 \pi r_{\rm i}^2 \rho_{\rm i} v_{{\rm K},{\rm i}}$, important for the detection of the decollimation of the jet. For simulations A - E at the time $\tau =500$, the axial mass fluxes (solid lines) at the outer axial boundary ( $z=80~r_{\rm i}$) are roughly comparable to the axial fluxes at the inner boundary ($z\sim 0$), the small difference being due to the decollimation of the jet. As the axial mass flux decreases, the radial flux can be seen to increase by roughly the same amount. In simulations B - E, with a negative or a flat radial slope (a) for the perturbation, the axial mass fluxes are little modified by the density structures and their average values do not change. The radial fluxes (dotted line) also do not change (in relation to simulation A).

When we use a positive slope for the perturbation in simulations F (a=1/2) and G (a=1), the axial mass fluxes decrease and the radial mass fluxes increase along the jet, indicating decollimation of the jet. In simulation G, the radial mass flux at the outer axial boundary ( $z=80~r_{\rm i}$) is roughly the same as the value of the axial mass flux at the inner boundary ($z\sim 0$). In these simulations, the axial mass fluxes are also modified by the density structures.

Figure 13 shows the sonic (dotted lines), Alfvénic (dashed-dotted lines) and magnetosonic Mach numbers (solid lines) for simulations A - D at $\tau =500$ in a cut at $r=5~r_{\rm i}$. Figure 14 shows these Mach numbers for simulations E - G. The sonic Mach numbers are computed, using the thermal pressure p in the sound speed.

	Figure 13: Simulations A - D: the sonic Mach numbers (dotted lines), Alfvénic Mach numbers (dashed-dotted lines) and magnetosonic Mach numbers (solid lines) along the z-axis at $\tau =500$ in a cut at $r=5~r_{\rm i}$. From simulation B to D, we show for each one the pair of panels: at the top we plot the sonic and magnetosonic Mach numbers; at the bottom we plot the Alfvénic Mach number.
Open with DEXTER

	Figure 14: Simulations E - G: the sonic Mach numbers (dotted lines), Alfvénic Mach numbers (dashed-dotted lines) and magnetosonic Mach numbers (solid lines) along the z-axis at $\tau =500$ in a cut at $r=5~r_{\rm i}$. We show for each of these simulations the pair of panels: at the top we plot the sonic and magnetosonic Mach numbers; at the bottom we plot the Alfvénic Mach number.
Open with DEXTER

The magnetosonic Mach number uses the effective sound speed which includes the thermal and Alfvénic turbulent pressure, $p_{\rm A}$. These figures show that the jet is already supersonic and super-Alfvénic near the base of the corona ($z\sim 0$) and super-magnetosonic at $z\sim 4~r_{\rm i}$. The presence of the density structures and the decollimation of the jet modify the Mach number distributions. For negative a, in simulations B - D, the sonic and magnetosonic Mach numbers show roughly the same average values. However, with null or positive values for a, these curves are modified and their average values decrease along the cut at $r=5~r_{\rm i}$ due to the decollimation of the jet. This scenario is more evident in simulation G (Fig. 14). One may observe that the sonic Mach numbers (dotted lines) have roughly the same character as the axial velocity distribution (Fig. 8: dotted lines). However, due to the decollimation, the sonic Mach numbers decrease in regions near the jet axis. They have structures which are out of phase (by roughly $90^\circ$) with respect to the density distribution (Fig. 8: solid line). The maximum values of Alfvénic Mach numbers (dashed-dotted line) are due to the decreased axial magnetic fields. Negative values correspond to a reversal of the axial magnetic field.

Beyond the Alfvén surface, where the Alfvénic Mach number is greater than unity, the toroidal field dominates the poloidal field. The dynamical importance of the magnetic field is evaluated using the quantity $\epsilon$, which is the ratio of the magnetic to the kinetic energy

$$\epsilon = \left ( B^2_{\phi} / 4 \pi \rho v^2_{\rm p} \right ).$$ (22)

The theory of an axisymmetric, stationary MHD outflow (e.g., Pelletier & Pudritz 1992) gives, from the induction equation, the value of the toroidal magnetic field at any point

$$\frac {B_\phi} {\rho} = \frac {\left ( v_\phi - r \Omega_0 \right )} {k},$$ (23)

where the constant $k = \rho v_{\rm p} / B_{\rm p}$ is the wind mass load (i.e., mass loss rate per unit time per unit magnetic flux that threads the outflowing gas flux density per unit of poloidal magnetic flux).

After some manipulation of Eqs. (23) and (22), we can evaluate the value of $\epsilon$ in the asymtotic limit $z \rightarrow \infty$ on any field line of radius $r_{\infty}$:

$$\epsilon_\infty = \left ( \rho v r^2 \right )_\infty \left (... ...\left[ \left( \Omega_\infty / \Omega_0 \right) - 1 \right ]^2.$$ (24)

From the steady state theory, when the jet conserves its material angular momentum far beyond the Alfvén surface, we have a low angular frequency $\Omega_{\infty}$ compared with the Keplerian angular frequency $\Omega_0$, which it initially had on the disc. Therefore, we ignore the ratio in the square brackets of Eq. (24). The first factor in the above equation is proportional to the total mass flux at infinity, a conserved quantity in this theory. We may, thus, replace it with the mass flux ( $\rho v_{\rm p,0}r^2$) at the surface of the disc and obtain

$$\epsilon_{\infty} \simeq \frac {1} {2^{3/2}} \left (\frac {B... ...0}} \right ) \frac {\zeta} {\left ( 1 + \zeta \right )^{3/2}},$$ (25)

where $\zeta = \left (B_{\phi,0}/4 \pi v_{\rm K,0} \right )/ k_0$.

We plot in Fig. 15 the dependence of $\epsilon _{\infty }$ (Eq. (25)) as a function of r and the radial dependence of $\epsilon$ (Eq. (22)), extracted from simulation A (with no perturbation) at 500 time units in a cut at $z=40.5~r_{\rm i}$.

	Figure 15: The quantity $\epsilon$ (solid line) (defined in Eq. (22)) is plotted along a cut at $z=40.5~r_{\rm i}$. The quantity $\epsilon _{\infty }$, obtainned from Eq. (25), is the dotted line. These two quantities are approximately equal (differing by 15%) at $r \sim 5~r_{\rm i}$.
Open with DEXTER

This is done, as well, for a cut at $r=5~r_{\rm i}$ (Fig. 16).

	Figure 16: The quantity $\epsilon$ (solid line) (defined in Eq. (22)) plotted along a cut at $r=5.0~r_{\rm i}$. The quantity $\epsilon _{\infty }$, obtainned from Eq. (25), is the dotted line. These two quantities are approximately equal (differing by 15%) from $z \sim 30~r_{\rm i}$ $r = 80~r_{\rm i}$.
Open with DEXTER

Observing Fig. 15, we see that at about $r=5~r_{\rm i}$, we obtain good agreement of $\epsilon$ (solid line) in the outer axial region of the jet ( $z=40.5~r_{\rm i}$, in this case), where $\epsilon _{\infty }$ (dotted line) obtained analytically from steady state theory. In a cut at $r=5~r_{\rm i}$ (Fig. 16), the agreement in the outer region from about $z=30~r_{\rm i}$ to $z=80~r_{\rm i}$ is clear. This result reinforces the steady state nature of our jet, for simulation A, in which we did not introduce perturbations in the disk.

Our results showed that episodic jets arise when perturbations are present. In order to understand the origin of the structures, we analyzed the individual terms of the radial component of the momentum equation (Eq. (6)). In the steady state, the momentum equation is

$$\left ( v_r \frac{\partial} {\partial r} + v_z \frac {\partia... ...tial} {\partial r} \left ( r^2 + z^2 \right )^{-1/2}}_{\rm II}$$

$$\hskip 2.55 cm\underbrace{+\; \frac{v_{\phi}^2} {r}}_{\rm III... ...c{\partial} {\partial r} \left ( r B_{\phi} \right )}_{\rm IV}$$

$$\hskip 2.55 cm \underbrace{-\; \frac{8 \pi} {\delta_{\rm i} ... ... i} \rho} B_z \frac{\partial B_r} {\partial z}} _{\rm VI}\cdot$$ (26)

The terms on the right side of Eq. (26) are as follows. The first term (I) is due to the effective pressure ( $p+p_{\rm A}$); the second term (II) is due to gravity (using $\Phi ^{\prime }=-1/\left( r^{\prime 2}+z^{^{\prime }2}\right) ^{1/2}$); the third term (III) is the centrifugal term; the fourth term (IV) is due to the toroidal magnetic field; the fifth term (V) is due to the pressure of the axial field and the sixth term (VI) is due to the axial and radial magnetic fields. Figure 17 shows the radial distribution of terms I-VI of Eq. (26) for simulation A (without perturbation) in a cut at $z=40.5~r_{\rm i}$.

Figure 17: A cut at $z=40.5~r_{\rm i}$ showing the radial distribution of each term in the steady state form of the radial momentum, Eq. (26). The solid line (S) denotes the sum of all the terms. Term I (dotted line) is due to the effective pressure ( $p+p_{\rm A}$), term IV (dot-dashed line) is due to the toroidal magnetic field and term V (dashed line) is due to the pressure of the axial field. On this scale, all the other terms make small contributions.

Open with DEXTER

Observing the sum of all of the terms (solid line: S), we see that there is a radial equilibrium of forces. For the region where $r \gtrsim 2.5~r_{\rm i}$, the pressure (dotted line: I) and the centrifugal force (dot-dashed line: III) tend to expand the jet radially. The jet is constrained in the direction of the axis by the toroidal magnetic field force (dot-dot-dashed line: IV). This behaviour holds for axial distances of $10~r_{\rm i}$ to $80~r_{\rm i}$. Thus, the toroidal component of the magnetic field acts to confine the jet along the z-axis. An equilibrium of this type is naturally subject to instabilities due to perturbations along the jet.

Chan & Henriksen (1980) studied the dynamical behaviour of axisymmetric jets in a "transverse self-similar model". They considered the cases in which the flow is confined and collimated by an external pressure or by pinching due to a toroidal magnetic field. In this "transverse self-similar model" the variables are separated in the r and z components, e.g., $B= \left ( b_r \frac {r} {R}, b_\phi \frac {r} {R}, b_z \right )$, where b_r, $b_\phi$, b_z and R are functions only of z (R is the scaling radius and is taken to be the radius of the jet boundary at z). The quantities b_r and $b_\phi$ are the values of B_rand $B_\phi$ on the jet boundary. Studying the magnetic pinching case in which the magnetic field is dynamically important ( $\epsilon_\infty > 0$, i.e., our case), with different parameters, they obtained jets which undergo quasi-periodic oscillations caused by the pinch instability due to the toroidal magnetic field. As explained by Chan & Henriksen (1980), the calculations showed that the pinch by the toroidal magnetic field is stopped by thermal and magnetic pressure and that the beam oscillates quasi-periodically. Moreover, they showed that when only magnetic piching occurs,

$$\epsilon_\infty = \frac {1} {4} \pi \alpha^2,$$ (27)

where $\alpha \equiv R_{\max} / \left (z_{\max} - z_{\min} \right )$, the maximum jet radius divided by the semiperiod. The period is measured from node to node or antinode to antinode. Empirically they found that in various cases, $\epsilon_\infty / \alpha^2$ is approximately constant, close to 2 or smaller. If Eq. (27) is strictly valid, the numerical constant should be less than unity. In our simulations, we obtained structures with periods of about $11~r_{\rm i}$ (or semiperiod of $\sim$5.5 $r_{\rm i}$). If we assume $R_{\max} = 7.5~r_{\rm i}$, considering the radial density distribution in Fig. 9 and $\epsilon_\infty \simeq 0.2$ from Fig. 15, we obtain the value $\epsilon_\infty / \alpha^2 \simeq 0.11$. This value is about 20 times smaller than the value ($\sim$2), obtained by Chan & Henriksen (1980). We note that they used $\epsilon _{\infty }$ at the few % level (e.g., $\epsilon_\infty \sim 2.3\%$), whereas we obtained $\epsilon_\infty = 20\%$, assuming $R_{\max} = 7.5~r_{\rm i}$.

In an interesting analytic work, Bogoyavlenskij (2000) reported some new exact solutions of the plasma equilibrium equations that model astrophysical jets in the comoving reference frame. The equilibrium solutions depend on arbitrary parameters which can be adjusted so that the topology of the magnetic field becomes complex. The fields (e.g., Fig. 1 in Bogoyavlenskij 2000) are a combination of nested cylindrical magnetic surfaces and nested toroidal surfaces. The innermost tori are circular magnetic axes that are linked to the neighboring closed magnetic field lines. Certain solutions are quasiperiodic in z, i.e., the magnetic field shows structures with similar, but not equal, periods. This implies that the magnetic field never repeats in the z direction, but can have a structure arbitrarily close to the initial data. Their "winding pattern" changes continuously with z and does not repeat itself.

5 Summary and conclusions

We made numerical 2.5-dimensional simulations of astrophysical jets from Keplerian disks, subject to a random perturbation in the velocity field. The simulations were made using different values for the parameter a, which defines the radial dependence of the perturbation ( $\alpha \; r^a$). In simulation A, no perturbation was introduced, while in simulations B - G, we used a=-3/2, -1, -1/2, 0, 1/2, 1, respectively.

For simulation A, where there is no perturbation, we obtained a steady state, highly collimated jet along the disk axis after a time $\tau =500$. The maximum axial velocity was of the order of the Keplerian velocity at the innermost radius of the accretion disk $r_{\rm i}$. In a cut at $r=5~r_{\rm i}$, the jet is already supersonic and super-Alfvénic near the base of the corona and super-magnetosonic at $z\sim 4~r_{\rm i}$.

The results common to simulations B - G are as follows:

1.

The formation of regularly spaced structures with a separation of roughly $11~r_{\rm i}$ along the jet axis, is observed,

2.

The axial mass fluxes are litle modified by the density structures, and

3.

The presence of the density structures modify the Mach number distibutions along the z-axis.

When we used a negative radial slope (a<0) for the perturbation (simulations B - D), we obtained the formation of a collimated jet where we observe the following:

a. The highest axial velocities are localized near the jet axis at about $r=5~r_{\rm i}$;

b. The axial mass fluxes at $z=80~r_{\rm i}$ are comparable to the axial fluxes at $z\sim 0$, their average values remaining roughly the same, as compared with the case for no perturbation (simulation A); and

c. The radial mass fluxes are similar to those of simulation A (where there was no perturbation).

The decollimation of the jet occurs when we use a positive radial power dependence (a > 0) for the perturbation. Simulations E and G show this effect where we observed that

a. The axial velocity increases in the regions far from the jet axis and decreases near the axis;

b. Decollimation modifies the sonic Mach number distributions with decreasing average values in the regions near the jet axis;

c. The region near the jet axis becomes less and less well-defined;

d. The axial mass flux decreases, while the radial mass flux increases along the jet axis (in simulation G with a=1, the radial mass flux at the outer axial boundary ( $z=80\,r_{\rm i}$) assumes the value of the axial mass flux at the inner boundary ($z\sim 0$)).

These simulations, introducing random perturbations in the accretion disk, produced a jet in which regularly spaced structures are present with periods of about $11\,r_{\rm i}$. Using this same numerical approuch, Ouyed & Pudritz (1999) obtained episodic jets introducing continuously the gas from the disk with velocities smaller than the $10^{-3} v_{{\rm K},{\rm i}}$ which we used. For injection velocities equal to this value, they obtained a highly collimated steady state jet with features similar to our jet, with no perturbation in the velocity field of the disk. When we introduced perturbations with specific features in our steady state jet, our results showed regular structures with separation distances ($\sim$11 $r_{\rm i}$) similar to the case when OP99 used a smaller velocity of $10^{-4}v_{{\rm K},{\rm i}}$. Therefore, introducing the perturbations we obtained the formation of strutuctures for the case in which that autors obtained a steady state jet.

We simulated a region very near the central object ( $20~r_{\rm i} \times 80~r_{\rm i}$). It is interesting to evaluate the spatial dimensions involved in our simulations. For a protostar of mass $0.5 ~M_{\odot}$, assuming that $r_{\rm i} = 0.04$ AU for the inner radius of the disk, $80~r_{\rm i}$ is equivalent to $\sim$10^-5 pc. For a black hole of $10^8~ M_{\odot}$ (the Schwartzschild radius is $\sim$1 AU), assuming $r_{\rm i} = 20.6$ AU, $80~r_{\rm i}$ is equivalent to $\sim$10^-2 pc. A distance of separation of $\sim$11 $r_{\rm i}$ for the structures corresponds, for the protostar case, to $\sim$10^-6 pc and for the black hole case, to $\sim$10^-3 pc. These scales are much smaller than the distances between the Herbig-Haro objects observed in the protostellar jets (about 10^-2 pc). They are also smaller than the distances observed between the emission knots, frequently found in AGN jets. We are, therefore, obtaining structures very near the formation region of the jets, involving scales not yet resolved by current instruments.

Acknowledgements

B.F.V. acknowledges the financial support of the Brazilian agency FAPESP under grant 95/99122. V.J.P. and R.O. would like to thank the Brazilian agencies CNPq and FAPESP for partial support. We would like to thank LCCA-USP for granting the use of the CRAY J-90 computer in our simulations and especially the staff for their invaluable help. R. Ouyed and especially C. Fendt contributed to this paper with encouraging and valuable discussions, for which we are grateful. The authors also would like to thank the project PRONEX (No. 41.96.0908.00) for partial support.

References

 
• Bell, A. R., & Lucek, S. G. 1995, MNRAS, 277, 1327 In the text NASA ADS  
• Bogoyavlenskij, O. I. 2000, Phys. Rev. Lett., 84, 9, 1914 In the text  
• Chan, K. L., & Henriksen, R. N. 1980, ApJ, 241, 534 In the text NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1997a, ApJ, 482, 712 (OP97a) NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1997b, ApJ, 484, 794 NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1999, MNRAS, 309, 233 (OP99) In the text NASA ADS  
• Königl, A., & Pudritz, R. 1999, Protostars and Planets IV, ed. V. Manning, A. Boss, & S. Russel (U. Arizona Press: Tucson) In the text  
• Pelletier, G., & Pudritz, R. E. 1992, ApJ, 394, 117 In the text NASA ADS  
• Shibata, K., & Uchida, Y. 1986, PASJ, 38, 631 In the text NASA ADS  
• Shu, F., Najita, H., Shang, Z., & Li, Z.-Y. 1999, Protostars and Planets IV, ed. V. Manning, A. Boss, & S. Russel (U. Arizona Press: Tucson) In the text  
• Stone, J. M., & Norman, M. L. 1992a, ApJS, 80, 753 NASA ADS  
• Stone, J. M., & Norman, M. L. 1992b, ApJS, 80, 791 NASA ADS  
• Stone, J. M., & Norman, M. L. 1994, ApJS, 433, 746 In the text  
• Ustyugova, G. V., Koldoba, A. V., Romanova, M. M., Chechetkin, V. M., & Lovelace, R. V. E. 1995, ApJ, 439, L39 In the text NASA ADS