A&A 377, 589-608 (2001)
DOI: 10.1051/0004-6361:20011145

Atomic T Tauri disk winds heated by ambipolar diffusion

I. Thermal structure

P. J. V. Garcia^1,2, - J. Ferreira³ - S. Cabrit⁴ - L. Binette⁵

1 - Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal
2 - CRAL/Observatoire de Lyon, CNRS UMR 5574, 9 avenue Charles André, 69561 St. Genis-Laval Cedex, France
3 - Laboratoire d'Astrophysique de l'Observatoire de Grenoble, BP 53, 38041 Grenoble Cedex, France
4 - Observatoire de Paris, DEMIRM, UMR 8540 du CNRS, 61 avenue de l'Observatoire, 75014 Paris, France
5 - Instituto de Astronomía, UNAM, Ap. 70-264, 04510 D. F., México

Received 21 March 2001 / Accepted 3 August 2001

Abstract
Motivated by recent subarcsecond resolution observations of jets from T Tauri stars, we extend the work of Safier (1993a,b) by computing the thermal and ionization structure of self-similar, magnetically-driven, atomic disk winds heated by ambipolar diffusion. Improvements over his work include: (1) new magnetized cold jet solutions consistent with the underlying accretion disk (Ferreira 1997); (2) a more accurate treatment of ionization and ion-neutral momentum exchange rates; and (3) predictions for spatially resolved forbidden line emission (maps, long-slit spectra, and line ratios), presented in a companion paper, Garcia et al. (2001b). As in Safier (1993a), we obtain jets with a temperature plateau around 10⁴ K, but ionization fractions are revised downward by a factor of 10-100. This is due to previous omission of thermal speeds in ion-neutral momentum-exchange rates and to different jet solutions. The physical origin of the hot temperature plateau is outlined. In particular we present three analytical criteria for the presence of a hot plateau, applicable to any given MHD wind solution where ambipolar diffusion and adiabatic expansion are the dominant heating and cooling terms. We finally show that, for solutions favored by observations, the jet thermal structure remains consistent with the usual approximations used for MHD jet calculations (thermalized, perfectly conducting, single hydromagnetic cold fluid calculations).

Key words: ISM: jets and outflows - stars: pre-main sequence - MHD - line: profiles - accretion disks

1 Introduction

Progresses in long slit differential astrometry techniques and high angular resolution imaging from Adaptive Optics and the Hubble Space Telescope have shown that the high velocity forbidden emission observed in Classical T Tauri Stars (CTTSs) is related to collimated (micro-)jets (e.g. Solf 1989; Solf & Böhm 1993; Ray et al. 1996; Hirth et al. 1997; Lavalley-Fouquet et al. 2000; Dougados et al. 2000; Bacciotti et al. 2000). Although outflow activity is known to decrease with age (Bontemps et al. 1996), CTTSs still harbor considerable activity (e.g. Mundt & Eislöffel 1998) and present the advantage of not being embedded. It is now commonly believed that such jets are magnetically self-confined, by a "hoop stress'' due to a non-vanishing poloidal current (Chan & Henriksen 1980; Heyvaerts & Norman 1989). The main reason lies in the need to produce highly supersonic unidirectional flows. Indeed, this requires an acceleration process that is closely related to the confining mechanism. The most promising models of jet production rely therefore on the presence of large scale magnetic fields, extracting energy and mass from a rotating object. However, we still do not know precisely what the jet driving sources are. Moreover, observed jets harbor time-dependent features, with time-scales ranging from tens to thousands of years. Such time-scales are much longer than those involving the protostar or the inner accretion disk. Therefore, although the possibility remains that jets have a non-stationary origin (e.g. Ouyed & Pudritz 1997; Goodson et al. 1999), only steady-state models will be addressed here.

Stationary stellar wind models have been developed (e.g. Sauty & Tsinganos 1994), however observed correlations between signatures of accretion and ejection clearly show that the disk is an essential ingredient in jet formation (Cohen et al. 1989; Cabrit et al. 1990; Hartigan et al. 1995). Therefore we expect accretion and ejection to be interdependent, through the action of magnetic fields. There are mainly two classes of stationary magnetized disk wind models, depending on the radial extent of the wind-producing region in the disk. In the first class (usually referred to as "disk winds''), a large scale magnetic field threads the disk on a large region (Blandford & Payne 1982; Wardle & Koenigl 1993; Ferreira & Pelletier 1993,1995; Li 1995; Li 1996; Ferreira 1997; Krasnopolsky et al. 1999; Casse & Ferreira 2000a,b; Vlahakis et al. 2000). Such a field is assumed to arise from both advection of interstellar magnetic field and local dynamo generation (Rekowski et al. 2000). In the second class of models (referred to as "X-winds''), only a tiny region around the disk inner edge produces a jet (Camenzind 1990; Shu et al. 1994; Shu et al. 1995; Shu et al. 1996; Lovelace et al. 1999). The magnetic field is assumed to originally come from the protostar itself, after some eruptive phase that linked the disk inner edge to the protostellar magnetosphere. Note that in both models, jets extract angular momentum and mass from the underlying portion of the disk. However, by construction, "disk-winds'' are produced from a large spread in radii, while "X-winds'' arise from a single annulus. Apart from distinct disk physics, the difference in size and geometry of the ejection regions should also introduce some observable jet features. Another scenario has been proposed, where the protostar produces a fast collimated jet surrounded by a slow uncollimated disk wind or disk corona (Kwan & Tademaru 1988,1995; Kwan 1997), but such a scenario still lacks detailed calculations.

So far, all disc-driven jet calculations used a "cold'' approximation, i.e. negligible thermal pressure gradients. Therefore, each magnetic surface is assumed either isothermal or adiabatic. But to test which class of models is at work in T Tauri stars, reliable observational predictions must be made and the thermal equilibrium needs then to be solved along the flow. Such a difficult task is still not addressed in a fully self-consistent way, namely by solving together the coupled dynamics and energy equations. Thus, no model has been able yet to predict the gas excitation needed to generate observational predictions.

One first possibility is to use a posteriori a simple parameterization for the temperature and ionization fraction evolution along the flow. This was done by Shang et al. (1998) and Cabrit et al. (1999) for X-winds and disk winds respectively. These approaches are able to predict the rough jet morphology, but do not provide reliable flux and line profile predictions, since the thermal structure lacks full physical consistency.

The second possibility is to solve the thermal evolution a posteriori, with the difficulty of identifying the heating sources (subject to the constraint of consistency with the underlying dynamical solution). Several heating sources are indeed possible: (1) planar shocks (e.g. Hartigan et al. 1987,1994); (2) oblique magnetic shocks in recollimating winds (Ouyed & Pudritz 1993,1994); (3) turbulent mixing layers (e.g. Binette et al. 1999); and (4) current dissipation by ion-neutral collisions, referred to as ambipolar diffusion heating (Safier 1993a,1993b). A further heating scenario (not yet explored in the context of MHD jets and only valid in some environments) is photoionization from OB stars (Reipurth et al. 1998; Raga et al. 2000; Bally & Reipurth 2001). Of all these previous mechanisms only ambipolar diffusion heating allows "minimal'' thermal solutions, in the sense that the same physical process - non-vanishing currents - is responsible for jet dynamics and heating. As a consequence no additional tunable parameter is invoked for the thermal description. Furthermore, Safier (1993b) was able to obtain fluxes and profiles in reasonable agreement with observations. In this paper, we extend the work of Safier (1993a,1993b) by (1) using magnetically-driven jet solutions self-consistently computed with the underlying accretion disk, and (2) a more accurate treatment of ionization using the Mappings Ic code and ion-neutral momentum exchange rates which include the thermal contribution. In a companion paper (Garcia et al. 2001, hereafter Paper II), we generate predictions for spatially resolved orbidden line emission maps, long-slit spectra, and line ratios.

This article is structured as follows: in Sect. 2 we introduce the dynamical structure of the disk wind under study, and present physical values of the density, velocity, magnetic field, and Lorentz force along streamlines; in Sect. 3 we describe the physical processes taken into account in the thermal evolution computations, whose results are presented and discussed in Sect. 4. Conclusions are presented in Sect. 5. Some important derivations, dust description and consistency checks of our calculations are presented in the appendices.

   
2 Dynamical structure

2.1 General properties

A precise disk wind theory must explain how much matter is deviated from radial to vertical motion, as well as the amount of energy and angular momentum carried away. This implies a thorough treatment of both the disk interior and its matching with the jets, namely to consider magnetized accretion-ejection structures (hereafter MAES). The only way to solve such an entangled problem is to take into account all dynamical terms, a task that can be properly done within a self-similar framework.

In this paper, we use the models of Ferreira (1997) describing steady-state, axisymmetric MAES under the following assumptions: (i) a large scale magnetic field of bipolar topology is threading a geometrically thin disk; (ii) its ionization is such that MHD applies (neutrals are well-coupled to the magnetic field); (iii) some active turbulence inside the disk produces anomalous diffusion allowing matter to cross the field lines. Two extra simplifying assumptions were used: (iv) jets are assumed to be cold, i.e. powered by the magnetic Lorentz force only (the centrifugal force is due to the Lorentz azimuthal torque), with isothermal magnetic surfaces (the midplane temperature varying as $T_0 \propto r^{-1}$) and (v) jets carry away all disk angular momentum. This last assumption has been removed only recently by Casse & Ferreira (2000a).

All solutions obtained so far display the same asymptotic behavior. After an opening of the jet radius leading to a very efficient acceleration of the plasma, the jet undergoes a refocusing towards the axis (recollimation). All self-similar solutions are then terminated, most probably producing a shock (Gomez de Castro & Pudritz 1993; Ouyed & Pudritz 1993). This systematic behavior could well be imposed by the self-similar geometry itself and not be a general result (Ferreira 1997). Nevertheless, such a shock would occur in the asymptotic region, far away from the disk. Thus, we can confidently use these solutions in the acceleration zone, where forbidden emission lines are believed to be produced (Hartigan et al. 1995).

  
2.2 Model parameters

The isothermal self-similar MAES considered here are described with three free dimensionless local parameters (see Ferreira 1997, for more details) and four global quantities:
(1) the disk aspect ratio

$$\varepsilon= \frac{h(\varpi)}{\varpi}$$ (1)

where $h(\varpi)$ is the vertical scale height at the cylindrical radius $\varpi$;
(2) the MHD turbulence parameter

$$\alpha_{\rm m}= \frac{\nu_{\rm m}}{V_{\rm A}h}$$ (2)

where $\nu_{\rm m}$ is the required turbulent magnetic diffusivity and $V_{\rm A}$ the Alfvén speed at the disk midplane; this diffusivity allows matter to cross field lines and therefore to accrete towards the central star. It also controls the amplitude of the toroidal field at the disk surface.
(3) the ejection index

$$\xi= \frac{{\rm d} \ln \dot M_{{\rm acc}} (\varpi)}{{\rm d} \ln \varpi}$$ (3)

which measures locally the ejection efficiency ($\xi=0$ in a standard accretion disk), but also affects the jet opening (a higher $\xi$translates in a lower opening);
(4) M_* the mass of the central protostar;
(5) $\varpi_{\rm i}$ the inner edge of the MAES;
(6) $\varpi_{\rm e}$ the outer edge of the MAES, a standard accretion disk lying at greater radii. This outer radius is formally imposed by the amount of open, large scale magnetic flux threading the disk and producing jets;
(7) $\dot{M}_{\rm acc}$, the disk accretion rate fueling the MAES and measured at $\varpi_{\rm e}$.

For our present study, we keep only $\xi$ and $\dot{M}_{\rm acc}$ as free parameters and fix the values of the other five as follows: The disk aspect ratio was measured by Burrows et al. (1996) for HH 30 as $\sim$0.1 so we fix $\varepsilon =0.1$. The MHD turbulence parameter is taken $\alpha _{{\rm m}}=1$ in order to have powerful jets (Ferreira 1997). The stellar mass is fixed at $M_* =0.5 \ M_{\odot}$, typical for T Tauri stars, and the inner radius of the MAES is set to $\varpi_{\rm i} =$ 0.07 AU (typical disk corotation radius for a 10 days rotation period): inside this region the magnetic field topology could be significantly affected by the stellar magnetosphere-disk interaction. The outer radius is kept at $\varpi_{{\rm e}}= 1 \rm {AU}$ for consistency with the one fluid approximation (Appendix C) and the atomic gas description. Regarding atomic consistency, Safier (1993a) solved the flow evolution assuming inicially all H bound in H₂. He found H₂ to completely dissociate at the wind base, for small $\varpi _0$. However, after a critical flow line footpoint H₂ would not completely dissociate, therefore affecting the thermal solution. This critical footpoint was at 3 AU for his MHD solution nearer our parameter space.

We note that our two free parameters are still bounded by observational constraints: mass conservation relates the ejection index $\xi$ to the accretion/ejection rates ratios,

$$2 \dot{M}_{{\rm J}}\simeq \xi \dot{M}_{{\rm acc}} \ln\frac{\varpi_{\rm e}}{\varpi_{\rm i}}\cdot$$ (4)

The observational estimates for the ratio of mass outflow rate by mass accretion rate are $\dot{M}_{{\rm J}}/\dot{M}_{{\rm acc}}\simeq 0.01$ (Hartigan et al. 1995). The uncertainties affecting these estimates can be up to a factor of 10 (Gullbring et al. 1998; Lavalley-Fouquet et al. 2000). The range of ejection indexes considered here (0.005-0.01) is kept compatible with Hartigan's canonical value. The accretion rates $\dot{M}_{\rm acc}$ are also kept free but inside the observed range of $10^{-8}~M_\odot ~{\rm yr}^{-1}$ to $10^{-5}~M_\odot~{\rm yr}^{-1}$ in T Tauri stars (Hartigan et al. 1995).

Table 1 provides a list of some disk and jet parameters. These local parameters were constrained by steady-state requirements, namely the smooth crossing of MHD critical points. Disk parameters are useful to give us a view of the physical conditions inside the disk. Thus, the required magnetic field B₀ at the disk midplane and at a radial distance $\varpi _0$ is

