© ESO, 2012

1. Introduction

Powerful supersonic bipolar jets are ubiquitous in young accreting stars, suggesting that jets may play a key role in extracting excess angular momentum from the accretion flow. The high ejection efficiency, and the universal jet collimation properties independently of circumstellar envelope density, both strongly favor a magneto-hydrodynamical (MHD) ejection and collimation process (see e.g. Konigl & Pudritz 2000; Shu et al. 2000; Cabrit 2007a). However, the main jet origin is still the subject of intense debate. Proposed options include a self-collimated, steady centrifugally-driven disk wind – either from the inner disk edge near corotation (“X-wind”) or from a more extended range of radii – a pressure-driven stellar wind, and unsteady ejections by reconnexions in the sheared stellar magnetosphere, all three processes probably coexisting to some degree in young stars (see Ferreira et al. 2006, and refs. therein).

Of the wide variety of models published, stationary solutions for “X-winds” and extended MHD disk-winds are of particular interest as they provide synthetic predictions in atomic lines that reproduce several aspects of jet collimation and acceleration in microjets from T Tauri stars (Shang et al. 1998; Shang et al. 2002; Cabrit et al. 1999; Garcia et al. 2001b; Pesenti et al. 2004; Pyo et al. 2003, 2006; Cabrit 2007b). Extended MHD disk winds have received further attention since the tentative detections of jet rotation signatures, that suggest an outer launch radius for the atomic jet¹ in T Tauri stars  ≃0.2−3 AU, well beyond the disk inner edge (Bacciotti et al. 2002; Anderson et al. 2003; Pesenti et al. 2004; Coffey et al. 2004, 2007). Evidence for depletion of refractory species in atomic jets also suggests ejection from relatively large disk radii, outside the dust sublimation radius (Nisini et al. 2005; Podio et al. 2006, 2009; Agra-Amboage et al. 2011). If they remove most of the angular momentum from the accretion flow, such extended MHD disk winds would have a strong impact on the conditions for planet formation and migration in disks (Combet & Ferreira 2008). Therefore it is important to find further observational tests of their presence and radial extent.

An important new test for proposed jet engine(s) is whether they can explain the molecular jet/wind counterparts recently detected within 1000 AU of young stars. Molecular jets were first discovered emanating from the so-called “Class 0” protostars at the earliest stage of star formation, in CO, SiO, H₂, and H₂O masers (Bachiller & Cernicharo 1990; Guilloteau et al. 1992; Gueth & Guilloteau 1999; McCaughrean et al. 1994; Claussen et al. 1998). These narrow molecular jets are clearly distinct from the “standard” swept-up molecular outflows in that they are much faster, with deprojected speeds of 60 − 150 km s^-1 as opposed to  ≤20 km s^-1 in the swept-up gas, much more highly collimated on-axis, and with specific chemical enhancements (Bachiller & Cernicharo 1990; Gueth & Guilloteau 1999; Santiago-García et al. 2009; Tafalla et al. 2010). Their collimation scale, ejection/accretion ratio, and variation of physical parameters with velocity are strikingly similar to those of atomic T Tauri jets, pointing to a similar origin (Cabrit et al. 2007; Lee et al. 2007a; Nisini et al. 2007). Counterparts in H₂, and more rarely in CO, were then identified towards more evolved infrared protostars of Class I (Cernicharo & Reipurth 1996; Davis et al. 2001b, 2002; Davis et al. 2006; Bally et al. 2007) and optically visible T Tauri stars of Class II (Takami et al. 2004, 2007; Herczeg et al. 2006; Pety et al. 2006; Beck et al. 2008). They consist of high-speed H₂ in the jet beam itself, and/or of slow wind/cavities at  ≃30 km s^-1 encompassing the fast ionic jet, on scales of a few 100 AU. The molecular emission tends to peak at lower velocities than atomic lines (Davis et al. 2003).

Three processes are generally invoked to produce molecular jets/winds in young stars: (i) entrainment of ambient molecular material; (ii) gas-phase reformation of H₂ in a dust-free atomic wind (e.g. stellar or magnetospheric); (iii) ejection/reformation of molecules in a dusty wind from the disk surface.

The entrainment hypothesis is intuitively appealing, but appears quite challenging. Turbulent mixing of ambient gas along the walls of a supersonic jet beam predicts too low columns of H₂ (Taylor & Raga 1995). Shock entrainment in narrow bowshock wings along the jet is potentially more efficient, but would require that ambient molecular gas can quickly mix-in and refill the cocoon behind preceding bowshocks (Davis et al. 2001a; Raga & Cabrit 1993). Swept-up molecular cavities driven into ambient gas by an atomic wide-angle wind have also been proposed in Class II sources where molecular counterparts are broader than the atomic jet (see e.g. Takami et al. 2007). However, the nature of such a wide-angle component remains unclear; the fast  ≃100 km s^-1 and powerful wide-angle wind predicted by the “X-wind” model, propagating into a flattened envelope, would sweep up a cavity of width to length  >1 at 200 AU from the source in less than 1000 years (Lee et al. 2001; Shang et al. 2006). This is difficult to reconcile with the moderate width/length  < 1 of molecular “cavities” in Class II sources of age  ≃10⁶ yrs, unless the wide-angle wind is much slower/weaker (e.g. a disk wind launched from large radii). In any case, an ambient swept-up cavity would not easily explain the narrow, high-velocity molecular jets seen in Class 0/I sources. The presence of blueshifted H₂ or CO up to 50 − 100 km s^-1 within 10 AU of Class II/Class I sources (Herczeg et al. 2006; and Herczeg, priv. comm.) also seems difficult to explain with entrainment.

Alternatively, it has been argued that molecular counterparts may trace the primary ejected material itself (Cernicharo & Reipurth 1996; Takami et al. 2004; Pety et al. 2006; Hirano et al. 2006). Studying the formation/survival of H₂ in protostellar winds is a difficult problem as it requires the combination of state of the art MHD wind models with chemical networks and realistic radiation fields and heating/cooling terms, few authors having addressed it self-consistently in the past.

Reformation of H₂ in dust-free stellar winds by ion reactions has been the most extensively studied so far. Spherical models with prescribed temperature laws showed that H₂ forms and survives only at high mass-flux rate and if temperature drops rapidly below 2000 K (Rawlings et al. 1988; Glassgold et al. 1989; Rawlings et al. 1993). Ruden et al. (1990) solved in parallel for the ionization, H₂ abundance, and thermal evolution, including ion-neutral “drag” heating. They found that less dense winds are hotter and more collisionally dissociated, so that the final H₂ abundance exceeds a few percents only at very high wind mass-fluxes  >10^-5   M_⊙ yr^-1. Such a model does not seem able to explain the presence of detectable H₂ in low-luminosity Class 0 − I − II jet sources, where total mass-flux rates are 1−3 orders of magnitude smaller (Lee et al. 2007a; Antoniucci et al. 2008; Hartigan et al. 1995). Far-ultraviolet excess from accretion further reduces the predicted H₂ (and SiO) abundance at low wind mass-flux, even though CO is efficiently formed (Glassgold et al. 1991). Alternatively, Raga et al. (2005) show that H₂ abundances of up to 10% could be reached behind dense internal shocks in a dust-free jet. However, the jet velocity variability amplitude should be less than 15 km s^-1, which is quite restrictive. The effect of magnetic cushioning and of far-ultraviolet (FUV) photons remains to be investigated.

Here we investigate the third possibility, namely whether molecular jet counterparts could be tracing dusty disk winds. We focus on self-similar, centrifugally-driven MHD disk winds. Pioneering calculations by Safier (1993) indicated that H₂ could survive collisional dissociation in the outer regions of such disk winds, beyond a minimum launch radius  ~0.5−3 AU depending on the adopted MHD solution. However, the rate of collisional dissociation may have been underestimated, as the ionization fraction was obtained assuming atomic gas. The faster dissociative recombinations occuring in molecular gas would enhance ion-neutral drag, increasing the wind temperature and the dissociation rate. The calculations also did not consider H₂ destruction by stellar FUV photons, coronal X-rays, or endothermic reactions with O and OH.

In the present paper, we explore this problem one step further by integrating in parallel for the detailed chemistry, ionization, and thermal state along MHD disk wind streamlines. We use a model developed for magnetized molecular interstellar shocks with an extended chemical network of 134 species (Flower & Pineau des Forêts 2003), adding irradiation by (attenuated) coronal X-ray photons and FUV photons from the accretion shock. Another difference with respect to Safier (1993) is the use of a slower and denser MHD accretion-ejection solution that better matches the rotation signatures, poloidal velocities, and ejection-accretion ratio in atomic T Tauri jets (Pesenti et al. 2004; Ferreira et al. 2006). The flow dynamics, thermo-chemistry, and treatment of FUV and X-ray radiation, are described in detail in Sect.  2, which may be skipped in a first reading. Our results are presented in Sect. 3 for a range of accretion rates and stellar masses representative of the Class 0, Class I, and most active Class II phases of solar-type stars. In Sect. 4, the results are compared to static irradiated disk models and to observational trends, and model limitations are discussed. Conclusions are summarized in Sect. 5.

2. Model description

2.1. Dynamical MHD solution

The adopted dynamical MHD model sets the total density, velocity, magnetic field and current in space. A large scale magnetic field of bipolar topology is assumed to thread the disk on a large radial extension. The wind model satisfies the full single-fluid MHD equations and has been obtained by assuming steady-state, axisymmetry, and a self-similar functional form for all physical variables. These models are an extension of the Blandford & Payne (1982) ideal MHD jet models in that they also solve for the full 3D dynamics of the underlying resistive Keplerian accretion disk (where accreting matter needs to cross field lines), including all dynamical terms in the vertical, azimuthal, and radial directions. In particular, the disk compression and (slight) sub-Keplerian rotation induced by the magnetic field, addressed in Shu et al. (2008), are accurately treated. The accretion-ejection transition, wind mass-loading, and slow point crossing, are thus all determined self-consistently, given a vertical profile of temperature and resistivity in the disk (see Ferreira 1997; Casse & Ferreira 2000, for details).

