4 The physical structure of the envelope

The envelope parameters, such as the density and temperature structures, are constrained in this Section using mainly the observed continuum emission.

4.1 Simple analysis

For optically thin dust emission one can derive the following expression for its intensity, in the Rayleigh-Jeans limit, as a function of the impact parameter b (Shirley et al. 2000)

$$\frac{I_{\nu}(b)}{I_{\nu}(0)} = \left (\frac{b}{b_0} \right )^{-\gamma} , \gamma=\alpha+\beta-1,$$ (5)

assuming that the density

$$n_{\rm d}(r) \propto r^{-\alpha}$$ (6)

and temperature

$$T_{\rm d}(r) \propto r^{-\beta}$$ (7)

follow power-law distributions. For IRAS 16293-2422 $\gamma \sim1.8$ is derived in the range $b \sim 15 \arcsec$-50$\arcsec$ for the 450 $\mu$m data. Assuming $\beta=0.4$ (cf., Doty & Leung 1994; Shirley et al. 2000; see also Fig. 5), $\alpha \sim2.4$ is obtained, consistent with the range of values determined by Shirley et al. for a sample of protostellar objects. As discussed by Shirley et al. there are a number of caveats when using this simplistic approach, e.g., the validity of applying the Rayleigh-Jeans approximation in the cool outer parts of the envelope and deviations of the dust temperature from a single power-law (see also Hogerheijde & Sandell 2000; Doty & Palotti 2002). Shirley et al. estimate the uncertainties in the derived $\alpha$ to be of the order $\pm$0.5 using this simple approach.

   
4.2 Power-law density model

In order to derive more reliable envelope parameters the dust radiative transfer model presented in Sect. 3.1 is used. The density of the dust (and gas) is assumed to follow a simple power-law

$$n_{\rm H_2}(r) = n_0 \left ( \frac{1000~{\rm AU}}{r} \right ) ^{\alpha},$$ (8)

where n₀ is the number density of H₂ at a distance of 1000 AU from the star. The models where the density structure is given by a single power-law will be referred to as static envelopes, since in the molecular excitation analysis no large-scale velocity field is included. The total amount of dust and its spatial distribution are simultaneously determined, i.e., $\tau _{100}$ (the optical depth at 100 $\mu$m), $\alpha$, and $r_{\rm e}/r_{\rm i}$ (the geometrical envelope thickness) are the adjustable parameters in our model. Scaling the model output to the source luminosity of 27 $L_{\odot}$ and a distance of 160 pc fixes the absolute scale and allows various physical parameters of the envelope around IRAS 16293-2422 to be determined.

To find the best fit model in the 3D parameter space the strategy adopted is to first analyze the model grid by applying the observed radial brightness distributions. As shown in Fig. 2, the 450 $\mu$m and 850 $\mu$m brightness distributions are sensitive to the slope of the density distribution and, to a lesser extent, the size of the dusty envelope. Changing the total amount of dust by changing $\tau _{100}$ has only a minor effect on the allowed values of $\alpha$ since normalized radial brightness distributions are used. Values of $\alpha$ in the range 1.5-1.9 are found to be acceptable, with a preferred value of 1.7. A typical accuracy of $\pm$0.2 in the derived value of $\alpha$ was also found by Jørgensen et al. (2002) when analyzing SCUBA images at 450 $\mu$m and 850 $\mu$m for a large sample of protostars. In general, the 450 $\mu$m data should provide a more reliable value of $\alpha$ because of the higher resolution, which provides a larger sensitivity to changes in the density structure.

The SED provides a good constraint on $\tau$ once the density profile is known (Fig. 2; see also Doty & Palotti 2002). For a density slope of 1.7 the optical depth at 100 $\mu$m is estimated to be approximately 4.5. At 450 $\mu$m and 850 $\mu$m the optical depths are $\sim$0.4 and $\sim$0.1 respectively. For $\alpha \geq1.7$ the size of the envelope is not well constrained which is not surprising given that only points in the brightness distributions out to 50$\arcsec$ ( $r_{\rm e}/r_{\rm i} \sim250$) are used. The circumstellar envelope will eventually merge with the more extended cloud material in which the object is embedded. The maximum outer radius of the envelope is fixed at the point where the dust temperature reaches 10 K. The envelope size $r_{\rm e}/r_{\rm i}$ is estimated to be 250 for the adopted $\alpha$ of 1.7.