$$B_0 = 0.3\ \zeta \left ( \frac{M_*}{M_{\odot}} \right)^\frac... ...1\,\rm {AU}} \right)^{\frac{\xi}{2} - \frac{5}{4}} \ \mbox{G}.$$ (5)

The global energy conservation of a cold MAES writes

$$P_{\rm acc} = 2P_{\rm jet}\ +\ 2P_{\rm rad}$$ (6)

where $P_{\rm acc}$ is the mechanical power liberated by the accretion flow, $P_{\rm jet}$ the total (kinetic, thermal and magnetic) power carried away by one jet and $P_{\rm rad}$ the luminosity radiated at one surface of the disk. For the solutions used, the accretion power is given by

$$P_{\rm acc} = \eta\ \frac{GM_* \dot M_{\rm acc}}{2\varpi_{\rm i}} 0.1\eta \left ( \frac{M_*}{M_{\odot}} \right) \left ( \frac{\dot M_{\rm acc}} {10^{-7}~M_{\odot}~\rm {yr}^{-1}} \right )$$ =

$$\times \left ( \frac{\varpi_{\rm i}}{0.07\,{\rm AU}} \right)^{-1}\ L_{\odot}$$ (7)

where $L_{\odot}$ is the solar luminosity and the efficiency factor $\eta= (\varpi_{\rm i}/\varpi_{\rm e})^\xi - (\varpi_{\rm i}/\varpi_{\rm e})$ depends on both the local ejection efficiency $\xi$ and the MAES radial extent. Typical values for our solutions are $\eta\simeq 0.9$. The ratio $P_{\rm jet}/P_{\rm rad}$ is given in Table 1. The jet parameters, mass load $\kappa$ and magnetic lever arm $\lambda$, have the same definition as in Blandford & Payne (1982). They are given here to allow a comparison with the solutions used in Safier's work.

 

 

Table 1: Isothermal MAES parameters. With $\varepsilon =0.1$ and $\alpha _{{\rm m}}=1$, the only parameter remaining free is $\xi$. Here, the magnetic lever arm $\lambda$, mass load $\kappa$ and initial jet opening angle $\theta _0$ are presented to ease comparison with Safier's models. However these parameters do not uniquely determine the MHD solution.

Solution	$\xi$	$\zeta$	$P_{\rm jet}/P_{\rm rad}$	$\kappa$	$\lambda$	$\theta_0 (\hbox{$^\circ$ })$
A	0.010	0.729	1.46	0.014	41.6	50.6
B	0.007	0.690	1.46	0.011	59.4	52.4
C	0.005	0.627	1.52	0.009	84.2	55.4

2.3 Physical quantities along streamlines

Figure 1: Several wind quantities along a streamline for model A (long-dashed line), B (solid) and C (dashed): jet nuclei density $\tilde n$, velocity, magnetic field, and Lorentz force. For the latter three vectors, poloidal components ( $v_{\rm p}$, $B_{\rm p}$, $(\vec J \times \vec B)_{\rm p}$) are plotted in black and toroidal components ($v_\phi$ $B_\phi$, $(\vec J \times \vec B)_\phi$) in red. The field line is anchored at $\varpi _0=0.1$ AU, around a $1~M_{\odot}$ protostar, with an accretion rate $\dot M_{\rm acc}= 10^{-6}~M_{\odot}~\rm{yr}^{-1}$.

Open with DEXTER

In order to obtain a solution for the MAES, a variable separation method has been used allowing to transform the set of coupled partially differential equations into a set of coupled ordinary differential equations (ODEs). Hence, the solution in the $(\varpi,z)$space is obtained by solving for a flow line and then scaling this solution to all space. Once a solution is found (for a given set of parameters in Sect. 2.2), the evolution of all wind quantities Q along any flow line is given by:

$$Q(\varpi,z) Q_0(\varpi_0) Q_\chi(\chi)$$ (8)

where the self-similar variable $\chi =z/\varpi _0$ measures the position along a streamline flowing along a magnetic surface anchored in the disk at $\varpi _0$. In particular, the flow line shape equation is given by $\varpi(z)= \varpi_0 \Psi(\chi)$, where the function $\Psi(\chi)$ is provided by solving the full dynamical problem. In Fig. 1 we plot the values of the jet nuclei density $\tilde n$ and the poloidal and toroidal components of jet velocity, magnetic field and Lorentz force. This is done for our 3 models, along a streamline with $\varpi _0$ = 0.1 AU, $M_{*}= 1 \ M_{\odot}$ and $\dot M_{\rm acc}= 10^{-6}~M_{\odot}~\rm{yr}^{-1}$. Values for other $\varpi _0$ and $\dot{M}_{\rm acc}$ can easily be deduced from these plots using the following scalings

 

$$\tilde{n} \propto \dot M_{\rm acc} M_\ast^{-\frac{1}{2}} \varpi_0^{\xi - \frac{3}{2}}$$

$$\propto M_\ast^\frac{1}{2} \varpi_0^{-\frac{1}{2}}$$ v

$$\propto \dot M_{\rm acc}^{\frac{1}{2}} M_\ast^\frac{1}{4} \varpi_0^{\frac{\xi}{2} - \frac{5}{4}}$$ B (9)

$$\propto \dot M_{\rm acc}^{\frac{1}{2}} M_\ast^\frac{1}{4} \varpi_0^{\frac{\xi}{2} - \frac{9}{4}}.$$ J

The terminal poloidal velocity is $v_{\rm {p},\infty} \simeq \sqrt{G M_\ast/ \xi \varpi_0}$ so that solutions with smaller opening angles also reach smaller terminal velocities, with higher terminal densities. The point where $v_\phi$ reaches a minimum is also the point where the jet reaches its maximum width (we call it recollimation point), before the jet starts to bend towards the axis. The numerical solution becomes unreliable as we move away from this point. The MHD solution is stopped at the super-Alfvénic point, which is reached nearer for higher $\xi$. An illustration of the resulting ($\varpi, z$) distribution of density $\tilde n$ and total velocity modulus for model A can be found in Fig. 1 of Cabrit et al. (1999).

3 Flow thermal and ionization processes

3.1 Main equations

Under stationarity, the thermal structure of an atomic (perfect) gas with density n and temperature T is given by the first law of thermodynamics:

 

$$P \vec{\nabla} \cdot \vec{v}\ +\ \vec{\nabla}\cdot U\vec{v} = \Gamma\ -\ \Lambda ,$$ (10)

where P= nkT is the gas pressure, $U=\frac{3}{2}nkT$ its internal energy density, $\vec{v}$ the total gas velocity and $\Gamma$ and $\Lambda$ are respectively the heating and cooling rates per unit volume. Since during most of the flow the ejected gas expands, we call the term $\Lambda_{{\rm adia}}= P \vec{\nabla} \cdot \vec{v}$the adiabatic cooling.

The gas considered here is composed of electrons, ions and neutrals of several atomic species, namely $n= n_{\rm e} + \overline{n_i} + \overline{n_n}$ where the overline stands for a sum over all present chemical elements. We then define the density of nuclei $\tilde{n}= \overline{n_i} + \overline{n_n}$ and the electron density $n_{\rm e}= f_{\rm e} \tilde{n}$. Correspondingly, the total velocity $\vec{v}$appearing in Eq. (10) must be understood as the barycentric velocity. As usual in one-fluid approximation, we suppose - and verify it in Sect. C.1 - all species well coupled (through collisions), so that they share the same temperature T. We also assume that no molecule formation occurs, so that mass conservation requires

 

$$\vec{\nabla}\cdot \tilde{n} \vec{v} = 0.$$ (11)

Under stationarity, the gas species ionization state evolves according to the rate equations,

 

$$\frac{{\rm D}f_{A^i}}{{\rm D}t} % \vec{v}\cdot \vec{\nabla} f_{A^i} = \frac{R_{A^i}}{\tilde n} ,$$ = (12)

$$\frac{{\rm D}f_{\rm e}}{{\rm D}t} % \vec{v}\cdot \vec{\nabla} f_{\rm e} = \sum_{A} \sum_{i=1}^{i \leqslant N_A} i \, \frac{R_{A^i}}{\tilde n} ,$$ = (13)

subjected to the elemental conservation constraint,

 

$$\sum_{i=0}^{i \leqslant N_A} R_{A^i} = 0 .$$ (14)

In the above equations $f_{A^i}= n_{A^i}/\tilde n$ is the population fraction of the chemical element A, i times ionized, R_Aⁱ the rate of change of that element and N_A is its maximum ionization state. Note that R_Aⁱ is a function of all the species densities n_Aⁱ, the temperature and the radiation field. In order to obtain the gas temperature and ionization state, we must solve energy equation coupled to the ionization evolution. This task requires to specify the ionization mechanisms as well as the gas heating ($\Gamma$) and cooling ($\Lambda$) processes (other than $\Lambda _{\rm adia}$).

3.2 Ionization evolution

3.2.1 The Mappings Ic code

We solve the gas ionization state (Eqs. (12) to (14)) using the Mappings Ic code - Binette et al. (1985); Binette & Robinson (1987); Ferruit et al. (1997). This code considers atomic gas composed by the chemical elements H, He, C, N, O, Ne, Fe, Mg, Si, S, Ca, Ar. We also added Na (whose ionization evolution is not solved by Mappings Ic), assuming it to be completely ionized in Na II. Hydrogen and Helium are treated as five level atoms.

The rate equations solved by Mappings Ic include photoionization, collisional ionization, secondary ionization due to energetic photoelectrons, charge exchange, recombination and dielectronic recombination. This is in contrast with Safier, who assumed a fixed ionization fraction for the heavy elements and solved only for the ionization evolution of H and He, considering two levels for H and only the ground level for He.

The adopted abundances are presented in Table 2. In contrast with Safier, we take into account heavy element depletion onto dust grains (see Sect. 3.2.3 and Appendix B) in the dusty region of the wind.

 

Table 2: Abundances by number of various elements with respect to hydrogen. The notation 3.55(-4) means $3.55\times 10^{-4}$. Column $Z_\odot$gives solar abundances, from Savage & Sembach (1996). Column $Z_{{\rm DL84}}$ elemental abundances locked in grains for a MRN-type dust distribution, from Draine & Lee (1984). Column $Z_{{\rm d}}$gives the abundances locked in grains in the diffuse clouds toward $\zeta$ Ophiuchi (Savage & Sembach 1996)), computed using the solar gas phase abundances from the same authors. The depleted abundances adopted here are $Z=Z_\odot -Z_{{\rm d}}$.

