MODELING MHD ACCRETION–EJECTION—FROM THE LAUNCHING AREA TO PROPAGATION SCALES

No physical scales are introduced in the equations above. The results of our simulations will be presented in non-dimensional units. We normalize all variables, namely $P, \rho, {\boldsymbol {V}},\ {\rm and}\ {\boldsymbol {B}}$, to their values at the inner disk radius R₀. Lengths are given in units of R₀, corresponding to inner disk radius. Velocities are given in units of V_{K, 0}, corresponding to the Keplerian speed at R₀. Thus, 2πT corresponds to one revolution at the inner disk radius. Densities are given in units of ρ₀, corresponding to R₀. Pressure is measured in $P_{0} = \epsilon ^2 \rho _{0} V_{0}^2$.

We thus may apply our scale-free simulations to a variety of jet sources. In the following, we show the physical scaling concerning three different object classes—brown dwarfs (BDs), YSOs, and AGNs. In order to properly scale the simulations, we vary the following masses for the central object, M = 0.05 M_☉ (BD), M = 1 M_☉ (YSO), M = 10⁸ M_☉ (AGN), and thus define a scale for the inner disk radius of

$$\begin{eqnarray} R_{\rm 0}& = & {\rm 0.1}\, {\rm AU}\quad {\rm (YSO)}\nonumber \\ & = & {\rm 0.01}\, {\rm AU}\quad {\rm (BD)}\nonumber \\ & = & 20\,{\rm AU} \left(\frac{R_{\rm 0}}{10 \,{R}_{\rm S}} \right) \left(\frac{M}{10^8\, {M}_{\odot }} \right) \quad {\rm (AGN)}, \end{eqnarray} \tag{ 6 }$$

where R_S = 2GM/c² is the Schwarzschild radius of the central black hole. For consistency with our non-relativistic approach, we require R₀ > 10R_S, implying that we cannot treat highly relativistic outflows. Table 1 summarizes the typical scales for leading physical variables. For more detailed scaling laws, we refer to Appendix A.

Table 1. Typical Parameter Scales for Different Sources

	YSOs	BDs	AGNs	Unit
R0	0.1	0.01	20	AU
M0	1	0.05	108	M☉
ρ0	10−10	10−13	10−12	g cm−3
V0	94	66	6.7 × 104	km s−1
B0	15	0.5	1000	G
T0	1.7	0.25	0.5	days
$\dot{M}_0$	3 × 10−5	2 × 10−10	10	M☉ yr−1
$\dot{J}_0$	3.0 × 1036	1.5 × 1030	3 × 1051	dyne cm
$\dot{E}_0$	1.9 × 1035	6.7 × 1029	2.6 × 1046	erg s−1

Note. Simulation results will be given in code units and can be scaled for astrophysical application.

Download table as:  ASCIITypeset image

We apply a numerical grid with equidistant spacing in the θ direction, but stretched cell sizes in the radial direction, considering ΔR = RΔθ. Our computational domain of a size R = [1, 1500R₀], θ = [0, π/2] is discretized with (N_R × N_θ) grid cells. We use a general resolution of N_θ = 128. In order to cover a factor 1500 in radius, we apply N_R = 600. This gives a resolution of 16 cells per disk height (2) in the general case. However, we have also performed a resolution study applying a resolution two times as high (or low), thus using 256 × 1200 (or 64 × 300) cells for the whole domain, or 35 (9) cells per disk height. We satisfy the Courant–Friedrichs–Lewy condition by using a CFL number of 0.4.

For the initial conditions, we follow a meanwhile standard setup, which has been applied in a number of previous publications (Zanni et al. 2007; Sheikhnezami et al. 2012; Fendt & Sheikhnezami 2013). The initial structure of the accretion disk is calculated as the solution of the steady state force equilibrium,

$$\begin{equation} \nabla P +\rho \nabla \Phi _{\rm g}- {\boldsymbol {J}} \times {\boldsymbol {B}} - \frac{1}{R} \rho V_{\rm \phi }^2 ({ {\boldsymbol {e}} }_R \sin \theta + { {\boldsymbol {e}} }_\theta \cos \theta) = 0. \end{equation} \tag{ 7 }$$

We solve this equation assuming radial self-similarity, i.e., assuming that all physical quantities X scale as a product of a power law in R with some function F(θ),

$$\begin{equation} X \equiv X_0 R^{\beta _X} F(\theta). \end{equation} \tag{ 8 }$$

Self-similarity requires in particular that the sound speed and the Alfvén speed scales as the Keplerian velocity, V_K∝r^−1/2, along the disk midplane. As a consequence, the power-law coefficients β_X are determined as follows: $\beta _{V_{\rm \phi }} = -1/2, \beta _{\rm P}= -5/2, \beta _\rho = -3/2$, and $\beta _{B_{\rm R}} = \beta _{B_{ \theta }} = \beta _{B_{ \phi }} = -5/4$.

An essential non-dimensional parameter governing the initial disk structure is the ratio between the isothermal sound speed $C_{\rm s}^{\rm T}= \sqrt{P/\rho }$ and the Keplerian velocity $V_{\rm K}= \sqrt{GM/r}$, evaluated at the disk midplane, $\epsilon = [ C_{\rm s}^{\rm T}/ V_{\rm K}]_{\,\theta =\pi /2}.$ This quantity determines the disk thermal scale height H_T = r. In our simulations, we generally initially assume a thin disk with = 0.1. Note that for the remainder of the paper, when discussing the dynamical properties of disk and outflow, we consider the adiabatic sound speed $C_{\rm s}=\sqrt{\gamma P/\rho }$. The geometrical disk height, namely, the region where the density and rotation significantly decrease, is about 2 (see Figure 1). We therefore define the geometrical disk height as H ≡ 2r.

Zoom In Zoom Out Reset image size

Following Zanni et al. (2007), our reference simulation is initialized only with a poloidal magnetic field, defined via the vector potential ${\boldsymbol {B}} = \nabla \times A {\boldsymbol {e}} _{\phi }$, with

$$\begin{equation} A = \frac{4}{3} B_{\rm p,0} r^{-1/4} \frac{m^{5/4}}{ (m^2 + ctg^2\theta)^{5/8} }. \end{equation} \tag{ 9 }$$

The parameter $B_{\rm p,0} = \epsilon \sqrt{2 \mu _{\rm 0} }$ determines the strength of the initial magnetic field, while the parameter m determines the degree of bending of the magnetic field lines. For m → ∞, the magnetic field is purely vertical. As we will see below, the long-term evolution of the disk–jet structure is insensitive to this parameter, since, due to advection and diffusion processes and the jet outflow, the magnetic field structure changes substantially over time. We therefore set m = 0.5 in general.

The strength of the magnetic field is governed by the magnetization parameter

$$\begin{equation} \mu = \frac{B_{\rm P}^2}{2P}\, {\,\theta =\frac{\pi }{2}}, \end{equation} \tag{ 10 }$$

the ratio between the poloidal magnetic field pressure and the thermal pressure, evaluated at the midplane, and is set to be constant with radius. As we will see below, the magnetic field distribution substantially changes over time, while the disk–jet dynamics is governed by the actual disk magnetization. Typically, the initial magnetization is μ₀ ≈ 0.01 in this simulations.

Outside the disk, the gas and pressure distribution is defined as hydrostatic "corona,"

$$\begin{equation} \rho _{\rm cor} = \rho _{\rm cor,0} R^{-1/(\gamma - 1) }, \quad P_{{\rm cor}} = \frac{\gamma -1}{\gamma } \rho _{\rm cor,0} R^{-\gamma /(\gamma - 1)}, \end{equation} \tag{ 11 }$$

where ρ_{cor, 0} = δρ_disk(R = 1, θ = π/2) with δ = 10⁻³.