Figure 2: $\chi ^2$ maps showing the sensitivity of the single power-law model to the adjustable parameters $\tau _{100}$, $\alpha$, and $\Delta R = r_{\rm e}/r_{\rm i}$ using the SED and the SCUBA maps as observational constraints. Contours are drawn at $\chi ^2_{{\rm min}}$+(2.3, 4.6, 6.2) indicating the 68% ("1$\sigma$''), 90%, and 95% ("2$\sigma$'') confidence levels, respectively.

The quality of the best fit model can be judged from the reduced $\chi ^2$ obtained from

$$\chi^2_{{\rm red}} = \frac{\chi^2_{{\rm min}}}{N-p},$$ (9)

where p is the number of adjustable parameters. It is found that the best fit model presented above has $\chi^2_{{\rm red}} =0.1$ when the SED is used as constraint and $\chi^2_{{\rm red}} =0.4$ for the 450 $\mu$m and 850 $\mu$m radial brightness distributions, respectively. In all, the observations are well reproduced. The best fit solution is presented in Table 4. The inner radius, $r_{{\rm i}}$, of the envelope is calculated to be $4.8\times10^{14}$ cm (32 AU) fixing the outer radius to $1.2\times10^{17}$ cm (8000 AU or 50$\arcsec$). Mundy et al. (1990) estimated the radial size of the envelope to be $\sim$4400 AU from observations of ammonia emission. The best fit model is presented in Fig. 3 overplotted on the observational constraints. The total mass of molecular hydrogen contained within the outer radius of 50$\arcsec$ is estimated to be 5.4 $M_{\odot}$. Blake et al. (1994) derived a H₂ mass content of $\sim$0.5-0.75 $M_{\odot}$ within a 10$\arcsec$ radius from an excitation analysis of the observed C¹⁷O molecular line emission assuming its abundance to be $3.8\times10^{-8}$. The present analysis gives a mass of $\sim$0.7 $M_{\odot}$ within the same radius, in excellent agreement. Thus, IRAS 16293-2422 appears to have one of the most massive envelopes of the known class 0 protostars (André et al. 2000; Jørgensen et al. 2002).

Although disk emission is typically only responsible for a small fraction of the total flux at sub-millimetre wavelengths it can contribute to the fluxes of the innermost points on the brightness profiles leading to a steeper inferred density profile. Tests in which the flux within a radius of one beam was reduced by 50% indicate that the best fit value of $\alpha$ is reduced by 0.1-0.2.

The temperature and density structures obtained from the best fit model are presented in Fig. 5. For comparison, the predicted temperature structure based upon an optically thin approximation (Chandler & Richer 2000) and scaled to the luminosity of IRAS 16293-2422 is shown for two different opacity laws. Both predict the dust temperature to follow a single power-law. Clearly, the temperature structure obtained from the detailed radiative transfer analysis is not well described by a power-law and has a significantly steeper gradient in the inner parts of the envelope where optical depth effects are important. The interstellar radiation field is potentially important for the temperature structure in the outer parts of the envelope. However, detailed modelling (S. Doty, priv. comm.) shows, for the envelope around IRAS 16293-2422, that this effect is small when assuming a typical interstellar radiation field. The difference in the dust temperature is $\sim$10-20% (a few K) at $r \ga 8\times10^{16}$ cm.

Figure 3: Best fit model, using a density structure described by a single power-law, compared with observed radial brightness distributions and the SED. In the $\chi ^2$-analysis of the radial brightness distributions only the data points separated by one full beam width were used, shown here with error bars. The error bars represent the rms scatter of the observations within each of the bins and are a combination of the noise as well as gradients and departure from spherical symmetry in the brightness map. The dotted line is the azimuthally averaged SCUBA beam at the time of the observation. In the SED panel the observations, represented by circles with error bars (open circles were not used in the $\chi ^2$-analysis), are overlayed with the output from the best fit model.

The dust opacities adopted in the previous analysis might not be appropriate in the inner hot parts of the envelope where the ices start to evaporate off the grains. Instead opacities for bare grains, without any ice mantles, should preferably be used. DUSTY is not set up to allow dust opacities with radial dependence so only the limiting cases with or without ice mantles can be compared. To test the sensitivity of the derived envelope parameters on the adopted set of dust opacities, the analysis is repeated using the bare grain opacities presented by Ossenkopf & Henning (1994) (Col. 2 in their Table 1; hereafter OH2). From the radial brightness distributions the same range of $\alpha$ and $\Delta R$ as when using OH5 is obtained. The optical depth, however, is significantly reduced by about a factor of two, and so is the total mass of the envelope. Adopting the model parameters as derived in Table 4, i.e., $\tau_{100} = 4.5$, $\alpha =1.7$, and $\Delta R =250$, but using the OH2 opacities increases the flux at 60 $\mu$m by $\sim$20% compared with the OH5 model. The 100 $\mu$m flux is not significantly affected whereas the 1.3 mm flux is roughly twice that of the OH5 flux. A model with dust properties varying with radius (OH2 in the inner warm parts and OH5 in the outer cool parts) would thus serve to improve the fit to the SED. Since the bulk of the mass (99.8%) is at low temperatures where the grains are coated with ice mantles, the model parameters obtained with the OH5 opacities were used in the chemical analysis.