$$\begin{tabular}{ >{\scriptsize }c >{\scriptsize }c >{\scripts... ...3($-$6) \\ Na &2.04($-$6)& &1.81($-$6)\\ \hline \end{tabular}$$

  
3.2.2 Photoionizing radiation field

For simplicity, the central source radiation field is described in exactly the same way as in Safier and we refer the reader to the expressions (C1-C10) presented in his Appendix C. This radiation field is diluted with distance but is also absorbed by intervening wind material ejected at smaller radii.

We treat the radiative transfer as a simple absorption of the diluted central source, namely

 

$$4 \pi J_\nu(r, \theta) = \frac{ L_\nu (\theta)}{4 \pi r^2} {\rm e}^{- \tau_\nu(r, \theta)},$$ (15)

where $J_\nu(r,\theta)$ is the local mean monochromatic intensity at a spherical radius r and angle with the disk vertical $\theta$, $L_\nu(\theta)$ is the emitted luminosity of both star and boundary layer and $\tau_\nu(r,\theta)$ the optical depth towards the central object.

We now address the question of optical depth. In our model, the flow is hollow, starting from a ring located at the inner disk radius $\varpi_{\rm i}$ and extending to the outer radius $\varpi_{\rm e}$. The jet inner boundary is therefore exposed to the central ionizing radiation, which produces then a small layer where hydrogen is completely photoionized. The width $\Delta r$ of this layer can be computed by equating the number of emitted H ionizing photons, $Q({\rm H}_0) = \frac{1}{2} \int_{\nu_0}^{\infty} L_\nu {\rm d}\nu / h \nu$, to the number of recombinations in this layer, $n_{\rm H}^2\alpha_{\rm B}(T)2 \pi r^2\Delta r$ for our geometry. We found that $\Delta r \ll r$, and thus assume that all photons capable of ionizing hydrogen are exhausted within this thin shell. Furthermore, there is presumably matter in the inner "hollow'' region, so the previous considerations are upper limits.

For the heavy elements, photoionization optical depths are negligible, due to the much smaller abundances, and are thus ignored. The opacity $\tau_\nu$ is assumed to be dominated by dust absorption (see Appendix B for details). Dust will influence the ionization structure at the base of the flow, where ionization is dominated by heavy elements.

To summarize, the adopted radiation field is a central source absorbed by dust, with a cutoff at and above the Hydrogen ionization frequency.

  
3.2.3 Dust properties and gas depletion

Safier showed that if dust exists inside the disk, then the wind drag will lift the dust thereby creating a dusty wind. Our wind shares the same property. We model the dust (Appendix B) as a mix of graphite and astronomical silicate, with a MRN size distribution and use for the dust optical properties the tabulated values of Draine & Lee (1984), Draine & Malhotra (1993), Laor & Draine (1993). For simplicity we assumed the dust to be stationary, in thermodynamic equilibrium with the central radiation field and averaged all dust quantities by the MRN size distribution.

In addition, we take into account depletion of heavy elements into the dust phase. This effect was not considered by Safier. In Table 2 we present the dust phase abundances needed to maintain the MRN distribution (Draine & Lee), and our adopted depleted abundances, taken from observations of diffuse clouds toward $\zeta$Ori (Savage & Sembach 1996). These are more realistic, although presenting less depletion of carbon than required by MRN. Depletion has only a small effect on the calculated wind thermal structure, but can be significant when comparing to observed line ratios based on depleted elements.

3.3 Heating and cooling mechanisms

3.3.1 MHD heating

The dissipation of electric currents $\vec J$ provides a local heating term per unit volume $\Gamma_{\rm MHD}= \vec J \cdot (\vec E + \frac{\vec v} {c} \times \vec B)$, where $\vec E$ and $\vec B$ are the electric and magnetic fields, $\vec{v}$ the fluid velocity and cthe speed of light. In a multi-component gas, with electrons and several ion and neutral species, the generalized Ohm's law writes

$$\vec E + \frac{\vec v}{c} \times\vec B = \overline{\eta} \vec{J} ... ...vec{\nabla} P_{\rm e}}{e n_{\rm e}} + \frac{\vec J \times\vec B}{e n_{\rm e}c},$$ (16)

where $\eta$ is the fluid electrical resistivity, $\rho$ and $\rho_{\rm n}$are the total and neutral mass density, $m_{\rm in}$ the reduced ion-neutral mass, $n_{\rm i}$ the ionic density and $\nu_{\rm in}$ the ion-neutral collision frequency. The overline stands for a sum over all chemical elements relevant to a given quantity. These last quantities depend on the gas ionization state, the temperature and the momentum exchange rate coefficients. The reader is referred to Appendix A for the expressions of these coefficients and the uncertainties affecting them.

The first term appearing in the right hand side of the generalized Ohm's law is the usual Ohm's term, while the second describes the ambipolar diffusion, the third is the electric field due to the electron pressure and the last is the Hall term. This last effect provides no net dissipation in contrary to the other three. It turns out that the dissipation due to the electronic pressure is quite negligible and has been therefore omitted (Appendix C). Thus, the MHD dissipation writes

 

$$\Gamma_{{\rm MHD}} \overline{\eta}J^{2} + \Big(\frac{ \overline{\rho\rm _n} } {\rho ... ...ec{J}\times \vec{B} \bigr\Vert^{2}}{\overline{m_{\rm in} n\rm _i \nu_{\rm in}}}$$ =

$$\equiv \Gamma_{{\rm Ohm}} + \Gamma_{{\rm drag}}.$$ (17)

The first term, Ohmic heating $\Gamma_{{\rm Ohm}}$, arises mainly from ion-electron drag. The second term is the ambipolar diffusion heating $\Gamma _{\rm drag}$ and is mainly due to ion-neutral drag. This last term is the dominant heating mechanism in Safier's disk wind models, as well as in ours.

An important difference with Safier is that we take into account thermal speeds in ion-neutral momentum exchange rate coefficients. This increases $\nu_{\rm in}$, and results in significantly smaller ionization fractions (Sect. 4.8).

3.3.2 Ionization/recombination cooling

Both collisional ionization cooling $\Lambda_{{\rm col}}$ and radiative recombination cooling $\Lambda_{{\rm rec}}$ effects are taken into account by Mappings Ic. These terms are given by,

$$\Lambda_{{\rm col}} \sum_{A} \sum_{i=0}^{i<N_A} R^{\rm c}_{A^i} I_{A^i}$$ = (18)

$$\Lambda_{{\rm rec}} \sum_{A} \sum_{i=0}^{i<N_A} n_{\rm e} n_{A^{i+1}} k T\beta_{\rm B}(A^i)$$ = (19)

where $R^{\rm c}_{A^i}$ is the collisional ionization rate, I_Aⁱ the ionization energy and $\beta_{\rm B}(A^i)$ is the case B radiative recombination rate to the ionization state Aⁱ.

These ionization/recombination effects, taken into account in part by Safier, are in general smaller than adiabatic and line cooling.

3.3.3 Photoionization heating and radiative cooling

Photoionization by the radiation field, not taken into account by Safier, provides an extra source of heating $\Gamma_{{\rm P}}$. This term, which is also computed by Mappings Ic, is given by

$$\Gamma_{{\rm P}} = \sum_A \sum_{i=0}^{i<N_A} n_{A^i} \int_{\nu_0^... ...\infty} \frac{4 \pi J_\nu}{h \nu} a_\nu^{A^i} (h\nu - h \nu_0^{A^i}) {\rm d}\nu$$ (20)

where $\nu_0^{A^i}$ is the threshold frequency for ionization of the chemical element A, i times ionized, $4\pi J_\nu$ is the radiation flux of the central source (described in Sect. 3.2.2) and $a_\nu^{A^i}$ the photoionization cross-section. We found it to be the dominant heating source at the base of the flow, at the inner radii and for high accretion rates.

Collisionally excited lines provide a very efficient way to cool the gas, thanks to an extensive set of resonance and inter-combination lines, as well as forbidden lines. This radiative cooling $\Lambda_{\rm rad}$ is computed by Mappings Ic by solving for each atom the local statistical equilibrium, and will allow us to compute emission maps and line profiles for comparison with observations (see Garcia et al. 2001). We include cooling by hydrogen lines, $\Lambda_{\rm rad}({\rm H})$, in particular H$\alpha$, which could not be computed by Safier (two-level atom description).

3.3.4 Other minor heating/cooling mechanisms

Several processes, also computed by Mappings Ic, appeared to be very small and not affecting the jet thermal structure. We just cite them here for completeness: free-free cooling and heating, two-photon continuum and Compton scattering.

We ignored thermal conduction, which could be relevant along flow lines, the magnetic field reducing the gas thermal conductibility in any other direction. Also ignored was gas cooling by dust grains and heating by cosmic rays. We checked a posteriori that these three terms indeed have a negligible contribution (see Appendix C.2).

3.4 Numerical resolution

In our study, flow thermodynamics are decoupled from the dynamics - cold jet approximation. The previous Eqs. (10) to (14) can then be integrated for a given flow pattern. The dynamical quantities (density, velocity and magnetic fields) are given by the cold MHD solutions presented in Sect. 2. For the steady-state, axisymmetric, self-similar MHD winds under study, any total derivative writes

$$\frac{{\rm D}}{{\rm D}t} = (\vec{v} \cdot \vec{\nabla}) = \frac{v_z}{\varpi_0} \frac{\rm d}{{\rm d} \chi}$$ (21)

where v_z is the vertical velocity and $\chi =z/\varpi _0$ is the self-similar variable that measures the position along a flow line anchored at $\varpi _0$. With this in hand, and the mass conservation condition for an atomic wind (Eq. (11)), the energy Eq. (10) becomes an ODE along the flow line:

 

$$\frac{{\rm d} \ln T}{{\rm d} \chi} \frac{2}{3} \frac{{\rm d}\ln\tilde{n}}{ {\rm d} \chi} - \frac{\rm... ...1+ f_{\rm e}) + \frac{\Gamma - \Lambda}{\frac{3}{2} n k T \frac{v_z}{\varpi_0}}$$ =

$$\simeq \frac{2}{3} \frac{{\rm d} \ln \tilde n}{{\rm d}\chi} \left ( 1 - ... ...\rm MHD}} + (\Gamma -\Lambda)_{{\rm Map}}} {\Lambda_{{\rm adia}}} \right )\cdot$$ (22)

The term ${\rm d}\ln(1+f_{\rm e})/{\rm d}\chi$ (due to variations in internal energy density) is in fact negligible and has not been implemented for computational simplicity (see Appendix C.2). The term $(\Gamma - \Lambda)_{{\rm Map}}= \Gamma_{{\rm P}} - \Lambda_{{\rm\rm rad}} - \Lambda_{{\rm col}} - \Lambda_{{\rm rec}}$ is provided by Mappings Ic. The wind thermal structure is computed by integrating along a flow line the energy Eq. (22) coupled to the set of electronic population equation (Eqs. (12) to (14)).

3.4.1 Initial conditions

The integration of the set Eqs. (12) to (14) and (22) along the flow is an initial value problem. Thus, some way to estimate the initial temperature and populations must be devised. All calculations start at the slow-magnetosonic (SM) point, which is roughly at two scale heights above the disk midplane (for the solutions used here).

To estimate the initial temperature, Safier equated the poloidal flow speed at the SM point to the sound speed. Although this estimation agrees with cold flow theory, it is inconsistent with the energy equation which is used further up in the jet. Our approach was then to compute the initial temperature and ionic populations assuming that

$$\frac{{\rm D} T}{{\rm D}t}\bigg\vert_{\chi=\chi_{\rm s}}=0 ~~~~~\mbox{and}~~~~~ \frac{{\rm D} n_{A^i}}{{\rm D}t}\bigg\vert_{\chi=\chi_{\rm s}} =0$$ (23)

is fulfilled at $\chi=\chi\rm _s$. We thus assume for convenience that, at the base of the jet, there is no strong variations neither in temperature nor in ionization fractions. Physically this means that the gas is in ionization equilibrium, at the obtained temperature, with the incoming radiation field. The temperature thus obtained is always smaller than that provided by the SM speed. This is due to the large opening of the magnetic surfaces, providing a dominant adiabatic cooling over all heating processes. For numerical reasons the minimum possible initial temperature was set to 50 K. It is noteworthy that the exact value of the initial temperature only affects the base of the flow, below a few thousand degrees (see Appendix C.3). These regions are too cold to contribute significantly to optical line emission, leaving observational predictions unaffected.

The initial populations are computed by Mappings Ic assuming ionization equilibrium with the incoming radiation field. However for high accretion rates and for the outer zones of the wind, dust opacity and inclination effects shield completely the ionizing radiation field. The temperature is too low for collisional ionization to be effective. The ionization fraction thus reaches our prescribed minimum - all Na is in the form Na II (Table 2) and all the other elements (computed by Mappings Ic) neutral. However, soon the gas flow gains height and the ionization field is strong enough such that the ionization is self-consistently computed by Mappings Ic.

3.4.2 Integration procedure

After obtaining an initial temperature and ionization state for the gas we proceed by integrating the system of equations. In practice the ionization evolution is computed by Mappings Ic and separately we integrate Eq. (22) with a Runge-Kutta type algorithm (Press et al. 1988). We maintain both the populations and Mappings Ic cooling/heating rates per $\tilde{n}^2$ fixed during each temperature integrating spatial step. After we call Mappings Ic to evolve the populations and rates, at the new temperature, during the time taken by the fluid to move the spatial step. This step is such, that the RK integration has a numerical accuracy of 10^-6and, that the newly computed temperature varies by less than a factor of 10^-4. Such a small variation in temperature allows us to assume constant rates and populations while solving the energy equation. We checked a few integrations by redoing them at half the step used and found that the error in the ionic fraction is <10^-3in the jet, and < 10^-2 in the recollimation zone; the temperature precision being roughly a few times better. This ensures an intrinsic numerical precision comfortably below the accuracy of the atomic data and the $\langle\sigma_{\rm {H}\,{\rm H}^+}v\rangle$ collision cross-sections which, coupled to abundance incertitudes, are the main limitating factors. Details on the actual numerical procedure used by Mappings Ic to compute the non-equilibrium gas evolution are given in Binette et al. (1985).

  
4 Thermal structure results

In this section we present the calculated thermal and ionization structure along wind flow lines, discuss the physical origin of the temperature plateau and its connection with the underlying MHD solution, discuss the effect of various key model parameters and finally compare our results with those found by Safier. The parameters spanned for the calculation of the thermal solutions are the wind ejection index $\xi$ describing the flow line geometry, the mass accretion rate $\dot{M}_{\rm acc}$ and the cylindrical radius $\varpi _0$ where the field is anchored in the disk.

4.1 Temperature evolution

In Fig. 2, solid curves present the out of ionization equilibrium evolution of temperature, electronic density, and proton fraction along flow lines with $\varpi _0=0.1$ and 1 AU, as a function of $\chi =z/\varpi _0$, for accretion rates ranging from 10^-8 to $10^{-5}~M_{\odot}$ yr^-1. For comparison purposes, dashed curves plot the same quantities calculated assuming ionization equilibrium at the local temperature and radiation field. For compactness we present only these detailed results for our model B, with an intermediate ejection index $\xi = 0.007$. We divide the flow in three regions: the base, the jet and the recollimation zone. These regions are separated by the Alfvén point and the recollimation point (where the axial distance reaches its maximum). We only present the initial part of the recollimation zone here, because the dynamical solution is less reliable further out, where gas pressure is increased by compression and may not be negligible anymore. Note that the recollimation zone was not yet reached over the scales of interest in the solutions used by Safier.

Figure 2: Several wind quantities versus $\chi =z/\varpi _0$ for model B. The out of ionization equilibrium calculations are the solid curves and, for comparison, the ionization equilibrium are the dashed ones. The vertical dotted lines mark the Alfvén point and recollimation point. Top: temperature, middle: electronic density $n_{\rm e}$, bottom: proton fraction $f_{\rm p}= n({\rm H}^+)/n_{\rm H}$. The accretion rate $\dot{M}_{\rm acc}$ increases in the direction of the arrow from 10-8 to $10^{-5}~M_{\odot}$ yr-1 in factors of 10.

Open with DEXTER

The gas temperature increases steeply at the wind base (after an initial cooling phase for high $\dot{M}_{\rm acc} \ge 10^{-6}~M_{\odot}$ yr^-1). It then stabilizes in a hot temperature plateau around $T \simeq$ 1-3 $\times 10^4$ K, before increasing again after the recollimation point through compressive heating. The plateau is reached further out for larger accretion rates and larger $\varpi _0$. Its temperature decreases with increasing $\dot{M}_{\rm acc}$. The temperature plateau and its behavior with $\dot{M}_{\rm acc}$ were first identified by Safier in his wind solutions. We will discuss in Sect. 4.4 why they represent a robust property of magnetically-driven disk winds heated by ambipolar diffusion.

4.2 Ionization and electronic density

The bottom panels of Fig. 2 plot the proton fraction $f_{\rm p}= n({\rm H}^+)/n_{\rm H}$ along the flow lines. It rises steeply with wind temperature through collisional ionization, reaching a value ${\simeq} 10^{-4}$ at the beginning of the temperature plateau. Beyond this point, it continues to increase but starts to "lag behind'' the ionization equilibrium calculations (dashed curves): the density decline in the expanding wind increases the ionization and recombination timescales. Eventually, for $\chi \gtrsim 100$, density is so low that these timescales become longer than the dynamical ones, and the proton fraction becomes completely "frozen-in'' at a constant level, typically a factor 2-3 below the value found in ionization equilibrium calculations (dashed curves).

