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).