BIPOLAR JETS LAUNCHED FROM MAGNETICALLY DIFFUSIVE ACCRETION DISKS. I. EJECTION EFFICIENCY VERSUS FIELD STRENGTH AND DIFFUSIVITY

For our numerical simulations, we apply the MHD code PLUTO 3.01 (Mignone et al. 2007), solving the conservative, time-dependent, resistive, inviscous MHD equations, namely, for the conservation of mass, momentum, and energy,

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \bm {v})=0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial (\rho \bm {v})}{\partial t} + \nabla \cdot \left[\bm {v} \rho \bm {v} - \frac{\bm {B} \bm {B}}{4\pi } \right] + \nabla \left[ P + \frac{B^2}{8\pi } \right] + \rho \nabla \Phi = 0, \end{equation} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{\partial e}{\partial t} + \nabla \cdot \left[ \left(e + P + \frac{B^2}{8\pi } \right)\bm {v} - \left(\bm {v} \cdot \bm {B} \right)\frac{\bm {B}}{4\pi } + (\bar{\bar{\eta }} : \bm {j}) \times \frac{\bm {B}}{4\pi } \right]\nonumber\\ &&\quad= - \Lambda _{\rm cool}. \end{eqnarray} \tag{ 3 }$$

Here, ρ is the mass density, $\bm {v}$ is the velocity, P is the thermal gas pressure, $\bm {B}$ is the magnetic field, and Φ = −GM/R is the gravitational potential of the central object of mass M, with the spherical radius $R = \sqrt{r^2 + z^2}$. In general, the magnetic diffusivity is defined as a tensor $\bar{\bar{\eta }}$ (see Section 2.5). The evolution of the magnetic field is described by the induction equation,

$$\begin{equation} \frac{\partial \bm {B}}{\partial t} - \nabla \times (\bm {v} \times \bm {B} - \bar{\bar{\eta }}: \bm {j}) = 0 \end{equation} \tag{ 4 }$$

with the electric current density $\bm {j}$ given by Ampère's law $\bm {j} = (\nabla \times \bm {B}) / 4\pi$. The cooling term Λ can be expressed in terms of Ohmic heating Λ = gΓ, with $\Gamma = (\bar{\bar{\eta }}: \bm {j}) \cdot \bm {j}$, and with g measuring the fraction of the magnetic energy that is radiated away instead of being dissipated locally. For simplicity, here we adopt g = 1. The gas pressure follows an equation of state P = (γ − 1)u with the polytropic index γ and the internal energy density u. The total energy density is

$$\begin{equation} e = \frac{P}{\gamma - 1} + \frac{\rho v^2}{2} + \frac{B^2}{8\pi } + \rho \Phi. \end{equation} \tag{ 5 }$$

Our simulations are performed in axisymmetry applying cylindrical coordinates. The CENO3 algorithm as third-order interpolation scheme is used for spatial integration (Del Zanna & Bucciantini 2002) together with a third-order Runge–Kutta scheme for time evolution and an HLL Riemann solver. For the magnetic field evolution, we apply the constrained transport method (FCT) ensuring solenodality $\nabla \cdot \bm {B} =0$.

Throughout the paper distances are expressed in units of the inner disk radius r_i, while p_{d, i} and ρ_{d, i} denote the disk pressure and density at this radius, respectively. The index "i" refers to a number value at the inner disk radius at z = 0 and time t = 0.

In Appendix A, we show for comparison the astrophysical scaling for YSO jets and AGN jets. Typically, r_i < 0.1 AU for YSO and r_i < 10 Schwarzschild radii for AGNs. Naturally, we cannot treat any relativistic effects of AGN jets with our non-relativistic setup. Velocities are measured in units of the Keplerian speed v_{K, i} at the inner disk radius. Time is measured in units of t_i = r_i/v_{K, i}, which can be related to the Keplerian orbital period τ_{K, i} = 2πt_i. Pressure is given in units of p_{d, i} = ²ρ_{d, i}v²_{K, i}. The magnetic field is measured in units of B_i = B_{z, i}. As usual we define the aspect ratio of the disk as the ratio of the isothermal sound speed to the Keplerian speed, both evaluated at disk mid-plane, ≡ c_s/v_K.⁵ For a more details see Appendix A.

We define initial conditions following a setup applied by other authors previously—a magnetically diffusive accretion disk is prescribed in sub-Keplerian rotation, above which a hydrostatic corona in pressure balance with the disk is located (Zanni et al. 2007; Murphy et al. 2010). The coronal density is chosen several orders of magnitude below the disk density, thus implying an entropy and a density jump from disk to corona. However, contrary to previous authors (Zanni et al. 2007; Tzeferacos et al. 2009; Murphy et al. 2010), we do not apply an initial vertical velocity profile in the disk and only prescribe a radial velocity profile. We have seen that for our long-term simulations the whole disk system adjusts to a new dynamical equilibrium which does not depend on the vertical profile of the initial velocity distribution.

2.3.1. Initial Disk Structure

We prescribe an initially geometrically thin disk with = H/r = 0.1 which is itself in vertical equilibrium between thermal pressure and gravity.

We follow the standard setup employed in a number of previous papers (Zanni et al. 2007; Murphy et al. 2010).

As initial disk density distribution we prescribe

$$\begin{equation} \rho _{\rm d} = \rho _{\rm d,i} \left(\frac{2}{5\epsilon ^2}\frac{r_{\rm i}}{r} \left[\frac{r}{R}- \left(1-\frac{5\epsilon ^2}{2}\right) \right] \right)^{3/2} \end{equation} \tag{ 6 }$$

(Murphy et al. 2010), while the initial disk pressure distribution follows:

$$\begin{equation} P_{\rm d} = P_{\rm d,i} \left(\frac{\rho _{\rm d,i}}{\rho _{\rm d}}\right)^{5/3}. \end{equation} \tag{ 7 }$$

The disk is set into slightly sub-Keplerian rotation accounting for the radial gas pressure gradient and advection,

$$\begin{equation} v_{\phi , \rm d}(r) = \sqrt{1-\frac{5\epsilon ^2}{2}}\,\sqrt{\frac{\it GM}{r}} \end{equation} \tag{ 8 }$$

following Murphy et al. (2010), but neglecting viscosity.

In our setup the initial poloidal velocity is imposed by hand following the prescription of Zanni et al. (2007):

$$\begin{equation} v_{r,\rm d}(r) = \epsilon ^2 \sqrt{2 \mu } \left(\frac{5}{4}+\frac{5}{3m^2}\right)\sqrt{\frac{\it G M}{r}} = \frac{r}{z}\,v_{z, \rm d}. \end{equation} \tag{ 9 }$$

Simulations with initially vanishing disk accretion resulted in the same asymptotic inflow–outflow evolution. Starting with a disk inflow as in Equation (9), the inflow–outflow structure is established on a shorter timescale.

2.3.2. Structure of the Coronal Region

Above the disk, we define an initially hydrostatic density and pressure stratification (the so-called corona),

$$\begin{equation} \rho _{\rm c}=\rho _{\rm a,i}\left(\frac{r_{\rm i}}{R}\right)^\frac{1}{\gamma -1}\!,\, P_{\rm c}=\rho _{\rm a,i}\frac{\gamma -1}{\gamma }\frac{\it GM}{r_{\rm i}}\left(\frac{r_{\rm i}}{R}\right)^\frac{\gamma }{\gamma -1}\!. \end{equation} \tag{ 10 }$$

The parameter δ ≡ ρ_c/ρ_d quantifies the initial density contrast between disk and corona. We adopt δ = 10⁻⁴.