The electron density ($n_{\rm e}$) evolution is shown in the middle panels of Fig. 2. In the jet region, where $f_{{\rm p}}$is roughly constant, the dominant decreasing pattern with $\chi$ is set by the wind density evolution as the gas speeds up and expands. Similarly, the rise in $n_{\rm e}$ in the recollimation zone is due to gas compression. A remarkable result is that, as long as ionization is dominated by hydrogen (i.e. $f_{{\rm p}} \gtrsim 10^{-4}$), $n_{\rm e}$ is not highly dependent of $\dot{M}_{\rm acc}$, increasing by a factor of 10 only over three orders in magnitude in accretion rate. This indicates a roughly inverse scaling of $f_{{\rm p}}$ with $\dot{M}_{\rm acc}$ (bottom panels of Fig. 2), a property already found by Safier which we will discuss in more detail later.

In regions at the wind base where $f_{\rm p} < 10^{-4}$, variations of $n_{\rm e}$ are linked to the detailed photoionization of heavy elements which are then the dominant electron donors. The respective contributions of various ionized heavy atoms to the electronic fraction $f_{\rm e}$ is illustrated in Fig. 3 for $\dot{M}_{\rm acc}=10^{-6}~ M_{\odot}$ yr^-1. While O II and N II are strongly coupled to hydrogen collisional ionization through charge exchange reactions, the other elements are dominated by photoionization. The sharp discontinuity in C II and Na II at the wind base for $\varpi _0=0.1$ AU is caused by the crossing of the dust sublimation surface by the streamline (see Appendix B). Inside the surface we are in the dust sublimation zone where heavy atoms are consequently not depleted onto grains and hence have a higher abundance. In contrast, for $\varpi _0=1$ AU, the flow starts already outside the sublimation radius, in a region well-shielded from the UV flux of the boundary-layer, where only Na is ionized. Extinction progressively decreases as material is lifted above the disk plane and sulfur, then carbon, also become completely photoionized.

	Figure 3: Ion abundances with respect to hydrogen ( $f_{A^i}= n_{A^i}/\tilde n$ and thus fAi depends also on the abundances) along the flow line versus $\chi \equiv z/\varpi _0$ for model B in out of ionization equilibrium, with $\dot M_{\rm acc}= 10^{-6}~M_{\odot}~\rm{yr}^{-1}$. The jump at $\chi \sim 0.5$ is due to depletion as the gas enters the sublimation surface.
Open with DEXTER

4.3 Heating and cooling processes

The heating and cooling terms along the streamlines for our out of equilibrium calculations are plotted in Fig. 4 for $\varpi _0=0.1$ and 1 AU, and for two values of $\dot{M}_{\rm acc}$= 10^-6 and $10^{-7}~M_{\odot}$ yr^-1.

Figure 4: Heating and cooling processes (in erg s-1 cm-3) along the flow line versus $\chi \equiv z/\varpi _0$ for model B. Top figures for $\dot M_{\rm acc}= 10^{-6}~M_{\odot}~\rm{yr}^{-1}$and bottom ones for $\dot{M}_{\rm acc}=10^{-7} ~M_\odot~{\rm yr}^{-1}$. Ambipolar heating and adiabatic cooling appear to be the dominant terms, although Hydrogen line cooling cannot be neglected for the inner streamlines.

Open with DEXTER

Before the recollimation point, the main cooling process throughout the flow is adiabatic cooling $\Lambda _{\rm adia}$, although Hydrogen line cooling $\Lambda_{\rm rad}({\rm H})$ is definitely not negligible. The main heating process is ambipolar diffusion $\Gamma _{\rm drag}$. The only exception occurs at the wind base for small $\varpi_0 \le$ 0.1 AU and large $\dot{M}_{\rm acc} \ge 10^{-6}~M_{\odot}$ yr^-1, where photoionization heating $\Gamma_{{\rm P}}$ initially dominates. Under such conditions, ambipolar diffusion heating is low due to the high ion density, which couples them to neutrals and reduces the drift responsible for drag heating. However, $\Gamma_{{\rm P}}$ decays very fast due to the combined effects of radiation dilution, dust opacity, depletion of heavy atoms in the dust phase, and the decrease in gas density. At the same time, the latter two effects make $\Gamma _{\rm drag}$ rise and become quickly the dominant heating term. In the recollimation zone, the main cooling process is hydrogen line cooling $\Lambda_{\rm rad}({\rm H})$, and the main heating term is compression heating ( $\Lambda _{\rm adia}$ is negative).

A striking result in Fig. 4, also found by Safier, is that a close match is quickly established along each streamline between $\Lambda _{\rm adia}$ and $\Gamma _{\rm drag}$, and is maintained until the recollimation region. The value of $\chi$ where this balance is established is also where the temperature plateau starts. We will demonstrate below why this is so for the class of MHD wind solutions considered here.

  
4.4 Physical origin of the temperature plateau

The existence of a hot temperature plateau where $\Lambda _{\rm adia}$exactly balances $\Gamma _{\rm drag}$ is the most remarkable and robust property of magnetically-driven disk winds heated by ambipolar diffusion. Furthermore, it occurs throughout several decades along the flow including the zone of the jet that current observations are able to spatially resolve.

In this section, we explore in detail which generic properties of our MHD solution allow a temperature plateau at ${\simeq} 10^4$ K to be reached, and why this equilibrium may not be reached for other MHD wind solutions.

4.4.1 Context

First, we note that the energy equation (Eq. (22)) in the region where drag heating and adiabatic cooling are the dominant terms (which includes the plateau region) can be cast in the simplified form:

 

$$\frac{{\rm d} \ln T}{{\rm d} \ln \chi} -\frac{2}{3} \frac{ {\rm d} \ln \tilde{n}}{ {\rm d} \ln \chi} \times \Big(\frac{\Gamma_{\rm drag}}{\Lambda_{\rm adia}}-1\Big)$$ =

$$\delta^{-1} \times \Big(\frac{G(\chi)}{F(T)}-1\Big),$$ = (24)

where:

 

$$\delta(\chi)^{-1} \equiv \Big(-\frac{2}{3}\frac{{\rm d} \ln \tilde{n}}{{\rm d}\ln\chi}\Big),$$ (25)

remains positive before recollimation and depends only on the MHD solution,

 

$$G(\chi)=-\frac{\frac{1}{c^2} \bigl\Vert \vec{J} \times \vec{B} \b... ...vec{\nabla})\tilde{n}} \propto \frac{M_{\star}^2}{\dot{M}_{{\rm acc}} \varpi_0}$$ (26)

is a positive function, before recollimation, that depends only on the MHD wind solution, and

$$F(T)={kT(1+f_{\rm e})}{\overline{m_{\rm in} f_{\rm i} f_{\rm n} \... ...m in}v\rangle}} \times {\Big( \frac{ \rho }{ \overline{\rho}_{\rm n} } \Big)^2}$$ (27)

also positive, depends only on the local temperature and ionization state of the gas. The functions G and F separate the contributions of the MHD dynamics and ionization processes in the final thermal solution.

Figure 5: Left: function F(T) in erg g cm3 s-1versus temperature assuming local ionization equilibrium and an ionization flux that ionizes only all Na and all C. Center: function $G(\chi )$ in erg g cm3 s-1 for models A, B, C (bottom to top), $\dot M_{\rm acc}= 10^{-6}~M_{\odot}~\rm{yr}^{-1}$ and $\varpi_0=0.1~{\rm AU}$. Right: temperature for model B from the complete calculations in ionization equilibrium (dashed), and assuming $T =T_{\ominus }$ as given by Eq. (30) ( dash-dotted). $\varpi _0=0.1$ AU and accretion rates $\dot{M}_{\rm acc}$ are 10-8 to $10^{-5}~M_{\odot}$ yr-1, from top to bottom.

Open with DEXTER

The function $\delta$ is roughly constant and around unity before recollimation, it diverges at the recollimation point and becomes negative after it. Throughout the plateau  $\delta \sim 1$.

The "wind function'' G is plotted in the center panel of Fig. 5 for our 3 solutions. It rises by 5 orders of magnitude at the wind base and then stabilizes in the jet region (until it diverges to infinity near the recollimation point). The physical reason for its behavior is better seen if we note that the main force driving the flow is the Lorentz force:

$$G\propto \Big\vert\Big\vert\frac{{\rm D}\vec{v}}{{\rm D}t}\Big\vert\Big\vert^2 \times \Big(\frac{{\rm D}\rho}{{\rm D}t}\Big)^{-1}.$$ (28)

For an expanding and accelerating flow $({\rm D}\rho/{\rm D}t)^{-1}$ is an increasing function. At the wind base the Lorentz force accelerates the gas thus causing a fast increase in G. Once the Alfvén point is reached, the acceleration is smaller and ${\rm D}\vec{v}/{\rm D}t$decreases. However this decrease is not so abrupt as in the case of a spherical wind, because the Lorentz force is still at work, both accelerating and collimating the flow. This collimation in turn reduces the rate of increase of $(D\rho/Dt)^{-1}$. The stabilization of G observed after the Alfvén point is thus closely linked to the jet dynamics.

The "ionization function'' F is in general a rising function of Tand is plotted in the left panel of Fig. 5 under the approximation of local ionization equilibrium. Two regimes are present: in the low temperature regime, $f_{\rm i}$  $\gg f_{\rm p}$ is dominated by the abundance of photoionized heavy elements and $F(T) \propto T\,f_{\rm i}$ increases linearly with T, for fixed $f_{\rm i}$. The effect of the UV flux in this region is to shift vertically F(T): for a low UV flux regime only Na is ionized and $f_{{\rm i}} \simeq f(Na{\sc ii})$; for a high UV flux regime were Carbon is fully ionized, $f_{{\rm i}} \simeq f(C{\sc ii})$. In the high temperature regime ( $T \geq 8000$ K) where hydrogen collisional ionization dominates, $f_{\rm i}$ $\simeq f_{\rm p}$, and $F(T) \propto T\,f_{\rm p}$ becomes a steeply rising function of temperature, until hydrogen is fully ionized around $T \simeq 2\times 10^4$ K. The following second rise in F(T) is due to Helium collisional ionization. As we go out of perfect local ionization equilibrium the effect is to decrease the slope of F(T) in the region where H ionization dominates. In the extreme situation of ionization freezing, F(T) becomes linear again as in the photoionized region: $F(T) \propto T f_{{\rm p,freezed}}$.

4.4.2 Conditions needed for a hot plateau

The plateau is simply a region where the temperature does not vary much,

 

$$\frac{{\rm d} \ln T}{{\rm d} \ln \chi} = \epsilon {\hspace{1cm }\rm with \hspace{1cm }} \vert\epsilon\vert\ll 1.$$ (29)

Naively, temperatures $T_\ominus(\chi)$ defined by,

 

$$F(T_\ominus)= G(\chi) \Longleftrightarrow \Gamma_{\rm {drag}} = \Lambda_{\rm {adia}},$$ (30)

will zero the right hand side of Eq. (24) and thus satisfy the plateau condition. This equality is the first constraint on the wind function G, because there must exist a temperature $T_\ominus$ such that the equality holds. However this condition is not sufficient. Indeed the above equality describes a curve $T_\ominus(\chi)$ which must be flat in order to satisfy the plateau condition (Eq. (29)). Therefore the requirement of a flat $T_\ominus$ translates in a second constraint on the variation of G with respect to F. This constraint is obtained by differentiating Eq. (30). We obtain

 

$$\frac{{\rm d} \ln G}{{\rm d} \ln \chi} = \frac{{\rm d} \ln F}{{\rm d} \ln T}\bigg\vert_{T = T_\ominus} \times \epsilon$$ (31)

and after using (Eq. (29)):

$$\Big\vert\Big\vert\frac{{\rm d} \ln G}{{\rm d} \ln \chi}\Big\vert... ...rm d} \ln F}{{\rm d} \ln T} \bigg\vert_{T = T_\ominus}(\chi)\Big\vert\Big\vert.$$ (32)

Thus only winds where the wind function G varies much slower than the ionization function F will produce a plateau. This is fulfilled for our models: below the Alfvén surface, G varies a lot, but collisional H ionization is sufficiently close to ionization equilibrium that F(T) still rises steeply around 10⁴ K (Fig. 5). For our numerical values of G, within our range of $\dot{M}_{\rm acc}$ and $\varpi _0$, we have $T_\ominus \simeq 10^4$ K and thus $\vert{\rm d} \ln G/{\rm d} \ln \chi\vert \ll \vert{\rm d} \ln F/{\rm d} \ln T\vert$. Further out, where ionization is frozen out, we have ${\rm d} \ln F/{\rm d} \ln T = 1$ (because $F\propto T f_{\rm p, freezed}$) but it turns out that in this region G is a slowly varying function of $\chi$, and thus we still have $\vert{\rm d} \ln G/{\rm d} \ln \chi\vert \ll \vert{\rm d} \ln F/{\rm d} \ln T\vert$.

Finally, a third constraint is that the flow must quickly reach the plateau solution $T(\chi) = T_\ominus(\chi)$ and tend to maintain this equilibrium. Let us assume that $T =T_{\ominus }$ is fulfilled at $\chi=\chi_0$, what will be the temperature at $\chi= \chi_0(1+x)$? Letting $T= T_\ominus(1 + \vartheta)$ and assuming $\epsilon \ll 1$, Eq. (24) gives us $\vartheta = \exp(- \alpha x)$, which provides an exponential convergence towards $T =T_{\ominus }$ as long as

$$\alpha= \frac{1}{\delta} \frac{{\rm d} \ln F}{{\rm d} \ln T}\bigg\vert_{T = T_\ominus} > 0 \ .$$ (33)

Note that $\alpha$ depends mainly on the MHD solution (independent of $\dot{M}_{\rm acc}$/$\varpi _0$) and that the steeper the function F, the faster the convergence. The above criterion is always fulfilled in the expanding region of our atomic wind solutions, where $\delta > 0$ and F increases with temperature. The physical reason for the convergence can be easily understood in the following way: it can be readily seen that if at a given point $T> T_\ominus(\chi)$, then $G(\chi)/F(T) < 1$ and the gas will cool (cf. Eq. (24) with $\delta > 0$). Conversely, if $T<T_\ominus(\chi)$, the gas will heat up. Thus, for $\delta > 0$, the fact that F(T) is a rising function introduces a feedback that brings and maintains the temperature close to its local equilibrium value $T_\ominus(\chi)$, and $\Lambda _{\rm adia}$ close to $\Gamma _{\rm drag}$.

