3 Line modeling

   
3.1 Origin of the FIR line emission

The molecular emission (H₂O, CO and OH) observed toward IRAS 4 can have at least three different origins: the two outflows powered by IRAS 4A and IRAS 4B, the PDR at the surface of the cloud, and the collapsing envelopes around the two protostars. The origin of the molecular line emission can be disentangled when the spatial distribution of the line emission and/or the line profiles are available. For example, lines arising in the envelope or in the molecular cloud have narrow profiles whereas lines arising from outflows show broadened profiles with extended wings. Unfortunately, in the case of the ISO-LWS observations, the relatively low spatial and spectral resolution do not allow to observationally disentangle the different components. However, the comparison between the central and the NE-red and SW-blue positions allows a first guess of the origin of the observed emission. The [OI] and [CII] lines have comparable line fluxes in the three observed positions. For this reason it is likely that the observed [OI] and [CII] emission is associated with the ambient diffuse gas, either emitted in the PDR or in the molecular cloud itself. On the contrary, only the lowest lying ( $J_{\rm up} \leq 17$) CO lines are detected on the NE-red position, while no H₂O or CO emission is clearly detected on the SW-blue position. On the other hand, the observations of the CO 3-2 line show that the high velocity gas (the fastest outflow component) peaks at the NE-red and SW-blue positions (Blake et al. 1995). The lack of water emission in these two outflow peak positions is not in favor of the hypothesis that the on-source water emission originates in the outflow. Although we cannot exclude a different origin and/or contamination for example from the densest parts of the outflow located in the ISO beam, in the following we explore the hypothesis of the envelope thermal emission and interpret the observed water line emission according to the CHT96 model. The first goal of our modeling is to verify that the thermal emission from the surrounding envelope can reproduce the water line observations, a necessary condition even though not sufficient to test this hypothesis. A following section will then address the possible origin of the CO and OH observed emission.

3.2 Model description

The CHT96 model computes in a self-consistent way the radiative transfer, thermal balance, and chemistry of the main gas coolants (i.e. O, CO and H₂O) across the envelope, in the inside-out framework (Shu 1977). Here we give a brief description of the main aspects of the model.

The initial state of the envelope is assumed to be an isothermal sphere in hydrostatic equilibrium, which density is given by:

$$n_{\rm H_{2}}(r) = \frac{a^{2}}{2 \pi \mu {\rm m}_{{\rm H}} G} r^{-2}$$ (1)

where a is the sound speed, $m_{{\rm H}}$ is the hydrogen mass, $\mu$ is the mean molecular mass, r the distance from the center and G the gravitational constant.

At t = 0 the equilibrium is perturbed and the collapse starts from inside out, propagating with the sound speed. The density in the inner collapsing region is given by the free-fall solution:

$$n_{{\rm H_{2}}}(r) = \frac{1}{4 \pi \mu m_{{\rm H}}} \frac{\dot{M}} {(2GM_{*})^{1/2}} r^{-3/2}$$ (2)

where M_* is the star mass, $\dot{M}$ is the accretion rate, related to the sound speed by:

$$\dot{M} = 0.975 \frac{a ^{3}}{G}\cdot$$ (3)

The free-fall velocity is given by:

$$v(r) = \left( \frac{2GM_{*}}{r} \right) ^{1/2}.$$ (4)

The gravitational energy is released as material falls at the core radius R_*, so that the luminosity of the protostar is:

$$L_{*} = \frac{GM_{*}\dot{M}}{R_{*}}\cdot$$ (5)

In the following L_* is the bolometric luminosity of IRAS 4, and we leave M_* and $\dot{M}$ as free parameters.

The radiative transfer in the envelope is solved in the escape probability approximation in presence of warm dust, following the Takahashi et al. 1983 formalism. The CHT96 model assumes that the initial chemical composition is that of a molecular cloud, and then it solves the time dependent equations for the chemical composition of 44 species, as the collapse evolves. H₂O, CO and O are of particular importance since they are the main coolants of the gas, and hence we study the chemistry of these species in detail. The CO molecule is very stable, and its abundance results constant across the envelope. H₂O is mainly formed by dissociative recombination of the H₃O⁺ in the cold outer envelope, while, at dust temperature above 100 K, icy grain mantles evaporate, injecting large amounts of water into the gas phase. When the gas temperature exceeds $\sim$250 K, the H₂O formation is dominated by the endothermic reactions O + H$_{2} \to$ OH + H followed by H₂ + OH $\to$ H₂O + H, which transform all the oxygen not locked in CO molecules into H₂O.