Figure 4: Best fit CO and CS models using a static envelope (dotted line) and a Shu-type collapsing model (full line), overlayed on the observed spectra (histograms) (see text for details). The calibration uncertainty in the intensity scale for the observed spectra is $\sim$15%.

In addition to the analysis of the continuum emission the observed molecular line emission is useful in constraining the physical properties of the envelope. Traditionally, CO and CS line emission have been extensively used for this purpose and are adopted here to test the validity of the best fit model obtained from the dust analysis. Using the radiative transfer code presented in Sect. 3.2, the total CO and CS abundances relative to H₂ obtained for the best fit model presented in Table 4 are ${\sim}4\times10^{-5}$ and ${\sim}3\times10^{-9}$, respectively. In the $\chi ^2$-analysis only the velocity-integrated intensities in the lines were used. These values are within a factor of about two of what is commonly derived for YSOs (van der Tak et al. 2000b). The abundances in combination with the quality of the fits (Fig. 4; see also Sect. 5), in particular the ratios among various transitions which are sensitive to the gas temperature and density, are reassuring and further strengthen the adopted physical model.

In the modelling of the molecular line emission the gas temperature is assumed to follow that of the dust. In models which self-consistently treat the energy balance the gas temperature is generally lower than that of the dust in the outer regions due to imperfect gas-grain coupling (Ceccarelli et al. 1996; Doty & Neufeld 1997; Ceccarelli et al. 2000a). To test the effects of a departure of the gas temperature from that of the dust due to gas-grain decoupling in the outer regions, the dust temperature was scaled by a constant factor. For $T_{{\rm gas}} \la0.7\times T_{\rm dust}$ the envelope becomes too cool to fit the observed line intensity ratios. Thus, the gas temperature appears to follow that of the dust within $\sim$30% in the region probed by the CO emission.

The static envelope model fails to explain the details of the individual spectra and the potential of the molecular emission to constrain the velocity fields will be investigated further in Sect. 4.3. The derived abundances are not very sensitive to the adopted value of $\alpha$, within the limits derived from the dust radiative transfer model. Similarly, increasing the outer radius by a factor of two only marginally affects the line intensities. The abundances derived for a wide variety of molecular species are further presented in Sect. 5.

Thus, the emerging picture from the analysis is that the envelope around IRAS 16293-2422 indeed has a region of dense and hot gas inside a radius of ${\sim}2\times10^{15}$ cm (150 AU, 1$\arcsec$), with temperatures decreasing to $\sim$10 K at $\sim$10¹⁷ cm (8000 AU).

Figure 5: Properties of the circumstellar envelope around IRAS 16293-2422 obtained from the dust modelling. Shown are the dust temperature (left) and density (middle) structures for best fit model (full line) assuming a single power-law distribution of the density. In the temperature panel, an optically thin prediction using an opacity law $\kappa _{\nu } \propto \nu ^{\beta }$ with $\beta = 1$ (dotted line; $T_{\rm d} \propto r^{-0.4}$) and $\beta = 2$ (dashed line; $T_{\rm d} \propto r^{-0.33}$), is also shown. In the density panel the best fit Shu-type collapsing envelope model is also shown (dash-dotted line). The velocity structure obtained from the Shu infall model is shown in the right panel.

Figure 6: $\chi ^2$-maps showing the sensitivity of the Shu-collapse model to the adjustable parameters in the modeling of the dust continuum emission (left) and the CS (CS and C34S transitions; centre) and CO (13CO, C17O, and C18O transitions; right) line emission compared to observations. Contours are drawn at $\chi^2_{{\rm min}}+(2.3,~4.6,~6.2)$ indicating the 68% ("1$\sigma$''), 90%, and 95% ("2$\sigma$'') confidence levels, respectively. The number of observational constraints used, N, are also shown. The quality of the best fit model can be estimated from the reduced chi-squared statistic $\chi ^{2}_{{\rm red}} = \chi ^{2}_{{\rm min}}/(N$-2).

   
4.3 Infall model