We conclude that three analytical criteria must be met by any MHD wind dominated by ambipolar diffusion heating and adiabatic cooling, in order to converge to a hot temperature plateau:
(1) Equilibrium: the wind function G must be such that $F(T)=G(\chi)$ is possible around $T \simeq 10^4$ K;
(2) Small temperature variation: the wind function $G(\chi )$ must vary slower than the ionization function F(T) such that $\vert{\rm d} \ln G/{\rm d} \ln \chi\vert \ll \vert{\rm d} \ln F/{\rm d} \ln T\vert$: (i) the wind must be in ionization equilibrium, or near it, in regions where G is a fast function of $\chi$; (ii) once we have ionization freezing, G must vary slowly, with $\vert {\rm d} \ln G / {\rm d}\ln \chi \vert \ll 1$;
(3) Convergence: $\alpha >0$, i.e. ${2 \over 3} ({\rm d} \ln \tilde{n} / {\rm d} \ln \chi) \times ({\rm d} \ln F / {\rm d} \ln T) \vert_{T = T_\ominus} < 0$, which is always verified for an atomic and expanding wind.

4.4.3 Comments on other MHD winds

Not all types of MHD wind solutions will verify our first criterion. Physically, the large values of $G(\chi )$ observed in our solutions indicates that there is still a non-negligible Lorentz force after the Alfvén surface. In this region (which we call the jet) the Lorentz force is dominated by its poloidal component which both collimates and accelerates the gas (Fig. 1). The gas acceleration translates in a further decrease in density contributing to a further increase in G. Models that provide most of the flow acceleration before the Alfvén surface might turn out to have a lower wind function $G(\chi )$, not numerically compatible with the steep portion of the ionization function F(T). These models would not establish a temperature plateau around 10⁴ K by ambipolar diffusion heating. They would either stabilize on a lower temperature plateau (on the linear part of F(T)) if our second criterion is verified, or continue to cool if G varies too fast for the second criterion to hold. This is the case in particular for the analytical wind models considered by Ruden et al. (1990), where the drag force was computed a posteriori from a prescribed velocity field. The G function for their parameter space (Table 3 of Ruden et al.) peaks at ${\sim}10^{-48}$ erg g cm³ s^-1 at ${\sim} 3~ R_{\star}$ and then rapidly decreases as $G \propto r^{-1}$ for higher radii. This translates into a cooling wind without a plateau.

4.5 Scalings of plateau parameters with $\dot{M}_{acc}$ and $\varpi _0$

Figure 6: Verification of the plateau scalings (Eq. (34)). Left: we plot the measured $T_\ominus f_{{\rm p},\ominus }$ versus $(\dot{M}_{{\rm acc}}\varpi_0)^{-1}$ for all models in out of ionization equilibrium. The evolution is linear except for the very edges of the $(\dot{M}_{{\rm acc}}\varpi_0)^{-1}$ domain. This is due to the failure of our assumptions: in the lower edge $f_{\rm e} < 10^{-4}$ and thus it isn't dominated by H; in the upper edge $f_{\rm e} > 0.1$ and thus the small ionization fraction approximation fails. The crosses are from Safier; because of a smaller momentum transfer rate coefficient they are quite above model C. Right: we plot the measured $f_{{\rm p},\ominus }$ versus $(\dot{M}_{{\rm acc}}\varpi_0)^{-1}$ for Model B in out of ionization equilibrium (solid) and ionization equilibrium (dashed). Note that all the variation is absorbed by $f_{{\rm p},\ominus }$ and thus $f_{{\rm p},\ominus} \propto (\dot{M}_{{\rm acc}}\,\varpi_0)^{-1}$. A straight line is also plot for comparison.

Open with DEXTER

The balance between drag heating and adiabatic cooling (Eq. (30)) can further be used to understand the scalings of the plateau temperature $T_\ominus$ and proton fraction $f_{{\rm p},\ominus }$ with the accretion rate $\dot{M}_{\rm acc}$ and flow line footpoint radius ($\varpi _0$). In the plateau region, ionization is intermediate, i.e., sufficiently high to be dominated by protons but with most of the Hydrogen neutral. Under these conditions we have $F(T)\propto T_{\ominus} f_{{\rm p},\ominus}$. On the other hand, self-similar disk wind models display $G(\chi) \propto (\dot{M}_{{\rm acc}}\varpi_0)^{-1}$ (Eq. (26)). Therefore, we expect

 

$$T_\ominus f_{\rm p} \propto \frac{1} {\dot{M}_{{\rm acc}}\,\varpi_0}\cdot$$ (34)

This behavior is verified in the left panel of Fig. 6 for our out-of-equilibrium results for the 3 wind solutions.

To predict how much of this scaling will be absorbed by $T_\ominus$and how much by $f_{{\rm p},\ominus }$, Safier considered the ionization equilibrium approximation: for the temperature range of the plateau ($T\sim10^4$ K) $f_{\rm p}\simeq$$f_{\rm i}$ is a very fast varying function of T that can be approximated as $f_{\rm p}\propto T^{a}$ with $a \gg 1$. One predicts that $T_{\ominus}\propto (\dot{M}_{{\rm acc}}\,\varpi_0)^{-1/(a+1)}$while $f_{{\rm p},\ominus}\propto (\dot{M}_{{\rm acc}}\,\varpi_0)^{-a/(a+1)} \simeq (\dot{M}_{{\rm acc}}\,\varpi_0)^{-1}$. Hence, the inverse scaling with ( $\dot{M}_{\rm acc}\varpi_0$) should be mostly absorbed by $f_{{\rm p},\ominus }$, while the plateau temperature is only weakly dependent on these parameters. This is verified in the right panel of Fig. 6, where $f_{{\rm p},\ominus }$ in ionization equilibrium is plotted as a function of $(\dot{M}_{{\rm acc}}\,\varpi_0)^{-1}$. The predicted scaling is indeed closely followed.

Let us now turn back to the actual out-of-equilibrium calculations. At the base of the flow we find that the wind evolves roughly in ionization equilibrium (see Fig. 2), however at a certain point the ionization fraction freezes at values that are near those of the ionization equilibrium zone. This effect implies that the ionization fraction should roughly scale as the ionization equilibrium values at the upper wind base. This in indeed observed in Fig. 6. This memory of the ionization equilibrium values by $f_{{\rm p}}$ (as observed in the solar wind by Owocki et al. 1983) is the reason why the scalings of $f_{{\rm p},\ominus }$ with ( $\dot{M}_{\rm acc}\varpi_0$) remain correct. We computed for our solutions the scalings and found for model B: $f_{{\rm p},\ominus}\propto\dot{M}_{{\rm acc}}^{-0.76}$, $f_{{\rm p},\ominus}\propto\varpi_0^{-0.83}$, $T_{\ominus}\propto\dot{M}_{{\rm acc}}^{-0.13}$ and no dependence of T on $\varpi _0$, confirming the memory effect on the ionization fraction only.

Finally, we note that for accretion rates in excess of a few times $10^{-5}~M_{\odot}$  yr^-1, the hot plateau should not be present anymore: Because of its inverse scaling with $\dot{M}_{\rm acc}$, the wind function G remains below 10^-47 erg g cm³ s^-1, and F = G occurs below 10⁴ K, on the linear low-temperature part of F(T) where $f_{{\rm i}} \simeq f(C{\sc ii})$(see Fig. 5). These colder jets will presumably be partly molecular. Interestingly, molecular jets have only been observed so far in embedded protostars with high accretion rates (e.g. Gueth & Guilloteau 1999).

4.6 Effect of the ejection index $\xi$

The importance of the underlying MHD solution is illustrated in Fig. 5. The ejection index $\xi$ is directly linked to the mass loaded in the jet ($\kappa$, see Table 1 and Ferreira 1997). Thus a higher $\xi$ translates in an stronger adiabatic cooling because more mass is being ejected. The ambipolar diffusion heating is less sensitive to the ejection index, because the density increase is balanced by a stronger magnetic field. Hence, the wind function $G(\chi )$ decreases with increasing $\xi$. As a result, the plateau temperature and ionization fraction also decrease (see Eq. (34)).

In Fig. 7 we summarize our results for the three models by plotting the plateau $f_{{\rm p},\ominus }$ versus $T_\ominus$, for several $\dot{M}_{\rm acc}$ (values of $\varpi _0$are connected together). In this plane, our MHD solutions lie in a well-defined "strip'' located below the ionization equilibrium curve, between the two dotted curves. For a given model, as $\dot{M}_{\rm acc}$ increases, the plateau ionization fraction and the temperature both decrease, as expected from the scalings discussed above, moving the model to the lower-left of the strip. Increasing the ejection index decreases $G(\chi )$, and it can be seen that this has a similar effect as increasing $\dot{M}_{\rm acc}$ (Eq. (34)).

Figure 7: Out of ionization equilibrium evolution of wind quantities in the plateau. Points are for all flow lines ( $\varpi _0=0.07\times 1.3^i$ AU and i=0,1,...,10) and accretion rates ( $\dot{M}_{\rm acc}=10^{-i}~M_\odot~{{\rm yr}}^{-1}$ with i=5,6,7,8). We plot, $f_{{\rm p},\ominus }$ versus $T_\ominus$ for all models C (red), B (black), A (blue) and in red the values expected from ionization equilibrium. The dashed/dotted lines are for models without depletion, the solid lines are for models with depletion. The accretion rates increase in the direction top right to bottom left. The thick solid curve traces $f_{{\rm p}}$ in ionization equilibrium, while the two dotted lines embrace the locus of our MHD solutions.

Open with DEXTER

4.7 Depletion effects on the thermal structure

In our calculations we take into account depletion of heavy species into the dust phase. We ran our model with and without depletion and found these effects to be minor. Changes are only found when $f_{{\rm p}} \lesssim 10^{-4}$. The temperature without depletion is slightly reduced (the higher ionization fraction reduces $\Gamma _{\rm drag}$) and as a consequence $f_{{\rm p}}$ is also smaller. Normally these changes affect only the wind base, as the temperature increases $f_{{\rm p}}$ dominates the ionization and we obtain the same results for the plateau zone. However for high accretion rates ( $\dot{M}_{\rm acc}=10^{-5}~M_\odot~{\rm yr}^{-1}$) in the outer wind zone (large $\varpi _0$) we still have $f_{{\rm p}} \lesssim 10^{-4}$ for the plateau and thus the temperature without depletion is reduced there.

  
4.8 Differences with the results of Safier

	Figure 8: $\langle\sigma_{\mbox{\tiny H H$^+$ }}v\rangle$ in cm3 s-1 using Draine expression (solid), using Geiss & Buergi expression (see Appendix A.2) (dash) and Draine expression but ignoring the thermal velocity (dot).
Open with DEXTER

The most striking difference between our results and those of Safier is an ionization fraction 10 to 100 times smaller. This difference is mainly due to both different $\langle\sigma_{{\rm H H}^+} v\rangle$ momentum transfer rate coefficients and dynamical MHD wind models.

The critical importance of the momentum transfer rate coefficient ( $\langle\sigma_{{\rm H H}^+} v\rangle$) for the plateau ionization fraction can be seen by repeating the reasoning in the previous section but including the momentum transfer rate coefficient in the scalings. We thus obtain

$$T_\ominus\,f_{{\rm p},\ominus} \propto \frac{1}{\langle\sigma_{{\rm H H}^+} v\rangle \dot{M}_{{\rm acc}} \varpi_0} \cdot$$ (35)

This shows that, because the freezing of the ionization fraction is correlated to the ionization fraction at the base of the wind (which is in ionization equilibrium), $f_{{\rm p}}$ will scale with the momentum transfer rate coefficient value. This means that if the momentum transfer rate coefficient is larger, there is a better coupling between ions and neutrals and hence a smaller drag heating. For the calculation of the $\langle\sigma_{{\rm H H}^+} v\rangle$, Safier ignored the contribution of the thermal velocity in the collisional relative velocity. This considerably reduces $\langle\sigma_{{\rm H H}^+} v\rangle$ and thus, increases $f_{{\rm p}}$. In Fig. 8 we plot the corresponding momentum transfer rate coefficient values. It can be seen that ignoring the thermal contribution to the momentum transfer rate coefficient decreases it typically by a factor of $\gtrsim$6. We also plot in this figure the value obtained by Geiss & Buergi (1986) illustrating the uncertainties in the momentum transfer rate coefficient (more on this in Appendix A.2).

5 Concluding remarks

We performed detailed non-equilibrium calculations of the thermal and ionization structure of atomic, self-similar magnetically-driven jets from keplerian accretion disks. Current dissipation in ion-neutral collisions - ambipolar diffusion -, was assumed as the major heating source. Improvements over the original work of (Safier 1993a,b) include: a) detailed dynamical models by Ferreira (1997) where the disk is self-consistently taken into account but each magnetic surface assumed isothermal; b) ionization evolution for all relevant "heavy atoms'' using Mappings Ic code; c) radiation cooling by hydrogen lines, recombination and photoionization heating using Mappings Ic code; d) H-H⁺ momentum exchange rates including thermal contributions; and e) more detailed dust description.

We obtain, as Safier, warm jets with a hot temperature plateau at $T \simeq 10^4$ K. Such a plateau is a robust property of the atomic disk winds considered here for accretion rates less than a few times $10^{-5}~M_{\odot}$ yr ^-1. It is a direct consequence of the characteristic behavior of the wind function $G(\chi )$ defined in Eq. (26): (i) $G(\chi )$ increases first and becomes larger than a certain value fixed by the minimum ionization fraction (see Fig. 5); (ii) $G(\chi )$ is flat whenever ionization freezing occurs (collimated jet region). More generally, we formulate three analytical criteria that must be met by any MHD wind dominated by ambipolar diffusion heating and adiabatic cooling in order to converge to a hot temperature plateau.

The scalings of ionization fractions and temperatures in the plateau with $\dot{M}_{\rm acc}$ and $\varpi _0$ found by Safier are recovered. However the ionizations fractions are 10 to 100 times smaller, due to larger H-H⁺ momentum exchange rates (which include the dominant thermal velocity contribution ignored by Safier) and to different MHD wind dynamics.

