4 Estimation of the CO column densities

Two different approaches were used to estimate the CO column density distribution across the cloud. Firstly, we used the traditional way to derive column densities directly from the observed lines by assuming that the excitation temperature for each transition is constant along the line of sight (Sect. 4.1). Secondly, starting from the physical cloud model (see Sect. 2) we used the Monte Carlo radiative transfer program developed by Juvela (1997) to simulate the observed profiles from the cloud. The physical quantities were derived by fitting the calculated spectra to the observed ones (Sect. 4.2). In the following, the two methods and the results are described in detail.

  
4.1 Line of sight homogeneity

4.1.1 Assumptions

The assumption of a constant excitation temperature implies that either the cloud is homogenous (the density and the kinetic temperature are constant) or that the cloud is isothermal and the transition in question is thermalised. Even if this simplistic assumption is valid, the derivation of the excitation temperatures of the observed isotopic CO lines involves the following difficulties: 1) The $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO/C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O abundance ratio X has been observed to change from cloud to cloud, and may also vary from the cloud surface to its interior parts. This depends on ¹³C fractionation and selective photodissociation (see e.g. Bally & Langer 1982; Smith & Adams 1984). 2) The populations of the rotational levels of C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O may deviate from LTE, i.e. from the situation where a single value of $T_{\rm ex}$describes the relative populations of all J-levels. The populations are controlled by the collisional excitation and selective photodissociation mechanisms, which depend on the rotational quantum number (Warin et al. 1996). According to Warin et al. (1996) the latter process causes that low-lying rotational levels are thermalised or overpopulated and higher levels are subthermally excited. In their model calculation for a dense dark cloud with a kinetic temperature of $T_{\rm kin}=10$ K (see their Figs. 6a-c), the excitation temperatures of both C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=1\mbox{--}0)$ and $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO $(J=1\mbox{--}0)$ lie close to , whereas $(J=2\mbox{--}1)$ and all higher transitions of C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O settle between 6 and 7 K.

In the derivation of the C¹⁸O column densities we have made the following assumptions:

1.

the C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=1\mbox{--}0)$ and $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO $(J=1\mbox{--}0)$ transitions have the same excitation temperature $(J=1\mbox{--}0)$, 

2.

the $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO/C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O abundance ratio X is constant throughout the cloud, and 

3.

all the higher transitions of C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O, i.e. C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=2\mbox{--}1, 3\mbox{--}2, \ldots)$ have the same excitation temperature $(J=2\mbox{--}1)$, which may differ from $(J=1\mbox{--}0)$. 

This method, based on the modelling results of Warin et al. (1996), has been used earlier by Harjunpää (2002). The possibility that Xcan be constant within a cloud is supported by the results of Harjunpää & Mattila (1996) and Anderson et al. (1999) in CrA, where the $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO/C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O and H $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO⁺/HC $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O⁺ abundance ratios have an almost invariable value of 10.

  
4.1.2 Equations and results

Figure 2: $T_{\rm R}^{13}$ vs.  $T_{\rm R}^{18}$ correlation plot. Intensities are taken from Gaussian fits to the spectra at the velocities of the C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O peaks. The solid line is the best fit of Eq. (3), using a non-weighted least-squares fit to all data points with ${\rm SNR}>3$ in both isotopes (asterisks). The curve implies a $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO/C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O abundance ratio X = 11.2, and $\ensuremath{T_{\rm ex}}\mbox{$(J=1\mbox{--}0)$ } = 8.0$ K. The dashed line is the fit to all data points (no SNR criterion, asterisks and crosses), yielding X = 8.2and $\ensuremath{T_{\rm ex}}\mbox{$(J=1\mbox{--}0)$ } = 8.4$ K. Box encompassed data points have velocity offsets >0.3  between the lines of the two isotopes and are not used in these fits. For comparison, the dotted line gives the low-opacity relation for the terrestrial isotopic ratio X=5.5.