Current theories of star formation state that protostars are formed from the gravitational collapse of cloud cores consisting of gas and dust (e.g., Shu et al. 1987). There is now growing evidence that low-mass protostars have parts of their circumstellar material in a state of collapse (e.g., Myers et al. 2000) and an obvious extension to the previous analysis, which uses a static envelope with a density distribution described by a single power-law, is to attempt to reproduce the results with an infall model. That the circumstellar envelope around IRAS 16293-2422 is in a state of collapse has been inferred previously from modelling of millimetre CS observations (Walker et al. 1986; Zhou 1995; Narayanan et al. 1998), although such conclusions based on low spatial resolution data have been questioned by Menten et al. (1987). Furthermore, Ceccarelli et al. (2000a,b) have suggested an infall model based upon a physical-chemical model. Here, the observed continuum and molecular line emission will be analyzed using the well known collapse model presented by Shu (1977).

In the Shu inside-out collapse model the self-similar solution is presented in terms of the dimensionless variable x = r/at, where ris the radial distance scale, and characterized by the isothermal speed of sound, a, and the time after onset of collapse, t. The location of the collapsing wave front at any instant t is described by $r_{\rm c} = at$. The density $\rho$ and velocity u have the form

$$\rho(r,t) = \frac{\alpha(x)}{4\pi Gt^2},\ u(r,t) = a v(x),$$ (10)

where G is the gravitational constant and $\alpha(x)$ and v(x) are tabulated by Shu (1977). In the static part of the envelope (x >1)

$$\alpha(x)=\frac{2}{x^2}, \ v(x)=0.$$ (11)

The asymptotic behavior as $x\rightarrow 0$ of the solution is

$$\alpha(x)=\left( \frac{m_0}{2x^3} \right)^{1/2},\ v(x)=-\left (\frac{2m_0}{x} \right)^{1/2},$$ (12)

where m₀ is equal to 0.975 for this particular solution. Given the limited resolution provided by the JCMT at 450 and 850 microns it is not surprising that the observations are well described by a single power-law with $\alpha$ in the range $\sim$1.5-2.0 which are the two extremes obtained from the Shu-model. The estimated value of $\alpha$ will depend on the location of the collapsing wavefront, $r_{\rm c}$. In comparison, Jørgensen et al. (2002) derived values of $\alpha$ for a large number of class 0 protostars and found values typically in the range $\sim1.3$-2.0 .

The input parameters to DUSTY are the same as for the single power-law models (see Table 4). In addition, the envelope size was fixed to a radius 5000 AU. Making the envelope larger will produce increasingly worse fits to the radial brightness distributions obtained from the SCUBA observations. The sensitivity of the two adjustable parameters a and $r_{\rm c} = at$ in the modelling is shown in Fig. 6 (left panel) where the observational constraints used are the SED and the SCUBA 450 $\mu$m radial brightness distribution. The best fit model is obtained using $a \sim0.6$ km s^-1 and $r_{\rm c} \sim 2 \times 10 ^{16}$ cm, putting the age at $\sim$10⁴ yr.

The mass accretion rate can be estimated from

$$\dot{M}=\frac{m_0a^3}{G},$$ (13)

where G is the gravitational constant. For the best fit model $\dot{M} \sim 5\times10^{-5}$ $M_{\odot}$ yr^-1 is obtained. Such a high accretion rate is able to account for the source luminosity. The total mass contained in the envelope within 5000 AU is $\sim$3.3 $M_{\odot}$, consistent with the estimate obtained from the static envelope model within the same radius.

The relative success of the dust modelling using a static envelope, with a single power-law to describe the density structure, makes it hard to discriminate between the two models in the present analysis. In Fig. 5 the density and velocity structures obtained from the best fit Shu model are presented and, for the density, compared with results from the static envelope model. The largest discrepancy occurs at small radial distances where the collapsing envelope model predicts about a factor of two to three lower densities. We stress the observational data set used in the dust modelling is not directly probing this region. At larger radii the model is better constrained and the two models agree well. It should be noted that the best fit single power-law model gives slightly better reduced $\chi ^2$ values for the combined set of observations. However, the molecular data provide further constraints since they have the potential to probe the large scale velocity field.