We performed detailed consistency checks for our solutions and found that local charge neutrality, gas thermalization, single fluid description and ideal MHD approximation are always verified by our solutions. However at low accretion rates, for the base of outer wind regions ( $\varpi_0\sim 1$ AU) and increasingly for higher $\xi$, single fluid calculations become questionable. For the kind of jets under study, a multi-component description is necessary for field lines anchored after $\varpi_0 > 1$ AU. So far, all jet calculations assumed either isothermal or adiabatic magnetic surfaces. But our thermal computations showed such an increase in jet temperature that thermal pressure gradients might become relevant in jet dynamics. We therefore checked the "cold'' fluid approximation by computing the ratio of the thermal pressure gradient to the Lorentz force, along ( $\beta _\parallel$) and perpendicular ( $\beta _\perp$) to a magnetic surface. Both ratios increase for lower accretion rates and outer wind regions. We found that for some solutions, thermal pressure gradients play indeed a role, however only at the wind base (possible acceleration) and in the recollimation zone (possible support against recollimation). Fortunately, (as will be seen in a companion paper, Paper II), the dynamical solutions which are found inconsistent are also those rejected on an observational ground. Therefore, it turns out that the models that best fit observations are indeed consistent.

Acknowledgements

P. J. V. G. acknowledges financial support from Fundação para a Ciência e Tecnologia by the PRAXIS XXI/BD/5780/95, PRAXIS XXI/BPD/20179/99 grants. The work of LB was supported by the CONACyT grant 32139-E. We thank the referee, Pedro N. Safier, for his helpful comments. We acknowledge fruitful discussions with Alex Raga, Eliana Pinho, Pierre Ferruit and Eric Thiébaut. We are grateful to Bruce Draine for communicating his grain opacity data, and to Pierre Ferruit for providing his C interface to the Mappings Ic routine. P. J. V. G. warmly thanks the Airi team and his adviser, Renaud Foy, for their constant support.

  
Appendix A: Multicomponent MHD equations

  
A.1 Single fluid description

Let us consider a fluid composed of $\alpha$ species of numerical density $n_\alpha$, mass $m_\alpha$, charge $q_\alpha$ and velocity $\vec v_\alpha$. All species are assumed to be coupled enough so that they have the same temperature T. To get a single fluid description, we then define

$$\rho \Sigma_\alpha n_\alpha m_\alpha$$ =

$$\rho \vec v \Sigma_\alpha m_\alpha n_\alpha \vec v_\alpha$$ =

$$\Sigma_\alpha n_\alpha k_{{\rm B}} T$$ P = (A.1)

$$\vec J \Sigma_\alpha n_\alpha q_\alpha \vec v_\alpha$$ =

as being the flow density $\rho$, velocity $\vec{v}$, pressure P and current density $\vec J$. We consider now a fluid composed of three species, namely electrons (e), ions (i) and neutrals (n). The equations of motion for each species are

   

$$\rho_{\rm n} \frac{D\vec{v}\rm _n}{Dt} -\vec{\nabla}P_{\rm n} - \rho_{\rm n} \vec{\nabla} \Phi_{\rm G} + \vec{F}_{\rm en} + \vec{F}_{\rm in}$$ = (A.2)

$$\rho_{\rm i} \frac{D\vec{v}\rm _i}{Dt} -\vec{\nabla}P_{\rm i} - \rho_{\rm i} \vec{\nabla} \Phi_{\rm G} +... ...\rm i} n_{\rm i} \left(\vec{E} +\frac{\vec{v}_{\rm i}}{c} \times \vec{B}\right)$$ = (A.3)

$$\rho_{\rm e} \frac{D\vec{v}_{\rm e}}{Dt} -\vec{\nabla}P_{\rm e} - \rho_{\rm e} \vec{\nabla} \Phi_{\rm G} +... ...\rm e}\,n_{\rm e} \left(\vec{E} +\frac{\vec{v}_{\rm e}}{c}\times \vec{B}\right)$$ = (A.4)

where $\Phi_{\rm G}$ is the gravitational potential and the collisional force of particles $\alpha$ on particles $\beta$ is given by $\vec{F}_{\alpha \beta} = m_{\alpha\beta} n_{\alpha} \nu_{\alpha\beta} (\vec{v}_\alpha - \vec{v}_\beta)$, $m_{\alpha\beta}= m_{\alpha} m_{\beta}/(m_{\alpha}+m_{\beta})$ being the reduced mass, $\nu_{\alpha\beta} = n_\beta \langle\sigma_{\alpha\beta} v\rangle$ the collisional frequency and $\langle\sigma_{\alpha\beta} v\rangle$ the averaged momentum transfer rate coefficient.

A single fluid dynamical description of several species is relevant whenever they are efficiently collisionally coupled, namely if they fulfill $\Vert\vec{v}_{\alpha=\rm e,n,i}- \vec{v} \Vert \ll \Vert \vec{v}\Vert$. Under this assumption and using Newton's principle ( $\Sigma_{\alpha, \beta} \vec F_{\alpha \beta}= \vec 0$), we get the usual MHD momentum conservation equation for one fluid

 

$$\rho \frac{D\vec{v}}{Dt} = -\vec{\nabla}P - \rho \vec{\nabla} \Phi_{\rm G} + \frac{1}{c}\vec{J} \times \vec{B}$$ (A.5)

by adding all equations for each specie. The Lorentz force acting on the mean flow is

 

$$\frac{1}{c}\vec{J} \times \vec{B} = (1 + X) (\vec F_{\rm in} + \vec F_{\rm en}) \ ,$$ (A.6)

where $X = \rho\rm _i/\rho\rm _n$. Even if the bulk of the flow is neutral, collisions with charged particles give rise to magnetic effects. In turn, the magnetic field is coupled to the flow by the currents generated there. This feedback is provided by the induction equation, which requires the knowledge of the local electric field $\vec E$. Its expression is obtained from the electrons momentum equation

$$\vec{E} + \frac{\vec{v}_{\rm e}}{c} \times \vec{B} \frac{1}{en_{\rm e}}\big(m_{\rm ie} n_{\rm i}\nu_{\rm ie} \vec{v}... ...nu_{\rm ne} \vec{v}_{\rm ne} \big) - \frac{\vec{\nabla} P_{\rm e}}{e n_{\rm e}}$$ = (A.7)

where $\vec{v}_{\alpha \beta} \equiv \vec v_\alpha - \vec v_\beta$ is the drift velocity between the two species. Due to their negligible contribution to the mass of the bulk flow, all terms involving the electrons inertia have been neglected (electrons quite instantly adjust themselves to the other forces).

All drift velocities can be easily obtained. The electron-ion drift velocity is directly provided by $\vec{v}_{\rm ie} = \vec J/e n_{\rm e}$. Using Eq. (A.6) and noting that $\rho_{\rm e} \nu_{\rm en} \ll \rho_{\rm i} \nu_{\rm in}$ we get the ion-neutral drift velocity

$$\vec{v}_{\rm in}=\frac{\frac{1}{c}\vec{J}\times\vec{B}}{(1+X) m_{... ..._{\rm en}}{m_{\rm in} n_{\rm i} \nu_{\rm in}} \frac{\vec{J}}{ e n_{\rm e}}\cdot$$ (A.8)

On the same line of thought, the electrons velocity is $\vec{v}_{\rm e} = \vec{v} -(\vec{v}-\vec{v}_{\rm e})$ where

$$\vec{v} - \vec{v}_{\rm e} \frac{\vec{v}_{\rm n} - \vec{v}_{\rm e}}{1+X} + \frac{X}{1+X}(\vec{v}_{\rm i} - \vec{v}_{\rm e})$$ = (A.9)

$$\simeq \frac{\vec J}{e n_{\rm e}} - \frac{1}{(1+X)^2} \frac{\frac{1}{c} \vec{J}\times\vec{B}}{m_{\rm in} n_{\rm i}\nu_{\rm in}}\cdot$$ (A.10)

Gathering these expressions for all drift velocities, we obtain the generalized Ohm's law

 

$$\vec{E} + \frac{\vec{v} }{c} \times \vec{B} \eta \vec{J} + \frac{\frac{1}{c}\vec{J}\times\vec{B}}{e n_{\rm e}} - \frac{\vec{\nabla} P_{\rm e}}{en_{\rm e}} - \frac{1}{(1+X)^2} \frac{\frac{1}{c^2}(\vec{J}\times\vec{B})\times\vec{B}}{m_{\rm in} n_{\rm i} \nu_{\rm in}}$$ (A.11)

where $\eta=(m_{\rm ne}n_{\rm n} \nu_{\rm ne} + m_{\rm ie} n_{\rm i}\nu_{\rm ie})/(en_{\rm e})^2$is the electrical resistivity due to collisions. The corresponding MHD heating rate writes

$$\Gamma_{\rm MHD}= \vec{J}\cdot\vec{E'} = \vec{J} \cdot \left (\vec{E} + \frac{\vec{v}}{c} \times \vec{B} \right )$$ (A.12)

where $\vec{E'}$ is the electrical field in the comoving frame. This expression leads to Eq. (17).

The generalization of this derivation for a mixture of several chemical elements has been done in a quite straightforward way. The bulk flow density becomes $\rho= \overline{\rho_{\rm i}} + \overline{\rho_{\rm n}} + \rho_{\rm e}$, where the overline stands for a sum over all elements (ions and neutrals), with $X=\overline{\rho_{\rm i}}/\overline{\rho_{\rm n}}$. The neutrals and ions velocities are means over all elements, $\langle \vec{v}_{\rm n,i}\rangle \equiv \sum_{\rm n,i} \rho_{\rm n,i} \vec{v}_{\rm n,i}/\overline{\rho_{\rm n,i}}$. The conductivity and collision terms are also sums over all elements, namely $\overline{\eta}=(\overline{m_{\rm ne}n_{\rm n}\nu_{\rm ne}} + \overline{m_{\rm ie} n_{\rm i}\nu_{\rm ie}})/(en_{\rm e})^2$ and $\overline{m_{\rm in} n_{\rm i} \nu_{\rm in}}$, and are computed using the expressions for the collision frequencies.

  
A.2 Momentum transfer rate coefficient

For ion-electron collisions we use the canonical from Schunk (1975), summed over all species:

$$\overline{m_{\rm ie} n_{\rm i} \nu_{\rm ie}}=m_{\rm e} n_{\rm e} ... ...B}} T_{\rm e}}{\pi m_{\rm e}}}\sum_i n_{\rm i} Z_{\rm i}^2 \ln{\Lambda_{\rm i}}$$ (A.13)

with the Coulomb factor $\Lambda_{\rm i}=(3/2 Z_{\rm i} e^3)\sqrt{(k_{{\rm B}}T_{\rm e})^3/ \pi n_{\rm e}}$.

For the collisions between electrons and neutrals we use the expression of Osterbrock (1961) for the collisional momentum transfer rate coefficient between a neutral and a charged particle, which corrects the classical one (e.g., Schunk 1975) for strong repulsive forces at close distances. Its expression is $\langle\sigma v\rangle_{\rm n,i-e} = 2.41 \pi e \sqrt{\alpha_{\rm n}/m_{\rm n,i-e}}$, where the polarizabilities $\alpha_{\rm n}$ used are also taken from Osterbrock. We thus obtain

$$\overline{m_{\rm en} n_{\rm n} \nu_{\rm ne}} = 2.41 \pi\, e \,n_{\rm e} \sum_{n={\rm H,He}} n_{\rm n} \sqrt{m_{\rm e} \alpha_{\rm n}}\,.$$ (A.14)

Finally, it is mainly the ion-neutral collision momentum transfer rate coefficient determines the ambipolar diffusion heating. It can be computed with the previous momentum transfer rate coefficient expression. However as noted by Draine (1980) the previous expression underestimates $\sigma$ at high velocities. Thus, as Draine, we take the "hard sphere'' value for the cross-section ( $\sigma_{\rm S} \simeq 10^{-15}$ cm²) whenever it is superior to the polarizability one. For intermediate to hight ionizations ( $f_{\mbox{\tiny H$^+$ }}\gtrsim10^{-4}$) the dominant ion-neutral collisions are those between H-H⁺. Charge exchange effects between these two species will amplify $\langle\sigma_{\mbox{\tiny H H$^+$ }}v\rangle$above the values expected by polarizability alone and thus it is computed separately (Eqs. (A.16) and (A.17)). We thus obtain for ion-neutral collisions

$$\overline{m_{\rm in} n_{\rm i} \nu_{\rm in}} = \frac{1}{2} m_... ..._{\rm S}\,\tilde{v}; \langle\sigma v\rangle_{{\rm He},\rm i})}$$ (A.15)

where $\tilde{v}=\sqrt{8k_{{\rm B}}T/\pi m_{\rm in} + v_{\rm in}^2}$. For $\tilde{v}<1000$ kms^-1.

The value of $\langle\sigma_{\mbox{\tiny H H$^+$ }}v\rangle$ which we used is given by Draine (1980),

 

$$\frac{\langle\sigma_{ \mbox{{\tiny H\,H$^+$ }}}v\rangle}{1 \,{\rm... ...10^{-9}\,\tilde{v}^{0.73}&\tilde{v}\geq 2 \:\mbox{km s}^{-1} \end{array}\right.$$ (A.16)

Safier (1993a) used the expression $\tilde{v}=v_{\rm in}$ which, as discussed in Sect. 4.8, results in a smaller momentum transfer rate coefficient. Geiss & Buergi (1986) computed another expression of the H-H⁺ momentum transfer rate coefficient, which provides

$$\langle\sigma_{ \mbox{{\tiny H\,H$^+$ }}}v\rangle\ = 1.12\ti... ...1+2\log{T_4})^2\big)\times10^{-9}\:\: {\rm cm}^3~{\rm s}^{-1}.$$ (A.17)

In Fig. 8 we compared both momentum transfer rate coefficients, they typically differ in 40%, which can be used as an estimate of their accuracy. It is thus the uncertainty in the H-H⁺momentum transfer rate coefficient that dominates the final intrinsic uncertainty of our calculations.

  
Appendix B: Dust implementation