A key free parameter of the models that is most amenable to observational constraints is ξ, the disk ejection efficiency, fixing the accretion rate as a function of the radius as Ṁ_acc ∝ r^ξ. The total mass outflow to disk accretion rate over a radial extension r_in to r_out is then 2Ṁ_j/Ṁ_acc = ξln(r_out/r_in). The parameter ξ is also linked to the jet physics. As long as wind magnetic torques dominate, it is inversely related to the Blandford & Payne magnetic lever arm parameter λ through (Ferreira 1997) (1)where r₀ is the midplane radius of the magnetic surface along which the matter flows and r_A is its cylindrical radius at the poloidal Alfvén point. The maximum asymptotic jet velocity and specific angular momentum along each streamline are given respectively by and (Blandford & Payne 1982), and are thus smaller for decreasing λ and increasing ejection efficiency ξ.

In Garcia et al. (2001a), we calculated the thermal-ionization structure for atomic MHD disk wind solutions with high lever arms (42 ≲ λ ≲ 84), obtained for a vertically isothermal disk profile (Ferreira 1997). Garcia et al. (2001b) presented synthetic forbidden line predictions and showed that the predicted line centroid velocities were on average too high when compared to observations, and that models with lower lever arms should be favored. This conclusion was later confirmed by the small rotation signatures observed in T Tauri atomic jets, which rule out values of λ ≥ 20 (Pesenti et al. 2004; Ferreira et al. 2006).

Therefore, in the present study we use a “slow” disk wind solution with a smaller lever arm parameter λ = 13.8, developed by Casse & Ferreira (2000). This particular value was chosen as best reproducing the possible rotation signatures reported across the outer layers of the DG Tau atomic jet, with an inferred outer launch radius r_out ≃ 3 AU (see Pesenti et al. 2004). The same model can also account for the highest velocity  ≃ 350 km s^-1 detected on the jet axis, if launching occurs down to the corotation radius r_in ≃ 0.07 AU (Agra-Amboage et al. 2011). Indeed, the predicted asymptotic speed in this model varies with anchor radius r₀ as (2)The corresponding ejection efficiency, ξ = 0.04, also appears compatible with the observed ejection to accretion ratio in the DG Tau atomic jet, given current uncertainties (Agra-Amboage et al. 2011). The resulting (r,z) distribution of density and streamlines is shown in Fig. 1. The density contours turn from horizontal to vertical at polar angles θ ≤ 18°, producing an apparent density collimation in good agreement with observed jet widths (Cabrit et al. 1999; Garcia et al. 2001b; Ray et al. 2007). This dense axial beam is surrounded by wider streamlines from larger r₀ (white curves in Fig. 1) that recollimate on larger scale and reach a lower terminal speed (see Eq. (2)). Such a steep transverse velocity decrease from axis to edge is observed across several resolved Class II jets including DG Tau and CW Tau (Lavalley-Fouquet et al. 2000; Bacciotti et al. 2000; Coffey et al. 2007) and seems in good agreement with the chosen MHD disk wind solution (Pesenti et al. 2004; Cabrit 2009).

“Slow” MHD solutions (λ < 20) that become super Alfvénic may be obtained by including a small entropy deposition at the disk upper layers z ≃ h (Casse & Ferreira 2000). The warm mass-loading region is denser and the ejection efficiency ξ is substantially enhanced, lowering the λ value (see Eq. (1)). Such additional heating was assumed to be powered by dissipation of MHD disk turbulence and parametrized as a tiny fraction f of the disk accretion power, with f = 8  ×  10^-4 in the solution chosen here (see Casse & Ferreira 2000, for a discussion of the effect of f on the wind dynamics). This new class of solutions, referred to as “warm” in Pesenti et al. (2004), is still “cold” in a dynamical sense as the initial thermal energy remains negligible compared to gravity, and wind acceleration is still mostly magnetic. Therefore, the wind temperature beyond the slow point does not affect the dynamics and we may recompute it a posteriori from the actual heating/cooling terms, without loss of self-consistency (see discussions in Sect. 4 and Appendix A).

Other input properties of the adopted solution are: a thermal (uncompressed) aspect ratio ε = h/r = C_s/V_K = 0.03 (with C_s the sound speed and V_K the Kepler speed), a turbulent resistivity in the disk midplane ν_m = α_mV_A,0h, with V_A,0 the Alfvén velocity and α_m = 2 (in order to ensure stationarity), and a turbulent viscosity ν_v = ν_m. Crossing of the slow magnetosonic point is obtained for a ratio of magnetic to thermal pressure in the disk midplane μ ≃ 0.33 (leading to quasi-Keplerian rotation; Shu et al. 2008). The computed second Blandford & Payne jet parameter is κ ≃ 0.1, and the field inclination at the disk surface is θ₀ ≃ 57.3°, with similar surface values of B_z, B_φ and B_r.

For a given self-similar MHD solution, all physical quantities at a given polar angle obey specific scalings with the stellar mass M_⋆, the disk accretion rate Ṁ_acc, and the anchoring radius r₀ of the magnetic surface in the midplane, as given in Eq. (9) of Garcia et al. (2001a). In particular, the midplane field in the adopted MHD solution is, (3)which follows the same scaling as the available magnetic field measurements in protostellar disks (Shu et al. 2007). The required field strengths thus appear plausible. In the following, we will vary M_⋆, Ṁ_acc, and r₀ in order to explore their effect on the molecular content and thermal state of the disk wind.

2.2. Thermo-chemical evolution

The thermo-chemical evolution along wind streamlines was implemented by adapting the 2003 version (Flower & Pineau des Forêts 2003) of a code constructed to calculate the steady-state structure of planar molecular MHD multi-fluid shocks in interstellar clouds (Flower et al. 1985). The ion-neutral drift, FUV field, and dust-attenuation were calculated following the methods developed by Garcia et al. (2001a) in the atomic disk wind case. Irradiation by stellar X-rays was also added, following the approach of Shang et al. (2002).

Given the high densities and low ion-neutral drift speeds in the wind (see Sect. 3) the same temperature and velocity is adopted here for all particles. The latter is prescribed by the underlying single-fluid MHD wind solution. We thus integrate numerically the following differential equations on mass density ρ, species number density n(a), and temperature T along the streamline, as a function of altitude z above the disk midplane: Here, v_z and (∇·v) are the bulk vertical flow speed and the (3D) flow divergence interpolated from the MHD wind solution, n_tot is the total number density of particles, T is the temperature, k_B is the Boltzmann constant, S and R are the rates of change in mass and number of particles, respectively, per unit volume, and Γ and Λ are the heating and cooling rates per unit volume. The equations on n(a) apply to each species a as well as to the individual populations of the first 49 levels of H₂ (up to an energy of 20 000 K) which are integrated in parallel with the other variables. The equations on ρ_f apply to each “fluid” f (neutral, positive, negative) and are used mainly for internal checking purposes, as the total mass density ρ(z) is prescribed by the MHD solution.

Cooling and heating mechanisms include

• radiative cooling by H₂ lines excited by collisions with H, H₂, He, and electrons (Le Bourlot et al. 1999);

• radiative cooling by CO, H₂O, and ¹³CO in the Large Velocity Gradient approximation (Neufeld & Kaufman 1993), and by OH and NH₃ in the low-density limit (Flower et al. 1985);

• atomic cooling by fine-structure and metastable lines of C, N, O, S, Si, C⁺, N⁺, O⁺, S⁺, Si⁺ (Flower et al. 2003) and Fe⁺ (Giannini et al. 2004);

• inelastic scattering of electrons on H and H₂ (Aggarwal et al. 1991; Hummer 1963; Rapp & Englander-Golden 1965);

• energy released by collisional ionization and dissociation and exo/endo-thermicity of chemical reactions (Flower et al. 1985);

• energy heat/loss through thermalization with grains (Tielens & Hollenbach 1985);

• ambipolar-diffusion heating by elastic scattering between the neutral fluid and charged ions and grains (Garcia et al. 2001a, see Sect. 2.3);

• ohmic heating arising from the drift between electrons and other fluids (Garcia et al. 2001a, see Sect. 2.3);

• photoelectric effect on grains (Bakes & Tielens 1994, Eq. (42)) irradiated by the (attenuated) FUV field of hot accretion spots (see Sect. 2.4);

• heating through cosmic-rays and (attenuated) coronal stellar X-rays (Dalgarno et al. (1999), see Sect. 2.5).

The reader is referred to the corresponding references for a discussion of the physical context and hypotheses involved in modelling each of the above processes.

2.2.1. Chemical network

The chemical network consists of 134 species, including atoms and molecules (either neutral or singly ionized) as well as their correspondents inside grain refractory cores, and on grain icy mantles. A representative polycyclic aromatic hydrocarbon (PAH) with 54 carbon atoms is also included, with a fractional abundance of 10^-6 per H nucleus. The total elemental abundances and their initial distribution among gas, grain cores, and icy mantles are taken from Tables 1 and 2 of Flower & Pineau des Forêts (2003).

The charge balance of grains and PAHs is treated as in Flower & Pineau des Forêts (2003). The number of grains per H nucleus and the mean grain cross section are determined from the abundance of depleted elements in grain cores, and the adopted grain size distribution (see Sect. 2.4.2). Most of the grains are charged, and assumed well-coupled to the charged fluid.