From the above equations and comments, the water line emission depends on the mass of the central object, the accretion rate and the abundance of H₂O in the outer envelope and in the warm region, where its abundance is dominated by the mantle evaporation. All these quantities directly enter into the H₂O line emission, and specifically into the determination of the H₂O column density. In fact, the accretion rate sets the density across the protostellar region (Eqs. (1) and (2)). The central mass of the protostar affects the velocity field, and hence indirectly the line opacity (Eq. (4)). This parameter also sets the density in the free-fall region (Eq. (2)) and therefore the gas column density in this region. The water emission also depends indirectly on the O and CO abundances, which enter in the thermal balance and hence in the gas temperature determination. Several recent studies (Baluteau et al. 1997; Caux et al. 1999b; Vastel et al. 2000; Lis et al. 2001) have shown that almost the totality of the oxygen in molecular clouds is in atomic form. Accordingly, we assume the oxygen abundance to be the standard interstellar value, i.e. $5 \times 10^{-4}$ with respect to H₂. With regard to the CO abundance, following Blake et al. 1995 we adopt a CO abundance of 10^-5 with respect to H₂, lower than the standard abundance as this molecule is believed to be depleted on IRAS 4. We will comment later on the influence of these parameters on our results. Finally, the water abundance in the cold and in the warm parts of the envelope are poorly known and are free parameters in our study.

To summarize, we applied the CHT96 model to IRAS 4, and to reproduce the observations we varied the four following parameters: the mass of the central object M_*, the accretion rate $\dot{M}$, the water abundance in the outer cold envelope  $X_{\rm out}$, and the water abundance in the region of mantle evaporation  $X_{\rm in}$. The principal limitation to the application of this model to IRAS 4 is that the ISO beam includes both IRAS 4A and IRAS 4B envelopes. As a first approximation, we assumed that the two envelopes contribute equally to the molecular emission. Finally we assumed that the two envelopes touch each other, namely they have a radius of 3000 AU (i.e. 30'' in diameter), in agreement with millimeter continuum observations (Motte & André 2001). In our computations we adopted a distance of 220 pc in agreement with JSD02 and a luminosity of 5.5  ${L}_{\odot}$ for each protostar, according to Sandell et al. (1991) when assuming such a distance.

3.3 Water line modeling results

In order to constrain the mass, accretion rate and water abundance across the envelope, we run several models varying the central mass from 0.3 to 0.8  ${M}_{\odot }$, the accretion rate from 10^-5 to 10^-4  ${M}_{\odot }$ $~{\rm yr}^{-1}$, the water abundance in the outer parts of the envelope  $X_{\rm out}$ between 10^-7 and 10^-6, and a water abundance in the inner parts of the envelope  $X_{\rm in}$ between 10^-6 and $2 \times 10^{-5}$ respectively. In the following we discuss the results of this modeling.

3.3.1 Accretion rate and water abundance in the outer parts of the envelope

One of the difficulties in constraining the central mass, accretion rate and water abundance in the envelope is that the water line intensity a priori depends on all the parameters. However, choosing appropriate lines can help constraining one parameter at once. Low lying H₂O lines are expected to rapidly become optically thick in the outer envelope, where they are easily excited. Hence these lines depend weakly on  $X_{\rm in}$ and M_* (which affect the line emission in the collapsing inner region). We therefore used the low-lying lines to constrain the other two parameters, namely $X_{\rm out}$ and $\dot{M}$. For this we minimized the $\chi^2=\frac{1}{N-1}\sum_{1}^{N}\frac{\left(\rm Observations - Model\right)^{2}}{\sigma^{2}}$ obtained considering only water lines having a  ${E_{\rm up}}$ lower than 142 cm^-1, and where $\sigma$ is the error associated with each line flux. Figure 2 shows the $\chi ^2$ surface as a function of $X_{\rm out}$ and $\dot{M}$ for M_* = 0.5  ${M}_{\odot }$. We obtained a similar plot for M_* = 0.3  ${M}_{\odot }$ and the result is basically the same. As suspected, the chosen lines constrain relatively efficiently the two parameters $X_{\rm out}$ and $\dot{M}$. The minimum $\chi ^2$ is obtained for a water abundance $X_{\rm out}$  $\sim 5 \times 10^{-7}$ and an accretion rate $\dot{M}$  $\sim 5 \times 10^{-5}$ ${M}_{\odot }$ $~{\rm yr}^{-1}$.