Although it is common to define the initial accretion velocity, balancing the imposed diffusivity V_R = ηJ_ϕ/B_θ, we find that this is not necessary in our set of parameters. It can even disturb the initial evolution of the disk accretion. Accretion requires corresponding torques to be sustained and since there is no initial poloidal current defined, B_ϕ = 0 (which takes time to build up from a weak poloidal field), a non-zero initial velocity will only lead to extra oscillations. Figure 2 illustrates this issue, showing two identical simulations with zero and non-zero accretion velocity. We found that the origin of the oscillations for both cases is the inner boundary. The first wave of accretion is somehow bounced outward by the inner boundary and results in an oscillatory pattern. We believe that the bouncing at the inner disk boundary can be caused by small inconsistencies between the boundary conditions and the intrinsic disk physics, a problem that may be present in other studies. Although they both result in the same final profile (a steady state has been reached only for small disk radii), the simulations with zero velocity profiles show fewer oscillations.

Zoom In Zoom Out Reset image size

We apply standard symmetry conditions along the rotational axis and the equatorial plane. Along the radial boundaries of the domain, we distinguish two different areas. That is (1) a disk boundary for θ > (π/2) − 2,¹ and (2) a coronal boundary for θ < (π/2) − 2, and consider different conditions along them (see Figure 3).

Zoom In Zoom Out Reset image size

Along the inner radial boundary for all simulations, we impose a constant slope for the poloidal component of the magnetic field