We consider 1143 reactions including neutral-neutral and ion-neutral reactions, recombination with electrons, charge exchange, cosmic ray induced desorption from grains, sputtering of grain icy mantles, and erosion of charged grain cores by impact of drifting heavy neutral species. Reaction rates between charged and neutral species are enhanced by ion-neutral drift, following the effective temperature prescription of Flower et al. (1985), with the drift speed v_in computed as described in Sect. 2.3.

H₂ collisional dissociation is treated level by level following Le Bourlot et al. (2002). Reformation on grains is computed with the sticking coefficient of Hollenbach & McKee (1979). Three-body reactions in the gas phase are not included, therefore we will limit ourselves to densities  ≤10¹² cm^-3.

Ionisation and dissociation reactions by far-ultraviolet (FUV) photons from accretion hot spots, and by stellar coronal X-rays, are included as described in Sects. 2.4 and 2.5, respectively.

2.2.2. Initial conditions and integration

The integration of the set of equations for temperature, mass and number of particles along the flow is an initial value problem. Thus initialization of temperature and initial populations have to be devised. As in Garcia et al. (2001a), we start all calculations from the wind slow magnetosonic point (located at 1.66h(r) for our adopted MHD solution), and assume that the temperature and species abundances there are at equilibrium with the local radiation field. These equilibrium values are obtained by performing a “steady-state” run over 10⁵ yrs where thermo-chemical equations are solved with the density held fixed. We then integrate the thermo-chemical equations along the flow using the DVODE integrator (Brown et al. 1989), until the recollimation point where the streamline reaches its maximum radial extension (radius of 40r₀ at z = 900r₀ for our chosen wind solution). We checked that the final temperature and H₂ abundance along the streamline do not depend sensitively on the initial equilibrium conditions, as long as the gas is fully molecular initially. This occurs because the ionisation (which controls the ion-neutral drift heating) adjusts rapidly at the dense wind base to the local radiation field, even though it becomes frozen-in at large distances.

2.3. Ambipolar diffusion coefficients

The ambipolar diffusion and ohmic heating terms require a detailed calculation of elastic momentum exchange rates and drift speeds. The ion-neutral drift speed is calculated from the MHD solution using the analytical formula derived in Appendix A of Garcia et al. (2001a), valid for ||v_in|| ≪ ||v||: where J is the electric current density, B is the magnetic field, c is the speed of light, q_e is the electron charge, n_e and n_i are the number densities of electrons and ions, X = ρ_i/ρ_n, and and are the total momentum transfer rate coefficients for ion-neutral and electron-neutral collisions, respectively.

The total ion-neutral momentum transfer rate coefficient is obtained by summing over the main neutrals (H, H₂, He) and over all ions, charged PAHs and charged grains: (9)where n(a) is the number density of species a, μ_ab is the reduced mass of particles a and b, and  ⟨ σv ⟩ _ab is their elastic collision rate coefficient; the latter is evaluated via the recent analytical fits to quantum-mechanical calculations provided in Table 2 of Pinto & Galli (2008a), when available, given as a function of the effective relative speed : For the H₂-H⁺ pair we use the updated fit in Pinto & Galli (2008b). For the rest of the pairs  ⟨ σv ⟩ _ab is taken as the maximum of the polarizability and hard sphere rate coefficients (Garcia et al. 2001a): (12)where α_b is the polarizability of H, He, or H₂ and the hard sphere cross section is taken as (13)Here denoted the average grain cross section over our size distribution. The total momentum transfer rate coefficient for collisions between electrons and neutrals is given by: (14)withthe  ⟨ σv ⟩ _ea also evaluated from the recent analytical fits of Pinto & Galli (2008a). The ion-electron momentum transfer rate entering ohmic heating is calculated with the classical formula of Schunk (1975; see Eq. (A.13) of Garcia et al. 2001a).

Figure 2 plots the momentum transfer rate coefficients (normalized by ) for collisions between the neutral fluid and the different charged fluids as a function of vertical distance above the disk midplane for a representative streamline. It may be seen that collisions with positive ions dominate the coupling, with an extra contribution from charged PAHs near the flow base.

The ambipolar drag heating and ohmic heating are given by (Garcia et al. 2001a): (15)with . Note that Γ_drag is proportional to the square of the Lorentz force per particle, and inversely proportional to the ionisation fraction in the wind.

2.4. FUV field

The FUV radiation will control the ion abundance (and therefore the ion-neutral drift) in the upper wind regions. It will also photodissociate molecules, and provide extra heating by collisions with warm dust and photoelectric effect on grains. Our treatment essentially follows the approach of Garcia et al. (2001a) with several updates, including an accretion spot geometry and a photochemical network.

2.4.1. Unattenuated radiation field

Instead of the optically thick equatorial boundary layer model adopted in Safier (1993) and Garcia et al. (2001a), we adopt here the paradigm of magnetospheric accretion, recently favored in T Tauri stars, where the inner disk is truncated at several stellar radii and accreting material flows along field lines in quasi free-fall, giving rise to shocked “hot spots” on the stellar surface (Bouvier et al. 2007). The accretion geometry in younger protostars is not yet constrained observationally, but their similar X-ray flare properties to T Tauri stars indicate an already well-developed stellar magnetosphere (Imanishi et al. 2003), hence we will adopt the same paradigm for simplicity.

The FUV radiation field thus comes from the central low-mass star of effective temperature T_⋆ = 4000 K, and from hot accretion spots of fixed temperature T_hs = 10   000 K, chosen to match the mean observed colour temperature of the FUV excess in T Tauri stars (Johns-Krull et al. 2000). Both radiation sources are treated as black bodies.

As in Garcia et al. (2001a), we neglect the scattering contribution to the FUV field². The (unattenuated) direct stellar flux at distance R from the star (with dimensions erg cm^-2 s^-1 Hz^-1) is calculated through the relation (16)where B_ν(T) denotes the specific intensity at frequency ν of a black body of temperature T, R_⋆ is the stellar radius, and δ_hs is the fraction of the total stellar surface  covered by hot spots (assumed uniformly distributed). The main difference with the boundary-layer model of Garcia et al. (2001a) is thus a fixed T_hs independent of Ṁ_acc, and a more isotropic radiation flux (due to the lack of projection effects or disk occultation).

In each model, δ_hs is determined by requiring that the hot spots radiate half of the accretion luminosity, i.e. . The value of one half is meant to be illustrative only: the actual fraction of accretion luminosity radiated in the accretion shock could exceed 80% for a large disk truncation radius  > 5 R_⋆, or decrease if a sizeable fraction of accretion energy is tapped to drive a stellar wind. None of these effects being well quantified, especially in Class 0 and Class I sources, we adopt 50% here for illustration and consistency with the boundary-layer model of Garcia et al. (2001a). We set R_⋆ = 3   R_⊙, the radius of a young accreting 0.5 M_⊙ star near the “birthline” with T_⋆ = 4000 K according to the models of Stahler (1988), noting that it remains a good approximation for a very young Class 0 (R_⋆  ≃  2   R_⊙ for a 0.1 M_⊙ protostar accreting at 5  ×  10^-6 M_⊙ yr^-1; cf. Stahler 1988). With the above assumptions, the modelled FUV hot spot continuum in the domain 1400 − 2000 Å is a factor 1.5−2.5 stronger than the accretion shock models of Gullbring et al. (2000) for DR Tau and DG Tau, using the same accretion rates (Ṁ_acc  =  3−5  ×  10^-7   M_⊙ yr^-1). On the other hand, if the FUV spectrum below 1400 Å were as flat as in the lower accretion stars BP Tau and TW Hya (Bergin et al. 2003), our blackbody model would underestimate the flux at 1000 Å by a factor 4−2. This remains sufficiently accurate for the present exploratory study.

Since we do not solve for radiative transfer, we neglect the extra FUV flux in the stellar Lyα line. Lyα pumping of highly excited levels of H₂ is occuring close to the star (Herczeg et al. 2002, 2006; Nomura & Millar 2005) but will be less important on the larger jet scales of interest here, especially since we focus on dense molecular jets where H₂ is well shielded.

2.4.2. Attenuation by dust

The stellar FUV field is attenuated mainly by dust in the disk wind along the line of sight to the star. The sublimation radius R_sub is calculated via Eq. (B.5) of Garcia et al. (2001a), considering a unique sublimation temperature T_sub = 1500 K. We used the dust optical properties tabulated by Draine & Lee (1984); Draine & Malhotra (1993); Laor & Draine (1993), assuming the standard “MRN” mixture of astronomical silicates and graphite and the grain radius distribution n(a) ∝ a^-3.5 proposed by Mathis et al. (1977), with grain radii in the interval [0.005, 0.25] μm. The sublimation surface for the Class 0/I/II models in Sect. 3 is plotted in Fig. 3 on top of inner disk wind streamlines.

The dust optical depth to the star at point of spherical coordinates (R,θ) is where is the dust extinction cross-section per H nucleus (integrated over the grain size distribution), and is the column density of hydrogen nuclei on dusty wind streamlines along the line of sight to the star. In our self-similar wind model, the wind density varies as R^-1.5f(θ), and  may then be integrated analytically as: (17)Here, n_H is the number density of H nuclei in all forms, n_H ≃ n(H) + 2n(H₂) + n(H⁺), and R_min(θ) is the spherical radius inside which there is no dust, i.e. either the sublimation radius, R_sub, or the radius of the innermost flow line at angle θ, R_in(θ), whichever is larger (see Fig. 3). The innermost wind streamline is assumed launched from r_in = 0.07 AU, a typical corotation radius in young low-mass stars.

The total dust geometrical cross section per H nucleus, calculated from our grain size distribution and the abundance of depleted elements in grain cores, is   ×  10^-21   cm², where  is the mean square grain radius; the ratio of A_V to N_H is then   ×  10^-22   cm².

2.4.3. Dust temperature

