Simulating protostellar jets simultaneously at launching and observational scales
We present the first 2.5-D MHD simulations of protostellar jets that include both the region in which the jet is launched magnetocentrifugally at scale lengths AU, and where the propagating jet is observed at scale lengths AU. These simulations, performed with the new AMR-MHD code AZEuS, reveal interesting relationships between conditions at the disc surface, such as the magnetic field strength, and direct observables such as proper motion, jet rotation, jet radius, and mass flux. By comparing these quantities with observed values, we present direct numerical evidence that the magnetocentrifugal launching mechanism is capable, by itself, of launching realistic protostellar jets.
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 discs by gravitational, magnetic, and centrifugal forces, has been extensively studied numerically (e.g., Uchida & Shibata 1985; Meier et al. 1997; Ouyed & Pudritz 1997a, b, 1999; Krasnopolsky et al. 1999; Vitorino et al. 2002; von Rekowski et al. 2003; Ouyed, Clarke, & Pudritz 2003; Porth & Fendt 2010; Staff et al. 2010). By treating the accretion disc as a boundary condition (e.g., Ustyugova et al. 1995), one can study jet dynamics independently of the disc (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 AU long (Bally, Reipurth, & Davis, 2007), and only recently have observations reached within 100 AU of the source (e.g., Hartigan, Edwards, & Pierson, 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-D) simulations of protostellar jets launched from the inner AU of a Keplerian disc, and follow the jet well into the observational domain (2500 AU). These calculations allows us to address the efficacy of the magnetocentrifugal mechanism, and to relate conditions near the disc 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 are among the best tested, documented, and most widely used astrophysical MHD codes available, though this is the first attempt to couple ZEUS-3D with AMR111ENZO, a hybrid N-body Eulerian code (O’Shea et al., 2004), links AMR with the hydrodynamical portion of ZEUS-2D.. We have implemented the block-based method of AMR detailed in Berger & Colella (1989) and Bell et al. (1994). Significant effort was spent minimising 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 Ramsey & Clarke (in preparation).
Observationally, the inner radius of a protostellar accretion disc, , is between 3–5 (Calvet et al., 2000) and, for a typical T Tauri star (, ), AU. Thus, following Ouyed & Pudritz (1997a), we initialise a hydrostatic, force-free atmosphere surrounding a protostar coupled to a rotating disc with 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, Clarke, & Pudritz, 2003).
We solve the equations of ideal MHD222AZEuS solves either the total or internal energy equation. We chose the latter because positive-definite pressures trump strict conservation of energy in these simulations; see Clarke (2010). () over a total domain of . To span the desired length scales, nine nested, static grids (refinement ratio 2) are initialised 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 and a resolution which we find sufficient to resolve the launching mechanism. Thus, the effective resolution for the entire domain is billion zones. The simulation highlighted in §3 was run to yr with an average time step in the finest grid of minutes and thus million time steps.
During the simulations, a thin region of low velocity and high poloidal magnetic field, , develops along the symmetry axis, the edge of which is defined by a large gradient in the toroidal magnetic field, . Insufficient resolution of 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 initialised in hydrostatic equilibrium (HSE; ). Because the LHS of the equation governing HSE,
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,
where and are the hydrostatic density and pressure given by:
Here, and are the initial density and pressure at and is assumed throughout the atmosphere at . 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 , thus modifying the first of equations (3) to:
If 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 . Here, we use .
The atmosphere is initialised with the force-free magnetic field used by Ouyed & Pudritz (1997a):
where is the vector potential, is the disc thickness (set to ), and is the magnetic field strength at , given by:
Here, and (plasma beta at ) are free parameters.
Finally, to ensure the declining density and magnetic field profiles do not fall below observational limits, we add floor values and (c.f. Bergin & Tafalla 2007, Vallée 2003) to equations (4) and (5). By imposing HSE and the adiabatic gas law at , a floor value on imposes effective floor values on and as well.
2.2 Boundary Conditions
In the accretion disc (, ), , the Keplerian speed, and is an “evaporation speed” at the disc surface. The disc and atmosphere are initially in pressure balance with a density contrast , while is initialised using equations (5).
Following Krasnopolsky et al. (1999), , and are held constant, , , (where is the induced electric field), , , , and . Since is sub-slow, these conditions are formally over-determined and should probably be allowed to float. Indeed, we allowed to be determined self-consistently in test simulations, and found only minor quantitative differences in the jet since the pressure gradient is only about 1% of the net Lorentz force at the disc surface. However, allowing to float in the boundary caused undue high temperatures in the disc, and thus small time steps. Therefore, the simulation proceeds more rapidly but otherwise virtually unchanged when is maintained at its initial value.
Inside (), we apply reflecting, conducting boundary conditions (). Thus, , and are reflected across , and magnetic boundary conditions are set according to , , and . At , and are evolved using the full MHD equations.
Finally, we use reflecting boundary conditions along the symmetry axis, and outflow conditions along the outermost and boundaries.
2.3 Scaling Relations
From equation (1) and the adiabatic gas law, one can show:
where is the spherical polar radius. From equations (6), (7), and the ideal gas law (, where is half a proton mass), we derive the following scaling relations to convert from unitless to physical quantities:
for . Note that is the only free parameter varied in this work.
3 Results for
Figure 1 depicts a jet with at yr from the highest resolution grid near the disc surface (bottom panel) to the coarsest grid in which the jet has reached a length of just under 2500 AU (top panel)333Time-lapse animations are available at http://www.ica.smu.ca/zeus3d/rc10/.. A few features worth noting include:
When (angle between and disc surface), Blandford & Payne (1982) show that cold gas near the disc is launched into a collimated outflow. Here, for all , but significant outflow is limited to inside the point where the slow surface intersects the disc ( AU = jet radius at the disc; second panel from top). Below , cold disc material has moved onto the grid and accelerated into the outflow. Above , the weak magnetic field has yet to drive enough disc material onto the grid to displace the hot atmosphere, and outflow is stifled. While gradually increases with time, the majority of mass flux originating from the disc is driven within (0.5 AU; bottom panel).
Jet material becomes super-fast () within a few AU of the disc, 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).
At large distances from the disc ( AU; top panel), the dynamics of the jet become dominated by , and the jet is led by an essentially ballistic, magnetic “nose-cone” with a Mach number of (e.g., Clarke, Norman, & Burns, 1986). Still, is a small fraction () of , consistent with Hartigan et al. (2007).
The knots dominating the bottom panel (c.f., Ouyed & Pudritz 1997b) are produced by the nearly harmonic oscillation of in , whereby fluctuates between and with a period . These oscillations result from the interplay between in-falling 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 AU of the symmetry axis, and gradually merge to form a continuous and narrow column of hot, magnetised material444The knots are resolved by 10–20 zones when they merge, and thus their merger is unlikely related to the ever-decreasing resolution of the nested grids. (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).
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.
4 Comparing simulations and observations
Table 1 summarises 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 , and made numerical measurements of the quantities in Table 1. Variation of other parameters (such as and ) is left to future work.
|proper motion (km s)||100 – 200 (500 max.)|
|rotational velocity (km s)||(5 – 25) 5|
|FWHM jet width (AU)||30 – 80 (at 200 AU)|
|mass-loss rate ()||0.01 – 1|
Note that is the initial value of the plasma beta at , and not the average in the jet. Indeed, Fig. 2a demonstrates that at very early time, , where , and where the average is taken over zones that exceed a certain threshold so that only out-flowing jet material is considered. Thus, the magnetic field within the jet is stronger than would suggest. Initially, is dictated by , but becomes dominated by within yr after launch. As time progresses, gradually increases but never rises above unity (at least for yr), even for . Still, one might speculate from Fig. 2a that with sufficient time, regardless of .
4.1 Proper motion
To understand this result physically, we begin with the magnetic forces:
(e.g., Ferreira 1997; Zanni et al. 2007) where are the gradients parallel and perpendicular to . For a given field line, a stronger at its “footprint” in the disc () generates a stronger which leads to stronger gradients in and thus, from equations (4.1), 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 , where is a weak function of the field strength at the footprint and thus of .
Following Spruit (1996), one can show that as a function of the “fast moment arm” (), the poloidal velocity at the fast point is:
since , the poloidal Alfvén speed at the fast point, is roughly proportional to . is the Keplerian speed at the footprint of the field line. We note that measured values of in our simulations vary as and agree with equation (14) to within 1% so long as the fluid is in approximate steady-state555Indeed, all four steady-state functions from Spruit (1996) remain constant in our simulations to within 5% along steady-state field lines, which we take as validation of our numerical methods..
After the poloidal force given by equations (4.1) decreases to 1% of its maximum value ( a few ), still follows a power law in with index and essentially unchanged from equation (14). Nearer the head of the jet where steady state is no longer valid, we find (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 persists.
Finally, (Fig. 2b and Table 2) is within of near the bow shock and maintains the same power-law dependence on . Thus, these jets are essentially ballistic, where the observed jet speed . In short, all measures of jet speed increase with , a trend that agrees with Anderson et al. (2005) who find for much less evolved jets, .
4.2 Toroidal velocity
Figure 2b and Table 2 show averaged over time and the jet volume for 100 AU as a function of . Like , asymptotes to a constant value. The region inside 100 AU is ignored because the torsion Alfvén wave at low has a non-negligible , is not part of the jet, and skews our results. By fitting a power law to these data, we find .
Unlike , we have not uncovered a rationale for this power law, yet it seems plausible one must exist given the tightness of fit. Eliminating from the power laws for and , we find that . To render this a useful observational tool, further work is needed to quantify the effects of other initial conditions such as and on both the power law index and the proportionality constant, as well as the effect our simplified disc model may have on conditions in the jet at observational length scales.
4.3 Jet radius and mass flux
The jet radius, , is defined by the contact discontinuity (steep temperature gradient in the second panel of Fig. 1) between shocked jet and shocked ambient material, which in turn is determined by where the radial jet ram pressure balances all external forces. Since ram pressure increases with , should increase with , just as observed in Table 2. At any given time, we find that varies with as a reasonable power law though, unlike or , the power index is not constant and decreases slowly in time, while itself increases in time, though at an ever-slowing rate.
The mass flux transported by the jet, , consists of material from both the disc and the atmosphere. Unlike previous simulations where jets are typically evolved long after the leading bow shock has left the grid, no part of any bow shock in our simulations reaches the boundary of the coarsest grid. Thus, each jet continues to entrain material from the atmosphere throughout the simulation at a rate that has a strong dependence on , as seen in Table 2. Indeed we find that varies with as a reasonable power law, with the power index decreasing slowly in time. As the atmosphere is depleted, the mass flux contribution from the disc (which, by design, is independent of ) becomes more important and the dependence of on diminishes.
We have presented the first MHD simulations of protostellar jets that start from a well-resolved launching region ( AU) and continue well into the observational domain (2500 AU). On the AU scale, each jet shows the characteristic and near steady-state knotty behaviour 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 () outflow led by a magnetically confined “nose-cone” (Clarke, Norman, & Burns, 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 (Fig. 2a), 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 to reach their asymptotic limits.
Our jet widths tend to be higher than those observed, particularly when one considers that the values for in Table 2 are at yr666Some simulations had not reached yr at the time of this writing., and that continues to grow in time (e.g., for the jet, AU by 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, Morse, & Raymond 1994), and thus temperatures in excess of 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 K, our mass fluxes drop by a factor of 10–100, in much better agreement with Table 1.
- Anderson et al. (2005) Anderson, J. M., Li, Z.-Y., Krasnopolsky, R., Blandford, R. D., 2005, ApJ, 630, 945.
- Bally, Reipurth, & Davis (2007) Bally, J., Reipurth, B., Davis, C. J., 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, K. Keil (Tucson: Univ. of Arizona Press), 215.
- Bell et al. (1994) Bell, J., Berger, M., Saltzman, J., Welcome, M., 1994, SIAM J. Sci. Comput., 15, 127.
- Berger & Colella (1989) Berger, M. J., Colella, P., 1989, JCoPh, 82, 64.
- Bergin & Tafalla (2007) Bergin, E.A., Tafalla, M., 2007, ARA&A, 45, 339.
- Blandford & Payne (1982) Blandford, R. D., Payne, D. G., 1982, MNRAS, 199, 883.
- Calvet et al. (2000) Calvet, N., Hartmann, L., Strom, S. E., 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, S. S. Russell (Tucson: Univ. of Arizona Press), 377.
- Clarke, Norman, & Burns (1986) Clarke, D. A., Norman, M. L., Burns, J. O., 1986, ApJ, 311, L63.
- Clarke (1996) Clarke, D. A., 1996, ApJ, 457, 291.
- Clarke (2010) Clarke, D. A., 2010, ApJS, 187, 119.
- Coffey et al. (2008) Coffey, D., Bacciotti, F., Podio, L., 2008, ApJ, 689, 1112.
- Dyson & Williams (1997) Dyson, J. E., Williams, D. A., 1997, The Physics of the interstellar medium (2nd ed.; Bristol: IOP Publishing).
- Ferreira (1997) Ferreira, J., 1997, A&A, 319, 340.
- Hartigan, Morse, & Raymond (1994) Hartigan, P., Morse, J. A., Raymond, J., 1994, ApJ, 436, 125.
- Hartigan, Edwards, & Pierson (2004) Hartigan, P., Edwards, S., Pierson, R., 2004, ApJ, 609, 261.
- Hartigan et al. (2007) Hartigan, P., Frank, A., Varniére, P., Blackman, E. G., 2007, ApJ, 661, 910.
- Krasnopolsky et al. (1999) Krasnopolsky, R., Li, Z.-Y., Blandford, R., 1999, ApJ, 526, 631.
- McKee & Ostriker (2007) McKee, C. F., Ostriker, E. C., 2007, ARA&A, 45, 565.
- Meier et al. (1997) Meier, D. L., Edgington, S., Godon, P., Payne, D. G., Lind, K. R., 1997, Nature, 388, 350.
- O’Shea et al. (2004) O’Shea, B. W., Bryan, G., Bordner, J., Norman, M. L., Abel, T., Harkness, R., Kritsuk, A., 2004, eprint (arXiv: astro-ph/0403044).
- Ouyed, Clarke, & Pudritz (2003) Ouyed, R., Clarke, D. A., Pudritz, R. E., 2003, ApJ, 582, 292.
- Ouyed & Pudritz (1997a) Ouyed, R., Pudritz, R. E., 1997a, ApJ, 482, 712.
- Ouyed & Pudritz (1997b) Ouyed, R., Pudritz, R. E., 1997b, ApJ, 484, 794.
- Ouyed & Pudritz (1999) Ouyed, R., Pudritz, R. E., 1999, MNRAS, 309, 233.
- Porth & Fendt (2010) Porth, O., Fendt. C., 2010, ApJ, 709, 1100.
- Pudritz et al. (2007) Pudritz, R. E., Ouyed, R., Fendt, C., Brandenburg, A., 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, K. Keil (Tucson: Univ. of Arizona Press), 277.
- Raga et al. (2002) Raga, A. C., et al., 2002, ApJ, 565, L29.
- Ray et al. (2007) Ray, T., Dougados, C., Bacciotti, F., Eislöffel, J., Chrysostomou, A., 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, K. Keil (Tucson: Univ. of Arizona Press), 231.
- Reipurth & Bally (2001) Reipurth, B., Bally, J., 2001, ARA&A, 39, 403.
- Spruit (1996) Spruit, H.C., 1996, in Evolutionary processes in binary starts, eds. R. A. M. J. Wijers, M. B. Davies, C. A. Tout, (Dordrecht: Kluwer academic publishers), 249.
- Staff et al. (2010) Staff, J. E., Niebergal, B. P., Ouyed, R., Pudritz, R. E., Cai, K., 2010, ApJ, 722, 1325.
- Uchida & Shibata (1985) Uchida, Y., Shibata, K., 1985, PASJ, 37, 515.
- Ustyugova et al. (1995) Ustyugova, G. V., Kolboda, A. V., Romanova, M. M., Chechetkin, V. M., Lovelace, R. V. E., 1995, ApJ, 439, 3.
- Vallée (2003) Vallée, J.P., 2007, NewAR, 47, 85.
- Vitorino et al. (2002) Vitorino, B. F., Jatenco-Pereira, V., Opher, R., 2002, A&A, 384, 329.
- von Rekowski et al. (2003) von Rekowski, B., Brandenburg, A., Dobler, W., Shukurov, A., 2003, A&A, 398, 825.
- Zanni et al. (2007) Zanni, C., Ferrari, A., Rosner, R., Bodo, G., Massaglia, S., 2007, A&A, 469, 811.