2.3.3. Magnetic Field Distribution

The initial magnetic field is prescribed by the magnetic flux function ψ following Zanni et al. (2007):

$$\begin{equation} \displaystyle \psi (r,z) = \frac{3}{4} B_{z,i} r_{\rm i}^2 \left(\frac{r}{r_{\rm i}}\right)^{3/4} \frac{m^{5/4}}{(m^2 + (z/r)^2)^{5/8} }. \end{equation} \tag{ 11 }$$

Here, B_{z, 0} measures the vertical field strength at (r = r_i, z = 0). The magnetic field components are calculated by rB_z = ∂Ψ/∂r and rB_r = ∂Ψ/∂z.

The computational domain spans a rectangular grid region applying a purely uniform spacing in the radial direction and a uniform-plus-stretched spacing vertically (Figure 2). We have run simulations in two different resolutions (see Table 1). For the low-resolution simulations, the domain extends to (96 × 288) inner disk radii r_i on a grid of (1500 × 3200) cells, resulting in a resolution of Δr = 0.064. For the high-resolution simulations, the resolution is increased to 0.025 at the cost of being limited to a smaller domain of (50 × 180)r_i. For all simulations presented here, a uniform grid is used for the magnetically diffusive disk.⁶

Zoom In Zoom Out Reset image size

For the typical disk model applied with = 0.1 the low resolution provides only 2 grid cells per disk scale height at the inner disk radius, while for somewhat larger radii, say at r = 5, we have about 10 cells per disk scale height. In the high-resolution runs, the disk resolution is increased by a factor 2.5. As discussed by Murphy et al. (2010), the low resolution barely resolves the jet-launching disk surface layer of the very inner disk, where steep gradients in density, pressure, and diffusivity are present. Here, numerical diffusivity is supposed to play a role and we have devoted a section to investigating this effect.

In order to follow the jet outflow over long distances and in order to provide a sufficient mass reservoir for disk accretion a large domain would be desirable. In order to resolve the wind-launching area, a high resolution would be required. However, as mentioned before, PLUTO in its current version does not allow for diffusive MHD simulations in a stretched grid. Furthermore, we experience that for large domains together with stretched grids (implying an elongated shape of the outer grid cells) the code had difficulties solving the conservative equation and the simulations crashed.

For the boundary conditions, axisymmetry on the rotation axis and equatorial symmetry for the disk mid-plane are imposed. At the upper z boundary, we use the standard PLUTO outflow boundary condition (zero gradient), but at the outer radial boundary condition, we apply a modified outflow condition as derived by Porth et al. (2011) in order to avoid artificial collimation. Most essential is the internal boundary enclosing the origin which we call a sink.

2.5.1. Internal Boundary—A Central Sink

2.5.2. Outer Boundary Conditions

A dissipative effect is required for steady-state accretion in order to allow matter to diffuse across the magnetic field threading the disk. Magnetic diffusivity allows the mass flux to cross the field lines and thus allows for accretion along the disk. We consider the magnetic diffusivity to be turbulent in nature, or "anomalous," and, thus, cannot be calculated self-consistently in our setup. The origin of the turbulent magnetic diffusivity is usually referred to being caused by the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998) in a moderately magnetized disk.⁸

Nevertheless, we argue that as we load turbulently diffusive disk plasma into the outflow, it is natural to assume that the outflowing material is initially turbulent as well. Therefore, we have the option of applying a diffusive scale height H_η = _ηr larger than the thermal scale height H = r of the disk, supposing that the turbulent disk material lifted into the outflow will remain turbulent for a few more scale lengths until turbulence decays. It is clear that the strong magnetization in the upper layers of the disk will prevent generation of turbulence by, e.g., the MRI, and will also quench the turbulence in the outflowing material. Without going into detail we may estimate the timescale for the decay of the turbulence potentially existing in the outflowing material as follows. The decay of (helical) MHD turbulence follows a power law E_M ∼ t⁻ⁿ, where E_M is the magnetic energy and the power-law index depends on further details but is in the range of n = 1/2–2/3 (Brandenburg & Nordlund 2011). Simulations indicate that MHD turbulence decays as fast as the hydrodynamic turbulence on a few eddy turnover times (Cho et al. 2002). With τ_η ⩾ H/c_S = r/c_S = r^1/2 (in code units), we find a timescale for the decay of MHD turbulence launched from the disk into the outflow of, e.g., τ_η(r = 5) ⩾ 3, thus in the range of half an orbital period. This must be compared to the kinematic timescale for the jet launching. If we consider the propagation of the initially slow disk wind over a few disk pressure scale heights τ_wind(r = 5) ≃ H/v_wind ≃ 0.5/0.05 = 10, we find that this time is of the order of a few eddy turnovers. Also, we will later see that jet launching is a rapid process with an outflow having been established after a few dynamical time steps.

Given the essential role of the magnetic diffusivity for wind/jet launching, we decided to investigate the launching process for different strengths and different scale heights of diffusivity. We have run simulations for a variety of combinations of parameter values, but for the sake of comparison, in this paper, we show results for diffusive scale height values, _η = 0.1, 0.2, 0.3, and 0.4.

Formally, we apply an α prescription (Shakura & Sunyaev 1973), similar to previous works (see, e.g., Zanni et al. 2007). We assume the diffusivity tensor to be diagonal with the non-zero components

$$\begin{equation} \eta _{\rm \phi \phi } \equiv \eta _{\rm p}, \quad \eta _{\rm rr}=\eta _{\rm zz} \equiv \eta _{\phi }. \end{equation} \tag{ 12 }$$

Here, we have defined a poloidal magnetic diffusivity η_p ≡ η_ϕϕ = η_{p, i}f(r, z) and the toroidal magnetic diffusivity η_ϕ ≡ η_rr = η_zz = η_{ϕ, i}f(r, z), respectively,⁹ with a function f(r,z) describing the diffusivity profile.

As is known from the literature, an anisotropic magnetic diffusivity is required to obtain stationary solutions (Wardle & Königl 1993; Ferreira & Pelletier 1995; Ferreira 1997). Most of the numerical simulations followed that approach and did find a stationary state from their simulations as well (see, e.g., Zanni et al. 2007; Tzeferacos et al. 2009; Murphy et al. 2010). Note that simulations of Casse & Keppens (2002, 2004) apply an isotropic diffusivity and also reach a steady state.

We define an anisotropy parameter by the ratio of the toroidal to poloidal diffusivity components, χ = η_{ϕ, i}/η_{p, i}. For the diffusivity function f several options can be considered. For the simulations in this paper, we apply

$$\begin{equation} f = v_{\rm A}(r,z=0) H \exp {\left(-2\frac{z^2}{H^2_{\eta }(r)}\right)}, \end{equation} \tag{ 13 }$$

where η_{p, i} and η_{ϕ, i} govern the strength of magnetic diffusivity. The Alfvén speed $v_{\rm A}(r,z=0) = {B_{\rm {z}}(r,z=0)}/{\sqrt{\rho (r,z=0)}}$ and the disk thermal scale height H(r) = c_s(r, z = 0)/Ω_K(r, z = 0) are both calculated along the mid-plane. The parameter _η measures the scale height of the disk magnetic diffusivity, similar to the hydrodynamic disk scale height parameter _η = H_η/r. For the present paper, we decided not to evolve the diffusivity profile in time. We find that since both the disk scale height H and the Alfvén speed v_A vary as the disk evolves, the strength of diffusivity may vary substantially from the initially prescribed value (up to a factor 10). This highly nonlinear feedback may disturb the progress of our simulation such that an artificially high diffusivity may be derived which will affect the numerical time stepping. A constant-in-time prescription of diffusivity simplifies our aim of disentangling the governing physical processes involved in jet launching.