Collisions with dust grains heated by the UV field can be an important gas heating/cooling term at the wind base (though unimportant further out). We assume that above the wind slow point (z ≥ 1.7   h), dust is optically thin to its own radiation. The dust temperature is then calculated at every step from the grain radiation equilibrium against incident (wind-attenuated) stellar photons: (18)where σ_B is the Stefan-Boltzmann constant, is the mean square grain radius, is the Planck-averaged grain emission efficiency (weighted by n(a)a²), and  is the dust absorption cross section per H nucleus defined earlier. Collisions with the gas are neglected in the grain thermal balance, as they compete with radiation only at densities  >10¹³   cm^-3 (Glassgold et al. 2004). We verified that when the FUV excess and the wind attenuation are negligible, our T_d values agree with detailed calculations at the top of the disk atmosphere by D’Alessio et al. (1999) for MRN dust and a 4000 K star (taking into account differences in R_⋆).

2.4.4. Photochemical network

A network of photoionisation and photodissociation reactions is implemented with rates in the form k_p = χ   γe ^− βA_V, where A_V is the visual extinction calculated above (Sect. 2.4.2), β,γ are constants depending on the reaction, taken from van Dishoeck (1988) and Roberge et al. (1991), and χ is the ratio at 1000 Å of the unattenuated FUV flux to the mean interstellar radiation field of Draine (1978), integrated over 4π (2  ×  10^-6   erg   s^-1   cm^-2   Å^-1). Such a “monochromatic” flux normalization at 1000 Å was chosen rather than the average ratio to the interstellar field over the interval 910 − 2066 Å, denoted G₀ in the literature (see e.g. Tielens & Hollenbach 1985), as photoreactions contributing to the main ionisation and dissociation reactions in our models (C, S, CH⁺, H₂, CO) occur in the narrow range 910 − 1200 Å.

Special cases of photodissociation reactions are the dissociation of H₂ and CO, which occur through line absorption. We thus need to consider not only the shielding caused by the dust particles, but also the additional self- and cross-shielding due to those species themselves. They were implemented in the model in an approximate way as described below.

2.4.5. H₂ photodissociation

We adopted an H₂ photodissociation rate coefficient per molecule of the form proposed by Draine & Bertoldi (1996): (19)where p₀(H₂) = 4.2  ×  10^-11 s^-1 is the unshielded rate in the Draine field, f_H2,d(A_V) = e ^− β_dA_V is the shielding factor by dust, with β_d = σ_d(1000   Å)/σ_d(5500   Å) = 6.3 for our dust mixture, and f_H2,H2(N_H2,b) is the H₂ self-shielding factor given by Eq. (37) in Draine & Bertoldi (1996), with the line width parameter b set here to the local sound speed c_s. The shielding column of H₂ molecules on the line of sight to the star, N_H2, depends on the H₂ abundance on all inner wind streamlines. In order to allow a local calculation in the present exploratory study, we assume a mean H₂ abundance on the line of sight to the star equal to half the local H₂ abundance at the current point of the streamline. This “local” assumption tends to overestimate self-shielding, but will not affect our results when photodissociation is not the dominant destruction mechanism. A conservative check of this hypothesis is made a posteriori (see Sect. 3).

2.4.6. CO photodissociation.

For CO, one needs to consider not only the self-shielding and the shielding by dust, but also the shielding by H₂ (through line overlap with CO). We have fitted the data of Table 11 from Lee et al. (1996) as a function of A_V and of the shielding columns N_CO and N_H2, and obtained the following analytical formulae for the shielding factors by dust, CO, and H₂ respectively: These analytical formulas are valid to within a factor of 2 up to N_CO ≃ 2  ×  10¹⁹ cm^-3, and over the whole range in N_H2 and A_V, as shown in Fig. 4.

The photodissociation rate coefficient per molecule is then (23)where p₀(CO) = 2  ×  10^-10 s^-1 is the unshielded rate for an isotropic Draine field (twice the value from Lee et al. 1996, who assume a one-sided cloud illumination). As done for H₂, the shielding column of CO molecules on the line of sight to the star, N_CO, is calculated assuming a mean CO abundance on inner disk streamlines of half the local value of the CO abundance at the current point.

2.5. X-rays

Hard coronal X-rays will be the dominant ionization process at the base of the dusty disk wind, where the stellar FUV flux is strongly extinguished by dense inner streamlines. Energetic secondary electrons generated by X-ray ionisation will also dissociate and heat the gas. Our adopted X-ray treatment and chemical network are described in the following.

2.5.1. X-ray spectrum and attenuation

The X-ray flux is modelled as a thermal spectrum with characteristic energy kT_X, F_X(E) ≡ dL_X(E)/dE = (L_X/kT_X)  ×  exp( − E/kT_X), The attenuated rate of X-ray energy deposition per H nucleus in the wind writes (24)where E₀ is the low energy cutoff of the incoming X-ray spectrum, σ_pe(E) is the photoelectric cross-section per H nucleus, and N_H(R,θ) is the total column of H nuclei to the star through inner wind streamlines (given by Eq. (17) with R_min(θ) = R_in(θ)). We calculate H_X analytically following Glassgold et al. (1997) and Shang et al. (2002), who adopted a low-energy cutoff E₀ = 0.1 keV and a power-law approximation to the cross-section, σ_pe(E) = 2.27  ×  10^-22   cm^-2   (keV/E)^p with p = 2.485. H_X then writes (25)where ξ₀ = E₀/kT_X, τ_X  ≡ σ_pe(kT_X)   N_H(R,θ) is the X-ray optical depth at kT_X, and J_p(τ_X) is an analytical fit given in Eq. (C3) of Shang et al. (2002), and plotted in Fig. 5. The fit is valid for , a condition verified in our models.

The effect of changing kT_X for a fixed L_X and fixed position (R,θ) (and therefore a fixed attenuating column N_H) may be visualised by rewriting H_X as (28)The function τ_XJ_p(τ_X), also plotted in Fig. 5, is almost flat over a broad range of τ_X (0.001 to 100). The regions where X-ray ionisation dominates in our models have values of τ_X ≃ 0.3 − 30 within this range for kT_X = 4 keV. Furthermore, the ionization fraction varies only with the square root of H_X. Varying kT_X in the range 2.5 − 15 keV (i.e. τ_X in the range 0.01 − 100) would thus not strongly modify our results.

2.5.2. X-ray ionization rates

X-ray ionization rates are calculated with the same simplifying assumptions as in Glassgold et al. (1997) and Shang et al. (2002): the mean energy to form an A ⁺  ion is approximated by its constant high-energy limit W(A ⁺ ), valid for primary photoelectrons with E > E₀  =  0.1 keV. Neglecting the small difference in energy between the absorbed X-ray photon and the primary photoelectron, and recalling that the absorbed energy H_X is per H nucleus, the X-ray ionization rate of species Aper atom is then given by (29)The limiting mean energies W(A ⁺ ) for X-ray ionization of H₂, H, and He in neutral gas mixtures with 10% of helium have been calculated by Dalgarno et al. (1999) who fitted them³ as linear functions of n(H)/n(H₂) (middle terms of Eqs. (30) − (32) below). We further simplify their expressions by taking 0.53 ≃ 0.5, 1.89 ≃ 2, and recalling that the number density of H nuclei n_H ≃ n(H) + 2n(H₂) at low ionization, and n(He) = 0.1n_H, which yields the right-hand side terms: Inserting the right-hand side expressions for W(A ⁺ ) in Eq. (29), the X-ray ionization rates per species A become independent of n_H and n(A) and may be written solely in terms of the parameter  that enters our chemical network: The total ionization rate per H nucleus, often denoted ζ_X in the literature, is related to  by regardless of the ratio H/H₂. The formation of H⁺ through dissociative ionization of H₂ is also included with a rate of /22 per molecule (Dalgarno et al. 1999). A constant cosmic ray ionization rate   ×  10^-17 s^-1 is added to the rate produced by stellar X-rays, but it plays a negligible role in our models.

2.5.3. X-ray induced dissociation

The rate of H₂ dissociation by X-ray induced electrons is taken as 0.5, following the results of Dalgarno et al. (1999) for primary electron energies  > 0.1 keV. The energetic electrons also collisionally excite H₂ to higher electronic levels, which then radiatively decay by emitting a flux of ultraviolet “secondary photons” able to dissociate other molecular species. We adopt the dissociation factors p_M calculated by Gredel et al. (1989) for various species, leading to a rate per unit volume of (Flower et al. 2007), where we adopt a grain albedo ω = 0.5. For CO, we adopt the high-abundance, high-temperature limit of p_CO = 28 for ω = 0.5 (Gredel et al. 1987). Photodetachment of electrons from PAHs and grains by secondary photons is also included (Flower & Pineau des Forêts 2003).

2.5.4. X-ray heating

The heating efficiency of the X-rays and cosmic rays is treated collectively. We use the results of Dalgarno et al. (1999), who studied the heating produced by high energy electrons interacting with a partly ionized gas mixture of H, H₂ and 10% of He. The heating efficiency η_h is calculated for ionization fractions  < 0.1 and includes the fraction of initial primary electron energy lost by elastic scattering with neutrals and electrons and by rotational excitation of H₂ (rapidly thermalized at our densities). The total heating rate per unit volume by X-rays and cosmic rays is given by: (38)where ϵ_H2 = 20 eV follows from the definition of  in Eq. (33), r = n(H₂)/n(H), and η_H2,η_H are respectively the heating efficiencies for the H₂-He and the H-He ionized mixtures, fitted as a function of the electron fraction x_e = n_e/n_H in the form (39)with η₀ = 0.055, c = 2.170, a = 0.366 for H₂-He, and η₀ = 0.117, c = 7.950, a = 0.678 for H-He (Dalgarno et al. 1999, Table 7).

3. Results