Figure 2: $\chi ^2$ surface as function of the water abundance in the outer envelope $X_{\rm out}$ and the mass accretion rate $\dot{M}$. The $\chi ^2$ has been obtained considering the lines with an upper level energy ${E_{\rm up}}$ lower than 142 cm-1 and for a central mass of 0.5  ${M}_{\odot }$.

3.3.2 Central mass and water abundance in the inner parts of the envelope

We then constrained the central mass M_* and the abundance in the innermost parts of the envelope  $X_{\rm in}$ using the high-lying lines. In fact to be excited, these lines require relatively high temperatures and densities which are likely to be reached only in the innermost parts of the envelope. Their intensities depend hence on the water abundance $X_{\rm in}$ in these parts and on the central mass M_*. The $\chi ^2$ surface as function of these two parameters is shown in Fig. 3, obtained considering the lines having an upper level energy ${E_{\rm up}}$ larger than 142 cm^-1, and assuming $\dot{M}$ = $5 \times 10^{-5}$ and $X_{\rm out}$ = $5 \times 10^{-7}$. In this case the $\chi ^2$has a minimum around M_* = 0.5  ${M}_{\odot }$ and  $X_{\rm in}$ = $5 \times 10^{-6}$. Specifically, if we adopt a constant H₂O abundance of $5 \times 10^{-7}$ ( $X_{\rm out}$) across the entire envelope, the model predicts intensities a factor between two and five lower than those observed (high lying lines, i.e. with  ${E_{\rm up}}$ $\geq 250$ K). In other words, the observed emission can only be accounted for if a jump in the water abundance is introduced when the dust temperature exceeds 100 K, the sublimation temperature of icy mantles. This jump has to be larger than about a factor 10.

Figure 3: $\chi ^2$ surface as function of the central mass M* and water abundance in the innermost parts of the envelope $X_{\rm in}$. The $\chi ^2$ has been obtained considering the lines with an upper level energy ${E_{\rm up}}$ larger than 142 cm-1 and assuming $\dot{M}$ = $5 \times 10^{-5}$ and $X_{\rm out}$ = $5 \times 10^{-7}$. Note that we did not include the 82 and 99.5 $\mu$m lines, which seems underestimated by our model (see text). We did not include the 113 $\mu$m line either, because of the blending with the CO $J_{\rm up} = 23$ line, which makes the estimate of the flux rather uncertain.

3.3.3 The best fit model

Assuming two identical envelopes, the best fit model is obtained with a central mass of 0.5  ${M}_{\odot }$, accreting at $5 \times 10^{-5}$ ${M}_{\odot }$ $~{\rm yr}^{-1}$ (Table 2).

 

 

Table 2: Best fit parameters and derived quantities for IRAS 4 and comparison with the values obtained towards IRAS 16293-2422 by Ceccarelli et al. (2000a).

Parameter	IRAS 4	IRAS 16293-2422
Mass ( ${M}_{\odot }$)	0.5	0.8
Accretion rate ( ${M}_{\odot }$ $~{\rm yr}^{-1}$)	$5 \times 10^{-5}$	$3 \times 10^{-5}$
Water abundance $X_{\rm out}$	$5 \times 10^{-7}$	$3 \times 10^{-7}$
Water abundance $X_{\rm in}$	$5 \times 10^{-6}$	$3 \times 10^{-6}$
Radius ( $T_{\rm dust} = 30$ K) (AU)	1500	4000
Radius ( $T_{\rm dust} = 100$ K) (AU)	80	150
Age (yr)	$1.0 \times 10^4$	$2.7 \times 10^4$

Assuming a constant accretion rate, this gives an age of 10 000 years, close to the dynamical age of the outflows. The abundance of water in the outer parts of the envelope is $5 \times 10^{-7}$ and it is enhanced by a factor 10 in the innermost regions of the envelope, where grain mantles evaporate. Figure 4 shows the ratio between the observed and predicted line fluxes as function of the upper level energy of the transition. The model reproduces reasonably well the observed water emission, with the exception of the lines at 99.5 $\mu$m and 82.0 $\mu$m that seems to be underestimated (by a factor 10) by our model.Since, on the contrary, lines in a comparable range of  ${E_{\rm up}}$ are well reproduced by our model, we think that this discrepancy is likely due to a rough baseline removal. Specifically, the estimate of both the 99.5 and 82.0 $\mu$m line fluxes may suffer of an incorrect baseline removal, as the lines lie on the top of a strong dust feature, which covers the 80-100 $\mu$m range (Ceccarelli et al. in preparation). Higher spectral resolution observations are required to confirm this explanation. Finally some unexplained discrepancy between the model and the observations may exist at the higher values of  ${E_{\rm up}}$.