In Paper II, presenting truly bipolar jet launching, we will discuss further models for the magnetic diffusivity.

We stress again the importance of providing the reader with the actual number values of the magnetic diffusivity applied.¹⁰ If we consider a η_{ϕ, i} ∼ 1 together with a maximum Alfvén speed in the disk of about 10⁻² (located at the inner disk radius), this gives a maximum disk diffusivity of about η(r, z)_i ≃ 0.01–0.1, through all the simulations. Figure 3 shows the distribution of the magnetic diffusivity. The diffusivity is mainly concentrated in the disk, increases with disk radius (consistent with a decreasing temperature or ionization degree with radius), and decreases with height resembling the transition from a cool disk to an ideal MHD outflow.

Zoom In Zoom Out Reset image size

The simulations are governed by a set of non-dimensional parameters for the following physical quantities.

• 1.  
  The magnetic field strength defined by plasma-β, and its initial geometry defined by the "bending" parameter m.
• 2.  
  The magnetic diffusivity, with the three parameters η_{ϕ, i}, η_{p, i}, and _η governing its strength, anisotropy, and the diffusivity scale height.
• 3.  
  The initial density contrast between disk and corona δ.
• 4.  
  The initial disk thermal scale height parameter = H/r.

An overview of these parameters is shown in Table 1 (first half) for the simulations presented in this paper. The second part of Table 1 shows the main dynamical quantities resulting from our simulations and will be discussed below. Table 2 shows similar numbers for the simulations of different diffusive scale heights.

Table 1. Input Parameters and Derived Dynamical Parameters of our Simulation Runs Following a Magnetic Diffusivity Profile Equation (13) with a Scale Height of the Disk Diffusivity _η = 0.4

	Case 1	Case 2	Case 3	Case 4	Case 5	Case 6	Case 7	Case 8	Case 9	Case 10
Δr	0.064	0.064	0.064	0.064	0.064	0.025	0.025	0.064	0.064	0.064
Δz	0.066	0.066	0.066	0.066	0.066	0.025	0.025	0.066	0.066	0.066
η	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4	0.4
ηp, i	0.03	0.01	0.15	0.09	0.03	0.03	0.03	0.03	0.03	0.03
ηϕ, i	0.09	0.03	0.09	0.27	0.01	0.09	0.09	0.09	0.09	0.09
χ	3	3	3/5	3	1/3	3	3	3	3	3
β	10	10	10	10	10	5000	10	50.0	250	500
rjet, z = 180	43.8	73.8	24.46	28.57	33.19	...	...	47.8	50.0	34.13
rjet, z = 60	25.01	38.41	15.37	17.72	25.02	21.69	19.6	27.54	24.18	17.22
rl	3.8	8.2	3.0	4.2	3.0	...	3.5	2.8	0.8	0.7
vjet, z = 170	0.4	0.23	0.85	0.56	0.8	0.47	0.42	0.56	0.49	0.27
vjet, z = 280	0.45	0.3	0.92	0.61	0.9	...	...	0.5	0.89	0.22
ζjet, 1a	2.37	0.013	5.28	5.1	2.58	0.89	4.6	1.56	0.32	2.08
$\dot{M}_z$	0.0005	6 × 10−6	0.0002	0.0002	0.0007	5 × 10−5	0.0008	0.0001	0.0004	0.0002
ζjet, 2	2.90	0.3	4.45	6.7	2.28	...	...	2.46	0.42	1.84
$\dot{M}_z$	0.001	0.0002	0.0003	0.0004	0.001	...	...	0.0003	0.0006	0.0004
$\dot{M}_{\rm acc}$b	0.015	0.022	0.0038	0.0035	0.022	0.001	0.013	0.005	...	...
$\dot{M}_{\rm ejec}$	0.008	0.016	0.001	0.001	0.007	0.005	0.009	0.004	...	...
${\dot{M}_{\rm ejec}}/{\dot{M}_{\rm acc}}$	0.5	0.7	0.2	0.2	0.31	...	0.6	...	...	...
ξ	0.3	0.5	0.09	0.09	0.1	...	0.39	...	...	...
$\dot{J}_{\rm kin}$	0.008	0.005	0.002	0.002	0.027	0.005	0.005	0.008	...	...
$\dot{J}_{\rm mag}$	0.03	0.01	0.0025	0.005	0.034	0.0	0.016	0.008	0.0	0.0
$\dot{J}_{\rm tot}$	0.038	0.015	0.004	0.007	0.06	0.005	0.021	0.017	...	...
Lkin	1 × 10−4	9 × 10−6	1.3 × 10−4	7.4 × 10−5	4 × 10−4	5.5 × 10−6	7.1 × 10−5	3.7 × 10−5	1.9 × 10−4	10−6

Notes. Displayed are the input parameters of grid resolution Δr, Δz, magnetic poloidal and toroidal diffusivity η_{p, i} and η_{ϕ, i}, the diffusivity anisotropy χ, and the plasma-beta parameter β. The following outflow parameters are derived from our numerical simulations: the asymptotic jet radius r_jet, calculated as mass flux weighted (see Equation (17)), the launching radius r_l, the typical speed of the outflow, v_jet, the average collimation degree ζ_jet (see Equation (18)), the mass accretion rate $\dot{M}_{\rm acc}$ and mass ejection rate $\dot{M}_{\rm ejec}$, the ejection index ξ, the diffusive scale height _η, the kinetic and magnetic angular momentum losses by the outflow, $\dot{J}_{\rm kin}$ and $\dot{J}_{\rm mag}$, respectively, as well as the total angular momentum loss $\dot{J}_{\rm tot}$, and its asymptotic kinetic luminosity $L_{\rm kin} = 0.5\,x \dot{M}_z v_{\rm jet,z}^2$. All parameters are given in code units. ^aThe collimation degree is measured from two different areas, ζ_{jet, 1} enclosing 0 ⩽ r ⩽ 40, 0 ⩽ z ⩽ 160, and ζ_{jet, 2} with the area 0 ⩽ r ⩽ 80 and 0 ⩽ z ⩽ 250. ^bThe mass flux values for the two simulation cases 9 and 10 with high β_i are omitted due to their highly perturbed behavior.

Download table as:  ASCIITypeset image

Table 2. Comparison Between Different Diffusive Scale Heights

	Case 1	Case 11	Case 12	Case 13
η	0.4	0.1	0.2	0.3
Δr	0.064	0.064	0.064	0.064
Δz	0.066	0.066	0.066	0.066
ηp, i	0.03	0.03	0.03	0.03
ηϕ, i	0.09	0.09	0.09	0.09
χ	3	3	3	3
β	10	10	10	10
rjet, z = 180	43.8	73.3	56.4	48.4
rjet, z = 60	25.01	45.8	39.2	34.4
rl	3.8	7.2	7.2	6.2
vjet, z = 170	0.46	0.35	0.33	0.41
vjet, z = 280	0.53	0.46	0.41	0.47
$\dot{M}_{\rm acc}$	0.015	0.024	0.022	0.018
$\dot{M}_{\rm ejec}$	0.005	0.015	0.015	0.01
$\dot{J} _{\rm kin}$	0.01	0.005	0.009	0.01
$\dot{J}_{\rm mag}$	0.033	0.015 ± 0.01	0.02 ± 0.003	0.028

Notes. Shown are some of the physical properties mentioned in Table 1 for different diffusive scale height runs which have been carried out up to t = 2000.

