SIMULATING PROTOSTELLAR JETS SIMULTANEOUSLY AT LAUNCHING AND OBSERVATIONAL SCALES

Jets and outflows from protostellar objects are fundamental aspects of the current star formation paradigm, and are observed anywhere star formation is ongoing. The mechanism proposed by Blandford & Payne (1982), in which jets are launched from accretion disks by gravitational, magnetic, and centrifugal forces, has been extensively studied numerically (e.g., Uchida & Shibata 1985; Meier et al. 1997; Ouyed & Pudritz 1997a, 1997b, 1999; Krasnopolsky et al. 1999; Vitorino et al. 2002; von Rekowski et al. 2003; Ouyed et al. 2003; Porth & Fendt 2010; Staff et al. 2010). By treating the accretion disk as a boundary condition (e.g., Ustyugova et al. 1995), one can study jet dynamics independently of the disk (e.g., Pudritz et al. 2007) though, in order to resolve the launching mechanism, numerical simulations have not followed the jet beyond 100 AU (e.g., Anderson et al. 2005).

In stark contrast, protostellar jets are ≳10⁴ AU long (Bally et al. 2007) and only recently have observations reached within 100 AU of the source (e.g., Hartigan et al. 2004; Coffey et al. 2008). This large-scale difference between observations and simulations makes direct comparisons difficult and, in this work, we aim to close this gap. We present axisymmetric (2.5-dimensional) simulations of protostellar jets launched from the inner AU of a Keplerian disk and follow the jet well into the observational domain (2500 AU). These calculations allow us to address the efficacy of the magnetocentrifugal mechanism and to relate conditions near the disk with directly observable properties of the jet.

The simulations presented herein are performed with an adaptive mesh refinement (AMR) version of ZEUS-3D (Clarke 1996, 2010) called AZEuS (Adaptive Zone Eulerian Scheme). The ZEUS-3D family of codes is among the best tested, documented, and most widely used astrophysical magnetohydrodynamic (MHD) codes available, though this is the first attempt to couple ZEUS-3D with AMR.¹ We have implemented the block-based method of AMR detailed in Berger & Colella (1989) and Bell et al. (1994). Significant effort was spent minimizing errors caused by passing waves across grid boundaries, which is of particular importance to this work. A full description of the code and the changes required for AMR on a fully staggered mesh will appear in J. P. Ramsey & D. A. Clarke (2011, in preparation).

Observationally, the inner radius of a protostellar accretion disk, r_i, is between 3 and 5 R_* (Calvet et al. 2000) and, for a typical T Tauri star (M = 0.5 M_☉, R_* = 2.5 R_☉), r_i = 0.05 AU. Thus, following Ouyed & Pudritz (1997a), we initialize a hydrostatic, force-free atmosphere surrounding a 0.5 M_☉ protostar coupled to a rotating disk with r_i = 0.05 AU. However, unlike Ouyed & Pudritz we use an adiabatic equation of state that conserves energy across shocks rather than an isentropic polytropic equation of state, as the distinction becomes important for supermagnetosonic flow (Ouyed et al. 2003).

We solve the equations of ideal MHD² (γ = 5/3) over a total domain of 4096 AU × 256 AU. To span the desired length scales, nine nested, static grids (refinement ratio 2) are initialized each with an aspect ratio of 4:1 (16:1 for the coarsest grid only) and bottom left corner at the origin. Our finest grid has a domain 4 AU × 1 AU and a resolution Δz = r_i/8 = 0.00625 AU which we find sufficient to resolve the launching mechanism. Thus, the effective resolution for the entire domain is >26 billion zones. The simulation highlighted in Section 3 was run to t = 100 yr with an average time step in the finest grid of ∼3 minutes and thus ∼18 million time steps.

During the simulations, a thin region of low velocity and high poloidal magnetic field, $B_{\rm p}=\sqrt{B_z^2+B_r^2}$, develops along the symmetry axis, the edge of which is defined by a large gradient in the toroidal magnetic field, ∂_rB_φ. Insufficient resolution of ∂_rB_φ can lead to numerical instabilities, and grids are added dynamically whenever this gradient is resolved by fewer than five zones.

2.1. The Atmosphere

The atmosphere is initialized in hydrostatic equilibrium (HSE; v_z = v_r = v_φ = 0). Because the left-hand side of the equation governing HSE,