$$\begin{equation} \varphi = 70^\circ \left(1+ \exp (- \frac{\theta - 45^\circ }{15^\circ } \right)^{-1}, \end{equation} \tag{ 12 }$$

where φ is the angle with respect to the unit vector ${ {\boldsymbol {e}} }_R$. The magnetic field direction is axial near the axis, θ = 0, while at the inner disk radius the inclination is 70° with respect to the disk surface. A smooth variation of the magnetic field direction is prescribed along the inner radial boundary. This is in concordance with Pelletier & Pudritz (1992), who showed that for a warm plasma the maximum angle with respect to the disk surface necessary to launch outflows is about 70°, and slightly larger than for a cold plasma (Blandford & Payne 1982).

The method of constraint transport requires the definition of only a tangential component, thus to prescribe B_θ along the innermost boundary, while the normal component B_R follows from solving $\nabla \cdot {\boldsymbol {B}} =0$. In order to implement the prescription of a constant magnetic field angle, we solve $\nabla \cdot {\boldsymbol {B}} =0$, taking into account the ratio of the cell-centered magnetic field components B_θ/B_R = −tan (φ). We start the integration from the axis (θ = 0), where B_θ = 0. Thus, by fixing the slope of the magnetic lines, we allow the magnetic field strength to vary.

Along the inner coronal boundary, we prescribe a weak inflow into the domain with V_P = 0.2. This is applied to stabilize the inner coronal region between the rotational axis and the disk jet, since the interaction between the current carrying, magnetized jet and zero-B_ϕ coronal region may lead to some extra acceleration of the coronal gas. As shown by Meliani et al. (2006), the pressure of such an inflow (e.g., a stellar wind) may influence the collimation of the jet, changing the shape of the innermost magnetic field lines.

In order to keep the influence of the dynamical pressure of the inflow similar during the whole evaluation (and also for different simulations), we set the density of this inflow with respect to the disk density at the inner disk radius. The density of the inflow corresponds to a hydrostatic corona ρ_infl = ρ_cor = ρ_disk|_midplane(t) · δ, where δ = 10⁻³. The inflow direction is aligned with the magnetic field direction. By choosing a denser inflow, we also increase the time step of our simulations by approximately three times, as the Alfvén speed in the coronal region lowers.

By varying the slope of the magnetic field φ along this inner corona in the range of 60–80 $\deg$, we found that it only slightly affects the slope of the innermost magnetic field lines. The global structure of the magnetic field is instead mainly governed by the diffusivity model. Since the inner boundary by design models the magnetic barrier of the star, we choose a rather steep slope in order to avoid the disk magnetic flux entering the coronal region.

Across the inner disk boundary (that is the accretion boundary), density and pressure are both extrapolated by power laws, applying ρR^−3/2 = const and PR^−5/2 = const, respectively. Both the toroidal magnetic field as well as the toroidal velocity components are set to vanish at the inner coronal boundary, B_ϕ = 0, V_ϕ = 0. For the inner disk boundary, we further apply the condition B_ϕ ∼ 1/r (J_θ = 0), and extrapolate the radial and the toroidal velocity by power laws, V_RR^−1/2 = const, and V_ϕR^−1/2 = const, respectively, while V_θ = 0.

For the inner disk boundary, only negative radial velocities are allowed, making the boundary to behave as a "sink," thus absorbing all material that is delivered by the accretion disk at the inner disk radius.

As the application of spherical coordinates provides an opportunity to use a much larger simulation domain compared to cylindrical coordinates, the outer boundary conditions have little influence on the evolution of the jet launched from the very inner disk. We therefore extrapolate ρ and P with the initial power laws and apply the standard PLUTO outflow conditions for V_R, V_θ, and V_ϕ at the outer boundary, thus, zero gradient conditions. We further require B_ϕ ∼ 1/r (J_θ = 0) for the toroidal magnetic field component, while a simple outflow condition is set for B_θ. Again, B_R is obtained from the $\nabla \cdot {\boldsymbol {B}} =0$.

For the radial velocity component, we distinguish between the coronal region, where we require positive velocities V_R ⩾ 0, and the disk region, where we enforce negative velocities V_R ⩽ 0.

As our application of a spherical geometry is new in this context, we summarize the boundary conditions in the Table 2.

Table 3. Comparison of Our Simulations with Simulations Performed by Other Authors

Reference	Cell/2		m	αm	χ	μ0	μact
Casse & Keppens (2004)	0.5	0.1		<1	1	≃1	⋅⋅⋅
Zanni et al. (2007)	2.5	0.1	0.35	0.1...1.0	1, 3	0.3	⋅⋅⋅
Tzeferacos et al. (2009)	2.5	0.1	0.4	0.1...1.0	3, 100	0.1–3.0	⋅⋅⋅
Sheikhnezami et al. (2012)	8.0	0.1	0.4	1.0	1/3, 3	0.002–0.1	⋅⋅⋅
This work, reference simulation	16	0.1	0.2–0.9	1.1–1.9	0.5	0.003–0.03	0.001–0.5
This work, resolution study	24.5	0.1	0.2–0.9	1.1–1.9	0.5	0.003–0.03	0.001–0.5

Note. The resolution is estimated for the inner disk (R = 1).

Download table as:  ASCIITypeset image

Accretion disks are considered to be highly turbulent, subject to the MRI in moderately magnetized disks (Balbus & Hawley 1991; Fromang 2013), and the Parker instability (Gressel 2010; Johansen & Levin 2008) for stronger magnetized disks. It is believed that when the magnetic field becomes sufficiently strong, the MRI modes become suppressed (Fromang 2013). On the other hand, a strong magnetic field may become buoyant, leading to the Parker instability. While the MRI is confined within the disk, the Parker instability operates closer to the surface of the disk where the toroidal magnetic field is stronger.

Extrapolating the results from a self-consistent, local treatment of turbulence to the mean field approach is not straightforward. In the local treatment, the extraction of angular momentum is due to both turbulence—operating on small scales—and torques by the mean magnetic field on large scales. Thus, the removal of angular momentum goes hand in hand with destroying the turbulent magnetic field or the effective magnetic diffusivity. In case of the mean field approach, there is no small-scale turbulence and, thus, no angular momentum removal by local turbulent motions. Here, the diffusivity plays only a role for leveling the magnetic field gradient, thus setting the overall structure of the magnetic field. Unfortunately, we lack complete knowledge of the disk turbulence, thus the connection between the mean magnetic field and the fluctuating part, or, in other words, the relation between the mean magnetic field and the effective torques, and the diffusivity and viscosity that turbulence provides. We believe that when moving from a local turbulence approach to the mean field approach, the relevance of the model should be approved by the relevance of the magnetic field distribution itself, and not by the diffusivity model. However, one should keep in mind that the magnetic field strength and structure of real accretion disks are also not known. Therefore, when considering any simulation results, the diffusivity model applied should always be taken into account.

A self-consistent study of the origin of the turbulence is beyond the scope of our paper. We therefore prescribe a certain model of the magnetic diffusivity. We apply an α-prescription (Shakura & Sunyaev 1973) for the magnetic diffusivity, implicitly assuming that the diffusivity has a turbulent origin. The diffusivity profile may extend up to one disk height above the disk surface (Figure 1).

Although we investigate different diffusivity models, all of them can be represented in a following form

$$\begin{equation} \eta _{\rm P}= {\alpha }_{\rm ssm}(\mu) C_{\rm s}\cdot H \cdot F_{\rm \eta } (z), \end{equation} \tag{ 13 }$$

where the vertical profile of the diffusivity is described by a function

$$\begin{equation*} F_{\rm \eta } (z) = \left\lbrace \begin{array}{@{}ll}1 & z\le H \\ \exp \left(- 2\left(\frac{z-H}{H} \right)^2 \right) & z > H, \\ \end{array} \right. \end{equation*}$$

confining the diffusivity to the disk region.

Although this parameterization of diffusivity is commonly used (except the profile function F_η(z)), there are no clear constraints upon the value α_ssm may take. As an example, King et al. (2007) discuss a magnitude of the turbulent α-parameter derived from observations and simulations, indicating observational values α_ssm ≃ 0.1..0.4. Numerical models with zero net magnetic field usually provide low numerical values α_ssm ≃ 0.01, reaching at most α_ssm ≃ 0.03 (Stone et al. 1996; Beckwith et al. 2011; Simon et al. 2012; Parkin & Bicknell 2013). On the other hand, numerical modeling of the MRI applying a non-zero net magnetic field (Bai & Stone 2013) indicates substantially higher values, α_ssm ≃ 0.08–1.0, with a corresponding magnetization μ = 10⁻⁴, 10⁻².

Obviously, different functions of α_ssm(μ) will lead to different evolution. We start from the well-known model for magnetic diffusivity applied by many authors previously (Casse & Keppens 2004; Zanni et al. 2007; Sheikhnezami et al. 2012),

$$\begin{equation} \eta _{\rm P}= {\alpha }_{\rm m}V_{\rm A}\cdot H \cdot F_{\rm \eta } (z) \end{equation} \tag{ 14 }$$

by applying ${\alpha }_{\rm ssm}= {\alpha }_{\rm m}\sqrt{2\mu }$, where $V_{\rm A}= B_{\rm P}/ \sqrt{\rho }$ is the Alfvén speed, and μ, C_s, and H are the magnetization, the adiabatic sound speed, and the local disk height, respectively, measured at the disk midplane. We evolve α_ssm and C_s in time, but for the sake of simplicity, we keep H and F_η(z) constant in time, thus equal to the initial distribution. The main reason is to avoid additional feedback, which favors the accretion instability (see below, or, e.g., Campbell 2009). Our test simulations evolving the disk height in time in fact indicate the rise of such instability earlier than in the case of a fixed-in-time disk diffusivity aspect ratio.

In general, the diffusivity tensor has diagonal non-zero components,

$$\begin{equation} \eta _{\rm \phi \phi } \equiv \eta _{\rm P}\quad \eta _{\rm RR}=\eta _{\rm \theta \theta } \equiv \eta _{\rm T}, \end{equation} \tag{ 15 }$$

where we denote η_P as the poloidal magnetic diffusivity, and η_T as the toroidal magnetic diffusivity, respectively.² The anisotropy parameter χ ≡ η_T/η_P quantifies the different strength of diffusivity in poloidal and toroidal directions.

In the literature, it is common to assume χ of the order of unity. Considering viscous disks, Casse & Ferreira (2000) showed that there is a theoretical limit for η_T, namely η_T > η_P. Highly resolved disk simulations indeed suggest χ ≃ 2...4 (Lesur & Longaretti 2009), implying that the toroidal field component is typically diffusing faster than the poloidal component.

The majority of simulations in the literature consider a magnetic field strength in equipartition with the gas pressure. However, also studying weakly magnetized disks, we find that there also exists an upper limit for the anisotropy parameter, above which the simulations show irregular behavior. On the other hand, it was pointed out by several authors (see, e.g., Zanni et al. 2007) that in case of a very low anisotropy parameter (thus, a weak toroidal diffusivity), the simulations might suffer from instabilities caused by strong pinching forces. Nonetheless, the existence of an upper limit for the anisotropy parameter was so far obscured by other processes.

Assuming a steady state and combining the poloidal component of the diffusion equation, M_R = α_ssmHJ_ϕ/B_P, with the relation for the Mach number $M_{R}= 2/\sqrt{\gamma } HJ_{ R}/B_{\rm P} \mu _{\rm act}$ (see below, Königl & Salmeron 2011), an interrelation between the toroidal and poloidal electric currents can be derived,

$$\begin{equation} \frac{J_{ \phi }}{J_{ R}} = \frac{\sqrt{2\mu }}{\sqrt{\gamma }{\alpha }_{\rm m}}. \end{equation} \tag{ 16 }$$

This relation has been proven to approximately hold for all of our simulations, thus indicating that a steady state has indeed been reached. Since the only free parameter in this relation is α_m, the choice of α_m governs the ratio of the electric current components.

As shown previously (Ferreira & Pelletier 1995; Ferreira 1997; Ferreira & Casse 2013), the toroidal component of the induction equation can be written as

$$\begin{equation} \eta _{\rm T}J_{ R}|_{\rm mid} = - R^2 \int _0^{2\epsilon } \mathbf {B_{\rm P}} \cdot \nabla \mathbf {\Omega } d\varphi - V_{\rm \theta }B_{ \phi } \end{equation} \tag{ 17 }$$

(here expressed for spherical coordinates), where J_R is computed at the disk midplane. This equation essentially states that the induction of the toroidal magnetic field component (from twisting the poloidal component) is balanced by the diffusion through the disk midplane and by the flux escaping through the disk surface. We emphasize that we do not neglect the V_θB_ϕ term, considering the assumption of a thin disk (Ferreira & Casse 2013), as we find it of key importance in our simulations, in particular in the regime of a moderately strong magnetic field, μ ⩾ 0.1.

Assuming that the induction of the magnetic field is primarily due to radial gradients, the radial component of the magnetic field can be approximated by a power law, B_R∝R^−5/4, and Equation (17) may be transformed into η_TJ_R|_mid ≃ B_RC_s − V_θB_ϕ, or

$$\begin{equation} ({\alpha }_{\rm ssm}\chi - M_{\theta }) J_{ R}\simeq \frac{B_{\rm R}}{H}, \end{equation} \tag{ 18 }$$

where

$$\begin{equation} M_{\theta }= - V_{\rm \theta }^{+}/C_{\rm s} \end{equation} \tag{ 19 }$$

is denoted as the ejection Mach number, where $V_{\rm \theta }^{+}$ is measured at the disk surface.

Using relation 16 between the poloidal and toroidal electric currents and defining the curvature part of the toroidal current $J_\phi ^{\rm curv} \equiv B_{\rm R}/H$, one may derive a constraint for the anisotropy parameter,

$$\begin{equation} \frac{J_\phi ^{\rm curv}}{J_{ \phi }} \simeq \left({\alpha }_{\rm m}\chi - \frac{M_{\theta }}{\sqrt{2\mu }}\right){\alpha }_{\rm m}\le 1. \end{equation} \tag{ 20 }$$

Any magnetic field geometry that is outwardly bent and has a decreasing field strength in the outward direction has to satisfy this relation, as J_ϕ consists of two positive terms, the gradient and the curvature. In cases where the vertical velocity term can be neglected (e.g., for a very weak magnetic field with μ ⩽ 0.02), the anisotropy parameter is $\chi < 1/{\alpha }_{\rm m}^2$, which, for our choice of α_m, is about 0.4. By probing the χ parameter space, we found that in order to obtain a stable accretion-outflow configuration for weakly magnetized disks, the χ should be in the range 0.3–0.7. We therefore decided to apply χ = 0.5 for all of our simulations. We note that in simulations applying an anisotropy parameter χ ⩾ 0.7, we faced the problem that the poloidal magnetic field lines were "moving" rapidly, such that the bending of the field lines along the disk midplane was actually inverted. This is a result of the combined effects of a strong outward diffusion and low torques at the midplane. On the other hand, in case of a rather low anisotropy parameter, the accretion is rapid and the jet does establish a steady behavior.

The reason why the commonly chosen anisotropy χ > 1 leads to a steady behavior is rather simple. As the disk magnetization grows during advection, the ejection Mach number grows as well (we find that M_θ∝6μ saturating at a level of 0.8; see below). Thus, in the case of a high disk magnetization—usually assumed in the literature—the above-mentioned upper limit for the anisotropy parameter is satisfied. In the case of weak magnetic field simulations, performed, for example, by Murphy et al. (2010), this limit is most likely satisfied by additional viscous torques.

In the Introduction we have already discussed the literature of accretion–ejection simulations. Here we want to explicitly emphasize specific details in which our simulations differ from previous works.

• 1.  
  A spherical grid has been applied, offering the opportunity of a much larger domain size as well as much higher resolution in the inner part of the disk. A new set of boundary conditions is used which is adapted to the spherical grid.
• 2.  
  We were able to explore a continuous range of simulation parameters. In particular, we were aiming to disentangle interrelations between the actual flow parameters, rather than an interrelation to the initial values.
• 3.  
  Altogether, our model setup allows very long-term simulations on a large grid—so far we have run simulations for approximately 30, 000 time units for a standard diffusivity model, and more than 150, 000 time units for our modified strong diffusivity model.
• 4.  
  We allow the inflow (into the coronal region) density to vary in time, thus keeping the ratio between the inflow and the disk densities δ the same as initially.

We have explored a vast range of the parameter space that covers the majority of simulations performed in the literature (see Table 3). Similar to all these papers, we assume a thin disk = 0.1. It is common to assume a magnetic diffusivity parameter α_m of about unity. In our case, we apply values α_m = 1.1...1.9. For the magnetic field bending parameter, we have chosen m = 0.5, which is slightly higher than the values usually adapted, m = 0.35...0.4. Although we finally show that m plays only a minor role, our choice is motivated by the fact that this value is more consistent with the inner boundary condition. The magnetic diffusivity anisotropy parameter χ is chosen to be smaller than unity, since it helps to keep sufficiently strong torques at the midplane and the bending of the magnetic field that supports launching.

Although we begin our simulation from a moderately weak initial magnetic field, the actual field strength in the inner disk at a certain radius may vary substantially over time. This allows us to study the interrelation between disk accretion and ejection physics and the actual magnetic field strength (thus the actual disk magnetization).

We define the disk surface as a surface where the radial velocity changes sign. In a steady state, the area where the radial velocity and the magnetic torque changes sign is almost identical.

Figure 5 shows the typical structure of the disk–jet system. Several, physically different regions can be distinguished—the inflow area, the jet acceleration area, the launching area, and the accretion domain. These are separated by colored lines. The white and black lines mark the disk surface that separates disk and corona regions. Two other lines separate the accretion, launching, and acceleration areas. We define the accretion region as the area where velocity is mainly radial, V_θ < 0.1V_R, and the acceleration region as the area where the flow velocity is parallel to the magnetic field, sin (angle(B, V)) < 0.1. We characterize the region in between as the launching area.

Zoom In Zoom Out Reset image size

In our simulations, the position of the disk surface as defined above remains about constant in time. This confirms our choice of the control volume (see below) and our choice to fix the diffusive scale height during the simulation.

According to our boundary conditions, we prescribe a weak inflow into the region between the inner disk radius and the rotational axis. This inflow provides the matter content as well as the pressure balance along the rotational axis. The astrophysical motivation can be the presence of a central stellar magnetic field or a stellar wind. In Figure 5, the inflow area is the area between the rotational axis and the red line that marks the magnetic field line rooted in the inner disk radius at the midplane. In all figures below, the magnetic field line closest to the axis always corresponds to the magnetic field line anchored at the inner disk radius.³

Here we briefly comment on the jet-launching mechanism in our reference simulation that is—as in previous simulations—the Blandford–Payne magneto-centrifugal driving.

As demonstrated above (see Figure 5), the magnetic torque rF_ϕ changes sign on the disk surface. It is negative in the disk and positive in the corona. Thus, the magnetic field configuration established extracts angular momentum from the disk. The angular momentum extraction relies on the induced toroidal magnetic field component, which plays a key role in transferring the angular momentum ∼B_rB_ϕ. Gaining angular momentum, the material that is loaded to the field lines from the accretion disk is pushed outward by the centrifugal force.

In order to illustrate the acceleration process, we show the magnitude of Lorentz force with respect to the thermal pressure and centrifugal forces. Figure 6 shows the contours where the perpendicular and parallel components of the Lorentz force are equal to the perpendicular and parallel components of the pressure and centrifugal forces, respectively. In the accretion disk, both the pressure and the centrifugal forces dominate the poloidal component of the Lorentz force. Below the disk surface, the Lorentz force (toroidal component) extracts the angular momentum from the disk.

Zoom In Zoom Out Reset image size

Since the Lorentz force increases along the outflow, it is worth checking the decomposed Lorentz force components F = ∇ × B × B in the directions parallel and perpendicular to the magnetic field (Ferreira 1997). The ratio between the toroidal and parallel components of the magnetic field,

$$\begin{equation} \frac{F_{||} }{F_{\rm \phi }} = - \frac{B_{ \phi }}{B_{\rm P}}, \end{equation} \tag{ 21 }$$

is shown in Figure 7. We see that the centrifugal force is stronger in the inner area of the disk rather than in the outer parts of the disk.

Zoom In Zoom Out Reset image size

At the sonic surface, the Lorentz force overcomes the pressure forces. From this point on, the main acceleration force is the centrifugal force. Further along the outflow—between the Alfvén surface and the fast surface—the Lorentz force becomes the main accelerating force.

We now explore the mass flux evolution, namely the accretion and ejection rates.

In the θ direction the control volume is enclosed by the disk surface (as defined above), S_S, and the disk midplane. The two other surfaces that enclose the control volume are marked by S₁ and S_R, and correspond to the vertical arcs at the innermost disk radius (R = 1) and at any other radius R. The control volume defined by these surfaces is denoted by V(R).

Thus, the disk mass enclosed by a radius R follows from

$$\begin{equation} M(R) = 2 \int _{V(R)}\rho dV, \end{equation} \tag{ 22 }$$

while the mass accretion rate at a certain radius R is

$$\begin{equation} \dot{M}_{\rm acc}(R) = - 2 \int _{S_1}^{S_R} \rho {\boldsymbol {V}} _p\cdot d {\boldsymbol {S}}, \end{equation} \tag{ 23 }$$

and the mass ejection rate is integrated along the disk surface,

$$\begin{equation} \dot{M}_{\rm ej}(R) = -2 \int _{S_S} \rho {\boldsymbol {V}} _p\cdot d {\boldsymbol {S}}. \end{equation} \tag{ 24 }$$

The mass accretion rate is defined as being positive if it increases the mass in the control volume. The mass ejection rate is defined as being positive if it decreases the mass of the control volume. The factor of two in front of the integrals takes into account the fact that only one hemisphere is treated. Note also a minus sign in front of the integrals.

It is common to introduce the ejection index ξ, which is based on the mass conservation law for a steady solution (Ferreira & Pelletier 1995). It basically measures the steepness of the radial profile of the accretion rate along the midplane. Setting the outer radius to r and the inner radius to unity, the ejection index interrelates ejection and accretion,

$$\begin{equation} \frac{\dot{M}_{\rm ej}}{ \dot{M}_{\rm acc}} = 1 - r^{-\xi }. \end{equation} \tag{ 25 }$$

We obtain the ejection index by a linear approximation of $\xi = -\log (1-\dot{M}_{\rm ej}/\dot{M}_{\rm acc})/\log (r)$ within r = [2, 10] The higher the ejection index, the higher the fraction of accreted matter being ejected within a given radius, and the less matter reaches the inner boundary. For our reference simulation, ξ ≃ 0.3 at T ≃ 10, 000.

Although the disk continuously loses mass (Figure 8), after dynamical times 1000–2000 the disk mass loss is much smaller that the corresponding ejection and accretion rates. We therefore state that the simulation evolves through a series of quasi-steady states. We find that a continuous disk mass loss is a typical feature of a simulation like our reference simulation. This is because the mass accretion from outside some outer disk radius is not able to sustain the mass that is lost by accretion and ejection within this radius. We also find that the standard diffusivity model typically leads to a magnetic field distribution in the disk that is almost constant in time (not in space). These two facts result in an increase of the disk magnetization in the inner disk, which in turn leads to more rapid accretion in the inner disk.

Zoom In Zoom Out Reset image size

Figure 8 shows the time evolution of the mass accretion and ejection rates. Note also the general decrease of the mass fluxes over time, which is a direct consequence of the decrease of the disk mass. The behavior and actual values of the mass fluxes are typical in the literature. One should note two distinct features of the mass fluxes. First, the higher integration volume, the higher the mass ejection rate. Second, in jet-launching disks, the mass accretion rate must increase with the radius. These plots also indicate that the evolution of the system can be seen as a consecutive evolution through a series of quasi-steady states.

Here we discuss simulations investigating the influence of the initial magnetic field bending parameter m. We have varied m from 0.2 (strongly inclined) to 0.9 (almost vertical).

Our main result is that although the simulations initially evolve slightly differently, in the long-term they are almost indistinguishable.

Figure 9 shows the time evolution of the ejection to accretion mass flux ratio. The fluxes are again computed for the control volume extending to R = 10. It seems to take a few 1000 dynamical time steps for the simulation to lose the memory of the initial magnetic field configuration, but at t = 10, 000 convergence has obviously been reached. This is also true for the fluxes of angular momentum and energy, and holds as well for the corresponding flux ratios.

Zoom In Zoom Out Reset image size

The reason why the simulations convergence into a single, specific configuration is the fact that it is mainly the diffusivity model that governs the evolution of the magnetic field evolution. When we begin the simulations with the same initial magnetic field strength at the midplane, this results in exactly the same magnetic diffusivity profile. Since we explore rather weak magnetic fields (weaker than the equipartition field), the underlying disk structure cannot be changed substantially by the Lorentz force. On the contrary, the magnetic field distribution adjusts itself in accordance with the diffusivity model, which has the same vertical profile ab initio.

The convergence we observe for these simulations, starting from an initial magnetic field with different bending, again confirms the reliability of our model in general.

Here, we briefly present example results of the resolution study. We have performed simulations with a grid resolution of (0.5, 0.75, 1.0, 1.5, 2.0) times our standard resolution of 128 cells per quadrant, corresponding to (64, 96, 128, 192, 256) cells per quadrant, or approximately (8, 12, 16, 24, 32) cells per disk height 2. Note that once the resolution in the θ direction and the radial extent of the disk is chosen, the resolution in the R direction is uniquely determined (see Section 2.1).

For the resolution study, all simulations were performed up to typically 10,000 time units. Figure 10 shows snapshots of some of these simulations. Essentially, we see an almost identical disk-outflow structure, indicating that numerical convergence has indeed been reached.

Zoom In Zoom Out Reset image size

Figure 11 shows the time evolution of the mass ejection to accretion ratio, again integrated throughout the control volume R < 10, for simulations of a different resolution. We note two particular issues. First, all curves bunch together at a mass ejection to accretion ratio of about 0.6, indicating convergence of the simulations (this is also true for other flux ratios). Second, the simulation with the highest resolution shows some intermittent behavior (see Figure 10). This might be related to the spatial reconstruction in non-Cartesian coordinates close to the symmetry axis or to the ability to resolve more detailed structures. We conclude that our simulations have converged for the launching region for the resolution chosen.

Zoom In Zoom Out Reset image size

Note that although the spherical grid is beneficial for a disk and outflow-launching studies, mainly due to the higher resolution of the inner disk, its application to the jet propagation further away from the jet source is limited because of the lack of resolution at larger radii.

Here we present the analysis of the angular momentum and energy transport in our simulations. It is common to explore the angular momentum and energy transport by means of their fluxes through a control volume. We define the accretion angular momentum flux $\dot{J}_{\rm acc}= \dot{J}_{\rm acc,kin}+ \dot{J}_{\rm acc,mag}$ as the sum of the kinetic and magnetic parts, keeping the same control volume as for the mass fluxes (see Appendix B).

Figure 13 shows the time evolution of the jet angular momentum flux for a number of simulations. We divide them into three different groups distinguished by their initial magnetization, μ₀ = 0.003, 0.01, 0.03. Different lines within each group represent simulations with different diffusivity parameters α_m. The jet angular momentum flux is calculated through the upper part of the control volume, thus, up to R = 10.

Zoom In Zoom Out Reset image size

As expected, the stronger the overall magnetic field strength is in the disk, the stronger torque it exerts, and thus the higher angular momentum fluxes we find.

Although the magnetic diffusivity parameter α_m differs within each group of lines (marked by different colors), the total torque measured for the corresponding simulations is about the same. This comes from the fact that the total torque is set by the global magnetic field. Thus, the evolution of the total torque is mainly set by the initial conditions.

We note that the simulations applying a strong initial magnetic field (represented by the upper bundle of curves in Figure 13) have been interrupted earlier compared to usual evolution times. For these cases, the inner part of the disk became highly magnetized. The strong magnetic field changes the inner disk structure such that our model for the magnetic diffusivity can no longer be applied. This is simply because the actual scale height of the disk significantly decreases and no longer coincides with the initial disk surface. We consider the case of strong magnetic fields as being beyond the scope of this paper. The underlying turbulence might be significantly suppressed as well.

In accordance with previous works (see, e.g., Zanni et al. 2007), we find that the ratio of angular momentum extracted by the jet to that provided by the disk accretion is always close to, but slightly larger than, unity, $\dot{J}_{\rm jet} / \dot{J}_{\rm acc}$ ≳ 1. The main reason is that the accretion rate in the outer part of the disk is too low to compete with a strong mass loss by the disk wind at these radii.

We define the accretion energy flux (accretion power) as the sum of the mechanical (kinetic and gravitational), magnetic, and thermal energy fluxes,

$$\begin{equation*} \dot{E}_{\rm acc}= \dot{E}_{\rm acc,kin}+ \dot{E}_{\rm acc,grv}+ \dot{E}_{\rm acc,mag}+ \dot{E}_{\rm acc,thm}, \end{equation*}$$

and, similarly, the jet energy flux (jet power) by

$$\begin{equation*} \dot{E}_{\rm jet}= \dot{E}_{\rm jet,kin}+ \dot{E}_{\rm jet,grv}+ \dot{E}_{\rm jet,mag}+ \dot{E}_{\rm jet,thm}. \end{equation*}$$

In contrast to the angular momentum flux, most of the energy flux is being released from the inner part of the disk. This makes the energy flux very sensitive to the conditions in the inner disk. Indeed, Figure 14 shows that the power liberated by the jet strongly depends on the diffusivity parameter, which is the main agent for governing the magnetic field strength. A weaker diffusivity parameter α_m leads to a higher magnetization, and thus a higher jet power. The same is true for the accretion power.

Zoom In Zoom Out Reset image size

We further find that the ratio of energy fluxes, namely the ratio of the jet to the accretion energy flux, is always close to, but slightly lower than, unity. This is also in accordance with Zanni et al. (2007).

In this section, we have provided some evidence that it is the actual magnetization in the disk that governs the fluxes of mass, energy, and angular momentum. In the next section, we show how exactly these fluxes are connected to magnetization.

In a steady state, diffusion and advection balance. Advection of the magnetic flux in principle increases the magnetic field strength, predominantly in the inner disk. On the contrary, the diffusion smooths out the magnetic field gradient. Therefore, the diffusivity model applied is a key ingredient for these processes, directly influencing the disk structure and evolution.

In the following analysis, we will not focus on the profile or the magnitude of the magnetic diffusivity, but concentrate on the resulting magnetization and its time evolution. As discussed above, by changing the magnetic diffusivity parameter α_m, we are able to explore how the actual disk magnetization influences various properties of the disk–jet system.

For each parameter run performed, we measure the actual physical variables in the disk–jet system, such as the time-dependent mean magnetization, accretion fluxes, jet fluxes, or the accretion Mach number. Naturally, the actual value of a certain property has evolved from the initial value during the simulation. With a mean value we denote the values averaged over a small area of the inner disk,

$$\begin{equation} X \equiv \langle X (r,z=0)\rangle. \end{equation} \tag{ 27 }$$

All mean quantities discussed are averaged over the inner disk midplane from R = 1.1 to R = 1.5. The choice of the averaging area is motivated as follows. First, in order to avoid any influence of the inner accretion boundary we have moved the inner integration radius about 10 grid cells away from it. Second, we are interested in examining the inner part of the disk, since it is the region where the magnetization is changing predominantly and it is the launching area for the most energetic part of the jet. Third, we want to avoid large magnetic gradients affecting our averaging area. Although most of jet energy is launched from regions broader than this small area, the area is seen as representative. We emphasize that the profiles of the jet power along the disk surface are similar for all simulations. In all cases we investigated, we find that the general behavior of the physical outflow or disk quantities with respect to underlying disk magnetization does not depend on the area where the averaging is performed(on neither the location nor the size).

Keeping all other parameters the same, we have carried out simulations varying the initial magnetization μ₀ and the strength of the magnetic diffusivity α_m. For all our simulations, starting with different initial magnetization, μ₀ = 0.003, μ₀ = 0.01, and μ₀ = 0.03, we find that the interrelation between the different jet or disk quantities and the disk magnetization is essentially the same. We therefore present only one group of simulations, namely that with μ₀ = 0.01. Although the initial magnetic field strength differs, we can already suspect at this point that it is the actual rather than the initial strength of the magnetic field in the disk that governs the disk accretion and ejection of the jet, and thus plays a major role in the launching process. This has not yet been discussed in the literature, as most publications parameterize their simulations by the initial parameters. An exception might be Sheikhnezami et al. (2012), who pointed out substantial changes in the disk plasma beta during the time evolution.

An interesting representation of the evolution of the main disk–jet quantities are (μ, X) plots, where X stands for the examining variable. Note that in these plots the time evolution is hidden.

4.2.1. Accretion Mach Number

As was shown by Königl & Salmeron (2011), there is a link between the mean accretion Mach number and the disk magnetization. In our terms, this relation can be expressed as

$$\begin{equation} {M}_{\rm R, act} = \frac{2}{\sqrt{\gamma }} q \mu _{\rm act}, \end{equation} \tag{ 28 }$$

where q is the magnetic shear,

$$\begin{equation} q = \frac{HJ_{ R}}{B_{\rm P}} = -\frac{B_\phi ^{\rm +}}{B_{\rm P}}, \end{equation} \tag{ 29 }$$

and where the plus sign denotes a variable estimated on the disk surface. Note that in our case there is no viscous contribution and the factor $1/\sqrt{\gamma }$ appears in the relation since the accretion Mach number is calculated using adiabatic sound speed.

Figure 15 shows that setting q to constant $q =2\sqrt{\gamma }$ (thus M_R = 4μ) is a good first approximation, especially for the strong magnetization cases. As we will see in the next section, in the case of a weak magnetic field, the magnetic shear q behaves far from being constant. Closer examination of the magnetic shear q reveals a presence of two different jet launching regimes.

Zoom In Zoom Out Reset image size

4.2.2. Magnetic Shear

A tight relation exists between the magnetic shear q and the ratio between the toroidal and poloidal magnetic field components. We define the magnetic shear q by the radial electric current at the disk midplane, since it does not require us to apply the notation of the disk surface, which, in case of a strong magnetic field, might change by up to 40%. Note that the magnetic shear is the first derivative of the accretion Mach number with respect to the magnetization. Therefore, it shows the growth rate of the local Mach number or steepness of the curve.

Figure 16 shows the magnetic shear with respect to the underlying inner disk magnetization. We see that the magnetic shear behaves in two different ways—in the case of low magnetization, μ ⩽ 0.03–0.05, the magnetic shear is substantially higher in comparison to the case of high magnetization, μ ⩾ 0.03–0.05. The explanation is straightforward: there is a turning point concerning the generation of the toroidal magnetic field versus flux losses through the disk surface (by the outflow). To understand this, one needs to set apart the generation processes of the toroidal magnetic field from the loss processes. The rate of the generation of the magnetic field in Keplerian disks is primarily set by the structure of the magnetic field and is rather constant in the case of a quasi-steady state On the other hand, the outflow speed (through the disk surface) is highly dependent on the actual disk magnetization. In the case of a weak magnetization, the outflow speed is rather small, which makes it possible to sustain a stronger magnetic shear. A strong disk magnetization results in a fast outflow, thus setting the maximum limit for the magnetic shear.

Zoom In Zoom Out Reset image size

4.2.3. Mass and Energy Flux

The magnetic shear has a great impact on outflow launching (Ferreira & Pelletier 1995). We confirm this finding by presenting the mass and energy ejection and accretion fluxes.

Figure 17 shows the ratio of the mass ejection rate, $\dot{M}_{\rm ej}(1.5) - \dot{M}_{\rm ej}(1.1)$, to the accretion rate, both averaged over the same area. Obviously, the ejection efficiency is higher for weaker magnetized disks.

Zoom In Zoom Out Reset image size

This is easy to understand considering Equation (21). In the case of a weak magnetic field, the strong magnetic shear (the high toroidal to poloidal magnetic field ratio) leads to faster poloidal acceleration, which is caused by the force component parallel to the magnetic field. This force also extracts the matter from the disk. In the case of a strong disk magnetization, the acceleration of matter is primarily supported by the centrifugal force.

We note that studying the ejection index (calculated within an area R = 2...10) with respect to the mean disk magnetization leads to very similar results, that is a saturation to values of about 0.35–0.40 in the case of high magnetization and a significant increase in the case of low magnetization.

Figure 18 shows the ratio of the energy ejection density to the average accretion energy, computed in the same way as for the mass fluxes. Compared to the mass fluxes, the energies show opposite behavior—the ejection to accretion power is an increasing function with magnetization. This is a highly important relation, since it relates two observables. Note that following our findings, there is not a fixed value for the ratio between the jet and accretion power. This should be considered when comparing observational results to theory.

Zoom In Zoom Out Reset image size

Essentially, this result shows the general importance of the magnetic energy flux compared to the mechanical energy. The mechanical energy flux is always negative, while the magnetic energy flux is positive. In the case of a strong magnetic field, our results are similar to Zanni et al. (2007), namely, that magnetic energy flux dominates the mechanical flux. We also see a saturation of the flux ratio in the case of a moderately strong magnetization. We find that in a weak magnetization case the energy flux ratio can be very small.

We also find that the ejection Mach number, Equation (19), increases almost linearly with magnetization.

Essentially, the general behavior of the mass and energy flux ratios does depend on the averaging area or its position.

The accretion power (see Appendix B, Equation (B3)) is mainly determined by the accretion rate at the inner disk radius. Assuming a typical scale height for the mass accretion as r, the accretion power can be estimated considering the magnetic shear and the actual magnetization of the disk,

$$\begin{equation} {\dot{E}_{\rm acc}} \simeq 0.06\,q\, \mu _{\rm act} \dot{E}_0, \end{equation} \tag{ 30 }$$

where $\dot{E}_0$ is the unit power (see Appendix A).

In the case of strongly magnetized disks, one can assume that the magnetic shear is approximately constant q ≈ 3, and this relation transforms into ${\dot{E}_{\rm acc}} \simeq 0.2 \mu _{\rm act} \dot{E}_0$. Note that this result connects two essential quantities—the accretion power, which manifests itself by the accretion luminosity, to the disk magnetization, which is intrinsically hidden from the observations.

The simple idea that the induction of the magnetic flux in steady state is compensated by the flux losses, both by diffusion and magnetic flux escape through the disk surface, becomes barely applicable in the case of a weak magnetic field. As discussed previously (see Section 2.5), in order to keep the magnetic field distribution properly bent, the magnetic diffusivity parameter α_m must be linked to the anisotropy parameter χ. Equation (20) was obtained considering the standard magnetic diffusivity model. A more general relation can be derived assuming a curvature of the magnetic field of about 0.5, which is the mean curvature of the initial field distribution (see Equation (9)),

$$\begin{equation} \left({\alpha }_{\rm ssm}\chi - M_{\theta }\right){\alpha }_{\rm ssm}\le \mu. \end{equation} \tag{ 31 }$$

Solving this inequality for α_ssm and assuming M_θ∝βμ, we find

$$\begin{equation} {\alpha }_{\rm ssm}\le \alpha _0 = \frac{\beta }{2} \mu _\chi + \sqrt{ \left(\frac{\beta }{2}\right)^2 \mu _\chi ^2 + \mu _\chi }, \end{equation} \tag{ 32 }$$

where μ_χ = μ/χ, and β ≈ 6 in our simulations. This relation shows that in order to keep the disk magnetic field properly bent, the α_ssm should behave differently in the two limits of magnetization: a linear relation to the magnetization in the case of a strong magnetic field, and proportional to the square root of magnetization in the case of a weak field. In the case of a strong deviation from this relation, the magnetic field structure will be substantially affected, resulting in a high field inclination (for α_ssm ≪ α₀), or a strong outward bending (for α_ssm ≫ α₀).

For Equation (32), we have implicitly assumed a linear relation between the ejection Mach number (at the disk surface) and the magnetization. In fact, our simulations approve such an interrelation. Figure 19 shows that for a moderately strong magnetic field, there is a linear relation between the ejection Mach number and the underlying disk magnetization. This behavior is also consistent with the ejection to accretion mass flux ratio we have discussed above.

Zoom In Zoom Out Reset image size

As a consequence of Equation (16), the α_ssm plays a direct role in determining the strength of the poloidal electric current with respect to the toroidal current. We find in our simulations that in order to sustain jets, the ratio of the poloidal to the toroidal current should be sufficiently high (about 15).

The difficulty in performing simulations of weakly magnetized outflows is that the specific torques may increase toward the disk surface area due to the low densities in that area. This will lead to a layered accretion along the disk surface, and thus much lower accretion along the midplane. Although this might be a relevant process in reality, our numerical setup is not suited for the treatment of such a configuration.

The standard, commonly used magnetic diffusivity model is parameterized by two constants, α_m and χ. In general, choosing a high anisotropy parameter χ implies a low diffusivity parameter α_m. Together, this will lead to a decrease in the poloidal electric current and an increase in the toroidal current (see Equation (16)). Thus, the resulting torque will not be sufficient to brake the dense matter along the midplane, and will lead to layered accretion. An anisotropy parameter lower than unity has proven to lead to a smoother time evolution, since it allows for stronger poloidal currents at the midplane, and thus a mass accretion that is developed over the full disk height. Another option to have χ > 1 would be to modify the vertical diffusivity profile such that it reaches the maximum not at the disk midplane, but at the disk surface. This would also help to develop a strong electric current in the disk midplane.

Here, we discuss another problem, which is common in the standard diffusivity model. As previously shown, most of the simulations suffer from the mass loss from the disk, leading to the increase of magnetization. However, the increase of magnetization further amplifies the mass loss. This is known as the accretion instability, first studied by Lubow et al. (1994), and later confirmed for more general cases (Campbell 2009). If, on the other hand, the diffusivity is too high (chosen by a high α_m parameter), the inevitable diffusion of the magnetic field will lead to the situation where a jet can no longer be sustained.

We emphasize that the main reason for the increase in magnetization is the mass loss, but not the actual magnetic field amplification. This is a direct consequence of the accretion instability, namely, a lack of sufficient feedback that could bring the accretion system back to a stable state. In other words, in order to run long-term simulations, one needs to apply a diffusivity model that provides a stronger feedback to the diffusivity profile than the standard choice $\propto \sqrt{2\mu }$.

The direct consequence of the accretion instability is that the magnetization increases toward the center. It is easy to show that the behavior of the magnetization has to be opposite. In accretion disks producing outflows, the mass accretion rate must naturally increase with radius. Assuming a radial self-similarity of the disk and taking into account that the accretion Mach number is linearly related to the magnetization and that $\rho \propto C_{\rm s}^3$, one derives

$$\begin{equation} \beta _{\mu } = \xi -2 - 4\beta _{C_{\rm s}}, \end{equation} \tag{ 33 }$$

where ξ is the ejection index and β_X represent the power-law index of a physical quantity X. Considering that magnetized disks are very efficient in producing outflows, ξ ≈ 0.2–0.4, one may expect β_μ to be positive (if $\beta _{C_{\rm s}} \approx -1/2$), thus, an increasing function with radius. However, the disk structure itself can be re-arranged such that $|{\beta _{C_{\rm s}}}| \le 1/2$, which eventually will satisfy the relation 33. This is indeed what we find.

One should, however, keep in mind that this equation is a rough estimate and might be subject to the different disk physics involved. If the magnetic torque is not the only supporter of the accretion, as in the case of the viscous simulations of Murphy et al. (2010), the above presented relation might be relaxed. However, similar analysis can be performed to set the limit of the magnetization with respect to other quantities.

In this section, we present a diffusivity model, which does not suffer from the accretion instability. Although the standard diffusivity model gave us the opportunity to probe a wider parameter space, it is not applicable for very long-term studies.

We emphasize that the transition from the a direct simulation of turbulence to the mean field approach, which lacks the small scales by design, is indeed subtle. So far, in the literature, the jet launching problem is addressed without considering the origin of the magnetic field by a dynamo. Therefore, the only way of amplifying the magnetic field is by advection (or stretching in case of a toroidal field). There is also no intrinsic angular momentum transport by the turbulence itself. The only term we need to model when applying a small-scale turbulence is the effective magnetic diffusivity. This might have surprising consequences. In order to suppress the turbulence, one should rather amplify the effective diffusivity—leading to a stronger decay of the magnetic field and resembling the quenching of diffusivity (or dynamo)—rather than decrease it, leading to stronger advection and thus an amplification of the magnetic field. The main motivation of our new model for the diffusivity is to consider stronger feedback by the disk magnetization. We overcome the accretion instability by assuming a stronger dependence of α_ssm on the magnetization,

$$\begin{equation} {{\alpha }_{\rm ssm}} = {\alpha }_{\rm m}\sqrt{2\mu _0} \left(\frac{\mu }{\mu _0}\right)^2, \end{equation} \tag{ 34 }$$

where we choose α_m = 1.55 and μ₀ = 0.01. Here we keep the previous overall form and constants, indicating that both models are the same at the magnetization μ₀. A quadratic dependence on μ was chosen in order to amplify the feedback. Choosing a power lower than two might have revealed other complications, for example, feedback too weak to work fast enough, and keep α_ssm under the constraint of Equation (32). We will refer to this diffusivity model as the strong diffusivity model.

Applying our strong diffusivity model enables us to overcome the accretion instability. As a result, we were able to perform our simulations for much longer times, reaching evolutionary time steps of t > 150, 000, which corresponds to approximately 25,000 revolutions at the inner disk orbit.

Figure 20 shows the typical computation domain and the initial dynamics of the system. As usual, the evolution starts with the propagation of the toroidal Alfvén wave, resulting in a propagating cocoon. At this point, the innermost area (R ⩽ 10) has reached a quasi-steady state, while the outer part has not even slightly moved from the initial state. Figure 21 presents the same simulation at a later state when a strong outflow has been developed. Notice the difference between the inner part, r ⩽ 200, and the outer part, r ⩾ 200. The inner part has already relaxed to a steady state, while the outer region shows rapid accretion and ejection patterns. This is a direct consequence of the new diffusivity model. The logic behind implementing the enhanced feedback is valid only when the accumulation of the flux is possible. The initial imbalance between advection and diffusion in the outer part leads to a rapid advection of the magnetic flux to the inner disk. As a result, the rapid accretion further leads to higher inclination angles of the magnetic field (smaller angle with respect to the disk surface). This results in a higher efficiency for the toroidal magnetic field induction, thus leading to even more rapid accretion and ejection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 22 shows the long-term evolution on a small sub-grid of the simulation with our strong diffusivity model. As can be seen, until time t = 10, 000, only a small fraction of the disk has dynamically evolved (up to R ≈ 50), while at later times the outer parts of the disk also reach a new dynamic state. A steady outflow is established from the whole disk surface (shown on this sub-grid) and reaches a super-fast magnetosonic speed. Note that the positions of the critical MHD surfaces are constant in time, which is a further signature of a steady state.

Zoom In Zoom Out Reset image size

The outflow reaches maximum velocities typically of the order of 100 km s⁻¹ for YSO, or 70,000 km s⁻¹ in case of AGNs. Concerning observationally relevant scales, our simulations compare to the following numbers. Our numerical grid is comparable to 150 AU for YSO, and 0.14 pc in the case of AGNs Physically more meaningful is the grid size where our simulation has reached a steady state. This is a size comparable to 25 AU for YSO, and 5000 AU in the case of AGNs, but can be extended by running the simulations longer. The dynamical timescale of 150,000 time units (or 25,000 disk orbits) corresponds to about 550 years in the case of YSO and about 200 years in the case of AGNs. Typical accretion rates of our simulations are 3 × 10⁻⁷ M_☉ yr⁻¹ for YSO and 0.1 M_☉ yr⁻¹ in the case of AGNs, but one has to keep in mind that these values depend not only on the intrinsic scaling of the central mass and inner disk radius, but also on the scaling of density. Therefore, we herewith present the most extended and longest MHD simulations of jet launching obtained so far—connecting the jet launching area close to the central object with the asymptotic domain, which is accessible by observations.

Although a spherical grid is computationally very efficient and may allow us to extend the computational domain to almost any radius, in reality its application for the jet-launching simulations is somewhat limited. There are two reasons for that. First, it takes obviously much longer time for the outer disk areas to evolve into a new dynamical steady state. Thus, the outer disk will remain close to the initial state of the simulation for quite some time. Second, a more severe drawback is the lack of resolution for the asymptotic jet. For example, for distances R > 500R₀ along the rotational axis, a jet radius of about r_jet ≈ 25 can be resolved only by about five grid cells (applying our typical resolution). We therefore restrict our computational domain for such grid size to about R_out ≃ 1000–2000R₀. The above-mentioned resolution issue is in fact one of the advantages for using cylindrical coordinates for jet formation simulations.

By design, the purpose of our strong diffusivity model was to avoid the accretion instability. As a consequence of this application, the magnetization profile does not decrease with radius (see Figure 23) Although both simulations (with standard and strong diffusivity model) start from the same initial disk magnetization, the disk evolution results in a quite different magnetization distribution. The standard diffusivity model (our reference simulation from above) results in a magnetization profile decreasing with radius. In contrast, for the strong diffusivity model, a rather flat magnetization profile emerges. In non-viscous simulations, assuming radial self-similarity, a flat (or not decreasing with radius) profile is essential for sustaining a continuous accretion flow at any given radius.

Zoom In Zoom Out Reset image size

As soon as a steady state is reached, the evolutionary track for this simulation is represented by a simple dot in all (μ, X) diagrams (at least from 1000 to 150,000 time units). The mean inner disk magnetization is μ ≈ 0.012. We find that this simulation fits to every relation presented above, such as mass and energy flux ratios, or magnetic shear, that were derived applying a standard diffusivity model. In other words, the aforementioned dots belong to the curves drawn in the (μ, X) diagrams.

There are several distinct features one can derive from the resulting magnetic field structure. Figure 24 shows the toroidal to poloidal magnetic field component ratio. Taking into account that the disk magnetization (calculated from poloidal component only) is uniform, three different regions can be distinguished. The first region is between the midplane and the disk surface where the toroidal magnetic field reaches its maximum and the torques change sign. The second region is between the disk surface and the Alfvén surface where the ratio of the field components is quite constant. The third region is beyond the Alfvén surface when the poloidal component of the magnetic field becomes weak enough to keep a rigid magnetic field structure and the toroidal component starts to dominate the poloidal one.

Zoom In Zoom Out Reset image size

In this subsection, we further explore the disk structure in a steady state. In Figure 25, we present the radial profiles of certain magnetohydrodynamical variables along the midplane. We show the profiles derived from our numerical simulations with their approximations by power laws, and compare them to the initial distribution. These radial profiles are obtained along the disk midplane, however, they also hold for at least one disk semi-height. The thetoidal profiles that normalized to the corresponding midplane value (not shown here) almost coincide with each other, indicating that the assumption of a self-similar disk is in fact reasonable, though different power indexes should be used.

Zoom In Zoom Out Reset image size

In particular, Figure 25 shows how the disk structure evolves from a certain initial power-law distribution into another power-law profile. We see distinct power-law profiles for radii up to R ⩽ 250. This corresponds to the area where the disk evolution has reached a steady state. For very small radii, R ⩽ 1.1, we also see a deviation from a power-law profile, which we consider as a boundary effect.

At time t = 150, 000, we find the following numerical values for the power-law coefficients β_X for the different variables at the midplane $X(r, \theta = \pi /2) \sim r^{\beta _X}$. The disk rotation remains Keplerian throughout the entire evolution, thus $\beta _{V_{\rm \phi }} = -1/2$. The radial profiles for density and gas pressure slightly change from their initial distribution. The density power-law index changes from β_ρ = −3/2 to β_ρ = −4/3, while the pressure changes from β_P = −5/2 to β_P = −20/9. For the accretion velocity, we find a profile of $\beta _{V_{R}} = -2/5$, and $\beta _{B_{ \theta }} = -10/9$ for the magnetic field. As a consequence, the profile for the magnetization remains about constant $\beta _{\mu } \sim 2\beta _{B_{ \theta }} / \beta _{\rm P} = (-20/9) / (20/9)$. The accretion velocity remains subsonic over the whole disk with an accretion Mach number of V_R/C_s ≃ 0.1.

Following Ferreira & Pelletier (1995) and considering the mass accretion $\dot{M}_{\rm acc}\sim R^2 \rho V_{R}$, it is easy to get the ejection index ξ = 0.26. This is in accordance with previous work Sheikhnezami et al. (2012).

Our strong diffusivity model, Equation (34), does not necessarily lead to a flat magnetization profile. The model does not directly force the magnetization to be uniformly distributed—in contrast, the profile is expected to be outwardly increasing. As we previously showed, the magnetization profile is linked to the ejection index and it has to be positive if the sound speed remains as initially distributed. What we found is that the disk hydrodynamics changes such that the magnetization of the disk remains flat, thus satisfying Equation (33).

Another way of reasoning is the following. In the case of a flat magnetization profile, the result of the simulation is no longer sensitive to the diffusivity model. In other words, it is possible to switch back from the new diffusivity model to the standard model when the radial magnetization profile is uniform (along the midplane). However, the diffusivity parameter α_m has to be correspondingly re-normalized. The only substantial deviation from a uniform profile is in the innermost disk, which might be influenced by the boundary condition. In fact, we performed such a test simulation, switching back from a strong diffusivity to a standard diffusivity model. As expected, the accretion instability begins to manifest itself similarly as before, leading to the typical magnetization profile (increasing toward the center). However, it takes a much longer time to substantially affect the outer parts of the disk.

The steady state solution achieved when using our strong diffusivity model perfectly fits the results obtained by the standard model (shown previously as a dot in the plots). This actually approves our understanding that the main agent in driving outflows is the actual magnetization, and that the magnetic diffusivity is only the mediator through which the magnetic field structure is being governed. A self-consistent treatment of turbulence is not yet feasible in the context of outflow launching. Therefore, what should be considered first is the resulting magnetic field strength and distribution, but not the diffusivity model itself.