Download table as:  ASCIITypeset image

Our reference simulation case 1 is carried out up to t = 5000 dynamical time, corresponding to 796 inner disk orbital periods. It is thus one of the longest simulations considering MHD jet launching. This huge time period corresponds, however, only to three rotations at a radius r = 40 r_i and to (only) one rotation at the outermost part of the disk r = 96 r_i. Consequently, while the inner part of the disk, and thus the outflow evolving from that part, has reached a steady-state situation, the outer disk corona is still dynamically evolving. We will therefore focus our discussion mainly on the inner disk areas.

The time evolution of density together with the magnetic field is shown in Figure 4. Note that we do not show field lines integrated from certain footpoint radii, but magnetic flux contours. This allows us to visualize the diffusion and advection of the magnetic field as a consequence of the disk evolution.

Zoom In Zoom Out Reset image size

The magnetic diffusivity distribution determines the coupling between the plasma and magnetic field, and, thus, affects the mass loading into the outflow. Since the high diffusivity is in the outer part of the disk, the coupling is weaker, and, thus, the mass loading less efficient. As a general result we observe a continuous and robust outflow launched from the inner part of the disk, expanding into a collimated jet.

Figure 5 shows the poloidal velocity distribution in jet and disk, and the normalized poloidal velocity vectors indicating the direction of the mass fluxes.

Zoom In Zoom Out Reset image size

The outflow is accelerated to superfast-magnetosonic speed (see the magnetosonic surfaces indicated in Figure 4)

A bow shock develops at the interface between the outflow and the surrounding hydrostatic corona. As the outflow propagates, the initial corona is swept out of the computational domain together with the bow shock. While the interaction with the ambient gas plays a role within the first evolutionary steps, the long-term evolution of the outflow—its acceleration and collimation—is purely determined by the force balance within the outflow and the physical conditions of the launching region.

As jets are launched from disks, the disk evolution is itself an essential part of the jet evolution and must be carefully considered. It is expected that the jet mass flux would somewhat correspond to the disk accretion rate (which would also depend on the actual mass of the disk).

Our simulations show that the disk evolution is complex, with smooth accretion phases followed by more turbulent stages. The time evolution of the accretion stream ρv_r is shown in Figure 6. We see that accretion starts first at small radii and fully establishes after t = 1000. Accretion shocks and vortices may destroy the quasi-steady-state situation. Close to the outer disk radius, mass is lost through the outer boundary. At late evolutionary stages we see a complex pattern in the ρv_r distribution indicating small shock and turbulent motion within the disk.

Zoom In Zoom Out Reset image size

Nevertheless, the negative mass fluxes in Figure 6, together with the negative velocity vectors in Figure 5 clearly indicate that accretion is established throughout the whole disk.

We now discuss the time evolution of the disk mass for our reference simulation. In order to measure the disk mass, we need to carefully define the control volume defining a corresponding "disk surface." Here, we consider as disk surface the location where the mass density at each radius falls bellow 10% of the mid-plane density (because of the vortex motion we could not use the negative velocity as a proxy for accretion). Measuring the mass contained in subsequent disk rings dM/dr = 2πrρdz, we see that most of the mass is indeed stored within the outer disk areas (Figure 7, top). We can identify two essential factors affecting the long-term disk evolution. First, the mass reservoir in the outer disk has decreased, mainly due to outwards mass loss across the outer disk boundary. Second, the mass of the inner jet-launching disk has also decreased, due to both outflow activity and accretion into the sink.

Zoom In Zoom Out Reset image size

Comparing the mass content of the disk at the initial and final stage, we find that the total disk mass decreases from M ≃ 520 at t = 0 to M ≃ 320 at t = 5000 (in code units). Integrating the mass loss by accretion and ejection into the jet (see Figure 7) over 5000 dynamical time steps, we find a total mass loss of M ≃ 75. The difference, i.e., M ≃ 125, is the mass which is lost from the outermost disk in vertical and radial directions.

In other words, a substantial amount of the disk mass is lost from the very outer part of the disk and does not directly influence the jet formation. The disk mass which is lost from the inner part is partly lost as disk wind/jet, partly accreted into the sink and partly replenished by accretion from outer radii.

The time evolution of the disk mass gives a similar picture (Figure 7, bottom). Integrated over the whole domain, the disk loses about 38% of its mass until t = 5000 in the reference simulation case 1. If we decrease the integration domain (outer radius r = 50), the inner disk loses less mass while part of its mass is being accreted from outer disk radii. In general this implies that up to time t = 5000, there is still sufficient mass to support the accretion process, but also that the disk mass shows a considerable decrease which may have changed the internal disk properties. Thus, for an even longer-lasting simulation one would have to invent another mass source for the disk accretion (e.g., by a physical mass inflow boundary condition at the outer disk radius properly taking into account angular momentum transfer, or a floor density distribution in the disk).

Along with the hydrodynamic disk evolution, the magnetic field distribution evolves, subject to the competing processes of advection and diffusion. This also changes the launching conditions, as the overall profiles of plasma-β, resp. magnetization, are affected. Figure 8 displays the time evolution of the magnetic flux surface Ψ = 0.1 for reference simulation case 1, initially rooted at the radius r = 2.0. We see that this flux surface is first moving (diffusing) outward "driven" by magnetic pressure gradient and tension, and then, once disk accretion has established again moves inward being advected with the mass flux. After about t = 3000 both processes balance and a quasi-steady-state situation in the system evolution is reached. Considering simply flux conservation, we may estimate the change in poloidal magnetic field strength and the corresponding change in magnetic diffusivity and plasma-β.

Zoom In Zoom Out Reset image size

Estimating an average field strength ${\bar{B}}_{\rm p}$ and a corresponding magnetic flux $\Psi \simeq {\bar{B}}_{\rm p} r^2$, flux conservation tells us that Ψ_{t = 1000} ≈ ψ_{t = 50}. Therefore, ${\bar{B}}_{\rm p,t=1000} \approx 10 \;B_{\rm p,t=50}$, since this flux surface (e.g., Ψ = 0.1) is now rooted at a different radius. With the increased poloidal magnetic field strength the Alfvén speed v_A ∼ B_p increases, and, thus, also the magnetic diffusivity parameterized as η ∼ v_A. Thus, we expect that in previous work (e.g., Zanni et al. 2007; Tzeferacos et al. 2009) the actual (time-dependent) value of the magnetic diffusivity may differ substantially from its initial value. We believe that this will impact the derived mass fluxes. A similar estimate can be made for the plasma-β or magnetization. Since β ∼ B⁻², the increase in plasma-β is by a factor 100, which is also expected to affect the jet formation severely. We will come back to this point later when we compare different parameter runs.

The main goal of this paper is to investigate what fraction of the mass accretion is loaded into the outflow, and how the mass loading is affected by various disk parameters. The derived mass ejection-to-accretion ratio can be an important ingredient for studies of AGN or YSO feedback, or another similar self-regulating outflow scenario.

Efficient accretion is feasible only if angular momentum is sufficiently removed. Since we do not consider physical viscosity in our simulations, angular momentum removal from the disk is mostly accomplished by the torque of the magnetized disk wind. Therefore, in our simulations disk accretion can only work if a disk wind has been established. In the following we discuss the inflow and outflow rates in our disk–jet simulations, concentrating first on the reference run (case 1).