In Fig. 4, spectra of CO and CS line emission as observed with the JCMT are presented. Lines which are optically thick, like CS, ¹³CO and C¹⁸O ( $J = 3\rightarrow 2$), show a distinct, narrow, absorption feature near the stellar velocity. This feature is due to effective self-absorption in the outer cool parts of the envelope. Also, the degree of the asymmetry in the line profiles increases with the optical depth in the lines. In the optically thin lines the self-absorption feature disappears and the lines are well described by a single Gaussian profile. The width of the self-absorption feature constrains the turbulent velocity to $\sim$0.3 km s^-1 in the outer envelope. For simplicity this value is adopted throughout the envelope. A turbulent velocity component varying with radius is beyond the scope of this article, see however Stark et al. (2002, in prep.).

The observed CS emission (including that from C³⁴S) is analyzed using the abundances derived from the static envelope model presented in Sect. 4.2 (see also Sect. 5). From only the integrated intensities it is possible to constrain both a and $r_{\rm c}$ as is shown in Fig. 6 (middle panel). The model spectra are presented in Fig. 4 together with the observations. The fit to the integrated intensities is worse than that obtained from the single power-law model presented in Sect. 4.2, however. In particular the C³⁴S ( $J =7\rightarrow 6$) line is poorly reproduced in the infall model. Using also the line profiles as constraints we find that a model where the collapsing wavefront is located at $\sim$ $1.0\times10^{17}$ cm and the value of a is $\sim$0.9 km s^-1 best reproduces the observations.

Analyzing the CO emission (¹³CO, C¹⁸O, and C¹⁷O) gives yet another set of estimates. At first the abundances derived from the static envelope models presented in Sect. 4.2 (see also Sect. 5) are adopted. From the analysis of the integrated intensities shown in Fig. 6 (right panel) and the line profiles presented in Fig. 4 a Shu-model with $a \sim1.0$ km s^-1 and $r_{\rm c} \sim 1.2 \times 10 ^{17}$ cm reproduces the CO observations well, in excellent agreement with the CS modelling. The derived dynamical age of the system is ${\sim}3\times10^{4}$ yr. In contrast with the CS modelling, the fit to the integrated CO intensities is equally good as obtained for the single power-law model. As discussed in Sect. 5, the CO abundance obtained from the static envelope model is about a factor of 2-3 lower than what is typically observed for interstellar gas. If instead the CO abundance relative to H₂ is assumed to be the "standard'' interstellar value of $1\times10^{-4}$and the standard isotopic ratios are assumed ([CO/¹³CO]  =60, [CO/C¹⁷O]  =2500, and [C¹⁸O/C¹⁷O]  =3.9) the estimate of a is $\sim$0.75 and $r_{\rm c} \sim 3 \times 10 ^{17}$ cm. However, the quality of the fit becomes worse in this case.

The results obtained from the best fit Shu-model are presented in Fig. 4 overlayed onto the observed CO and CS spectra. The integrated intensities and, to some extent, line profiles can be modelled with the spherically symmetric infall solution. The details of the spectra will, however, intricately depend on the adopted geometry, velocity fields, and chemical gradients. The velocity field, in particular position-velocity maps of the source, form a stringent test of the dynamical models. For IRAS 16293-2422 it will likely be necessary to include a rotational component to the velocity field (Menten et al. 1987; Zhou 1995; Narayanan et al. 1998).

Other estimates of the infall radius, based on analysis of molecular line emission, range between about 5 and $15\times10^{16}$ cm (Walker et al. 1986; Zhou 1995; Narayanan et al. 1998; Ceccarelli et al. 2000a) in excellent agreement with the values obtained here from analysis of CO and CS emission. However, there appears to be a discrepancy between the dust and molecular line analysis, possibly reflecting the fact that the simple Shu-collapsing core model is not fully adequate to describe the state of the infalling material and/or that some of the CO and CS emission is associated with the outflow and surrounding cloud. Moreover, gradients in the molecular abundances will affect the parameters derived.

  Table 5: Derived abundances using a constant molecular abundance f_X relative to H₂ throughout the envelope.

$^{\rm a}$ The reduced $\chi^2$ of the best fit model. Based on the $\chi^2$-analysis the derived abundances are accurate to about 20-30% (1$\sigma$) when the fit is good, i.e., $\chi^2_{{\rm red}} \sim1$.
$^{\rm b}$ The number N of independent observational constraints used in the modelling.
$^{\rm c}$ Abundance derived by Blake et al. (1994) and van Dishoeck et al. (1995) using a simple excitation analysis.
$^{\rm d}$ Lowest energy of the upper level involved in the transitions used as constraints in the modelling.
$^{\rm e}$ Highest energy of the upper level involved in the transitions used as constraints in the modelling.
${\dagger}$ Derived from optically thin isotope(s) assuming standard isotopic ratio(s).