$$\begin{equation} \nabla {p}+\rho \nabla \phi =0, \end{equation} \tag{ 1 }$$

is not a perfect gradient, differencing it directly on a staggered-mesh can commit sufficient truncation error to render the atmosphere numerically unstable. Thus, we replace ∇ϕ with the corresponding poloidal gravitational acceleration vector,

$$\begin{equation} \vec{g}=-\frac{1}{\rho _{\rm h}}\nabla {p_{\rm h}}, \end{equation} \tag{ 2 }$$

where ρ_h and p_h are the hydrostatic density and pressure given by

$$\begin{equation} \rho _{\rm h}=\rho _{\rm i}\left(\frac{r_{\rm i}}{\sqrt{r^2+z^2}}\right)^{\frac{1}{\gamma -1}}\quad {\rm and}\quad \quad {p}_{\rm h}=\frac{p_{\rm i}}{\gamma }\left(\frac{\rho _{\rm h}}{\rho _{\rm i}}\right)^{\gamma }. \end{equation} \tag{ 3 }$$

Here, ρ_i and p_i are the initial density and pressure at r_i and p ∝ ρ^γ is assumed throughout the atmosphere at t = 0. In this way, differencing Equation (1) maintains HSE to within machine round-off error indefinitely.

However, Equations (3) as given are singular at the origin where truncation errors are significant regardless of resolution. These errors can launch a supersonic, narrow jet from the origin destroying the integrity of the simulation. To overcome this problem, we replace the point mass at the origin with a uniform sphere of the same mass and a radius R₀, thus modifying the first of Equations (3) to

$$\begin{equation} \left(\frac{\rho _{\rm h}}{\rho _{\rm i}}\right)^{\gamma -1}=\left\lbrace {\begin{array}{@{}ll@{}}\displaystyle \frac{r_{\rm i}}{\sqrt{r^2+z^2}},&{r^2+z^2\ge {R_0}^2};\\ \ms {\displaystyle \frac{r_{\rm i}}{R_0}\,\frac{3{R_0}^2-r^2-z^2}{2{R_0}^2},}&{r^2+z^2<{R_0}^2}. \end{array}}\right. \end{equation} \tag{ 4 }$$

If R₀ is sufficiently resolved (e.g., four zones), the numerical jet is eliminated. The resulting "smoothed potential" is superior to a "softened potential" since the former has no measurable effects beyond R₀. Here, we use R₀ = r_i.

The atmosphere is initialized with the force-free magnetic field used by Ouyed & Pudritz (1997a):

$$\begin{equation} \begin{array}{rcl}A_{\varphi }&\!\!\!=\!\!\!&{\displaystyle \frac{B_{\rm i}}{\sqrt{2-\sqrt{2}}}\ \frac{\sqrt{r^2+(z+z_{\rm d})^2}-(z+z_{\rm d})}{r}};\\ [18pt] B_z&\!\!\!=\!\!\!&{\displaystyle \frac{1}{r}\frac{\partial (rA_{\varphi })}{\partial {r}},\quad \quad {B}_r=-\frac{\partial {A}_{\varphi }}{\partial {z}}},\quad \quad {B}_{\varphi }=0, \end{array} \end{equation} \tag{ 5 }$$

where A_φ is the vector potential, z_d is the disk thickness (set to r_i), and B_i is the magnetic field strength at r_i, given by

$$\begin{equation} B_{\rm i}=\sqrt{\frac{8\pi {p}_{\rm i}}{\beta _{\rm i}}}. \end{equation} \tag{ 6 }$$

Here, p_i and β_i (plasma beta at r_i) are free parameters.

Finally, to ensure that the declining density and magnetic field profiles do not fall below observational limits, we add floor values ρ_floor ∼ 10⁻⁶ρ_i and B_z,floor ∼ 10⁻⁵B_i (cf. Bergin & Tafalla 2007; Vallée 2003) to Equations (4) and (5). By imposing HSE and the adiabatic gas law at t = 0, a floor value on ρ imposes effective floor values on $\vec{g}$ and p as well.

2.2. Boundary Conditions

2.3. Scaling Relations

From Equation (1) and the adiabatic gas law, one can show