In order to calculate the integrated properties of inflow and outflow, such as mass flux or angular momentum flux, we need to carefully select the integration domain.¹¹ The derived fluxes depend strongly on the choice of the integration boundaries—thus, on how we distinguish between material belonging to accretion or ejection. In reality, there is no disk surface, but a smooth transition between accretion and ejection. The initial setup of a thin disk with aspect ratio = 0.1 dynamically evolves in time, and so does the disk surface.

Our control volume to measure the disk accretion and ejection rates is defined by an axisymmetric sector enclosed by two surfaces perpendicular to the equatorial plane located at r = r₁ and r = r₂, and two other surfaces being the disk mid-plane, and a surface s₁ which is close by and parallel to the initial disk surface (see Figure 9). We usually adopt r₁ = 1.0.

Zoom In Zoom Out Reset image size

The accretion rate is then calculated as

$$\begin{equation} \dot{M}_{\rm acc}(r) = - 2 \pi r \int _0^{a H} \rho v_r dz, \end{equation} \tag{ 14 }$$

where the parameter a controls the scale height of the integration and H = H(r, t = 0) = r is the initial thermal scale height of the disk. Considering the evolution of the large-scale disk velocity (see Figure 6), we have chosen a = 1.6 for our reference simulation. Similarly, we calculate the ejected mass flux as

$$\begin{equation} \dot{M}_{\rm ejec}(r)\vert _{\rm S} = \int _{r_1}^{r} \rho \bm {v}_{\rm p}\cdot d\bm {A}_{\rm S}, \end{equation} \tag{ 15 }$$

where the integration is performed along the inclined surface area S from r₁ to r₂ considering the area element $d\bm {A}_{\rm s}$.

The time evolution of the accretion and ejection rate for the reference simulation is shown in Figure 10. We have calculated the accretion and ejection rates for different sizes of the integration domain (thus with different outer radius r₂). We observe that for small radii, the accretion rate saturates at an approximately constant value. For the outer parts of the disk, however, a longer dynamical time is required for the disk to evolve into a new dynamical steady state. This is visible in the time evolution of the accretion rates—for increasing radii r₂ = 3, 5, and 10, the time when the plateau phase is reached, is increasing to t = 300, 800, and 1500. Beyond the plateau phase, the mass fluxes slightly decrease, most probably due to the overall decrease of the disk mass itself, and the subsequent change in the internal disk dynamics. The accretion rates at large radii are larger than those at smaller radii—the mass difference is ejected as outflow. For the control volume with larger r₂ the ejection rate increases, which is simply a consequence of the increased integration area.

Zoom In Zoom Out Reset image size

The last panel in Figure 10 shows the ejection-to-accretion rate $\dot{M}_{\rm ejec}/\dot{M}_{\rm acc}$. Again, depending on the size of the control volume this fraction changes. Although the outflow quickly evolves from the disk surface, and extends soon to large radii, several orbital periods are required to establish full accretion. At earlier times and for large radii, the accretion process is not fully established, resulting in somewhat arbitrarily high ejection-to-accretion ratios even above unity. For the inner part of the disk within radii <10, about 60% of the accreting material is diverted into the outflow for our reference simulation. This is a large fraction and similar to simulations in the literature, but definitely more than derived in steady-state models (Pelletier & Pudritz 1992; Ferreira & Pelletier 1995). For the other cases we have investigated, also smaller rates were obtained (see below).

Applying a radially self-similar approximation of the MHD equations, Ferreira & Pelletier (1995) have introduced a so-called ejection index ξ,

$$\begin{equation} \frac{\dot{M}_{\rm ejec}}{\dot{M}_{\rm acc}} = 1 - \left(\frac{r_{\rm i}}{r_{\rm e}} \right)^\xi \end{equation} \tag{ 16 }$$

essentially resulting from local mass conservation (where r_i = r₁ and r_e = r₂ in our notation). With self-similar solutions, Ferreira (1997) constrained the ejection index to 0.04 < ξ < 0.08. In comparison, for our reference run we find both larger numerical values and also a wider range for the ejection index, 0.1 < ξ < 0.5. Table 1 provides an overview of the ejection indexes we measured. We can apply our reference run to different astrophysical sources (see Appendix A), thus for a stellar jet with the central mass ≈1 M_☉ and ρ_i ≈ 10⁻¹⁰, the accretion rate is ≈10⁻⁶ M_☉ yr⁻¹. This value for AGNs with the central mass ≈10⁸M_☉ and ρ_i ≈ 10⁻¹² is about ≈7 M_☉ yr⁻¹.

In order to investigate the launching, the acceleration and the collimation processes of the outflow, we now consider the forces acting in the disk-outflow system.

To identify the forces for acceleration and collimation explicitly, it is preferable to project them along or perpendicular to a certain flux surface, respectively. By comparing these projected force components, we may disentangle the dominant driving and collimation mechanism. Figure 11 shows the force components along a flux surface (field line) rooted at (5, 0), and also the critical MHD surfaces—the slow-magnetosonic, the Alfvén, and the fast-magnetosonic surface. The forces involved in driving the outflow are the centrifugal force F_C, the Lorentz force F_L, the gas pressure gradient F_P, and gravity F_G. Among them, the pressure gradient and the centrifugal force are de-collimating while gravity and Lorentz force have collimating components. The pressure gradient, centrifugal force, and Lorentz force also contribute to accelerate the outflow. In agreement with the literature, we see that upstream of the slow-magnetosonic point, the acceleration is mainly by the gas pressure gradient. For the superslow outflow, the Lorentz force and the centrifugal force dominate. For the superfast flow, the Lorentz force plays a major role. In the case of the collimating forces the situation is different. Before the Alfvén point, the de-collimating centrifugal force dominates. Beyond the Alfvén point the Lorentz force has a comparable contribution, and finally after the fast-magnetosonic point, the Lorentz force dominates and collimates the outflow. The Lorentz force behavior is a sign of the electric current distribution in the disk/jet system. It compresses the disk while collimating the corona.

Zoom In Zoom Out Reset image size

The forces involved in launching the outflow are mainly the vertical gas pressure gradient, which counteracts the tension forces of the poloidal field, and gravity (see Figure 11). The thermal pressure gradient is always positive, and thus supports launching. It increases from the disk mid-plane to the surface and then decreases in the corona. The tension term is mostly negative, thus compressing the disk, and does not support launching. The magnetic pressure (toroidal component and total) are both positive above the disk, and, thus, may support the launching process. The toroidal component of the Lorentz force F_{ϕ, L} provides the magnetic torque braking the disk, and loading the outflow with both angular momentum and energy (not shown). In order to brake the disk, the torque must be negative in the disk and change sign at the disk surface (Ferreira 1997). This is confirmed by our simulation.

Figure 12 shows for comparison the three main force components (top) and the net force (bottom). We see that in the disk the vertical force components almost balance, while above the disk a considerable net force remains which launches and accelerates the outflow. In the disk, the (positive) gas pressure is almost balanced by the (negative) Lorentz force. However, the gas pressure slightly prevails and counteracts gravity. Essentially the figure shows that launching is a process which happens at the disk surface. The outflow material is not lifted from the mid-plane into the disk wind. It is the disk material accreting along the disk surface that is loaded into the outflow. The net (specific) force which is responsible for launching and initial acceleration is about 10% of the value of the single force components.

Zoom In Zoom Out Reset image size

The radius of the asymptotic jet and its opening angle can be measured by the observations rather easily. Therefore, it is interesting to provide a comparison from the simulations. Simple energy conservation arguments tell us that jets with the observed kinetic energy must originate from a disk area very close to the central object. This is a fact which seems to hold for all astrophysical jet sources. We will later consider the question of how these asymptotic properties are related to the conditions of the jet-launching process. Here, we first provide a clear numerical definition of these properties and apply them to our reference simulation. In order to make a quantitative comparison, we need to define a "radius of the outflow." We suggest a definition of an axial mass-flux-weighted jet radius, measured at a certain distance z = z_m from the source (see also Porth & Fendt 2010),