To investigate how the chemical content of protostellar disk winds evolves in time with the decline of accretion rate, we considered three sets of parameters representative of the main evolutionary stages of young low-mass stars with bright jets:

• a very young Class 0 protostar (Ṁ_acc = 5  ×  10^-6   M_⊙ yr^-1, M_⋆ = 0.1   M_⊙), where accretion is high and the star has not yet reached its final mass; the chosen parameter values describe e.g. the heavily embedded sub-millimeter exciting source of the HH 212 jet (Lee et al. 2006);

• a Class I source (Ṁ_acc = 10^-6   M_⊙ yr^-1, M_⋆ = 0.5   M_⊙), where the star has accumulated most of its final mass, but residual infall and accretion still proceed at a fast pace; the chosen parameters describe e.g. the infrared exciting source of HH26 (Antoniucci et al. 2008);

• an active Class II star (Ṁ_acc = 10^-7   M_⊙ yr^-1, M_⋆ = 0.5   M_⊙), representative of optically visible T Tauri stars with bright optical microjets, e.g. DG Tau or RW Aur (Gullbring et al. 2000).

With our choice of parameters, and the self-similar scalings of the MHD solution, the wind density drops by a factor 10 from one Class to the next. Our adopted stellar radius of 3 R_⊙ gives the same accretion luminosity of 5   L_⊙ in the Class 0 and Class I models, and 10 times smaller in the Class II, leading to a hot spot coverage fraction δ_hs of 3.2% and 0.32% respectively (cf. Sect. 2.4.1).

We adopt a thermal X-ray spectrum of characteristic energy kT_X = 4 keV and luminosity L_X = 10³⁰ erg/s, typical of the “hard” time-variable component in solar-mass young stars and jet-driving protostars (Imanishi et al. 2003; Güdel et al. 2007). L_X is a time-averaged value of 30% of the median flare luminosity in solar-mass YSOs, based on a typical flare interval of 4 − 6 days (Wolk et al. 2005). The exact value of kT_X has no major effect on our results (see Sect. 2.5.1).

Note that the flow crossing timescale out to the recollimation point is very short, only 50 yrs for a streamline launched at 1 AU in the Class I model. The chemical composition of the disk wind will thus deviate substantially from published “static” disk atmosphere models such as, e.g. Glassgold et al. (2004); Nomura & Millar (2005). The presence of adiabatic cooling and drag heating also introduce major differences in the thermal structure. These effects are discussed in Sect. 4.

We first analyse in detail the results obtained for the 3 evolutionary classes, for a representative anchor⁴ radius r₀ = 1 AU, and for a streamline launched just outside sublimation, r₀ ≃ R_sub. The effect of increasing r₀ is discussed later on, in Sect. 3.2.

3.1. Results for r₀ = 1 AU and r₀ = R_sub for the 3 classes

Figure 6 shows the main physical parameters integrated as a function of z along the 1 AU streamline for the Class 0/I/II models: temperature, ion-neutral drift speed, ionisation fraction, H₂ and CO abundance. Also plotted are the (prescribed) flow density and bulk speed, and the calculated A_V to the central star, attenuated X-ray ionisation rate of H₂, and “effective” FUV field at 1000 Å in Draine units (e^−3A_V being the attenuation factor for Carbon photoionisation in our chemical network). In Fig. 7, we show for comparison the results obtained for streamlines with an anchor radius just beyond the sublimation radius, computed assuming no H₂ or CO self-shielding (i.e. assuming no molecules on streamlines launched inside R_sub).

3.1.1. Radiation field and ionization structure

The ionization fraction is a key parameter governing the temperature profile, through the drag heating term. It is controlled by attenuation of the radiation field through inner wind streamlines. The behavior is qualitatively similar along the streamlines launched from 1 AU and from R_sub.

The A_V towards the star is highest at the wind base, where the line of sight crosses the densest wind layers (see Fig. 1) and drops steadily as material climbs along the streamline (Figs. 6, 7c). This causes a steep rise in the effective attenuated radiation field, until the 1/R² dilution factor takes over (Figs. 6, 7e,f). Combined with the wind density fall-off (Figs. 6, 7b), this generates a global increase in ionisation fraction along the streamlines out to 30 − 100 AU, until recombination sets in (Figs. 6, 7g). In our models, ionization is dominated by stellar X-rays for A_V > 3 mag, and by FUV photons further out.

As the wind density drops by a factor 10 from Class 0 to Class I to Class II (Figs. 6, 7b), the ionisation fraction, X(i⁺) ≡ n(ions)/n_H, increases by a factor 3 to 10. Indeed, the smaller attenuating column through the inner wind (Figs. 6, 7c) leads to higher X-ray and FUV ionization rates. At the same time, the lower wind density reduces recombination. Both effects work in the same direction.

Figure 8 shows the main ionization contributors along the 1 AU streamline for the Class 0, Class I, and Class II models. From these results, it can be seen that the main contributors depend on the overall ionization degree, X(i ⁺ ):

• when X(i⁺) ≤ 10^-6:

  the main contributors to the negative charge are the PAHs.The main positive carriers are molecular ions, with the mostabundant being HSO⁺ and H₃O⁺; and H (a direct product of X-ray ionization) also contributing in the Class II.

• when X(i⁺) > 10^-6:

  the main contributors to the negative charge are the electrons. The main positive carriers are the atomic ions S⁺, C⁺, and H⁺, which recombine more slowly with electrons than molecular ions do.

An interesting result of our calculations is the relatively high abundances of CH⁺, SH⁺, and H⁺ of 10^-6−10^-4 reached along the Class I and Class II streamlines (Fig. 8). CH⁺ and SH⁺ are formed by endothermic reaction of C⁺ and S⁺ with H₂, when the wind temperature exceeds about 3000 K. H⁺ forms by X-ray ionisation and charge exchange, but also by photodissociation of CH⁺, with carbon playing the role of a catalyst through the reaction chain: Through this repeated formation cycle, H₂ is partly converted into H⁺, helping to offset H⁺ destruction by charge exchange. The flow crossing timescale is also too short for significant H⁺ radiative recombination. This example clearly illustrates the out-of-equilibrium and hybrid nature of the chemistry in MHD disk winds, combining elements from both C-shocks (strong heating by ion-neutral drift) and photodissociation regions (abundant C⁺ and S⁺).

3.1.2. Drift speed, heating/cooling processes, and temperature behavior

Thanks to X-ray and UV ionization, we find that the drift speed remains relatively low along the streamlines, at 10% of the bulk poloidal speed or less, as seen Figs. 6, 7d. This validates a posteriori the single fluid approximation made in computing the MHD dynamical solution. It also keeps the drag heating and gas temperature to a moderate value. Note that the ratio of drift speed to poloidal speed scales as , so that the single fluid approximation remains fulfilled in the Class 0 jet (small M_⋆, large Ṁ_acc) despite its low ionization.

The main heating and cooling terms for r₀ = 1 AU are plotted in Fig. 9. Collisions with warm dust can be important near the wind base when densities exceed a few 10⁹ cm^-3, but ambipolar “drag” heating by ion-neutral collisions, Γ_drag, quickly prevails as density drops. The predominance of ambipolar diffusion over other sources of heating (shocks excepted) is a widespread characteristic of steady MHD-driven protostellar winds from low-mass sources, and results from the strong accelerating J  ×  B force combined with a low wind ionisation (see Ruden et al. 1990; Safier 1993; Garcia et al. 2001a). Cooling is dominated by adiabatic expansion out to z  ≃  1 AU, and by radiative cooling further out, due to the high molecular abundances. The main coolant is H₂ in the Class I and Class II models, and H₂O in the Class 0 model.

The resulting behavior of temperature along the 1 AU streamlines is qualitatively similar regardless of the accretion rate (Fig. 6a): after initial cooling due to adiabatic expansion, drag heating takes over and heats up the gas to 1000−5000 K. This leads to a sharp increase in molecular cooling (denoted Λ_rad thereafter) which (assisted by collisional dissociation cooling in the Class I/II models) eventually overcomes Γ_drag, leading to a temperature turnover followed by a slow decline where both terms remain in close balance.

The balance between Γ_drag and Λ_rad on the 1 AU wind streamline is maintained thanks to the powerful thermostatic behavior of molecular cooling, which is highly sensitive to temperature. Any slight excess/deficit of drag heating leads to an increase/decrease in radiative cooling until the two terms match again. In the outer wind, radiative cooling approaches the low-density limit (i.e. ), therefore the asymptotic jet temperature profile is determined by the condition . The latter term slowly declines with distance in our models, leading to the slowly declining temperature. Furthermore, the self-similar scalings for J,B,n_H in Eq. (9) of Garcia et al. (2001a) yield a scaling for , explaining why we find much cooler jets in the Class 0 case; final temperatures are only  ~500−700 K, compared to 2000 K in the Class I, and 3000 K in the Class II models.

The streamlines from r₀ = R_sub in the Class I and II cases show a slightly different asymptotic temperature behavior, with a flatter isothermal “plateau” at large distances (Fig. 7a). Molecules are heavily photodissociated on these two streamlines, so that ambipolar diffusion heating is balanced by adiabatic cooling alone. The flat temperature “plateau” resulting from this balance between Γ_drag and Λ_adia is a well known typical property of atomic self-similar MHD disk winds (Safier 1993; Garcia et al. 2001a). The lower atomic plateau temperature  ≃ 4000 K found here, compared to  ≃10⁴ K in Garcia et al. (2001a), results from our different adopted MHD wind solution, which is typically 2 times slower and 5 times denser (see Sect. 2.1).