$$\begin{equation} c_{\rm s}^2=\gamma \frac{p}{\rho }=\left(\gamma -1\right)\frac{GM_{\rm *}}{R}=\left(\gamma -1\right)v_{\rm K}^2, \end{equation} \tag{ 7 }$$

where R is the spherical polar radius. From Equations (6) and (7), and the ideal gas law (p = ρkT/〈m〉, where 〈m〉 is half a proton mass), we derive the following scaling relations:

$$\begin{eqnarray} &\displaystyle p_{\rm i}=(160 \;{\rm \;dyne \;cm^{-2}})\left(\frac{\beta _{\rm i}}{40}\right)\left(\frac{B_{\rm i}}{10\;{\rm G}}\right)^2;& \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \displaystyle \frac{\rho _{\rm i}}{\langle {m}\rangle }&=&(5.4\times 10^{12}\;{\rm cm}^{-3})\,\left(\frac{\beta _{\rm i}}{40}\right)\,\left(\frac{B_{\rm i}}{10\;{\rm G}}\right)^2\left(\frac{r_{\rm i}}{0.05\;{\rm AU}}\right)\nonumber\\ &&\times \left(\frac{0.5\;M_{\odot }}{M_{\rm *}}\right); \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &\displaystyle T_{\rm i}=(2.2\times 10^5\;{\rm K})\,\left(\frac{0.05\;{\rm AU}}{r_{\rm i}}\right)\left(\frac{M_{\rm *}}{0.5\;M_{\odot }}\right);& \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &\displaystyle c_{\rm s,i}=(77\;{\rm km\;s}^{-1}) \left(\frac{0.05\;{\rm AU}}{r_{\rm i}}\right)^{1/2}\left(\frac{M_{\rm *}}{0.5\;M_{\odot }}\right)^{1/2};& \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &\displaystyle \tau _{\rm i}=\frac{r_{\rm i}}{c_{\rm s,i}}=(9.7\times 10^4\;{\rm s}) \left(\frac{r_{\rm i}}{0.05\;{\rm AU}}\right)^{3/2}\left(\frac{0.5\;M_{\odot }}{M_{\rm *}}\right)^{1/2},&\nonumber\\ \end{eqnarray} \tag{ 12 }$$

for γ = 5/3. Note that β_i is the only free parameter varied in this work.

Figure 1 depicts a jet with β_i = 40 at t ≃ 100 yr from the highest resolution grid near the disk surface (bottom panel) to the coarsest grid in which the jet has reached a length of just under 2500 AU (top panel).³ A few features worth noting include the following.

• 1.  
  When θ < 60° (angle between $\vec{B}_{\rm p}$ and disk surface), Blandford & Payne (1982) show that cold gas near the disk is launched into a collimated outflow. Here, θ < 60° for all r>r_i, but significant outflow is limited to inside the point where the slow surface intersects the disk (r_j,d ∼ 30 AU = jet radius at the disk; the second panel from top). Below r_j,d, cold disk material has moved onto the grid and accelerated into the outflow. Above r_j,d, the weak magnetic field has yet to drive enough disk material onto the grid to displace the hot atmosphere, and outflow is stifled. While r_j,d gradually increases with time, the majority of mass flux originating from the disk is driven within r_i < r < 10 r_i (0.5 AU; bottom panel).
• 2.  
  Jet material becomes superfast (M_f ≲ 5) within a few AU of the disk, and the boundary between jet and entrained ambient material is defined by a steep temperature gradient (contact discontinuity; second panel). Portions of the original atmosphere, which remain virtually stationary throughout the simulation, are still visible above and ahead of the bow shock (top panel).
• 3.  
  At large distances from the disk (≳500 AU; top panel), the dynamics of the jet become dominated by B_φ, and the jet is led by an essentially ballistic, magnetic "nose-cone" with a Mach number of ∼10 (e.g., Clarke et al. 1986). Still, B_φ is a small fraction (10⁻³) of B_i, consistent with Hartigan et al. (2007).
• 4.  
  The knots dominating the bottom panel (cf., Ouyed & Pudritz 1997b) are produced by the nearly harmonic oscillation of $\vec{B}_{\rm p}$ in r_i < r < 2 r_i, whereby θ fluctuates between 55° and 65° with a period ∼30 τ_i. These oscillations result from the interplay between infalling material along the symmetry axis, and under/over pressurisation near the central mass. The knots are denser and hotter than their surroundings, and bound by magnetic field loops. They occupy a region within ∼2 AU of the symmetry axis and gradually merge to form a continuous and narrow column of hot, magnetized material⁴ (third panel). As such, they are unlikely to be the origin of the much larger-scale knots observed in some jets (e.g., HH111; Raga et al. 2002).