$$\begin{equation} r_{\rm jet}\vert _{z_{\rm m}} = \frac{\big[\int _0^{r_{\rm m}} r \rho v_{\rm z} dr\big]_{z_{\rm m}}}{\big[\int _0^{r_{\rm m}} \rho v_{\rm z} dr\big]_{z_{\rm m}}}. \end{equation} \tag{ 17 }$$

This radius measures the bulk of the mass flux contained in the jet at a specific distance z_m from the mid-plane. In order to derive the jet-launching area, we now adopt that flux surface which follows the bulk of the mass flux, thus passing through the point (r_jet(z_m), z_m), and trace the same flux surface back to the disk surface where it is rooted. This footpoint radius defines the launching area of the outflow. For reference run case 1, we find a mass-flux-weighted asymptotic jet radius for the high-velocity component of r_jet ≃ 44 at z_m = 180, corresponding to ≃ 4 AU for a protostellar scaling applying r_i = 0.1 AU. For the jet-launching radius, we measure r_l ∼ 0.4 AU in our reference run. Similarly, the launching area for the extended low-velocity component of the DG Tauri wind was measured by Anderson et al. (2003) as extending from ∼0.3 to 4 AU from the star. They invented a method relying on the MHD energy and angular momentum conservation along the jet. By comparing the observed kinetic energy and jet rotation they inferred the necessary disk rotation and launching area in the wind-launching region.

The asymptotic jet opening angle is the other jet characteristic. One way to measure the jet opening angle is by the inclination angle of the flux surface defined by the asymptotic jet radius as discussed above. Another measure of the collimation was suggested by Fendt (2006) who assigned an average collimation degree ζ of the outflow by comparing the vertical and radial mass fluxes in the outflow (applying a proper normalization per surface area of a cylinder of height z = z_m and radius r = r_m),

$$\begin{equation} \zeta _{\rm jet} = 2 \frac{z_{\rm m}}{r_{\rm m}} \frac{ \big[\int _0^{r_{\rm m}} r \rho v_{\rm z} dr \big]_{z_{\rm m}} }{ \big[\int _0^{z_{\rm m}} r_{\rm m} \rho v_{\rm r} dz \big]_{r_{\rm m}} }, \end{equation} \tag{ 18 }$$

where only positive poloidal velocities are considered. Outflow collimation would simply imply that ζ_jet > 1. For our reference run, we find ζ_jet = 2.37 implying that about twice as much mass is propagating along the jet axis than away from the jet axis. For comparison, the flux surface which encloses the bulk of the jet mass flux has an (half) opening angle of about 10° at (r, z) = (44, 180), but collimates slightly more further downstream. We note that this definition is related also to the concentration of mass flux across the jet, and, thus, provides somewhat different information than the opening angle of the field lines. For example, a cylindrical jet with zero degree opening angle may have a narrow or a broad radial density or mass flux profile. With our definition the narrow mass flux would be interpreted as more collimated.

Before we further investigate the physical effects, we discuss the results of our resolution study. Numerical diffusivity will add to the physical diffusivity, as it is a natural consequence of the finite-difference scheme applied in the PLUTO code. This effect could be particularly important in the jet-launching regime close to the disk surface where strong gradients in density, pressure, or magnetic diffusivity exist, and which may not be resolved. Murphy et al. (2010) have claimed that jet launching could be in fact present just due to numerical diffusivity.

In order to estimate the impact of numerical diffusivity we have repeated our reference simulation case 1 with a three times higher resolution (case 7), but with the same physical diffusivity and anisotropy parameter as for the reference simulation. We find that the leading disk and outflow properties are similar to the reference case. The accretion rate is slightly decreased from $\dot{M}_{\rm acc} = 0.015$ (case 1) to $\dot{M}_{\rm acc} = 0.013$ (case 7), while the ejection rate is slightly increased from $\dot{M}_{\rm ejec} = 0.008$ to $\dot{M}_{\rm ejec} = 0.009$ (see Table 1). Due to the smaller computational domain for the high-resolution simulation, we cannot compare the asymptotic jet radii (at z = 180); however, we can compare the radius of the bulk mass flux similar to Equation (17) for lower altitudes (at z = 60). The simulation with higher resolution seems to result in a slightly more collimated jet, with a jet radius of r_jet = 19.6 compared to r_jet = 25 for the reference simulation. This results in a similarly smaller jet-launching radius. The maximum jet velocities are just the same v_jet = 0.5.

Figure 14 shows the time evolution of the mass fluxes. We see that the accretion rates and ejection rates for case 1 and case 7 saturate at the same level. It appears that the high-resolution simulation needs more time to establish an outflow, while the accretion evolves similarly in both runs. One may see indication for a slightly larger accretion rate for the low-resolution run, which would fit into the picture that the disk material can more easily diffuse across the field lines due to the numerical diffusivity.

The same picture holds for the angular momentum fluxes (see Figure 15), which are very similar for the kinematic part and slightly offset for the magnetic part.

There is a common agreement in the literature that jet formation requires a certain amount of magnetic flux to be present in the jet-launching regime. On the other hand, the maximum magnetic flux which can be supported by the accretion disk is limited by the disk equipartition field strength (we will neglect the question of the origin of the magnetic field in this paper).

We have also studied the impact of the magnetic field strength on jet launching, governed by the plasma-β parameter. We note that we start all simulations with the same magnetic field profile; however, due to diffusion and advection in the disk, the field distribution may change substantially.

The magnetic field strength determines the amount of magnetic energy which is available for jet acceleration and is directly interrelated with the length of the lever arm of the magnetic torque on the disk (e.g., Pelletier & Pudritz 1992).

Due to the larger magnetic torque in case of a strong field, the disk angular momentum could be removed more efficiently. In order to investigate these effects quantitatively, we will compare simulations with different plasma-β, such as case 1 with initial β = 10, case 8 with β = 50, case 9 with β = 250, and the weak field case 6 with β = 5000. Note that the plasma-β is, however, a space- and time-dependent function of the simulation, β = β(r, z, t) and not only a single parameter we prescribe initially at the inner disk radius. In the high plasma-β regime, the flow of matter controls the dynamical structure of the system. For all cases of a weak magnetic field, the disk–jet system evolution seems highly perturbed—visible in the global structure of the outflow (Figure 13, cases 6, 8, and 9).

Note that the regimes with high plasma-β are known to be dominated by the internal, turbulent torque, which is not taken into account in our simulations by an α-viscosity.

In particular, for cases 9 and 10, the mass flux evolution is exceptionally perturbed (so we decided not to include them in all plots and just show the magnetic angular momentum flux). Figure 15 proves that with increasing plasma-β, less angular momentum is removed from the disk, and, subsequently, the accretion rate is reduced. This is particularly visible when we compare the simulation with β = 50 with the reference run with β = 10. In the case of the weakest magnetic field, thus for the simulations with β = 250 and β = 5000, the magnetic angular momentum removal is close to negligible. For the weak field cases, the accretion rate decreases, and in cases 9, 10, and 6 no efficient accretion is observed. In summary, we confirm the hypothesis that efficient magnetocentrifugal jet driving requires a strong magnetic flux (i.e., a low plasma-β), together with a large enough magnetic torque in order to produce a powerful jet.