Figure B.1: Dust sublimation surfaces geometry for the adopted radiation field.

As shown by Safier if there is dust in the disk, the wind is powerful enough to drag it along. Thus disk winds are dusty winds. Dust is important for the wind thermal structure mainly as an opacity source affecting the photoionisation heating at the wind base. To compute the dust opacity we need a description of its size distribution, its wavelength dependent absorption cross-section and the inner dust sublimation surface. In the inner flow zones and for high accretion rates the strong stellar and boundary layer flux will sublimate the dust, creating a dust free inner cavity (see Fig. B.1). Results on the evolution of dust in accretion disks by Schmitt et al. (1997) show that at the disk surface the initial dust distribution isn't much affected by coagulation and sedimentation effects. Thus we assume a MRN dust distribution (Mathis et al. 1977; Draine & Lee 1984):

$${\rm d}n_{\rm i}=n_{\rm H}\,A_i\,a^{-3.5}\,{\rm d}a$$ (B.1)

where ${\rm d}n_{\rm i}$ is the number of particles of type i ("astronomical silicate'' - Sil or graphite - C) with sizes in $[a,a+{\rm d}a]$, and $0.005~\mbox{$\mu$ m} \leq a \leq 0.25 ~\mbox{$\mu$ m}$, $A_{\rm Sil}=10^{-25.11}$ cm^2.5 H^-1 and $A_{\rm C}=10^{-25.16}$cm^2.5 H^-1. We then proceed by averaging all relevant grain quantities function of size and species (F_i(a)) by the size/species distribution,

 

$$\langle F_i(a)\langle_a = \int_{a_{\rm min}}^{a_{\rm max}} \sum_{i={\rm Sil,C}} F_i(a) \, \frac{{\rm d}n_{\rm i}}{N\rm _T}\cdot$$ (B.2)

In order to compute the sublimation radius, some description of the dust temperature must be made. For simplicity, we assume the dust to be in thermodynamic equilibrium with the radiation field, the dominant dust heating mechanism. In our case, the central source radiation field will dominate throughout the jet, except probably in the recollimation zone, where the strong gas emission overcomes the central diluted field. However in this region dust is no longer relevant for the gas thermodynamics and we will therefore only consider dust heating by the central source. The dust temperature $T_{\rm gr}$ for a grain of size a is obtained by equating the absorbed to the emitted radiation (e.g., Tielens & Hollenbach 1985),

 

$$4 \pi a^2 \langle Q_a\rangle(T_{\rm gr}) \sigma T_{\rm gr}^4 = \pi a^2 \,\int_0^\infty Q_a^{\rm abs}(\nu) 4\pi J_\nu {\rm d}\nu$$ (B.3)

where a is the grain size $\langle Q_a\rangle(T_{\rm gr})$ is the Planck-averaged emissivity (Draine & Lee 1984; Laor & Draine 1993; Draine & Malhotra 1993), $\sigma$ is the Stefan-Boltzmann constant, $Q_a^{\rm abs}(\nu)$ is the dust absorption efficiency and $4\pi J_\nu$ is the central source radiation flux at the grain position given by Eq. (15). Averaging out the previous equation by the size/species distribution (Eq. (B.2)) we obtain,

$$4 \langle Q_a^{\rm em}(T_{\rm gr})\rangle \sigma T_{\rm gr}^4 = \... ...angle Q_a^{\rm abs}\rangle(\nu) F(\nu) {\rm e}^{-\tau_\nu(r,\theta)} {\rm d}\nu$$ (B.4)

where we describe the central source flux by $F(\nu)$ which is attenuated only by the dust opacity $\tau$. For simplicity $F(\nu)$ is taken as exactly the same as in Safier, i.e. a classical boundary layer (Bertout et al. 1988). The sublimation radius is obtained from the previous expression by noting that at its position $\tau_\nu=0$,

$$\frac{\langle r_{\rm sub}(\theta)\rangle}{R_\ast}=\sqrt{\frac{g_\... ...ngle T_{\rm bl}^4 }{4 \langle Q_a^{\rm em}(T_{\rm sub})\rangle T_{\rm sub}^4 }}$$ (B.5)

where $T_{\rm bl}$ and $T_{\ast}$ are the boundary layer and star temperatures, $R_\ast$ the stellar radius and $g_{\rm bl}(\theta)$/ $g_\ast(\theta)$ are the $\theta$ dependent terms of the radiation field (given in Bertout et al. and Safier). We assume a dust sublimation temperature $T_{\rm sub}$ of 1500 K.

With the dust sublimation radius in hand we can now proceed to compute the dust optical depth defined as,

$$\tau_\nu(r,\theta) = \int_{r_{\rm in}(\theta)}^{r} \kappa_\n... ...}^{r} n_{\rm H} (r',\theta) \overline{\sigma(\nu)_a} {\rm d}r'$$ (B.6)

where $r_{\rm in}(\theta)$ is the radius inside which there is no dust. This radius is given by the inner flow line $r_{\varpi_{\rm i}}(\theta)$ and by the sublimation radius $\langle r_{\rm sub}(\theta)\rangle$ (see Fig. B.1) such that $r_{\rm in} (\theta) = {\rm max} (\langle r_{\rm sub}\rangle;r_{\varpi_{\rm i}})$. The dust absorption cross-section ( $\overline{\sigma(\nu)_a}$) is,

$$\overline{\sigma(\nu)_a}=\int_{a_{\rm min}}^{a_{\rm max}} \pi a^2... ...nu) A_{\rm Sil} + Q^{\rm abs}_{\rm C}(a,\nu) A_{\rm C} \big] a^{-3.5} {\rm d}a.$$ (B.7)

Using the self-similarity of $n_{\rm H} (r,\theta)$ we can integrate Eq. (B.6) to obtain,

$$\tau_\nu(r,\theta) =n_{\rm H}(r,\theta)\,\overline{\sigma(\nu)} \, 2r\Big(\sqrt{\frac{r}{\langle r_{\rm in}(\theta)\rangle}}-1\Big)$$ (B.8)

which was used in Eq. (15). We note that at large distances from the source, the optical depth converges to a finite value, proportional to $\dot{M}_{\rm acc}$ and whose $\theta$variation is function of the self similar wind solution and central source radiation field. Thus for high accretion rates, although the central source radiation hardens, the outer zones of the wind base are less photoionized than for smaller accretion rates.

  
Appendix C: Consistency checks

  
C.1 Dynamical assumptions

Figure C.1: Top left: we plot the relevant drift speeds normalized to the fluid poloidal velocity. The worst case of one fluid approximation violation is obtained for model C, $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=1$ AU. Top right: we plot the $\tau _{\alpha \beta }$ versus dynamical $\tau _{\rm dyn}$time-scales (s) versus $\chi$, normalized to the Keplerian period at the line footpoint (for a 1 $M_\odot$ star). We plot only the longer time-scales ( $\tau _{\rm in} = \tau _{\rm ni}/f_i$ and $\tau _{\rm en} = \tau _{\rm ne}/f_{\rm e}$). The worst case is obtained again for model C, $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=1$ AU. Bottom left: ideal MHD tests for the worst situation (model C, $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=1$ AU). As expected from our heated winds, the ambipolar diffusion term is the dominant one. Bottom right: ratios of the thermal pressure gradient to the Lorentz force versus $\chi$, along ( $\beta _\parallel$, solid) and across ( $\beta _\perp$, dash) a magnetic surface anchored at $\varpi _0$. The worst case for $\beta _\parallel$ is obtained for model A, $\varpi _0=1$ AU, $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$, the best for model A, $\varpi _0=0.1$ AU and $\dot{M}_{\rm acc}=10^{-5}~M_\odot~{\rm yr}^{-1}$. With respect to $\beta _\perp$, the worst case is for model C, $\varpi _0=1$ AU and $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$, while the best is for model A, $\varpi _0=0.1$ AU and $\dot{M}_{\rm acc}=10^{-5}~M_\odot~{\rm yr}^{-1}$. Although definitely not negligible in some models, those compatible with observations do not show important deviations from the "cold'' jet approximation.

First, local charge neutrality is always achieved. For example, we achieve a maximum Debye length of $r_{{\rm D}}\sim10^5 {\rm cm}$ at the outer radius of the recollimation zone (model C, lowest $\dot{M}_{\rm acc}$).

Second, single fluid approximation requires that relative velocity drifts of all species ($\alpha=$ ions, electrons, neutrals) $\Vert\vec{v} - \vec{v}_{\alpha}\Vert / \Vert\vec{v}\Vert$ are smaller than unity. These drifts are higher for lower accretion rates and at the outer wind base (due to the decrease in density and velocity, see Eq. (9)). In Fig. C.1 we present the worst case for the drift velocities, showing that our jets can be indeed approximated by single fluid calculations.

We assumed gas thermalization, which is achieved only if collisional time-scales between species $\tau_{\alpha\beta} = 1/\nu_{\alpha \beta}$ are much smaller than the dynamical time-scale $\tau_{{\rm dyn}} = \varpi_0 ({\rm d} v_z/{\rm d} \chi)^{-1}$. In the collision network considered here, the longer time-scales involve collisions with neutrals. However, even in the worst situation (see Fig. C.1), after the wind base they remain comfortably below the above dynamical time-scale. Our dynamical jet solutions were derived within the ideal MHD framework. This assumption requires that all terms in the right hand side of the generalized Ohm's law (Eq. (A.11)) are negligible when compared to the electromotive field $\vec{v} \times \vec{B}/c$. We consider Ohm's term $\Vert\eta \vec{J} \Vert$, Hall's effect $\Vert\vec{J}\times\vec{B}\Vert/c\,e n_{\rm e}$ and the ambipolar diffusion term $(\frac{\overline{\rho_{\rm n}}}{\rho})^2 \Vert(\vec{J} \times\vec{B}) \times\vec{B}\Vert/c^2 \overline{m_{\rm in} n_{\rm i} \nu_{\rm in}}$ (effects due to the electronic pressure gradient are small compared to the Lorentz force - Hall's term -). In Fig. C.1 we present the worst case for our ideal MHD checks. We find that deviations from ideal MHD remain negligible, despite the presence of ambipolar diffusion. As expected, this is the dominant diffusion process in our (non turbulent) MHD jets. Ambipolar diffusion is larger for low accretion rates and at the outer wind base (because the ratio of the ambipolar to the electromotive term scales as $(\dot{M}_{{\rm acc}}\,f_i)^{-1}$).

The worst case for the previous three tests is, as expected, for the model that attains the lowest density: Model C, with the lowest accretion rate ( $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$) and at the outer edge footpoint ( $\varpi_0=1~{\rm AU}$).

The dynamical jet evolution was calculated under the additional assumption of negligible thermal pressure gradient (cold jets). Since it is the gradient that provides a force, one should not just measure (along one field line) the relative importance of the gas pressure to the magnetic pressure (usual $\beta= P/(B^2/8\pi)$parameter). Instead, we compare the thermal pressure gradient to the Lorentz force, along ( $\beta _\parallel$) and perpendicular ( $\beta _\perp$) to the flow, namely

$$\beta_\parallel \frac{\nabla_\parallel P}{F_\parallel}= c \frac{\vec{v}_{\rm p} \cdot \nabla P}{ \vec{v}_{\rm p} \cdot (\vec J \times \vec B)}$$ = (C.1)

$$\beta_\perp \frac{\nabla_\perp P}{F_\perp} = c \frac{\nabla a \cdot \nabla P}{ \nabla a \cdot (\vec J \times \vec B)}\cdot$$ = (C.2)

Here $a(\varpi, z)$ is the poloidal magnetic flux function, hence $\nabla a$ is perpendicular to a magnetic surface. High values of $\beta _\parallel$ imply that the thermal pressure gradient plays a role in gas acceleration, whereas high values of $\beta _\perp$show that it affects the gas collimation.

In Fig. C.1 we plot the worst case of cold fluid violation and best case of cold fluid validity. Again the worst case appears at lower accretion rates and in outer wind zones. It can be seen that high values of $\beta _\perp$ and $\beta _\parallel$ can be attained, hinting at the importance of gas heating on jet dynamics (providing both enthalpy at the base of the jet and/or pressure support against recollimation further out). We underline that models inconsistent with the cold fluid approximation are those found to have the largest difficulty in meeting the observations (Paper II). Conversely, models that better reproduce observations also fulfill the cold fluid approximation. For those models, the thermal pressure gradient appears to be fairly negligible with respect to the Lorentz force.

   
C.2 Thermal assumptions

Figure C.2: Ignored heating/cooling terms (in erg s-1 cm-3). Left: we compare the cooling term $- \frac{3}{2}k \tilde{n} T Df_{\rm e}/Dt$ (dashed) with adiabatic cooling $\Lambda _{\rm adia}$(solid) for the worst case (model C, $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=0.1$ AU). Right: grain heating/cooling rate $\vert\Gamma _{\rm gr}\vert$ (dashed) compared with ambipolar diffusion heating $\Gamma _{\rm drag}$ (solid) for the worst case (model A, $\dot{M}_{\rm acc}=10^{-5}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=1$ AU). $\Gamma _{\rm gr}$ is only significant at the very base of the wind, and will not affect the thermal state further out in the jet.

Finally, we check that all ignored heating/cooling processes are not relevant when compared to adiabatic cooling and ambipolar diffusion heating.

The first ignored process is the term $- \frac{3}{2}k \tilde{n} T {\rm D}f_{\rm e}/{\rm D}t$. This term decreases for increasing accretion rate and $\varpi _0$ due to the lower ionizations found in these regions. It is plotted in Fig. C.2 for the worst case (model C, $\dot{M}_{\rm acc}=10^{-8}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=0.1$ AU). There, it reaches at most $\sim$13% of $\Lambda _{\rm adia}$. Typical values for higher accretion rates are only $\simeq$0.1% of $\Lambda _{\rm adia}$.

Next we consider heating/cooling of the gas by collision with dust grains, given by Hollenbach & McKee (1979):

$$\Gamma_{{\rm gr}} = n \langle n_{{\rm gr}} \sigma_{{\rm gr}} \ran... ...k_{{\rm B}} T }{\pi m_{{\rm H}} } } f(2k_{{\rm B}} T_{\rm gr} - 2k_{{\rm B}} T)$$ (C.3)