It is also noteworthy that, for both values of r₀, the initial temperature rise is much slower in the Class 0 jet; as shown in Garcia et al. (2001a) (Eqs. (24) and (26)), the slope of the initial temperature rise, where Λ_adia dominates the cooling, is determined by the ratio of Γ_drag to Λ_adia, which scales as in self-similar disk winds for a fixed ionisation fraction. This scaling is smaller in our Class 0 than in our Class I case by a factor 5³ = 125, explaining the delay to heat the gas in the Class 0 jet despite a 10 times lower ionization fraction than in the Class I jet.

3.1.3. H₂ survival

The calculated H₂ abundances along streamlines launched from 1 AU and from R_sub are plotted in Figs. 6, 7h, for the 3 classes of sources. Below, we first analyse the results in detail for the Class I case, and then discuss the Class II and Class 0 cases.

Class I case: the streamline launched just outside R_sub = 0.2 AU suffers heavy H₂ destruction (see Fig. 7h). Self-shielding is not operative (we assumed no molecules on inner, dust-free streamlines) and the FUV field is very strong. Hence most H₂ molecules are quickly photodissociated at z ≃ 1 AU. Balance with reformation on dust yields a final H₂ abundance  ≃10^-4.

In contrast, half of the H₂ molecules are found to survive along the 1 AU streamline (see Fig. 6h). This may be understood by examining the timescales of the main formation/destruction processes of H₂ along this streamline, plotted in Fig. 10. The flow timescale up to the current point, t_f, is also plotted for comparison and is seen to be very short (50 yrs out to 1000 AU). With our assumption of a mean H₂ abundance of half the local value, the self-shielded H₂ photodissociation timescale (red curve in Fig. 10) is now much longer than the flow timescale. As a conservative check of this conclusion, we also plot in Fig. 10 (thick pink curve with symbols) the minimum H₂ photodissociation timescale, where N(H₂) and the corresponding self-shielding factor are calculated assuming an H₂ abundance on inner streamlines equal to its minimum value, as obtained on the streamline r₀ = R_sub. Even in this overly pessimistic case, the minimum H₂ photodissociation timescale remains an order of magnitude longer than the flow time at all z (Fig. 10), confirming that H₂ will escape photodissociation on the 1 AU streamline in the Class I model.

Only two processes are found to be fast enough to affect the H₂ abundance on the 1 AU streamline: the first to occur is endothermic neutral-neutral reactions with O and OH to form H₂O and H atoms (denoted collectively as O + 2H₂ = H₂O + 2H in Fig. 10): The second process is collisional dissociation by H and H₂. The timescale shortens by 6 orders of magnitude as temperature increases from 2000 K to 4000 K, and approaches the flow time around z ≃ 2 AU, so that half of H₂ is destroyed by 10 AU. Further out, collisional dissociation shuts off due to both the temperature decline and the density fall-off; the H₂ abundance thus remains “frozen-in” at half its initial value (Fig. 6h). We conclude that streamlines will turn from mostly atomic to mostly molecular around r₀ ≃ 1 AU for our Class I model.

Active Class II case: the H₂ abundance shows a similar behavior to the Class I case. The streamline from R_sub suffers heavy photodissociation, with a final abundance of  ≃ 2  ×  10^-5 (Fig. 7h), while half of the molecules survive on the 1 AU streamline (Fig. 6h). However, the minimum photodissociation timescale on the 1 AU streamline, computed with the same method as above, is now shorter than the flow time around z ≃ 30 AU; therefore, our local approximation to self-shielding might significantly underestimate H₂ photodestruction on the Class II streamline launched at 1 AU. Our H₂ abundance in this particular case should thus be viewed conservatively as an upper limit, pending more detailed non-local calculations.

Class 0 case: here, wind densities are so high that dust screening against photodissociation is very effective, even on the streamline just outside R_sub. In addition, temperature is so low that oxydation reactions and collisional dissociation are slower than the flow crossing time, except briefly on the R_sub streamline. Therefore, the Class 0 disk wind retains 90 − 100% of its full initial H₂ content on dusty streamlines launched from r₀ ≥ R_sub = 0.2 AU.

3.1.4. Main O, C, and S-bearing species along the 1 AU streamline

The fractional abundances of the main O- and C-bearing neutral species with respect to H nuclei are plotted in Fig. 11 along the 1 AU streamline for the 3 classes of sources. S and SO are also plotted, as they were identified recently as tracers of molecular jets in Class 0 sources (Dionatos et al. 2009, 2010; Tafalla et al. 2010).

CO:

we find that no CO destruction occurs on the Class 0 streamline from 1 AU, thanks to the very efficient dust shielding provided by dense inner wind layers. In contrast, CO is partly destroyed by FUV photodissociation along the Class I/II streamlines. Once the FUV field drops down, CO reforms mainly through the reaction HCO ⁺  + e ⁻  → CO + H but the timescale is longer than the flow time and the final abundance per H nucleus is only 5  ×  10^-5 in the Class I jet, and 10^-5 in the Class II jet. These latter values should be considered conservatively as upper limits, as CO self-shielding could be significantly overestimated by our approximate “local” treatment (see Sect. 3.2 below). We thus predict an abundance ratio of H₂/CO  ≥  0.5  ×  10⁴ in all cases.

H₂O:

when the wind temperature exceeds  ≃ 400 K, oxygen not locked in CO is converted to H₂O through the endothermic reactions of O and OH with H₂ in Eqs. (43), (44). Together with H₃O⁺ recombination, this reaction is efficient enough to balance the water destruction processes active in Class I/II jets (FUV photodissociation, photoionization, reaction with H). The asymptotic abundance of water is thus similar and quite high on the 1 AU streamline for all 3 classes⁵, at  ≃ 5−10  ×  10^-5.

C, O, S atoms, OH and SO:

OH and atomic C and O have low abundances along the Class 0 streamline, where CO and H₂O are well shielded. In contrast, high abundances of these species  ≃ 0.5−1  ×  10^-4 are predicted in the outer regions of the Class I and Class II streamlines, following CO and H₂O partial photodissociation (see Fig. 11). Concerning sulfur, it is mainly in the form of atomic S beyond 100 AU along the Class 0 streamline, and is mainly photoionized into S⁺ in the Class I/II cases, where the FUV field is more intense (see Fig. 8). The asymptotic abundance of SO is substantial in the Class 0 only, at a level of 3  ×  10^-7. These predicted characteristics for Class 0 jets are in line with the relative abundance of SO to CO  ≃  2  ×  10^-3 reported by Tafalla et al. (2010), and with the mass-flux traced by atomic sulfur lines being possibly as large as that inferred from CO (Dionatos et al. 2010).

3.2. Effect of increasing launch radius in the Class I model

In Fig. 12 we illustrate in the Class I case the effect of the launch radius on the wind chemistry and temperature. We plot a range of r₀ of 1−9 AU corresponding to a range in final flow speed of 90 km s^-1 to 30 km s^-1. The streamline from R_sub = 0.2 AU is also shown for comparison. The main effect is that the flow has an “onion-like” thermal-chemical structure, with streamlines launched from larger radii having higher H₂ abundance and lower temperature and ionization. The behavior with r₀ is discussed in more detail below.

We do not present results for the interval 0.2 < r₀ < 1 AU in the Class I jet, as the minimum H₂ photodissociation timescale there becomes similar to or less than the flow time (see discussion in Sect. 3.1.3). Accurate H₂ abundances on these inner streamlines await a detailed non-local treatment of H₂ shielding, not taken into account in the present preliminary approach. For the same reason, we do not present results at r₀ > 1 AU for the Class II jet, but we note that the effect of increasing r₀ will be qualitatively the same as in the Class I case (i.e. higher H₂ abundance and lower temperature). We also do not present results at r₀ > 1 AU in the Class 0 case, since most H₂ and CO molecules survive already at a launch radius just beyond the sublimation radius R_sub = 0.2 AU (see Sect. 3.1.3).

3.2.1. Ionization and drift speed at various r₀

Wind streamlines are initially less ionized at larger r₀ (Figs. 12f,g), due to the larger attenuation of X-rays and FUV photons by intervening streamlines (Fig. 12c). At large distances, they recombine to an asymptotic fractional ionization of a few 10^-5, dominated by S⁺.

The drift speed does not exceed about 10 km s^-1 and remains less than 30% of the bulk flow speed out to r₀ = 9 AU (Fig. 12d). Therefore, magnetocentrifugal forces appear able to lift molecules out to disk radii of 9 AU for the typical parameters of Class I sources, in the MHD solution investigated here.

3.2.2. Temperature and molecule survival at various r₀

The gas temperature at larger launch radii follows a similar behavior to that found at r₀ = 1 AU (i.e. a gradual rise on a scale z ≃ r₀, followed by a shallow decrease), but the temperature peak is shifted to larger z (i.e. lower n_H) and reaches a smaller peak value (see Fig. 12a). Therefore, H₂ chemical oxydation and collisional dissociation are less efficient than for r₀ = 1 AU. As a result, the H₂ abundance increases with r₀ and more than 90% of the molecules survive for launch radii  ≥ 3 AU (Fig. 12h). Photodissociation of H₂, which was already too slow for r₀ = 1 AU, is further reduced by the additional screening from intervening wind streamlines and plays no role.

In contrast, we find that photodissociation remains the major destruction process for CO out to r₀  =  9 AU. The drop in CO abundance towards the inner r₀ = R_sub streamline largely exceeds the factor 2 assumed in our “local” computation of the CO self-screening column (see Fig. 12i). The plotted CO abundances for r₀  >  R_sub should thus be considered as upper limits only, and a detailed non-local treatment of self-shielding will be needed to obtain more accurate values.

4. Discussion

4.1. Comparison with static disk atmosphere models