Simulation case 6 applies a very weak magnetic field with β = 5000. Figure 16 compares the two simulations: case 6 with the standard setup case 7 in high resolution. No smoothly structured outflow is obtained for case 6, but a highly turbulent outflow with non-negligible outflow speed and mass flux. Due to the weak poloidal field, the outflow in case 6 is super-Alfvénic right from the launching point, thus magnetocentrifugal acceleration cannot play a role. Interestingly, the size of the turbulent features increases with distance from the origin. Also the poloidal magnetic field is highly tangled (Figure 13). The mass flux we measure for case 6 is about $\dot{M}_{\rm ejec} \simeq 0.001$ with outflow velocities of v_jet ≃ 0.4. While the velocities are comparable to case 7 (or case 1 with lower resolution), the ejected mass flux is substantially lower (factor two). Accretion in case 6 is weak (the smallest of all simulations), consistent with the low angular momentum losses (the smallest of all simulations). The accretion rate is even smaller than the ejection rate (factor five), and we may call such a disk an ejection disk instead of an accretion disk.

Zoom In Zoom Out Reset image size

The question is: what is launching and accelerating such a turbulent, high plasma-β outflow? Our interpretation is that the initial acceleration to superescape speed is by the toroidal magnetic pressure gradient (induced by the differential rotation between the rotating disk and the static corona). When the bow shock has left the domain, this acceleration process decays, and the remaining outflow acceleration is due to weak Lorentz force (as we are in the super-Alfvénic regime). The vertical mass flux for case 6 (very weak field) is one order of magnitude below the mass flux in case 7 (with strong field β = 10), measured at the same altitude (see Table 1). In summary, launching conditions as in case 6 with β = 5000 do not produce a strong jet flow, but a light disk wind of superescape speed.

We point out that the physical regime of acceleration does not only depend on the field strength, but also on the mass flux. Magnetocentrifugal effects are more evident for outflows with low-mass flux, while in heavy jets the magnetic pressure gradient may play a substantial role (see, e.g., Anderson et al. 2005).

The relative field strength (measured by the plasma-β) is a function of space and evolves in time. The initial field and gas pressure distribution, denoted by the β_i, is changed by the dynamical evolution including advection and diffusion of the disk magnetic field. Figure 17 shows the evolution of the plasma-β for case 1. It is clearly visible that the plasma-β along the disk surface changes substantially. The inner disk (for r < 15) reaches a steady state with a β somewhat higher than the initial value. Later, around t = 3000 when the disk mass decreases, the gas pressure decreases as well, resulting in a lower plasma-β in the disk corona. The disk corona therefore remains sufficiently magnetized and can support a fast jet. In case 1 the magnetic field structure reaches a steady state balanced by field diffusion advection.

Zoom In Zoom Out Reset image size

In simulation case 7 with the same parameter setup but a higher resolution, a steady state is reached as well. However, even for late time steps the magnetic flux continues to diffuse outward. This is somewhat surprising since the numerical diffusion is smaller and cannot be the reason (Figure 18, top). For case 8 with the higher initial plasma-β, advection seems to dominate outward diffusion of flux. The poloidal field distribution reaches a steady state with a magnetic flux concentration close to the inner edge of the disk (Figure 18, bottom). Thus, this setup seems favorable for jet launching as well, although the initial plasma-β was high.

Zoom In Zoom Out Reset image size

As an extreme example, for simulation case 9 with the initial β = 250 the plasma-β is considerably changed during the disk evolution. Case 9 shows very weak accretion (see the discussion above), so outward diffusion of magnetic flux dominates advection. Outward diffusion of flux implies a decrease of magnetic field strength. As a result the disk magnetization decreases and the plasma-β reaches number values above 10⁵. Therefore, this disk is not able to launch strong jets. So far no steady state is achieved.

The interrelation between the magnetic field strength and outflow launching has been discussed by other authors (Tzeferacos et al. 2009; Murphy et al. 2010), indicating that regions with weak field are not able to generate an outflow, and that both the collimation degree and the ejection rate increase with stronger field. Murphy et al. (2010) have concluded that jet launching for cases of weak magnetization may be artificially supported by numerical diffusivity within the disk surface layer, which should heat the gas, producing additional gas pressure. They suggest that ejection is possible as the magnetization reaches unity at the disk surface due to the steep density decrease. No ejection is reported when the mid-plane magnetization becomes too small. Nevertheless, the asymptotic jet velocity remained too low to explain the observed jet's speed. From our simulations, we find that even for weak magnetization in the disk the disk corona is sufficiently magnetized for jet launching and, depending on the efficiency of mass loading, fast jets could be driven (for comparison, see Figure 17 for case 1). The mass loading depends on resistivity, and we will discuss this aspect in the next section.

We finally consider the possibility to observe the MRI in our simulations. This is interesting since we are dealing with magnetic fields of various strength in combination with a differentially rotating system. To generate the MRI in numerical simulations, two conditions are essential—a large enough (but not too large) plasma-β in the disk as to initiate instability, and also a grid of high enough resolution in order to avoid damping of small wavelengths. The wavelength of the most unstable MRI mode is given by λ_MRI = 2πv_{A, z}/Ω_K for a Keplerian rotation Ω_K (see, e.g., Romanova et al. 2011). When we calculate this wavelength for a simulation with β = 10, we obtain λ_MRI ≃ 0.057, implying that we cannot resolve the MRI with our setup as the grid size is about half of the MRI wavelength.

We point out that the magnetic diffusivity applied in our simulations also contributes to suppressing the generation of MRI turbulence.

The strength and anisotropy of the magnetic diffusivity—controlled by the coefficients η_{ϕ, i}, η_{p, i}, and χ (see Table 1)—are essential parameters governing the coupling between the magnetic field and the plasma. Efficient magnetocentrifugal driving requires strong coupling between the magnetic field and the rotating disk material, thus a rather low diffusivity. A similar argument holds for the launching by magnetic pressure gradient. On the other hand, mass loading from accretion to ejection requires a certain degree of diffusivity in order to transport mass across the field lines.

So far, no general model for the magnetic diffusivity distribution in disk–jet structures is available. In many simulations considering magnetocentrifugal-driven outflows, the outflow is governed by ideal MHD, while the diffusivity is confined to the disk. In order to investigate the impact of diffusivity on jet launching, we have therefore performed simulations with varying strength of the diffusivity components η_{ϕ, i} and η_{p, i}.

We first compare simulations with the same anisotropy parameter χ = 3 as reference run case 1. In general, the strength of magnetic diffusivity governs the disk accretion rate. Our reference run case 1 reaches a quasi-steady state, establishing a balance between inward advection and outward diffusion of magnetic flux. This happens after t = 1000 and is disturbed again in the late stages of the dynamical evolution, most probably as a result of the change in the overall disk dynamics due to the decreasing disk mass. A change in diffusivity will also affect this balance. We now compare case 4 with higher and case 2 with lower diffusivity. In the less diffusive case, advection dominates and we obtain a higher accretion rate. We find that not only the accretion rate, but also the ejection rates are higher for disks with lower diffusivity. This trend is shown in Figure 19, where we display the accretion rate and also the ejection rate as a function of the poloidal diffusivity.

Zoom In Zoom Out Reset image size

Two physical processes affect the disk dynamics—advection and diffusion. In a less diffusive disk, the plasma is stronger coupled to the field and advection of magnetic flux dominates. Consequently, the magnetic flux surfaces are located further in the inner part of the disk, "rotate" faster, and a stronger B_ϕ is induced. The stronger toroidal field may lead to a stronger acceleration of the outflow, by stronger Lorentz forces (either magnetic tension of pressure forces). In addition, by comparing the vertical profiles of net launching forces at t = 1000 for a less diffusive disk (case 2), and reference run (case 1), we find a larger net force in case 2 by which implies that more disk material can be lifted into the outflow.