where $\langle n_{{\rm gr}} \sigma_{{\rm gr}}\rangle$ is computed from the adopted MRN distribution, and f=0.16 is the sticking parameter that takes into account charge and accommodation effects for a warm gas (Hollenbach & McKee 1979). With these values the grain heating/cooling becomes,

$$\Gamma_{{\rm gr}} = 4.78 \times 10^{-34} n^2 (T_{{\rm gr}}-T) \hspace{.5cm} \mbox{erg s$^{-1}$\space cm$^{-3}$ }.$$ (C.4)

This term increases with increasing $\varpi _0$ and accretion rate. In Fig. C.2 we compare it (in absolute value) with $\Gamma _{\rm drag}$ in the case where its contribution is the most important (model A, $\dot{M}_{\rm acc}=10^{-5}~M_\odot~\mbox{yr}^{-1}$ and $\varpi _0=1$ AU). $\Gamma _{\rm gr}$ is initially positive (dust hotter than the gas), but changes sign at $\chi \simeq 0.4$, where the gas becomes hotter than the dust, becoming an effective cooling term. It is only significant at the very base of the wind, where it exceeds the ambipolar diffusion heating term by a factor $\simeq$3.5. However, this effect will not have important consequences in terms of observational predictions: we will show in next section that the thermal state in the hot plateau (where forbidden line emission is excited) is not sensitive to the initial temperature. Furthermore, the outer streamlines at $\varpi_0 \simeq 1$ AU contribute much less to the line emission than the inner ones. At lower accretion rates ${\le} 10^{-6}~M_\odot~\mbox{yr}^{-1}$, $\Gamma _{\rm gr}$ is always $\le$10% of $\Gamma _{\rm drag}$.

Heating due to cosmic rays, which could be important in the outer tenuous zones of the wind is (Spitzer & Tomasko 1968),

$$\Gamma_{\rm cr} = n_{\rm H} \zeta \Delta E = 1.9\times 10^{-28} \, n_{\rm H} \mbox{ erg s$^{-1}$\space cm$^{-3}$ }$$ (C.5)

where $\zeta$ is the ionization rate which we took as $\zeta=3.5\times 10^{-17}$ s^-1 (Webber 1998) and $\Delta E = 3.4$ eV is the average thermal energy transmitted to the gas by each ionization (Spitzer & Tomasko 1968). This effect is at most ${\sim} 3.6\times 10^{-6}$ times $\Gamma _{\rm drag}$ for model A, $\dot{M}_{\rm acc}=10^{-5}~M_\odot~{\rm yr}^{-1}$ and $\varpi _0=0.1$ AU. Finally the thermal conductivity along magnetic field lines was computed with the Spitzer conductivity for a fully ionized gas (Lang 1999) and is irrelevant (at most ${\sim}10^{-6}$ of the adiabatic cooling term) the maximum being achieved at the recollimation zone where the physical validity of our MHD solutions ends. It should be pointed out that Nowak & Ulmschneider (1977) compute the thermal conductivity for a partially ionized mixture in ionization equilibrium and found that for low temperatures $T\sim10^{3-4}$K the Spitzer expression underestimates the conductivity by a factor 10². However this is still too small to be important.

   
C.3 Dependence on initial conditions

Figure C.3: Effect of initial temperature on the thermal evolution for model B, $\dot M_{\rm acc}= 10^{-6}~M_{\odot}~\rm{yr}^{-1}$. The several initial temperatures are in dashed 50 K, 100 K, 500 K, 1000 K, 2000 K and 3000 K. In solid we plot the solution obtained by our standard initial conditions. Note that the almost vertical evolution of the temperature for very low initial temperatures is not an artifact. The field line anchored at $\varpi _0=0.1$ AU crosses the sublimation surface at $\chi \simeq 0.5$.

Formally, our temperature integration is an initial value problem. In the absence of a self-consistent description of the disc thermodynamics, there is some freedom in the initial temperature determination. It is therefore crucial to check that the subsequent thermal evolution of the wind does not depend critically on the adopted initial value.

Safier obtained the initial temperature by assuming the poloidal velocity at the slow magnetosonic point ( $v_{\rm p,s}$) to be the sound speed for adiabatic perturbations $T_{\rm s} = \mu m_{\rm H} v_{\rm p,s}^2/ \gamma k_{{\rm B}}$. Here, we have chosen to compute the initial temperature assuming local thermal equilibrium ${\rm D}T/{\rm D}t=0$. Our method produces lower initial temperatures than Safier due to adiabatic cooling.

For high accretion rates $\dot{M}_{\rm acc}\,\geq 10^{-6}~M_\odot~{\rm yr}^{-1}$ our initial temperature versus $\varpi _0$ has a minimum at the beginning of the dusty zone: inside the sublimation cavity, the thermal equilibrium is between photoionization heating and adiabatic cooling. Just beyond the dust sublimation radius, photoionization heating is strongly reduced, but the ionization fraction is still too high for efficient drag heating, resulting in a low initial equilibrium temperature.

The initial ionization fraction is similarly determined by assuming local ionization equilibrium ${\rm D}f_A^{\rm i}/{\rm D}t = 0$ for all elements. It decreases with $\varpi _0$. For $\dot{M}_{\rm acc}=10^{-5}~M_\odot~{\rm yr}^{-1}$ and $\varpi_0 \geq 0.8$ AU, the initial ionization fraction is set to a minimum value by assuming that Na is fully ionized, which is somewhat arbitrary. However, as gas is lifted up above the disk plane, the dust opacity decreases and the gas heats up, so that ionization becomes dominated by other photoionized heavy species and by protons, all computed self consistently.

In order to check that our results do not depend on the initial temperature, we have run model B for a broad range of initial temperatures. As shown in Fig. C.3), we find that the thermal and ionization evolution quickly becomes insensitive to the initial temperature. If we start with a temperature lower than the local isothermal condition, the dominant adiabatic cooling is strongly reduced, and the gas strongly heats up, quickly converging to our nominal curve. If we start with a higher initial temperature, adiabatic cooling is stronger, and we have the characteristic dip in the temperature found by Safier. Our choice of initial temperature has the advantage of reducing this dip, which is somewhat artificial (see Fig. C.3). In either case, we conclude that our results are robust with respect to the choice of initial temperature. In particular, the distance at which the hot plateau is reached, which has a crucial effect on line profile predictions, is unaffected.

References

 
• Bacciotti, F., Mundt, R., Ray, T. P., et al. 2000, ApJL, 537, L49 In the text NASA ADS  
• Bally, J., & Reipurth, B. 2001, ApJ, 546, 299 In the text NASA ADS  
• Bertout, C., Basri, G., & Bouvier, J. 1988, ApJ, 330, 350 In the text NASA ADS  
• Binette, L., Cabrit, S., Raga, A., & Cantó, J. 1999, A&A, 346, 260 In the text NASA ADS  
• Binette, L., Dopita, M. A., & Tuohy, I. R. 1985, ApJ, 297, 476 In the text NASA ADS  
• Binette, L., & Robinson, A. 1987, A&A, 177, 11 In the text NASA ADS  
• Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 In the text NASA ADS  
• Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858 In the text NASA ADS  
• Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437 In the text NASA ADS  
• Cabrit, S., Edwards, S., Strom, S. E., & Strom, K. M. 1990, ApJ, 354, 687 In the text NASA ADS  
• Cabrit, S., Ferreira, J., & Raga, A. C. 1999, A&A, 343, L61 In the text NASA ADS  
• Camenzind, M. 1990, Rev. Mod. Astrophys., 3, ed. G. Klare (Springer-Verlag, Berlin) In the text  
• Casse, F., & Ferreira, J. 2000a, A&A, 353, 1115 In the text NASA ADS  
• Casse, F., & Ferreira, J. 2000b, A&A, 361, 1178 In the text NASA ADS  
• Chan, K. L., & Henriksen, R. N. 1980, ApJ, 241, 534 In the text NASA ADS  
• Cohen, M., Emerson, J. P., & Beichman, C. A. 1989, ApJ, 339, 455 In the text NASA ADS  
• Dougados, C., Cabrit, S., Lavalley, C., & Ménard, F. 2000, A&A, 357, L61 In the text NASA ADS  
• Draine, B. T. 1980, ApJ, 241, 1021 In the text NASA ADS  
• Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 In the text NASA ADS  
• Draine, B. T., & Malhotra, S. 1993, ApJ, 414, 632 In the text NASA ADS  
• Ferreira, J. 1997, A&A, 319, 340 In the text NASA ADS  
• Ferreira, J., & Pelletier, G. 1993, A&A, 276, 625 In the text NASA ADS  
• Ferreira, J., & Pelletier, G. 1995, A&A, 295, 807 In the text NASA ADS  
• Ferruit, P., Binette, L., Sutherland, R. S., & Pecontal, E. 1997, A&A, 322, 73 In the text NASA ADS  
• Garcia, P. J. V., Cabrit, S., Ferreira, J., & Binette, L. 2001, A&A, 377, 609, Paper II In the text  
• Geiss, J., & Buergi, A. 1986, A&A, 159, 1 In the text NASA ADS  
• Gomez de Castro, A. I., & Pudritz, R. E. 1993, ApJ, 409, 748 In the text NASA ADS  
• Goodson, A. P., Böhm, K., & Winglee, R. M. 1999, ApJ, 524, 142 In the text NASA ADS  
• Gullbring, E., Hartmann, L., Briceno, C., & Calvet, N. 1998, ApJ, 492, 323 In the text NASA ADS  
• Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736 In the text NASA ADS  
• Hartigan, P., Morse, J. A., & Raymond, J. 1994, ApJ, 436, 125 In the text NASA ADS  
• Hartigan, P., Raymond, J., & Hartmann, L. 1987, ApJ, 316, 323 In the text NASA ADS  
• Heyvaerts, J., & Norman, C. 1989, ApJ, 347, 1055 In the text NASA ADS  
• Hirth, G. A., Mundt, R., & Solf, J. 1997, A&AS, 126, 437 In the text NASA ADS  
• Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 In the text NASA ADS  
• Krasnopolsky, R., Li, Z., & Blandford, R. 1999, ApJ, 526, 631 In the text NASA ADS  
• Kwan, J. 1997, ApJ, 489, 284 In the text NASA ADS  
• Kwan, J., & Tademaru, E. 1988, ApJL, 332, L41 In the text NASA ADS  
• Kwan, J., & Tademaru, E. 1995, ApJ, 454, 382 In the text NASA ADS  
• Lang, K. R. 1999, Astrophysical Formulae, vol. I (Springer Verlag) In the text  
• Laor, A., & Draine, B. T. 1993, ApJ, 402, 441 In the text NASA ADS  
• Lavalley-Fouquet, C., Cabrit, S., & Dougados, C. 2000, A&A, 356, L41 In the text NASA ADS  
• Li, Z. 1995, ApJ, 444, 848 In the text NASA ADS  
• Li, Z. 1996, ApJ, 465, 855 In the text NASA ADS  
• Lovelace, R. V. E., Romanova, M. M., & Bisnovatyi-Kogan, G. S. 1999, ApJ, 514, 368 In the text NASA ADS  
• Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 In the text NASA ADS  
• Mundt, R., & Eislöffel, J. 1998, AJ, 116, 860 In the text NASA ADS  
• Nowak, T., & Ulmschneider, P. 1977, A&A, 60, 413 In the text NASA ADS  
• Osterbrock, D. E. 1961, ApJL, 134, 270 In the text NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1993, ApJ, 419, 255 In the text NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1994, ApJ, 423, 753 In the text NASA ADS  
• Ouyed, R., & Pudritz, R. E. 1997, ApJ, 484, 794 In the text NASA ADS  
• Owocki, S. P., Holzer, T. E., & Hundhausen, A. J. 1983, ApJ, 275, 354 NASA ADS  
• Raga, A., López-Martín, L., Binette, L., et al. 2000, MNRAS, 314, 681 In the text NASA ADS  
• Ray, T. P., Mundt, R., Dyson, J. E., Falle, S. A. E. G., & Raga, A. C. 1996, ApJL, 468, L103 In the text NASA ADS  
• Reipurth, B., Bally, J., Fesen, R. A., & Devine, D. 1998, Nature, 396, 343 In the text NASA ADS  
• Rekowski, M. V., Rüdiger, G., & Elstner, D. 2000, A&A, 353, 813 In the text NASA ADS  
• Ruden, S. P., Glassgold, A. E., & Shu, F. H. 1990, ApJ, 361, 546 In the text NASA ADS  
• Safier, P. N. 1993a, ApJ, 408, 115 In the text NASA ADS  
• Safier, P. N. 1993b, ApJ, 408, 148 In the text NASA ADS  
• Sauty, C., & Tsinganos, K. 1994, A&A, 287, 893 In the text NASA ADS  
• Savage, B. D., & Sembach, K. R. 1996, ARA&A, 34, 279 In the text NASA ADS  
• Schmitt, W., Henning, T., & Mucha, R. 1997, A&A, 325, 569 In the text NASA ADS  
• Schunk, R. W. 1975, Planet. Space Sci., 23, 437 In the text NASA ADS  
• Shang, H., Shu, F. H., & Glassgold, A. E. 1998, ApJL, 493, L91 In the text NASA ADS  
• Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 In the text NASA ADS  
• Shu, F. H., Najita, J., Ostriker, E. C., & Shang, H. 1995, ApJL, 455, L155 In the text NASA ADS  
• Shu, F. H., Najita, J., Ostriker, E. C., & Shang, H. 1996, ApJL, 459, L43 In the text NASA ADS  
• Solf, J. 1989, in Low Mass Star Formation and Pre-main Sequence Objects, 399 In the text  
• Solf, J., & Böhm, K. H. 1993, ApJL, 410, L31 In the text NASA ADS  
• Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971 In the text NASA ADS  
• Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722 In the text NASA ADS  
• Vlahakis, N., Tsinganos, K., Sauty, C., & Trussoni, E. 2000, MNRAS, 318, 417 In the text NASA ADS  
• Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218 In the text NASA ADS  
• Webber, W. R. 1998, ApJ, 506, 329 In the text NASA ADS