It is instructive to compare our results with thermo-chemical calculations of irradiated static molecular disks at a radius of 1 AU. Glassgold et al. (2004) focussed on X-ray disk irradiation, with a similar X-ray flux to ours (L_X = 2  ×  10³⁰ erg s^-1, kT_X =  1 keV). Nomura & Millar (2005) studied FUV irradiation by the star with standard ISM dust properties as adopted here, and a mean unshielded flux at a distance of 140 AU of G₀ ≃ 1000 times the mean interstellar flux, averaged over the interval 910 − 2000 Å (see their Fig. 4), similar to our Class I case (see Fig. 6f).

We find that the thermo-chemical structure of MHD disk winds launched from 1 AU in Class I sources markedly differs from irradiated static disks at the same radius: H₂ survives to much greater altitudes and the gas temperature rises more gradually. We briefly discuss and explain these differences below.

4.1.1. H₂ survival and wind shielding

In static disk models at 1 AU, H₂ is destroyed by X-rays or FUV photons above z ≃ 0.1 AU (Glassgold et al. 2004; Nomura & Millar 2005). In contrast, on the 1 AU streamline of the Class I disk wind, half of the gas remains molecular beyond z  ≃  10 AU and until the end of the streamline (Fig. 6h). The survival of H₂ against photodissociation arises from several effects:

• (i)

  The first key effect is the enhanced screening provided by dense inner wind streamlines compared to a hydrostatic flared disk geometry. As shown in Fig. 1, the self-similar MHD disk wind solution chosen here produces density contours that are horizontal just above the disk surface. In contrast, density contours in a hydrostatic disk are strongly flared out, allowing a more direct illumination of the disk surface by the central star (see Fig. 1 of Nomura & Millar 2005). This difference is illustrated by the large A_V values towards the central star in Fig. 6c; in the Class I streamline at 1 AU, the “effective” dust-attenuated FUV field at z ≃ 0.3 AU is χexp ^− 3A_V ≃ 2 compared to G₀ ≃ 10⁷ in the static disk of Nomura & Millar (2005) at the same (r,z). The peak X-ray ionisation rate is   ×  10^-11 s^-1, instead of 6  ×  10^-9 s^-1 at the top of the disk atmosphere in Glassgold et al. (2004).

• (ii)

  The second key factor is the very short flow timescale; only 50 years from the slow point out to 1000 AU for the 1 AU streamline. Static disk models assume chemical equilibrium, so that H₂ disappears wherever destruction reactions are faster than reformation. In disk winds, H₂ can still survive this situation as long as the destruction time is longer than the flow time (see Fig. 10).

4.1.2. Gas temperature profile

In static disk models at 1 AU, FUV photoelectric heating or X-ray heating induce a sharp increase in gas temperature up to 1500 − 5000 K within a few scale heights, at z/r₀ ≃ 0.15 − 0.2 (Nomura & Millar 2005; Glassgold et al. 2004). In our disk wind model, the rise in temperature is much more gradual, with temperatures of a few thousand K reached only around z/r₀ ≃ 1. This difference stems both from the greater attenuation of Xrays and FUV photons by inner wind streamlines, and from the powerful adiabatic expansion cooling experienced by the wind after launching. The latter term largely exceeds FUV and X-ray heating and is only slightly offset by drag heating (see Fig. 9), yielding a slow net heating rate. A related consequence is that temperatures above 2000 K are reached only in wind regions of low density so that H₂ collisional dissociation and chemical oxydation is limited.

4.2. Limitations of the model and planned improvements

The results presented here are only a preliminary attempt at addressing the highly complex problem of the coupled ionisation, chemical, thermal, and dynamical structure of MHD disk winds from young stars. Several limitations have been introduced that we discuss below with our planned improvements:

• (i)

  We ignored feed-back of the thermal structure on thedynamics. This is in fact well justified because thermal pressuregradients are already negligible above the slow point where ourintegration starts, and are only dynamically important in themass-loading region located deeper in the disk, aroundz ≃ h (see Casse & Ferreira 2000). Therefore, any deviations between our computed wind temperature and the prescription adopted by Casse & Ferreira (2000) to compute the dynamical solution do not alter the self-consistency of the model. Nevertheless, we show in Appendix A that the two temperatures fall close to each other near the wind base. The existence of such “slow” warm MHD solutions might then be a natural outcome of the heating mechanisms present near the surface of an irradiated molecular disk. Thermal calculations extending into the optically thick and resistive mass-loading region are under way to test this conjecture (Garcia et al., in prep.).

• (ii)

  We also ignored feed-back of the ion-neutral drift on the flow dynamics; we showed that this single fluid approximation is well fulfilled at launch radii  ≤ 9 AU but tends to degrade at larger r₀ or small Ṁ_acc, due to the scaling of v_in/V_p ∝ X(i ⁺ )^-1Ṁ_acc^-1. Therefore, ion-neutral drift may be one of the factors (together with the available magnetic flux) limiting the outer radial extension of molecular MHD disk winds in young stars. The drift will also decrease rotation signatures in neutral jet tracers at large r₀ compared to the single fluid prediction. Multi-fluid MHD calculations would be useful to further quantify this effect.

• (iii)

  We adopted a simplified treatment of shielding against photodissociation. An upcoming improvement of our model will be the inclusion of a non-local treatment of the self-screening, in order to obtain more accurate predictions of CO and/or H₂ abundances in the Class I and Class II jets. Shielding against photodissociation also depends on the assumed composition of inner wind streamlines launched between r_in = 0.07 AU and R_sub. As done in Ruden et al. (1990) and Garcia et al. (2001a), we assumed that all wind material with R(θ) > R_sub contributes to the A_V, i.e. that dust reforms efficiently on these dense streamlines after they exit the sublimation surface. This may be too optimistic. On the other hand, when computing the minimum H₂ abundance at r₀ = R_sub, we assumed that no molecules are present on more inner streamlines to provide self-shielding, which is pessimistic. Indeed, CO overtone infrared line profiles in Class II disks indicate the presence of rotating CO down to r_in ≃ 0.04 AU (Najita et al. 2003). Another planned improvement in our model will thus be to perform thermal-chemical modelling of dust-free streamlines launched inside R_sub, including H₂ gas phase formation processes through ion chemistry and three-body reactions, and advection from outer radii.

• (iv)

  The dust size distribution and dust to gas ratio were assumed to be the same as in the ISM. This is probably adequate for the young Class 0 and Class I, where the disk is young and still being fed by an infalling envelope, but there is evidence for both dust growth and dust settling starting to occur during the Class II phase, where the envelope has dissipated. This will decrease the efficiency of dust screening against photodissociation. On the other hand, grain growth may lead to a smaller dust sublimation radius, so that more molecular streamlines can contribute to self-shielding. Knowledge of the dust size stratification in turbulent magnetized disks will be needed to assess the net effect on H₂ survival in MHD disk winds from evolved Class II stars.

• (v)

  The present exploratory work considered only one particular steady MHD disk wind solution, with a magnetic lever arm parameter λ ≃ 14 best reproducing current tentative rotation signatures in atomic jets from Class II (see Sect. 2.1). If jet rotation has been overestimated, slower models with smaller values of λ would have to be considered. These denser winds (higher ejection efficiency ξ) would be better self-screened, and thus less subject to photodissociation. On the other hand, the effect on H₂ collisional dissociation is difficult to predict without detailed temperature calculations: the smaller accelerating Lorentz force per particle will tend to decrease ambipolar heating, while the lower ionization (due to enhanced screening) will tend to increase it (Eq. (15)). Such solutions are being developed and will be explored in forthcoming papers.

4.3. Predicted observational trends and future tests

Our preliminary results on the thermo-chemistry of a dusty centrifugal MHD disk wind exhibit general properties that seem promising to explain several observed trends in the molecular counterparts of stellar jets. They are summarized below, and future observational tests are outlined.

Concerning low-mass Class 0 sources, the model predicts that dusty flow streamlines will keep most of their H₂ and CO content at least down to the sublimation radius of  ≃0.2 AU, i.e. up to flow speeds of  ≃100 km s^-1 for M_⋆ = 0.1   M_⊙. This seems promising to explain the frequent detection of collimated fast H₂ and CO in Class 0 jets, as well as the high ratio of H₂ to CO  ≃  10⁴ recently estimated in one of them, HH211 (Dionatos et al. 2010). Furthermore, tentative rotation signatures reported so far in Class 0 jets (Lee et al. 2007b, 2008) seem consistent with steady MHD centrifugal launching from the expected range of disk radii r₀  ≃  0.15−0.6 AU, when the most likely poloidal speeds are adopted⁶. Finally, if the MHD disk wind operates over a radial extent from r_in = 0.2 AU to r_out of 1−9 AU, the ratio of molecular mass-flux to accretion rate would be ξlog (r_out/r_in)  ≃  0.06−0.15, in the same range as observations (Lee et al. 2007a).

As the accretion rate drops and the stellar mass increases, our modelling indicates that the MHD disk wind is more irradiated (less screening against FUV radiation) and hotter (self-similar scaling of Γ_drag with M_⋆ and Ṁ_acc). Therefore, the molecular H₂ zone moves to larger launch radii, around r₀  ≥  1 AU in our Class I source, with a typical terminal speed of  ≃100 km s^-1 (cf. Eq. (2)). The molecular region will move even further out in Class II sources. This predicted trend is in line with possible rotation signatures reported so far in molecular Class I/II jets (HH26 and CB26), which do appear to suggest larger magneto-centrifugal launch radii  ≃ 3 − 10 AU than in Class 0 jets, although statistics are admittedly still limited (Chrysostomou et al. 2008; Launhardt et al. 2009; Cabrit 2009). This prediction also seems promising to explain the trend towards slow, wide molecular counterparts in older Class II jet sources. The striking example of the broad slow H₂ wind encompassing the HL Tau and DG Tau atomic jets, and the slow CO conical flow surrounding the atomic jet in HH 30, both suggest a “hollow” molecular wind structure consistent with this idea (Takami et al. 2007; Beck et al. 2008; Pety et al. 2006).