Figure 4: Ratio between the line fluxes predicted by our best fit model and observed ones as function of the upper level energy of the transition ${E_{\rm up}}$. Triangles represent ortho H2O transitions, and squares represent para H2O transitions. Note that the model assumes an ortho to para ratio equal to 3.

In the figure we also report different symbols for the ortho and para water transitions respectively. In our model we assumed that this ratio is equal to 3. The comparison between the observations and predictions is consistent with this assumption. Plots of the predicted intensity of various lines as function of the radius are reported in Fig. 5. Finally, the envelope model predicts an intensity of $1.8 \times 10^{-13}$ erg s^-1 cm^-2 and a linewidth of $\sim$1 km s^-1 for the ground water line at 557 GHz, equivalent to an antenna temperature of 30 mK in the SWAS beam ($\sim$4'). SWAS detected a $T_{\rm a}^* \sim 100$ mK line, which linewidth is $\sim$18 km s^-1, and self-absorbed at the rest velocity (Bergin et al. 2002). The observed line is undoubtedly dominated by the outflow emission, and only a small fraction of its intensity can be attributed to the envelope emission, in agreement with our model predictions.

Figure 5: Predicted intensity of various lines as function of the radius. The upper panel shows the water lines at 179 $\mu$m (solid line), 108 $\mu$m (dotted line), 75 $\mu$m (dashed line) and 82 $\mu$m (dash-dotted line). The lower panel shows the CO $J_{\rm up} = 3$ (solid line), CO $J_{\rm up} = 14$ (dashed line) and C18O $J_{\rm up} = 3$ (dash-dotted line).

3.4 [OI], [CII], OH and CO emission

[OI] 63 $\mu$m and [CII] 157 $\mu$m emission is widespread, and probably associated with the cloud. A plausible explanation is that the two lines are emitted in the PDR resulting from the UV and/or X-ray illumination of this cloud from the several young stars that it harbors. The comparison of the observed fluxes with the model by Kaufman et al. (1999) suggests a PDR with a density of about 10⁴ cm^-3 and a incident FUV of $\sim$5 G₀. This PDR would account for the total observed flux of the [CII] 157 $\mu$m line and OI. The parameters we derive are in agreement with those quoted by Molinari et al. (2000), who studied the region around SVS13.

The thermal emission from the envelope predicts no C⁺ emission, of course, as no source of ionization is considered in the CHT96 model. The atomic oxygen, on the contrary, is present all along the envelope and it is predicted to emit $1.8 \times 10^{-12}$ erg s^-1 cm^-2. This is similar to the observed [OI] 63 $\mu$m flux. The fact that we do not see any [OI] 63 $\mu$m enhancement towards the source with respect to the surroundings can be explained if IRAS 4 is well embedded in the parental cloud. Being the ground transition, the [OI] 63 $\mu$m line is relatively easily optically thick, and an emission from an embedded source can be totally absorbed by the foreground material (Poglitsch et al. 1996; Baluteau et al. 1997; Caux et al. 1999b; Vastel et al. 2000).

Finally, our model predicts OH and CO $J_{\rm up}\geq 14$ line fluxes more than ten times lower than those observed. Note that the FIR CO lines predicted by the CHT96 model are optically thick and not sensitive to the adopted abundance, and therefore increasing the CO abundance would not change the result. An extra heating mechanism is evidently responsible for the excitation of the FIR CO lines observed in the central position. Shocks have been invoked in the literature (e.g. Ceccarelli et al. 1998; Nisini et al. 1999; Giannini et al. 2001), but this hypothesis has its own drawbacks and flaws (see Introduction). Another possibility is that the FIR CO lines are emitted in a superheated layer of gas at the surface of a flaring disk, as seen in the protostar El 29 (Ceccarelli et al. 2002), and/or at the inner interface of the envelope itself (Ceccarelli, Hollenbach, Tielens et al. in preparation). In the first case (disk surface), the gas is "super-heated'' because the grain absorptivity in the visible exceeds the grain emissivity in the infrared (e.g. Chiang & Goldreich 1997). In the latter case (envelope inner interface) the extra heating is provided by the X-ray photons from the central source. A full discussion about the FIR CO emission origin is beyond the scope of this article and we refer the interested reader to the above mentioned articles.