Figure 3: C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O column densities (contours) overlaid on the POSS-II red plate (image). The optical image has been smoothed to 10 resolution, and the transfer function is chosen to emphasise qualitatively the decrease in diffuse background brightness towards the centre of B68. The column densities are derived from our 2-transition observations of C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O, assuming a constant $(J=1\mbox{--}0)$ excitation temperature of 8 K (see text and Fig. 2). The angular resolution is 50 . Only observed positions with an uncertainty $\Delta\ensuremath{\!\:\!}\left(N(\mbox{C$\ensuremath{\:\!} ^{18}\ensuremath{\!\... ...emath{\:\!} {}N(\mbox{C$\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!} $ O}) < 30\%$ have been used for this map. Contours are drawn from 1.5 to 5.5 $\ensuremath{\;\!}\times\ensuremath{\;\!}$1014  in steps of 0.5 $\ensuremath{\;\!}\times\ensuremath{\;\!}$1014  , using alternating line thicknesses. The outer contours closely follow the optical outline of the globule. In the centre, the CO map is rather flat, while the diffuse light suggests an even further density increase. This behaviour is most likely to be due to CO depletion (see also Fig. 4).

The assumptions 1 and 2 lead to the following formula for the radiation temperatures at a certain velocity:

$$\frac{T_{\rm R}^{13}}{T_\infty} \; = \; 1 - \left[1 - \frac{T_{\rm R}^{18}}{T_\infty} \right]^X \; ,$$ (3)

where $T_{\rm R}^{13}$ and $T_{\rm R}^{18}$ are the radiation temperatures of the $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO $(J=1\mbox{--}0)$ and C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=1\mbox{--}0)$ lines, respectively; $T_\infty$ is defined by

$$T_\infty \equiv \frac{h\nu}{k} ~ \left[ f_\nu(T_{\rm ex}) - f_\nu(T_{\rm bg}) \right] \; ,$$ (4)

and the function $f_\nu(T)$ is defined by

$$f_\nu(T) \equiv \frac{1}{{\rm e}^{h\nu/kT} - 1} \; ,$$ (5)

where $\nu$ is the transition frequency, h is the Planck constant, and k is the Boltzmann constant. In Fig. 2 we have plotted the observed $\ensuremath{\:\!} ^{13}\ensuremath{\!\;\!}$CO and C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=1\mbox{--}0)$ radiation temperatures at the velocities of the C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O peaks. The best fit of Eq. (3) to the data with X=11.2 and $T_\infty=3.4$ K is presented as a solid curve. The obtained value for $T_\infty$ implies that $\ensuremath{T_{\rm ex}}\mbox{$(J=1\mbox{--}0)$ } = 8.0$ K . To give an impression of the accuracy of the determined parameters, we also show a second fit (dashed line), which also uses data points of low signal-to-noise ratio (SNR). The value for X (now 8.2) is substantially lower, but $(J=1\mbox{--}0)$ (now 8.4 K) differs by only 5%.

Using the derived value for $(J=1\mbox{--}0)$, the excitation temperature $(J=2\mbox{--}1)$ can be solved from the following equation for the observed C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=2\mbox{--}1)$/ $(J=1\mbox{--}0)$ integrated intensity ratio:

$$\frac{\int T_{\rm R}(2\mbox{--}1) \ensuremath{{\rm d}} v}{\int T_... ...)}{f_{\nu_{10}}(T_{\rm ex}(1\ensuremath{\!\:\!} -\ensuremath{\!\:\!}0))}} \cdot$$ (6)

For this calculation, the two C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O data sets were convolved to a common 50 Gaussian beam. The total C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O column density was calculated by using assumption 3. The derived distribution of the C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O column density is shown in Fig. 3. As can be seen in this figure, the column density distribution is relatively flat. The maximum is 5.5 $\ensuremath{\;\!}\times\ensuremath{\;\!}$10¹⁴  , and the familiar figure of the cloud is outlined by the level 3.0 $\ensuremath{\;\!}\times\ensuremath{\;\!}$10¹⁴  . In addition, the excitation temperature shows little variation in the central part of the cloud. The average $(J=2\mbox{--}1)$ above the column density level 2.5 $\ensuremath{\;\!}\times\ensuremath{\;\!}$10¹⁴  is 6.2 K and the sample standard deviation is 0.6 K.

Figure 4: a) C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O column density profile of B68. Crosses mark the same data points presented in the contour map of Fig. 3 (i.e. 50 resolution, uncertainty <30%). The centre of the globule is at $\alpha(2000)=\mbox{17\ensuremath{^{\rm h}} 22\ensuremath{^{\rm m}} 40\ensuremath{^{\rm s}} }$, $\delta=\mbox{$-$ 23\ensuremath{^\circ} 49\ensuremath{^\prime} 48\ensuremath{^{\prime\prime}} }$. The solid line results from averaging over rings at intervals of 10 . For comparison, the long-dashed lines give the BES profile with $\xi _{\rm max}$ = 7.0 and angular radius 106 ; black: convolved to our 50 resolution, grey: unconvolved. The normalised theoretical profile is arbitrarily scaled in y-direction with a factor of 30.3. b) Assuming a depletion law according to Eq. (7), the depleted and smoothed profile is given as a black dashed line. The parameter Z in Eq. (7) and the scaling in y-direction were determined in a -fit to match the measured profile (solid line); individual measurements (crosses in a)) have been omitted for visibility. The scaling corresponds to a no-depletion curve as given by the grey dashed line (unsmoothed). While the dashed lines show the best-fitting depletion model (Z=180 yields the lowest ), the other lines show profiles for other values of Z; dash-dot: Z=20; dash-dot-dot-dot: Z=0. In all cases, the grey lines show the unconvolved total (gas+dust) C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O column densities.

  
4.1.3 Column density profile