A third promising trend is that gas temperature on molecular streamlines is predicted to increase from T  ≃  700 K in Class 0 jets to T  ≃  2000−3000 K in Class I/II jets. This is in good agreement with temperatures  ≃300−1000 K recently inferred in Class 0 jets from H₂ pure rotational lines studied with Spitzer (Dionatos et al. 2009, 2010), and  ≃2500 K in Class I/II jets from H₂ ro-vibrational lines (Takami et al. 2006; Beck et al. 2008). As rovibrational lines might have a greater contribution from shock-heated gas, H₂ temperature estimates in Class II jets from pure rotational lines would be useful to confirm this result.

While these trends are encouraging, more specific tests need to be carried out to validate the applicability of centrifugal MHD disk winds to molecular protostellar jets. Closer comparison to observations await more detailed calculations of self-shielding for the Class I/II, and computations of synthetic maps and spectra in various tracers, taking into account beam smearing and NLTE effects (Yvart et al., in prep.). Several aspects promise to be particularly discriminant:

A first important test would be a confirmation of rotation in molecular jets, with a resolved spatial pattern consistent with centrifugal acceleration from the disk. Current measurements are affected by beam smearing, which can significantly lower the observed rotation signature compared to the actual underlying rotation curve (Pesenti et al. 2004). This effect could be important in the HH212 and HH211 Class 0 jets in Orion-Perseus, whose widths  ≤0.2′′ are smaller than the best resolution of current submm interferometers (Cabrit et al. 2007; Lee et al. 2007b). The lack of resolution also makes it more difficult to recognize possible contamination of rotation signatures by jet precession, orbital motions, or shock asymmetries (see e.g. Cerqueira et al. 2006; Correia et al. 2009, for discussions of possible artefacts). The ALMA interferometer will be essential to progress on this issue.

A second prediction of the MHD disk wind model considered here is that, beyond the dust sublimation radius, the wind base is substantially shielded from stellar photons and the gas heats up only gradually through ambipolar diffusion (see Sect. 3.1.2). CO rovibrational lines, excited around 1000 K, would thus be formed higher up in the atmosphere (at z  ≃  1 AU) than in static disks heated mostly by irradiation. Comparison of predicted CO rovibrational line profiles with observed ones in Class I/II sources in the near-infrared (Najita et al. 2003; Bast et al. 2011) may thus provide an interesting test. The high predicted abundance of H₂O  >  10^-5 in warm MHD disk winds ejected from  ≃ 1 AU, at all evolutionary stages (see Fig. 11) is also an important characteristic that may be testable with Herschel or infrared observations (e.g., Kristensen et al. 2010; Pontoppidan et al. 2010).

Finally, we note that the wind thermo-chemical properties may be locally modified by internal shock waves forming in the jet due to time-variability or instabilities; which molecules are destroyed or reformed will depend on whether the shock is of C (“continuous”) or J (“jump”) type, which in turn depends crucially on the ionization fraction, magnetic field intensity, and H₂ fraction in the preshock gas (Le Bourlot et al. 2002; Flower & Pineau des Forêts 2003). The present calculations provide for the first time the appropriate shock “initial conditions” in an MHD disk wind, setting the stage to calculate the predicted thermo-chemistry in internal shocks, and to compare it with shock observations. Although such developments are beyond the scope of the present paper, we note that SiO and CH₃OH, observed to be enhanced by 2 to 4 orders of magnitude in Class 0 jet knots (Bachiller et al. 1991; Cabrit et al. 2007; Tafalla et al. 2010), could in principle be released in a dusty disk wind following shock erosion or vaporization of grains (Gusdorf et al. 2008; Guillet et al. 2009, 2011).

5. Conclusions

We have investigated the non-equilibrium thermal-chemical structure of dusty streamlines in a “slow” self-similar MHD disk wind compatible with current observational constraints in atomic T Tauri microjets. We considered a range of accretion rates and stellar masses representative of the Class 0, Class I, and early Class II stages of low-mass star formation, and probing a decrease of 2 orders of magnitude in wind density. A detailed chemical network, a complete set of heating/cooling terms, and the effect of irradiation by stellar X-rays and FUV photons, were considered. Our main conclusions are the following:

• The MHD disk wind has an “onion-like” thermal-chemical structure, with temperature, ionization, and radiation field all decreasing as the launch radius of the streamline increases. The wind is sufficiently ionized to accelerate neutrals out to disk radii of  ≃9 AU, but is sufficiently self-screened and cool to remain molecular on streamlines launched beyond some minimum radius r₀. For the MHD solution explored here, this radius is r₀  ≃  0.2 AU (sublimation radius) for Class 0 parameters, r₀  ≃  1 AU for the Class I, and r₀  >  1 AU for the Class II.

• Key elements for the survival of H₂ in MHD disk winds, as opposed to static irradiated disk atmospheres, are: the short flow crossing timescales (50−100 yrs on the 1 AU streamline), the efficient shielding provided by inner wind streamlines against stellar FUV and X-ray photons, and the strong adiabatic cooling which delays gas heating and limits collisional and chemical destruction.

• Balance between ion-neutral drag and molecular line cooling establishes asymptotic temperatures  ≃ 700−2000−3000 K on the 1 AU streamline, increasing with protostellar stage. Temperatures at the wind slow magnetosonic point are within a factor 2 of the surface heating prescribed in the MHD solution, suggesting that the assumed high wind mass-loading might be naturally achievable (although thermal calculations deeper into the disk are needed to test this conjecture).

• The chemistry of dusty irradiated MHD disk winds is of a complex and hybrid nature, combining the out-of-equilibrium “warm” chemistry of C-shocks heated by ambipolar diffusion, and the photoprocesses present in PDRs and XDRs. Illustrations are the efficient gas-phase synthesis of H₂O at all stages, and the substantial abundances reached by SH⁺, CH⁺, and H⁺ in the hotter and more ionized Class I/II. In these latter cases, OH, O, and C also reach levels similar to H₂O while CO is strongly photodissociated.

• The model predictions appear promising to explain several observational trends in molecular wind counterparts, namely: the presence in low-mass Class 0 sources of highly collimated molecular jets with typical speed up to 100 km s^-1, low temperature  ≃700 K, and estimated magneto-centrifugal launch radii  ≃0.2−0.6 AU; as well as the trends for higher H₂ temperatures  ≃2500 K, larger magneto-centrifugal launch radii, and lower collimation and speed as the source evolves from Class 0 to Class I to Class II.

• More discriminant tests of the applicability of MHD disk winds to molecular jets/winds in young stars require synthetic maps and line profiles in H₂, CO, and H₂O, including a non-local treatment of self-shielding and proper account for beam dilution and NLTE excitation. Such calculations will be presented in following papers in this series, and compared with infrared and submm observations. The effect of internal shocks, in particular the enhancement of species believed to be released from grains, such as SiO and CH₃OH, will also be explored.

---

Acknowledgments

We are grateful to G. Herczeg and P. Lesaffre for useful suggestions and discussions, and to the anonymous referee for many insightful comments that helped to improve the paper. D. Panoglou, S. Cabrit, and P. Garcia wish to acknowledge financial and travel support through the Marie Curie Research Training Network JETSET (Jet Simulations, Experiments and Theory) under contract MRTN-CT-2004-005592. This research has made use of NASA’s Astrophysics Data System.

References

Appendix A: Comparison with the assumed heating in the dynamical solution

The MHD solution adopted in the present work has been obtained with a prescribed self-similar enthalpy injection function peaking near the disk surface, representing a tiny fraction f = 8  ×  10^-4 of the accretion power (Casse & Ferreira 2000). This small surface heating is sufficient to obtain denser, slower disk winds with mass fluxes and poloidal and rotation speeds consistent with observations (see Sect. 2.1), but the nature of the heating process is not yet identified.

In this context, it is interesting to compare the thermal pressure assumed in the underlying MHD accretion-ejection solution with that obtained from our thermal-chemical calculations where the various known heating, cooling and chemical processes are taken into account. Since the gas density ρ is prescribed by the solution, the relevant quantity to compare thermal pressures is just the ratio T/μ (with μ the mean molecular weight). This comparison is shown in Fig. A.1 along a 1 AU streamline for our Class 0, Class I, and Class II models.

The value of T_dyn/μ_dyn in the dynamical solution first undergoes an increase from the disk surface at z = h = 0.03   r₀ to about z_t = 10   h and then an adiabatic decrease. Although the initial rise is not exactly followed by our calculations, it is interesting that the difference between the two curves is less than a factor 2 near the slow magnetosonic point (our starting integration point for the thermo-chemistry). This suggests that the warm gas surface temperatures required to load dense, slow MHD disk winds compatible with observational constraints might be naturally achievable under a wide range of accretion rates, especially as MHD “turbulent/wave” heating will take place in the deeper resistive disk layers, and wind mass-loading could be further enhanced via upward motions induced by MRI turbulence (Suzuki & Inutsuka 2009).

The discrepancy between computed and assumed T/μ can exceed a factor 10 at z_t = 10   h, and a factor 100 beyond. This discrepancy has however no relevant impact on the dynamical behavior of the solution. Indeed, at z_t, the ratio of plasma pressure to magnetic pressure in the dynamical solution is already down to only a few 10^-3, so even a large discrepancy in thermal pressure will not affect the wind dynamics. Although the thermal pressure gradient did play an important role in the mass loading (done within the resistive MHD disk layers), it plays no role in the jet acceleration and collimation itself which is mainly due to magnetic forces. Thermal pressure gradients will eventually start affecting the jet dynamics after the jet refocusses towards the axis. For this reason we have not extended our thermal-chemical calculations beyond the recollimation point of the MHD solution.

All Figures