Zoom In Zoom Out Reset image size

Further details of this and other simulations of protostellar jets are left to a future paper, and we focus here on a few properties directly comparable with observations.

Table 1 summarizes a few observational characteristics of protostellar jets. To connect these attributes to conditions in the launching region, we have performed a small parameter survey in β_i and made numerical measurements of the quantities in Table 1. Variation of other parameters (such as ζ and ρ_i) is left to future work.

Note that β_i is the initial value of the plasma beta at r_i, and not the average β in the jet. Indeed, Figure 2(a) demonstrates that at very early time, 〈β〉 = 8π〈p〉/〈B²〉 ≲ β_i/5, where B² = B²_p + B²_φ, and where the average is taken over zones that exceed a certain threshold v_z so that only outflowing jet material is considered. Thus, the magnetic field within the jet is stronger than β_i would suggest. Initially, 〈β〉 is dictated by B_p, but becomes dominated by B_φ within ≲10 yr after launch. As time progresses, 〈β〉 gradually increases but never rises above unity (at least for t < 100 yr), even for β_i ≫ 1. Still, one might speculate from Figure 2(a) that with sufficient time, 〈β〉 → 1 regardless of β_i.

Zoom In Zoom Out Reset image size

4.1. Proper Motion

For t ≳ 10 yr, the velocity of the tip of the jet, v_jet, attains a steady value and, from Figure 2(b) and Table 2, we find v_jet ∝ B^0.44±0.01_i.

Table 2. Simulation "Observables" v_jet and 〈v_φ〉 are Asymptotic Values While r_jet and $\dot{M}_{\rm jet}$ are Measured at z = 200 AU and t = 20 yr

βi	160	40	10	2.5	1.0	0.4	0.1	α
Bi (G)	5	10	20	40	63.2	100	200
vjet (km s−1)	84	125	161	230	270	330	460	0.44 ± 0.01
〈vφ〉 (km s−1)	2.6	3.0	6.2	10.1	13.1	18.4	31	0.66 ± 0.01
2 rjet (AU)	21	40	60	85	94	104	130	0.35 ± 0.04
$\dot{M}_{\rm jet}$ (10−6 M☉ yr−1)	0.44	1.9	2.8	4.2	6.9	10.1	17.9	0.92 ± 0.09

Note. Uncertainties in α are from the fitting procedure.

Download table as:  ASCIITypeset image

To understand this result physically, we begin with the magnetic forces:

$$\begin{eqnarray} F_{\parallel } & = -\frac{B_\varphi }{r} \nabla _{\parallel }(rB_\varphi); \nonumber \\ F_\varphi & = \frac{B_{\rm p}}{r} \nabla _{\parallel }(rB_\varphi); \\ F_{\perp } & = -\frac{B_\varphi }{r} \nabla _{\perp }(rB_\varphi) + J_\varphi B_{\rm p} \nonumber \end{eqnarray} \tag{ 13 }$$

(e.g., Ferreira 1997; Zanni et al. 2007), where ∇_∥ and ∇_⊥ are the gradients parallel and perpendicular to $\vec{B}_{\rm p}$, respectively. For a given field line, a stronger B_p at its "footprint" in the disk (r = r₀) generates a stronger B_φ which leads to stronger gradients in rB_φ and thus, from Equations (13), greater magnetic forces to accelerate the flow. In practice, we find that most of the acceleration occurs before the fast point (and not the Alfvén point) located at r = r_f, where r_f is a weak function of the field strength at the footprint and thus of B_i.

Following Spruit (1996), one can show that as a function of the "fast moment arm" (ξ ≡ r_f/r₀), the poloidal velocity at the fast point is

$$\begin{equation} v_{\rm p,f}=\sqrt{a_{\rm p,f}\,v_{\rm K,0}}\left(\xi ^2+\frac{2}{\xi }-3\right)^{1/4}\propto \sqrt{B_{\rm i}}, \end{equation} \tag{ 14 }$$