This is exactly what we observe in the mass fluxes—the ejected mass flux in case 2 is two times the ejected mass flux in case 1. Note that a stronger B_ϕ will impose stronger pinching forces ∼(B_ϕ · ∇)B_ϕ on the disk and the ejected material, thus opposing the loading process (as the toroidal tension has a component along the field line decelerating the outflow meterial). However, the magnetic pressure gradient of the toroidal field ∼∇B²_ϕ has an outward directed component along the field line, thus supporting ejection. By comparing the vertical profiles of both terms for the two cases with β = 10 and β = 50, we find that (1) the toroidal magnetic field pressure gradient component along the field line is always dominating the pinching force of the toroidal field along the field line, and that (2) for higher plasma-β the (∇B²_ϕ)_z is larger than for the lower plasma-β case. Thus, in our simulations the toroidal field is supporting ejection.

In summary, although a certain magnitude of diffusivity is required for mass loading, the lower diffusivity allows for enhanced mass loading of accreting material into the outflow. Consequently, less diffusive disks tend to have higher ejection rates. This result is in agreement with the previous literature (Zanni et al. 2007).

Our simulations confirm that the strength of magnetic diffusivity does affect the asymptotic speed of the outflow. Figure 20 shows the typical jet velocity versus the poloidal diffusivity. Comparing different simulation runs (see Table 1), we find that the outflow from the less diffusive disk (case 2) remains slower, while the outflows formed from more diffusive disks (cases 3 and 4) are accelerated to higher speed. For case 2 with three times lower magnetic diffusivity, the asymptotic speed is reduced by ≈30%, while for case 4 with three times higher magnetic diffusivity the increase is about ≈30% (see Table 1). This is reasonable from a physical point of view, since in a weakly diffusive disk more material is diverted into the outflow. Thus, with the same amount of magnetic flux available, correspondingly less magnetic energy could be transferred per outflow mass unit. The more massive outflows are only accelerated to lower speed if launched with similar magnetization. Alternatively, a higher magnetic diffusivity, resulting in a weaker mass load, leads to a faster outflow.

Zoom In Zoom Out Reset image size

When comparing our results with the previous literature, one should keep in mind that we have applied a mass flux weighted velocity which we consider as the typical velocity of the bulk mass flux. The velocity value is generally lower than the maximum speed we measure in the outflow and which is mostly given in the literature. While the typical speed ranges from about 0.5 to 0.9 (in code units), the maximum speed in the outflow ranges from 1.2 to 1.8 inner disk Keplerian velocities and is comparable to other results in the literature. Furthermore, we have detected a slight time evolution in the velocity (see also the mass flux evolution shown in Figure 14). In our reference run case 1, the maximum asymptotic speed varies from v_{p, max} = 1.8 at t = 500 to v_{p, max} = 1.2 for t > 1000. For comparison, Zanni et al. (2007) and Tzeferacos et al. (2009) give a maximum speed of v_{p, max} ≃ 1.5–4.5 at t = 400.

Next, we investigate a possible impact of anisotropy in the diffusivity tensor, parameterized by χ. We compare case 1 and case 5 with the same polodial diffusivity η_p = 0.03, but with case 5 having a lower diffusivity in the toroidal direction. By comparing their leading properties (see Table 1), we find them to have quite similar properties. In particular, for both cases (with the same poloidal diffusivity) the ejection rates are also similar. However, case 5 with less toroidal diffusivity has a higher accretion rate.

Anisotropy of magnetic diffusivity also impacts the rotation of the outflow. It is believed that jets are rotating and the rotation is basically driven by the underlying disk rotation and the coupling between matter and field (Bacciotti et al. 2002; Anderson et al. 2003; Coffey et al. 2004; Fendt 2011). For a high toroidal diffusivity, the coupling is weak, and thus the acceleration in toroidal direction is also weak. With our simulations we can confirm this concept. Figure 21 shows the rotational velocity for the lower part of the jet. The outflow resulting from simulation case 1 with a nine times higher toroidal diffusivity shows a 50% lower rotation rate than for case 5. Similarly, the outflows in simulation case 4 with a three times higher toroidal diffusivity shows a 50% lower rotation rate than case 3.¹³

Zoom In Zoom Out Reset image size

We observe a close correlation between the accretion rate $\dot{M}_{\rm acc}$ and angular momentum flux $(\dot{J}_{\rm kin}+\dot{J}_{\rm mag})$ from the disk. In a system with higher angular momentum removal, a higher accretion rate is observed (Figure 15), with the accretion rate in case 5 being larger than for case 1, and for case 1 larger than for case 3. This confirms the common belief that in order to obtain higher accretion rates, a more efficient angular momentum removal is required.

To discuss the interrelation between collimation and diffusivity, we again apply a collimation degree by comparing the directed mass fluxes (see Equation (18)). As a general result we find that outflows launched from disks of higher poloidal diffusivity η_p are more collimated. This interrelation is displayed in Figure 22.

Zoom In Zoom Out Reset image size

We see that case 2 with the magnetic diffusivity of η_{p, i} = 0.01—maybe below a critical value—evolves differently from the others and has both a very low degree of collimation and an exceptionally large mass-flux-weighted jet radius. The inner part of this outflow appears rather hollow (Figure 13).

The comparison of the (mass flux weighted) jet radii for these runs provide a similar picture. We find an interrelation such that for increasing magnetic diffusivity the jet radius decreases,¹⁴ such that r_jet(case 3) ⩽ r_jet(case 4) < r_jet(case 5) < r_jet(case 1) < r_jet(case 2).

We have argued above that the disk material which is launched to form an outflow is likely to be turbulent, so we may expect a magnetically diffusive disk wind above the disk surface. Thus, since mass loading requires poloidal magnetic diffusivity, the jet-launching area can be extended into higher altitudes above the disk surface. Furthermore, the jet-launching area is extended into a domain where the plasma-β is comparatively low. Taking into account both effects, one may therefore expect (1) to launch jets from disks which would themselves be insufficiently magnetized for jet driving, and also (2) to launch outflows which have a rather low-mass flux and subsequently higher terminal velocity. In order to investigate these effects, we have therefore prescribed different height scales H_η for the magnetic diffusivity in our simulations (see Equation (13)). The physical properties measured for different diffusive scale heights are shown in Table 2. In general, we find different stages in the temporal evolution of the mass flux (accretion and ejection) in the disk–jet system (Figure 23). For small diffusive scale heights H_η, a turbulent pattern and perturbations appear in the outflow resulting in a more filamentary outflow. It seems that these perturbations have larger amplitudes in case of a smaller diffusive scale height. The amplitudes decrease when the diffusive scale height decreases. For _η = 0.4 they are damped completely. We believe that these perturbations are generated by a physical process within the launching region and that they are simply smoothed out by the magnetic diffusivity. A further study is needed to clarify this issue.

Zoom In Zoom Out Reset image size

For increasing scale height the accretion rate and ejection rate decrease, confirming the correlation derived above between the increasing magnitude of diffusivity and decreasing mass fluxes. In other words, an increasing diffusivity scale height simply corresponds to a higher diffusivity with the consequences as discussed above. Interestingly, all our simulations for different diffusivity scale heights converge to the same kinetic angular momentum flux in the outflow. On the other hand, we find a trend of decreasing magnetic angular momentum flux with decreasing diffusivity scale height.