We have plotted the radial distribution of the C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O column density in Fig. 4a. Also shown in this figure is the BES column density profile deduced from extinction measurements (see Sect. 2 and Alves et al. 2001). After convolving the BES column densities with a 2-dimensional 50 (full width at half maximum) Gauss function, we can exclude a constant CO abundance in the gas phase. The difference between the measured and the smoothed profile strongly suggests CO depletion in the inner part of B68.

To quantify the observed depletion, we assume a depletion law of the form

$$f_{\rm d} \equiv \frac{\ensuremath{\widehat{Y}_{\rm C^{18}O}}... ...18}O}} } = 1+Z~\left(\frac{n}{\ensuremath{n_{\rm c}} }\right),$$ (7)

where the depletion factor $f_{\rm d}$ is an increasing function of gas density; we write $\ensuremath{Y_{\rm C^{18}O}}\equiv{}n(\ensuremath{\mbox{C$\ensuremath{\:\!}^{18}\ensuremath{\!\;\!}$ O},{\rm gas}} )/n(\mbox{H$_2$ })$ for the fractional C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O abundance and designate $\ensuremath{\widehat{Y}_{\rm C^{18}O}}\equiv{}n(\ensuremath{\mbox{C$\ensuremath{\:\!}^{18}\ensuremath{\!\;\!}$ O},{\rm gas+dust}} )/n(\mbox{H$_2$ })$ as the no-depletion fractional abundance. Equation (7) follows from the assumption that the fractional CO abundance is governed by accretion onto dust grains and desorption processes, which return the molecules into the gas phase. The accretion rate in an isothermal cloud is $A~n(\ensuremath{\mbox{C$\ensuremath{\:\!}^{18}\ensuremath{\!\;\!}$ O},{\rm gas}} )~n(\mbox{H$_2$ })$, while the desorption rate can be written in the form $B~n(\ensuremath{\mbox{C$\ensuremath{\:\!}^{18}\ensuremath{\!\;\!}$ O},{\rm dust}} )$. The adsorption and desorption coefficients A and B depend primarily on grain properties, temperature and the flux of heating particles, as discussed in Sect. 6.1. In a steady state condition the relation $A~\ensuremath{n_{\rm H_2}} ^2~\ensuremath{Y_{\rm C^{18}O}} =B~\ensuremath{n_{\rm H_2}} ~(\ensuremath{\widehat{Y}_{\rm C^{18}O}} -\ensuremath{Y_{\rm C^{18}O}} )$ holds, and hence

$$f_{\rm d} = 1+\left(\frac{A}{B}\right) n(\mbox{H$_2$ })~.$$ (8)

Without knowing what actually is, we can determine Z from Fig. 4, calculating later $A/B=Z/(\ensuremath{n_{\rm c}}\cdot{}n(\mbox{H$_2$ })/n)$; we assume $n(\mbox{H$_2$ })/n$ to be 5/6 (Sect. 2).