since a_p,f, the poloidal Alfvén speed at the fast point, is roughly proportional to B_i. $v_{\rm K,0}=\sqrt{GM_{\rm *}/r_0}$ is the Keplerian speed at the footprint of the field line. We note that measured values of v_p,f in our simulations vary as B^0.5_i and agree with Equation (14) to within 1% so long as the fluid is in approximate steady state.⁵

After the poloidal force given by Equations (13) decreases to 1% of its maximum value (≳a few r_f), v_p still follows a power law in B_i with index 0.52 ± 0.04 and essentially unchanged from Equation (14). Nearer the head of the jet where steady state is no longer valid, we find 〈v_p〉 ∝ B^0.45±0.02_i (where the momentum-weighted average is taken across the jet radius), only slightly shallower than Equation (14). Thus, while the conditions in the jet have changed, some memory of the steady-state conditions at r_f persists.

Finally, v_jet (Figure 2(b) and Table 2) is within ∼10% of 〈v_p〉 near the bow shock and maintains the same power-law dependence on B_i. Thus, these jets are essentially ballistic, where the observed jet speed v_jet ∝ B^0.44±0.01_i. In short, all measures of jet speed increase with B_i, a trend that agrees with Anderson et al. (2005) who find for much less evolved jets, v_p ∝ B^1/3_i.

4.2. Toroidal Velocity

4.3. Jet Radius and Mass Flux

We have presented the first MHD simulations of protostellar jets that start from a well-resolved launching region (Δz_min = 0.00625 AU) and continue well into the observational domain (2500 AU). On the AU scale, each jet shows the characteristic and near steady-state knotty behavior first reported by Ouyed & Pudritz (1997b), though the origin of our knots is quite different. On the 1000 AU scale, each jet develops into a ballistic, supersonic (8 ≲ M ≲ 11) outflow led by a magnetically confined "nose-cone" (Clarke et al. 1986) and a narrow bow shock, consistent with what is normally observed.

On comparing Tables 1 and 2, our simulations comfortably contain virtually all observed protostellar jets on these four important quantities. We note that these tables would not have been in agreement had we stopped the jet at, say, 100 AU and measured these values then. It is only because our jets have evolved over five orders of magnitude in length scale that we can state with some confidence that the magnetocentrifugal launching mechanism is, by itself, capable of producing jets with the observed proper motion, rotational velocity, radius, and mass outflow rate. Indeed, our jets are still very young, having evolved to only 100 yr, and allowing them to evolve over an additional one or two orders of magnitude in time may still be useful. For example, it would be interesting to know whether 〈β〉 rises above unity for any of the jets (Figure 2(a)), and thus enter into a hydrodynamically dominated regime. It would also be interesting to see how long it takes for the power laws in jet radius and mass flux as a function of B_i to reach their asymptotic limits.

Our jet widths tend to be higher than those observed, particularly when one considers that the values for r_jet in Table 2 are at t = 20 yr,⁶ and that r_jet continues to grow in time (e.g., for the β_i = 40 jet, 2 r_jet ∼ 100 AU by t = 100 yr). As our jet radii mark the locations of the contact discontinuity while observed radii mark hot, emitting regions, our widths should be considered upper limits. That our values contain all observed jet widths is a success of these simulations.

Similarly, our numerical mass fluxes are higher than observed values by at least an order of magnitude. Since observed mass-loss rates account only for emitting material (e.g., in forbidden lines; Hartigan et al. 1994), and thus temperatures in excess of 10⁴ K (Dyson & Williams 1997, p. 104), our mass fluxes are necessarily upper limits as well. Indeed, if we measure our mass fluxes near the jet tip (instead of at 200 AU for Table 2) and restrict the integration to fluid above 10⁴ K, our mass fluxes drop by a factor of 10–100, in much better agreement with Table 1.

We thank the referee for timely and helpful comments on the manuscript, Marsha Berger for her AMR subroutines, and Sasha Men'shchikov for early work on AZEuS. Use of MPFIT by C. B. Markwardt and JETGET by J. Staff, M. A. S. G. Jørgenson, and R. Ouyed is acknowledged. This work is supported by NSERC. Computing resources were provided by ACEnet which is funded by CFI, ACOA, and the provinces of Nova Scotia, Newfoundland & Labrador, and New Brunswick.