Starting from the BES density profile, we first calculated a depleted column density distribution for various values of Z, then smoothed the results to our 50 resolution. Finally, we used the -method to scale the theoretical depleted profile to the measured profile, which corresponds to determining $\ensuremath{N_{\rm c}} (\ensuremath{\mbox{C$\ensuremath{\:\!}^{18}\ensuremath{\!\;\!}$ O},{\rm gas+dust}} )$. Figure 4b shows the profile with the lowest , which has Z=180 and $\ensuremath{N_{\rm c}} (\ensuremath{\mbox{C$\ensuremath{\:\!}^{18}\ensuremath{\... ...$ O},{\rm gas+dust}} )=\mbox{${3.9}\times10^{16}$ }~\ensuremath{{\rm cm}^{-2}}$. This profile reproduces the observed profile from $r=20\ensuremath{^{\prime\prime}}$ to $r=100\ensuremath{^{\prime\prime}}$ almost perfectly. For small values of Z the fits quickly become worse, with increasing by 15% for Z=50 and by 100% for Z=20. Larger values of Z cannot be excluded from the depletion analysis alone, because increases by only a few percent for $Z=180\rightarrow\infty$. This is so because for very large Z the scaling (which determines the gas+dust abundance) can compensate any futher increase in Z. Only by requiring that does not rise into physically unrealistic regimes we can give an upper limit for Z, which will be discussed in Sect. 6.

  
4.2 Monte Carlo simulations

We constructed a series of isothermal model clouds with kinetic temperatures between 6 and 16 K and distances between 40 and 300 pc, all with a BES-like density distribution ( $\ensuremath{\xi_{\rm max}} =7.0$ and $\theta_{\rm R} = 106\ensuremath{^{\prime\prime}}$). The radiative transfer problem was solved with Monte Carlo methods (Juvela 1997). The radiation field was simulated with a large number of model photons resulting from background radiation and emission from within the cloud, where a microturbulent velocity field was assumed. The simulations were used to determine the radiation field in each of the 40 spherical shells into which the model was divided and to derive new estimates for the level populations of the molecules. The whole procedure was iterated until the relative change in the level populations were $\sim$10^-4 between successive iterations.

The observed C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O spectra were averaged over rings at intervals of 20 at distances 0 to 100 from the selected centre position RA 17 19 37.${\rm ^s}$1, Dec -23 46 59 (1950.0). The effective resolution of the averaged C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=1\mbox{--}0)$ and $(J=2\mbox{--}1)$ spectra are 60 and 30 respectively. Corresponding spectra were calculated from the models and averaged with a Gaussian beam to the resolution of the observations. The correspondence between the observed spectra and the model was measured with a value summed channel by channel over all spectra. Averaged spectra at different distances from the centre position were all given equal weight in the fit.

For each pair of T and D, two parameters were optimised: the fractional C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O abundance and the turbulent line width. The best fit was obtained for $\ensuremath{T_{\rm kin}} =6$ K and D=80 pc with $\ensuremath{Y_{\rm C^{18}O}} =\mbox{${6.9}\times10^{-8}$ }$. However, this solution still produces far too little intensity at large offsets, the observed intensity at 100 being more than twice the model prediction.

We further studied a series of models assuming a non-constant C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O fractional abundance, $\ensuremath{Y_{\rm C^{18}O}}\equiv\ensuremath{Y_{\rm C^{18}O}} (n)$, in order to account for CO depletion in the cloud centre. The density dependency was assumed to be $\ensuremath{Y_{\rm C^{18}O}} = \ensuremath{\widehat{Y}_{\rm C^{18}O}} /(1+Z~n/\ensuremath{n_{\rm c}}$) , as suggested in the homogenuous model analysis (Sect. 4.1.3). For the depletion parameter Z we tried a large range of values between 4 and 10⁴.

As in Sect. 4.1.3, we find Z not to be well constrained, lying somewhere between 20 and 300. In contrast to our earlier analysis, we now do find an upper limit for Z, just as we can give a lower limit. The reason for this is that in the homogenuous approach, high densities n, which go along with a high Z, do not influence the excitation conditions ( being determined earlier), while now the total density n and not only $n(\mbox{C$\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!} $ O})$ is taken into account.

For both Z=50 and Z=200, the best fitting values for the distance and the temperature are 70 pc and 7 K respectively. Some other specific combinations of Z, D and T within the ranges Z=20-200, D=70-120 pc and T=6-8 K result in almost equally good fits. The Monte Carlo simulations do not favour a particular model, as the values are not significantly different. They do, however, favour the models including depletion as compared to the Z=0 ones. For Z in the given range the values are a factor of 1.6 lower than in the no-depletion models, and the spectra are well fitted from the centre out to the edge of the cloud. Figure 5 shows the correspondence between the model and the observed spectra.

Figure 5: Comparison of observed and modelled spectra. The observed spectra (histograms) are averages over rings at distances between 0 and 100 from the centre. The modelled spectra (smooth lines) are generated using Monte Carlo simulations as described in the text. The C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=1\mbox{--}0)$ transition is shown in the left column, C $\ensuremath{\:\!} ^{18}\ensuremath{\!\;\!}$O $(J=2\mbox{--}1)$ in the right column. Each row stands for one angular distance from the centre. The model behind the simulation presented in this figure is an isothermal cloud with BES-like density structure at 70 pc distance with a kinetic temperature of 7 K and CO depletion according to Eq. (7) with Z=200.