C18O DEPLETION IN STARLESS CORES IN TAURUS

Amanda Brady Ford; Yancy L. Shirley

doi:10.1088/0004-637X/728/2/144

1. INTRODUCTION

Understanding the transformation of a starless core into a protostellar object is the first step in understanding the overall process of low-mass star formation. A key component in understanding low-mass star formation is determining the evolutionary states of starless cores. The competing models for star formation—ambipolar diffusion and gravo-turbulent fragmentation—each have very different timescales (see, e.g., Shu et al. 1987; Klessen & Ballesteros-Paredes 2004; Mac Low & Klessen 2004; Evans 1999; Hartmann et al. 2001).

Ambipolar diffusion posits that collapse of a starless core is slowed by neutrals slipping relative to ions and the decay of magnetic fields (Shu et al. 1987). The ambipolar diffusion timescale is approximately 5 × 10⁶ yr, assuming an ionization fraction of 10⁻⁷. Gravo-turbulent fragmentation proposes that large-scale turbulence causes fragmentation of individual cores, which then collapse on a free-fall timescale. For H₂ densities of 5 × 10⁵ cm³, that implies a timescale of 5 × 10⁴ yr.

One approach is to use chemical evolution, or the process by which the chemistry of a core changes over time, as another piece of information to distinguish between these theories. In the cold (⩽ 10 K), well-shielded interiors of dense starless cores, significant freeze-out (depletion) of molecular species such as CO, CS, CCS, and HCO⁺ is observed (e.g., Tafalla et al. 2004a). The more evolved cores are expected to have more freeze-out, or depletion (Tafalla et al. 2004a). For example, Crapsi et al. (2007) found that in L1544, 93% of the C¹⁷O has frozen out in the center of this very centrally condensed core (central density of 8 × 10⁵ cm⁻³). Another example is the work of Redman et al. (2002) toward L1689B (central density of 2.8 × 10⁵ cm⁻³; Bacmann et al. 2000) where 90% of CO is depleted after 4 × 10⁴ years, or half the free-fall timescale at that density. However, this is a lower limit, as they assume a sticking coefficient of unity, and a constant density over time (which, as they point out, is unreasonable—the density would increase over time).

The state-of-the-art theoretical studies are chemodynamical models that use a chemical network at each time step during the collapse of a starless cores to calculate abundances. For instance, the models of Lee et al. (2004) show significant chemical depletion happening at roughly 10⁶ yr. It is therefore theoretically possible that the depletion profile can be used to distinguish between ambipolar diffusion and gravo-turbulent fragmentation.

In this paper, we map the J = 2–1 transition of C¹⁸O toward a sample of eight starless cores in the nearby Taurus molecular cloud. By comparing observed integrated intensities with intensities obtained from radiative transfer modeling, we place constraints on the amount of depletion in terms of how much CO is depleted (fraction of depletion f_d), and where this occurs (the depletion radius r_d). Our work is unique because we use a larger sample set (eight cores), all located in the same molecular cloud and are then able to compare the amount of depletion in these cores to evolutionary parameters, such as central density. We do not assume a constant T_g = 10 K profile for our models, as has been done in previous work, and investigate degeneracies in choice of temperature versus choice in abundance. We also investigate possible effects of the gas temperature and the dust temperature not being equal. Unless otherwise specified by T_dust, temperatures referenced are the kinetic gas temperature, T_g.

This paper is organized as follows. In Section 2 we discuss our sample selection and motivation. We also describe our observations and modeling methodology, including the implementation of the RATRAN (Hogerheijde & van der Tak 2000) radiative transfer code. In Section 3 we test our modeling procedure, including potential degeneracies of various parameters, and present results for each individual core. In Section 4 we present our analysis and in Section 5 we give our conclusions.

2. MOTIVATION AND PROCEDURES

2.1. Sample Selection and Motivation

Our sample consists of cores in Taurus, at approximately 140 pc (Loinard et al. 2005), so they are easy to compare at the same spatial scale. We chose this particular set because they span over two orders of magnitude in central density (see Table 1). Figure 1 shows observations of the dust continuum emission at 850 μm taken with SCUBA, as well as C¹⁸O 2–1 integrated intensity maps. These cores have at least some rough central dust continuum concentration.

**Figure 1.** Left: 850 μm dust continuum emission, from SCUBA data, normalized (arbitrary units). Right: C¹⁸O integrated intensity. Beam size indicated by solid circles. Contours are percentage of the peak for the 850 μm continuum, and value indicated for the C¹⁸O-integrated intensity.
Download figure:
Standard image High-resolution image

Table 1. Characteristics of these Cores

Core	R.A. (J2000)	Decl.	v_LSR	Central Density	Tinner	Touter	$\sigma _{I({\rm C}^{18}O\, 2-1)}$
			(km s⁻¹)	(cm³)	(K)	(K)	(K km s⁻¹)
TMC-2	04:32:45.5	+24:25:07.7	+6.30	2 × 10⁴	10.7	10.9	0.180
L1498	04:10:53.0	+25:10:08.6	+7.78	2 × 10⁴	10.2	10.3	0.194
L1512	05:04:08.6	+32:43:24.6	+7.10	8 × 10⁴	8.2	8.6	0.157
L1489	04:04:47.6	+26:19:17.9	+6.85	1 × 10⁵	11.7	12.4	0.166
L1517B	04:55:18.3	+30:37:49.8	+5.78	2 × 10⁵	8.7	9.4	0.200
L1521E	04:29:14.9	+26:13:56.6	+6.90	3 × 10⁵	7.6	8.5	0.187
L1495A-S	04:18:39.9	+28:23:15.9	+7.29	4 × 10⁵	9.7	12.8	0.187
L1544	05:04:17.2	+25:10:43.7	+7.18	8 × 10⁵	8.5	11.4	0.166

Download table as: ASCII Typeset image

Given a constant temperature profile (see Section 3.2 for how this changes if the temperature varies) and a constant [C¹⁸O]/[H₂], one would expect to observe more intense C¹⁸O emission in the center than that at outer radii. For many of these cores, that is not what is observed, as seen in Figure 1. Given the high central densities of these cores, the C¹⁸O could be depleted by freeze-out onto dust grains. The greater the density, the more freeze-out one would expect, so this freeze-out would occur preferentially in the centers of the cores (Goldsmith 2001). We shall model this effect with radiative transfer.

2.2. Observations

We took observations of the J = 2–1 transition of C¹⁸O (219.56036 GHz) during 2007 February 5–8, using the sideband-separating ALMA 1 mm (Band 6) prototype receiver on the 10 m Heinrich Hertz Submillimeter Telescope (HHT). The beam size (full width at half-maximum, FWHM) at the observed frequencies is 34''. The main beam has been measured to be Gaussian to detectable levels using maps of Mars (J. Bieging 2008, private communication). For the ALMA prototype receiver, the side lobes are insignificant, as those are down by at least 20 dB. The C¹⁸O line was tuned to the lower sideband (LSB) with a 6.0 GHz intermediate frequency (IF). Receiver noise temperatures were between 90 and 100 K with system temperatures <250 K. Typical sideband rejections were greater than 13 dB and were ignored in the subsequent calibration. The ALMA 1 mm prototype receiver is both sideband separating and dual polarization resulting in 4 IF outputs (Vpol LSB, Hpol LSB, Vpol upper sideband (USB), and Hpol USB) to the spectrometers. Unfortunately, we were only able to use one polarization (Vpol) of the LSB since there was only one high spectral resolution back end (Δν < 100 kHz) available in 2007 February. A Chirp Transform Spectrometer with 32 kHz spectral resolution was used for all C¹⁸O observations. HHT uses the standard Chopper Wheel Calibration Method (Penzias & Burrus 1973). The data were baselined and calibrated using the CLASS reduction package. Mapping observations were obtained on a 10'' grid covering 2' × 2'. Each position was observed for 2 minutes of integration time in position-switching mode. OFF positions were checked for C¹⁸O emission. Calibration scans of Saturn were made every 2 hr after pointing and focusing the telescope. The final reduced spectra were placed on the T_mb = T*_A/η_mb scale with an average main beam efficiency of η_mb = 0.74.

RMS errors for the observed integrated intensity points generally ranged from 0.1 to 0.2, the rms error of the observed position with the largest σ_I for each core is given in Table 1. Errors are defined as

$\begin{equation} \sigma _{I} = \sqrt{\Delta V_{{\rm line}} \Delta V_{{\rm channel}}} \times \sigma _T, \end{equation} \tag{ 1 }$

where ΔV_line is the velocity interval over the whole line, from line wing to line wing. ΔV_channel is the channel width, of 0.157 km s⁻¹, and σ_T is the baseline rms error, assuming a linear baseline.

2.3. Modeling

A radiative transfer calculation is appropriate for this work, as computing a column density from observations requires the assumption of a single excitation temperature throughout the core, which is incorrect. Also, because we are in non-LTE, the kinetic temperature is not equal to the excitation temperature along the entire line of sight. Therefore, we use the one-dimensional version of the Monte Carlo radiative transfer code, RATRAN, developed by Hogerheijde & van der Tak (2000), to generate radial profiles of C¹⁸O emission. RATRAN uses Accelerated Lambda Iteration to efficiently reach convergence between the level populations and calculations of J_ν, the mean intensity of the radiation field.

RATRAN inputs for our models include the following.

1.
C¹⁸O molecular collision rate file from the LAMBDA database (see Schöier et al. 2005).
2.
Dust grain properties. We use OH5, see Ossenkopf & Henning (1994), Table 1, Column 5. These opacities typically give the best fit in starless core dust continuum radiative transfer models (see Shirley et al. 2005).
3.
Radial points. We used 100 radial shells, from 20 AU to 30,000 AU, with equal logarithmic spacing.
4.
H₂ density profiles. The density profiles are constrained by the Shirley & O'Malia dust continuum radiative transfer models (see also Shirley et al. 2005 for methodology). The models assume isothermal, static Bonnor–Ebert (B-E) spheres with various central densities (Bonnor 1956; Ebert 1955) which are derived under assumptions of pressure-bounded hydrostatic equilibrium. See Figure 2 for our density profiles. While the B-E spheres are isothermal, their resulting temperature profile from the dust continuum radiative transfer calculation is non-isothermal; however, Evans et al. (2001) found that the correction to n(r) for the non-isothermal case is negligible for low-mass starless cores. See Shirley et al. (2005) for further discussion. Our central densities are not identical to those used in other published work, specifically in Young et al. (2004), because Shirley & O'Malia use a finer central density grid in the modeling. We also assume an H₂ ortho–para (o:p) ratio of 3:1, as in Flower et al. (2006), and find that our results do not depend strongly on o:p ratio.Running these density profiles through RATRAN requires choosing a cutoff radius, beyond which we assume H₂ (and, by extension C¹⁸O) has its standard interstellar medium (ISM) value (see below). The cutoff radius is necessary for modeling; it is not physically motivated. In reality, the edge of a core is difficult to define, and there is more likely to be a gradual decline in density versus a sharp discontinuity between the core and surrounding ISM, we investigated how sensitive our results were to this choice. This was especially important given that work by Onishi et al. (1996) found large amounts of C¹⁸O in the ambient Taurus molecular cloud, not necessarily associated with starless cores. We chose 30,000 AU as a reasonable outer radius for our models, consistent with the 0.1 pc typical core size given by Bergin & Tafalla (2007). For the C¹⁸O observations, we mapped only to ±60'' (corresponding to 8400 AU at 140 pc) from the dust continuum peak, as the dust maps were restricted by chopping (θ_chop/2 = 60''). Figure 3 shows that most of the contribution to the total column density occurs in the center parts of the cloud, suggesting our choice of outer radius is reasonable.
5.
C¹⁸O abundance profiles. We assume a step function, wherein C¹⁸O is depleted by some fraction f_d at some radius r_d. We define
$\begin{equation} {f}_{d}\equiv \frac{([{\rm C}^{18}\rm O]/[{\rm H}_{2}])_{\rm ISM}}{([{\rm C}^{18}\rm O]/[{\rm H}_{2}])_{\rm interior\ to\ {r}_{d}}}. \end{equation} \tag{ 2 }$
We vary f_d from 1.0 (no depletion) to 10³ in steps of 2ⁿ, with n varying from 0 to 10. We vary r_d from 1250 AU to 30,000 AU in steps of 1250 AU.Exterior to r_d, [C¹⁸O]/[H₂] has its standard ISM abundance, which we define in this work as 1.6 × 10⁻⁷, hereafter referred to as the "canonical abundance." This is slightly lower than the work by Frerking et al. (1982), their 1.7 × 10⁻⁷ versus our 1.6 × 10⁻⁷; we adjust it to be consistent with the work by Wilson & Rood (1994), who find a slightly lower ¹⁸O/¹⁷O ratio. Using χ²_r fitting, we find r_d and f_d that best match both the observed spectra for the central region and the observed integrated intensity of C¹⁸O emission across the core. Previous work on two of our cores, L1498 and L1517B, by Tafalla et al. (2004a) has shown that the abundance drop is quite sharp and therefore that a step function is reasonable (see also Lee et al. 2003). See Figure 4 for our C¹⁸O density profiles.
6.
Temperature profile of C¹⁸O. The temperature profiles are self-consistently calculated from dust continuum radiative transfer models of the cores. The dust continuum emission from the eight cores in this study was modeled as part of a larger sample of starless cores that were observed with SCUBA (Y. Shirley & K. O'Malia 2011, in preparation). A one-dimensional dust continuum radiative transfer code, CSDUST3, (Egan et al. 1988) is used to match submillimeter intensity profiles and spectral energy distributions. The heating in the dust continuum models is supplied by the Interstellar Radiation Field (ISRF) which is allowed to vary in overall strength and the degree of extinction (see Y. Shirley & K. O'Malia 2011, in preparation). For our purposes in RATRAN, we approximate the T(r) profile as the same for the gas kinetic temperature and the dust temperature. Work by Young et al. (2004) has indicated that this assumption holds true for densities above 10⁴ in L1498, above 10³ in L1512, and above 5 × 10⁴ in L1544. Discussion of the implications of these assumptions can be found in Section 3.2. Figure 5 shows the temperature profile used in each core.
7.
Velocity profile. We use a static velocity field with turbulent and thermal line broadening described by the Doppler b parameter = 0.2, since that quantity fits the observed line widths. Late stages of collapse would produce line widths broader than those observed. Work by Myers (2005) has indicated that the shape of the B-E sphere density profiles does not change substantially even in the collapsing case, until late in the collapse history.

**Figure 2.** H₂ number density vs. radius profiles used for each of the cores.
Download figure:
Standard image High-resolution image

**Figure 3.** We explore possible effects of choice of R_outer. We integrate the column density from r = 0 (the center) to r = R_outer (the outer radius, which we vary from zero to our model cutoff of 30,000 AU), then divide that value by the integrated column density from 0 to 30,000 AU, and plot that ratio as a function of R_outer. We do this for three model cores with densities of 2 × 10⁴, 3 × 10⁵, and 8 × 10⁵ cm⁻³, representing the full span of our sample's central density.
Download figure:
Standard image High-resolution image

**Figure 4.** Density of C¹⁸O vs. radius for representative combinations of *r_d* (in AU) and *f_d*. Overplotted dashed line is the density of H₂, just shifted down by the canonical abundance, [C¹⁸O]/[H₂] = 1.6 × 10⁻⁷, so as to fit on the same plot.
Download figure:
Standard image High-resolution image

**Figure 5.** Temperature vs. radius profiles from dust continuum radiative transfer models used for each of the cores.
Download figure:
Standard image High-resolution image

After the model has been run through RATRAN, we use the Miriad image analysis package (Sault et al. 1995) to convolve the outputs with our 34'' beam size and integrate over velocity to determine the integrated intensity in K km s⁻¹.

Before we began any of our modeling, we first attempt to replicate previously published results, as a check on our procedures for running the radiative transfer code. We were graciously provided density and temperature profiles by Mario Tafalla from work on L1521E (Tafalla & Santiago 2004). See Figure 6 for those profiles, as compared to our inputs for the same core. We were able to replicate their published best-fit curves through our analysis, even though they used a different radiative transfer code (Bernes 1979). This indicates that any differences in our results are based on differences in real physical inputs, such as density and temperature, and not differences in modeling procedures. L1521E results are further discussed in Section 3.3.

**Figure 6.** For L1521E, comparison of temperature and density profiles used in Tafalla & Santiago (2004) to those in this work. The bottom panels compare observed and model values of integrated intensity (K km s⁻¹) of the C¹⁸O J = 2–1 for both this work and theirs. In the bottom left panel, blue crosses are their observed points. The solid black line was made by using their temperature and density profiles, which we ran through RATRAN and our associated analysis. It appears to be identical to their best-fit curve found in their Figure 3. In the bottom right panel, the black points are our observed points. The red (upper) line is for no depletion, with our temperature and density profiles (twice the canonical abundance), run through RATRAN. The green (lower) line is our overall best fit with *r_d* = 2500 AU and *f_d* = 64. For further discussion, see Section 3.3.
Download figure:
Standard image High-resolution image

3. DISCUSSION

We begin our modeling by using model temperature and density profiles as discussed above and then varying r_d and f_d, as in previous work to find best fits. r_d and f_d are considered free parameters and are fully explored as such. As these parameters represent the amount of depletion in each core, which are then used to make chemical evolution timescale arguments, it is important to understand how sensitive the results are to the assumptions that go into their calculation. We explore how much these calculated values of r_d and f_d would change if the RATRAN inputs changed. Of the seven RATRAN inputs described in Section 2.3, we explore factors influencing the canonical abundance and temperature profiles of C¹⁸O. Inputs for molecular collision rates, dust grain properties, and velocity profiles are already as well constrained as they are likely to be for now, and there is good agreement in the literature, so we do not explore those further. H₂ density profiles are well constrained by dust continuum modeling, and choice of cutoff radius is of little consequence to the model column densities.

We choose to explore the canonical abundance of C¹⁸O because there is a large amount of disagreement in the literature over this value, and as we ourselves have no good way of independently verifying its correct value. We choose also to explore the sensitivity to the temperature profile, because while we self-consistently calculate a varying profile, other work assumes a constant profile. Furthermore, there are physically motivated reasons to believe that choosing different values of canonical abundance and/or temperature than described in the previous section would yield different best-fit r_d and f_d. The point of this analysis is not to find an absolute best fit for this multi-degree parameter space, since input data are too limited to make such a finding meaningful. The point is to see how much caution is needed—if we use r_d and f_d as proxies for chemical evolution, and if these proxies are extremely sensitive to our input assumptions, then our work and future work need to be mindful of this.

3.1. Variations in the ISM C¹⁸O Abundance

As discussed in Section 2.3, we assume that, in the absence of freeze-out, the ratio of [C¹⁸O]/[H₂] is fixed for all of our cores at the canonical abundance of 1.6 × 10⁻⁷. However, not everyone agrees that this is the correct value to use. A comparison of published work in the field shows that values of the canonical [C¹⁸O]/[H₂] abundance differ by a factor of seven, from 0.7 × 10⁻⁷ (Tafalla & Santiago 2004) to 1.5 × 10⁻⁷ (Tafalla et al. 2004a) to 4.8 × 10⁻⁷ (Lee et al. 2003). It may be that more than one of these values is correct, as the canonical abundance could even vary across the Taurus region. This may be a reasonable possibility in Taurus, given the large-scale structure of filaments, cavities, and rings and subsequent variations in the radiation field found by Goldsmith et al. (2008).

Because of the wide range of published values for this particular input parameter, we choose to explore it in further detail by varying the canonical abundance and seeing how that affected f_d and r_d. A full exploration of all possible values of canonical abundance is beyond the scope of this work, but we felt it necessary to at least gain a quick understanding of the nature of the impact. We ran models with all of the same input parameters described in Section 2.3, but for half the canonical abundance, as well as twice the canonical abundance. In changing the canonical abundance, some of the cores were best fit by entirely different values of f_d and r_d than for an abundance ratio of 1.6 × 10⁻⁷ (see the next section for further discussion). To our knowledge, even this cursory exploration of variances in canonical abundance has not been done before. Further constraints on canonical abundance are needed to be confident that f_d and r_d values generating lowest χ²_r for choice of one canonical abundance have physical meaning. More details can be found in Section 3.3.

3.2. Sensitivity to Temperature Profiles

We also explored how sensitive best-fit f_d and r_d values were to choice of temperature values. We started this investigation by using all of the same input parameters as before (see Section 2), but changing the temperature profiles. We explored both temperature profiles that vary and those remaining constant across the core. Our objective was not, in this way, to constrain the temperature profiles of these cores but to understand how sensitive results are to this parameter. We believe we are able to constrain well the temperature profiles from dust continuum radiative transfer modeling, described in Section 2.3. Additionally, a varying temperature profile, as used in our work, makes physical sense: the areas near the center of the core, being shielded from the ISRF, should be colder than those on the outside (Zucconi et al. 2000). Dust continuum modeling of our sample of cores suggests that while the temperature variation is only a few Kelvin, this means that the outer radius of the core is up to 30% warmer than the center (Figure 5).

However, much of the other work in this field has assumed that the temperature remains constant throughout the core, with a typical value of T = 10 K for all radii. Given that our work does not make this assumption, we explore how temperature profiles could affect which f_d and r_d best fit the observational data. It is, of course, the combination of temperature and density that produce the observed intensity. If a core is cooler than assumed, then one would require a higher density to match the observations and vice versa. This would lead one to believe a core is less depleted than if one started with a higher temperature.

We ran a representative set of models for B-E spheres with central densities of 3 × 10⁵ and 2.5 × 10⁴ (see Figures 7 and 8). These figures show the C¹⁸O 2–1 intensity versus radius for various combinations of r_d and f_d, using various temperature profiles. We plotted three constant temperature profiles: T = 7 K, T = 10 K, and T = 13 K, as well as a temperature profile that varied with radius between 8 and 13 K. For Figure 7, note how close the red line (temp = 10 K) is to the black line (the varying temperature profile). Since the varying temperature profile goes from 8–12 K, it seems one can simply take the average of those numbers and use that as the single temperature. Indeed, the shapes of all the curves are fairly similar. This is consistent with work by Tafalla et al. (2004a), who found that modeling with isothermal profiles can well reproduce the data. However, the exact choice of that isothermal temperature is important. From Figures 7 and 8, one can see that not all cores can be modeled at T = 10 K. One would need further information, such as self-consistent temperature modeling from dust emission, as in this work, to know which isothermal temperature to use. Figure 5 shows the temperature profiles we adopt for each core in the rest of this work.

**Figure 7.** Intensity vs. radius for various combinations of *r_d* and *f_d*. This is for a theoretical B-E sphere with an example density of 3 × 10⁵ cm⁻³. The bottom panel shows the temperature profiles used in the nine panels above, with matching line styles.
Download figure:
Standard image High-resolution image

**Figure 8.** Intensity vs. radius for various combinations of *r_d* and *f_d*. This is for a theoretical B-E sphere with a central density of 2.5 × 10⁴ cm⁻³. The bottom panel shows the temperature profiles used in the nine panels above, with matching line styles.
Download figure:
Standard image High-resolution image

Another common assumption is that the gas and dust temperatures are well coupled (we make that assumption in our work here). At high densities this is true, as collisions between the gas and dust temperatures result in equal kinetic gas and dust temperatures. At low densities, however, as in those found on the outside of the core, this is not necessarily the case. Indeed, recent work by Young et al. (2004) has shown that gas and dust temperatures can vary substantially in those conditions. However, precisely because T_d ≠ T_gas only on the outside of the core, in the lowest density areas that contribute least to the column density (see earlier section), we do not expect this to be a large effect. We were graciously provided gas and dust temperature profiles by Young et al for three cores for which they did detailed energetics calculations: L1498, L1512, and 1544. We used our same model density profiles but their calculated gas and dust temperature. See Section 3.3 for further discussion.

We will see in Section 3.3 that in many cases we find different f_d and r_d values for many of the cores than previously published. Based on the above analysis, we believe some of this discrepancy is likely due to our temperature profiles being different from those used in other work. As another quick check, we simultaneously varied both the canonical abundance and the temperature for a small slice of parameter space. We found that there can also be a degeneracy in the abundance and the temperature used, as in Figure 9, which worsens the problem. One can well fit the data by assuming a higher canonical abundance and a lower temperature, or a lower canonical abundance with a higher temperature.

**Figure 9.** Illustration of modeling degeneracies present in L1521E, as an example. Note that our data points (black) can be fit well with either our temperature profile, which varies with radius and double the canonical abundance of C¹⁸O, or with a constant temperature profile of T = 10 K and the regular abundance. Model points exactly match coordinates of observed points.
Download figure:
Standard image High-resolution image

3.3. Summary of Findings for Each Core

Figure 10 plots the C¹⁸O-integrated intensity χ²_r for depletion profiles with combinations of f_d and r_d. Due to the high scatter in the observed intensity points, we azimuthally averaged our observed points before comparing them to our one-dimensional model. We divided the data into 10'' annuli around the center core, and averaged the observed and model points over that range. We calculated a χ²_r as follows:

$\begin{equation} \chi ^{2}_r\equiv \sum \frac{ ({I}_{{\rm obs}}- {I}_{{\rm mod}})^{2}}{\sigma _{{\rm obs}}*({\rm number \,of \,annuli - 1}) }. \end{equation} \tag{ 3 }$

**Figure 10.** χ²_r maps of each of our eight sources, assuming a canonical abundance, 1.6 × 10⁻⁷. Observed and model points have been azimuthally averaged in 10'' annuli. The black regions represent χ²_r less than minimum(χ²_r)+1. For TMC2, L1521E, L1495AS, and L1544, no combination of *r_d* and *f_d* gives a χ²_r less than 1. Note that a wide range of *f_d* gives very similar χ²_r values, showing this value is poorly constrained. Blue diamonds represent the location of the lowest χ²_r, note they are in the bottom left corner for TMC2 and L1521E.
Download figure:
Standard image High-resolution image

We define I_obs as the mean intensity of each of the observed points, for each 10'' annuli. We define I_mod as the intensity of the middle model point in the same annuli. We define σ_obs as the variance in the intensity of the observed points for each annuli. For annuli with only one observed points (for some of the cores, there was only one observed data point in the 0''–10'' annuli), we define the variance as simply the error on that one point. Errors for these individual points are defined in Equation (1). We have six annuli, going from 0'' to 60''.

For TMC-2, L1521E, L1495A-S, and L1544, no combination of f_d and r_d with a canonical abundance gave a χ²_r< 1. However, as seen in Figure 11, increasing the canonical abundance by a factor of two was able to give χ²_r< 1 values under 1, except for L1544 (discussed below). L1512 and L1517B had combinations of f_d and r_d that gave χ²_r< 1 for both the canonical abundance and twice the canonical abundance, however using twice the canonical abundance ultimately gave lower χ²_r's. None of the cores had lower χ²_r when the canonical abundance was halved, so we do not favor this possibility. In addition to these χ²_r maps, we also calculate χ²_r of the line profiles for the central 10'' regions of these cores, as shown in Figure 12. Figure 13 shows observations of the central regions with best-fit models (see Table 2). Both χ²_r calculations have significant overlap between the best-fit values of f_d and r_d. While all of these models use the temperature profile calculated from dust continuum radiative transfer models, Figure 14 plots the integrated intensity χ²_r for three cores using the temperature profiles from the energetics calculations of Young et al. (2004), where T_d ≠ T_gas for all radii. This gives very similar results to Figure 10, consistent with the conclusion that using uncoupled gas and dust temperatures does not produce a substantial effect on the C¹⁸O J = 2–1 emission.

**Figure 12.** χ²_r of the central regions, showing the models that best match the central spectrum, assuming canonical abundance. The contours (starting with the minimum χ²_r value and increasing by 1) are similar in shape to Figure 10, but with higher χ²_r, as there is more noise in the line. L1495AS is not well fit due to the flat line shape, see Figure 13. Blue diamonds represent the location of the lowest χ²_r, note they are in the bottom left corner for TMC2 and L1521E. Note this figure is not azimuthally averaged since it is focused only on the central region of each core.
Download figure:
Standard image High-resolution image

**Figure 13.** Central spectra of each core, for both the observed values (black solid) and the model value (red dotted) with the lowest χ²_r for the whole region. See Table 2.
Download figure:
Standard image High-resolution image

**Figure 14.** χ²_r maps of cores for which work by Young et al. (2004) found that gas and dust temperatures are not the same. Data and model points have been azimuthally averaged in 10'' bins; contour levels start with the minimum χ²_r value. Blue diamonds represent the location of the lowest χ²_r.
Download figure:
Standard image High-resolution image

Table 2. Best Fit: r_d, f_d, and Abundances (as a Multiple of the Canonical Abundance)

Core	r_d	f_d	Abundance	χ²_r
	(AU)
TMC-2	5000	512	2	0.094
L1498	7500	16	1	0.063
L1512	10000	256	2	0.13
L1489	7500	1000	1	0.48
L1517B	10000	32	2	0.087
L1521E	2500	64	2	0.18
L1495A-S	7500	4	2	0.42
L1544	7500	16	2	1.5

Download table as: ASCII Typeset image

Inspection of Figures 10 and 11 shows that there is some degeneracy between exactly where the depletion occurs and by exactly how much. In many cases it is easier to constrain r_d than f_d. These figures also show the difficulty in constraining f_d, as a range of f_d values all give χ²_r < 1. Figure 15 shows the observations for each core, as well as models with various χ²_r values. Despite this degeneracy, it is clear that L1521E has low depletion while cores such as L1512 and L1517B, for example, have higher depletion. The best-fit f_d and r_d values for each source are listed in Table 2. Details for each core are given below.

1.
TMC2. This is a difficult core to model, since, based on dust emission (see Figure 1), it is not as centrally condensed, and is non-azimuthally symmetric. In the case of canonical abundance, it is best matched by a model profile with no depletion. However, even this profile has a fairly high χ²_r. In order to get the lowest χ²_r, we had to double the C¹⁸O ISM abundance. The lowest χ²_r for TMC-2 is a depleted model with twice the canonical abundance.
2.
L1498. We find this core is well modeled with canonical abundance of C¹⁸O with significant depletion. We find similar depletion to work by Willacy et al. (1998), which sets a lower limit on the depletion fraction at 8. Work by Tafalla et al. (2004a) found C¹⁸O to be entirely depleted (f_d > 1000) inside a radius of 9940 AU. We find a r_d of 7500 AU but only a factor of 16 depletion. However, as shown in Figure 10, we find that our findings and those of Tafalla et al. (2004a) to have very similar χ²_r. The slight discrepancy may be due to our azimuthally averaging our data, while Tafalla et al. (2004a) did not.
3.
L1512. The best match for this core is twice the canonical abundance of C¹⁸O with significant depletion. This is consistent with work by Lee et al. (2003). They give a greater r_d than this work (15,000 AU instead of this work's 10,000 AU) but a lower f_d (25 versus this work's 256). As in the case of L1498, Figure 11 shows both findings to have similar χ²_r's, and could be explained by our azimuthally averaging our data, while Lee et al. (2003) did not.
4.
L1489. This core is well modeled with canonical abundance of C¹⁸O with significant depletion.
5.
L1517B. Our results for L1517B show significant depletion, as does work by Tafalla et al. (2004b). Their work, however, found a best fit (to their observations) with f_d > 1000 and r_d = 11,620 AU, which is more depletion than our work (see Table 2), perhaps because they use a T = 10 K temperature profile, higher than ours.
6.
L1521E. Our best fit for this core is twice our canonical abundance with r_d = 2500 AU and f_d = 64. This represents little or no depletion, consistent with work by Tafalla & Santiago (2004), and suggests an early evolutionary state, consistent with work done by Hirota et al. (2002). Our exact findings for r_d and f_d, are not identical to those in this work, however, as Tafalla & Santiago (2004) find L1521E best fit by no depletion at all. However, Tafalla & Santiago (2004) use a different density profile and a different canonical abundance of [C¹⁸O]/[H₂] (ours is 1.6 × 10⁻⁷, versus 0.7 × 10⁻⁷), and assumes a constant temperature T = 10 K at all radii (Figure 6). Additionally, Tafalla & Santiago (2004) observational values for the C¹⁸O 2–1-integrated intensity are 30% higher than ours (Figure 6). Discussions with Tafalla have indicated that some of this discrepancy may be the result of different beam sizes, as the FWHM of his observations is 12'' versus our 34'', and systematic calibration errors (the IRAM observations were double sideband, while the HHT observations were single sideband).
7.
L1495A-S. Very high χ²_r, not well matched by any r_d and f_d with canonical abundance. Similar to TMC-2 in that it is best matched by a profile with no depletion if one assumes a canonical abundance. The lowest χ²_r value, however, comes from doubling the abundance and assuming significant depletion.
8.
L1544. Our results show this core to be best fit by using twice the canonical abundance but still significantly depleted, which is consistent with the conclusion of Lee et al. (2004). They find an f_d of 25 at an r_d of 9000 AU. This is more depletion than our findings, f_d of 16 at an r_d of 7500 AU, however they are using a higher canonical abundance, 4.8 × 10⁻⁷. Caselli et al. (1999) found this core to have an abundance drop by a factor of approximately 10 at a radius of 6500 AU. Our Figure 11 shows this result to have a higher χ²_r than our best fits, but still fairly consistent with our results. Aikawa et al. (2001) find a central hole of 4000 AU, inside of which no CO is present. This is more depletion than our findings, but their work assumes a higher H₂ number density than ours, 1.5 × 10⁶ cm⁻³. Our best-fit value still has a χ²_r greater than 1, as seen in Table 2. This is likely due to L1544 being not round, but inspection of Figure 15 shows that we were able to find a combination of f_d and r_d that, at least by eye, appears to be a good fit.

**Figure 15.** Observational data (in black), as well as the combination of *r_d* and *f_d* that minimized the reduced χ²_r (in blue). See Table 2. Errors in observed points are as defined in Equation (1). TMC2, L1521E, L1495A-S, L1517B, and L1544 are fitted with an abundance twice canonical. Plotted in black are the combination of *r_d* and *f_d* that give the highest χ²_r, less than minimum(χ²_r)+1. Black and blue lines are similar in most cases, illustrating the difficulty in constraining these quantities. See contour plots in Figures 10 and 11. This shows the difficulty in constraining *f_d*. For L1544, only one *r_d* and *f_d* combination gives χ²_r < min(χ²_r)+1.
Download figure:
Standard image High-resolution image

4. ANALYSIS

We find that all cores but one, L1521E, are best fit by models that include significant depletion. Our conclusion that all but one of these cores is significantly depleted holds even taking into account reasonable variances in outer radius, temperature, and canonical abundance. In no case did we find varying any of these parameters substantially affected the results. This work is unique in that we have tested a more complete parameter space of r_d, f_d, and C¹⁸O canonical abundance than previously published work. While outer radius does not impact our results, choices of temperature profile and canonical abundance can. Even fixing these parameters, we cannot, however, robustly distinguish between a factor of 10 depletion and a factor of 1000 in f_d, as the χ²_r are very similar (e.g., Figures 15, 10, and 11). In contrast, the depletion radius, r_d, for some of our cores is more strongly constrained than f_d; however, r_d values for many of these cores are separated by little more than a beam width, and four of the cores have the same value for r_d. Work by Lee et al. (2003) also found f_d hard to constrain.

At the distance of Taurus, information in the HHT maps contains information on scales larger than θ_mb/2 = 2100 AU. Greater spatial resolution would better constrain r_d. Improved resolution may even help constrain f_d. For instance, future observations with the 50 m Large Millimeter Telescope operating at 1.3 mm would increase our resolution by a factor of five. The difficulty in determining the absolute value of f_d limits our ability to calculate an accurate chemical depletion timescale. The rate of freeze-out, or the chemical depletion timescale of CO, is given in Rawlings et al. (1992) and Redman et al. (2002). The version in Redman et al. (2002) gives ${\dot{n}}_{\rm CO}\propto {\it n}_{\rm CO}$ . The solution of this differential equation implies that the amount of gas-phase CO decreases exponentially with time for a constant grain temperature, n_CO ∝ e^−t. The inability to tell the difference between f_d = 10 versus f_d = 1000 depletion translates into a factor of four difference in the timescales. As pointed out in Redman et al. (2002), for a core with H₂ density of 2.8 × 10⁵ (near the median central density of cores modeled by Y. Shirley & K. O'Malia 2011, in preparation), 90% of the CO is depleted at 43,000 years, or approximately half the central density free-fall timescale. But if the CO is actually depleted by 99%, this brings it to the free-fall timescale, and 99.9% depletion would be longer than the free-fall timescale. The difference between f_d = 10 versus f_d = 1000 could mean the difference between a core being much older than the free-fall timescale, hence disfavoring gravo-turbulent fragmentation as the dominant core evolution scenario, or being younger than the free-fall timescale, and consistent with gravo-turbulent collapse timescales.

Better constraints on the amount of depletion will be important for constraining ambipolar diffusion timescales, which depend on ionization fraction. Caselli & Walmsley (2001) found that freeze-out of C¹⁸O can be accompanied by depletion of metal ions, significantly lowering the ionization fraction. In the case of L1544, they find that ionization fraction can be 10⁻⁹, which would put the ambipolar diffusion timescale on the order of the free-fall timescale. They argue that better constraints on amount of CO depletion could place constraints on the ionization fraction, so these issues could be coupled. This timescale could be shorter than the chemical evolution timescale (depending on the amount of depletion, as discussed earlier), which would imply ambipolar diffusion is not the dominant process.

In Figure 16, we compare the best χ²_r model parameters to the physical parameters of the cores (central density, central temperature, as well as density and temperature at best-fit r_d values). The plotted points indicate the best-fit f_d and r_d (minimum χ²_r), as given in Table 2. Error bars show the range of f_d and r_d that correspond to χ²_r < min(χ²_r)+1. Figure 16 shows that it is much more difficult to constrain f_d than r_d. We find no strong correlations for any of these parameters with r_d, and none for f_d, although large error bars on f_d make teasing out trends difficult. One would expect the chemical evolution to track the dynamical evolution; as a core contracts, more CO would freeze out. Our work suggests something more complicated, as the cores with greater central density are no more likely to be significantly depleted than those with low density. Warmer cores are also no less likely to be depleted. For instance, L1521E has one of the lowest dust temperature profiles yet has the least depletion.

Our findings are consistent with the conclusion that these cores may be evolving at different dynamical rates. There are several possible explanations for why cores in Taurus could be evolving at different rates. For example, that the magnetic field strengths are different across Taurus. Since the magnetic pressure ∝B², small changes in magnetic field would lead to large changes in pressure support. This would lengthen the collapse timescale until the cores become supercritical. Unfortunately, the magnetic field strengths in Taurus starless cores are poorly constrained and the observed variation from dust polarization measurements is over an order of magnitude (e.g., ≈10 μG in L1498 to >100 μG in L1544; Kirk et al. 2006; Crutcher et al. 2004).

Rates of evolution may also be affected by turbulence. For instance, the external pressure is an important factor in determining the gravitational stability of the core Bonnor (1956). While the starless cores themselves have very narrow line widths that are not dominated by turbulence, the turbulence outside the cores may play an important role in the bounding pressure. The transition from thermal to turbulent mediums can be quite sudden in dense cores (Pineda et al. 2010). The amount of turbulent pressure varies across Taurus (Goldsmith et al. 2008). This external pressure is in addition to the weight of cloud material surrounding the cores (Lada et al. 2008). Observations of the surrounding turbulent environment (e.g., ¹³CO maps; Goldsmith et al. 2008) compared with observations of the transition between the coherent core and more turbulent surrounding medium (e.g., NH₃ maps with the new 7-beam receiver at the Green Bank Telescope) are needed to elucidate the importance of the bounding pressure and the role of turbulence on core evolution. Further information on the effects of magnetic fields and turbulence specific to Taurus can be found in Lee et al. (2009) and Goldsmith et al. (2008).

It is also possible that freeze-out is not as simple as we may expect, as desorption processes may be important. Roberts et al. (2007) investigated the effects of desorption mechanisms such as X-rays heating, direct impact of cosmic rays, UV irradiation, and exothermic reactions on dust grain surfaces (e.g., formation of H₂). They dismiss X-ray effects, as there are not many X-rays in these types of molecular clouds, but find that the other three mechanisms can potentially be significant. Calculating the magnitude of these effects, however, depends on poorly constrained parameters, such as the composition of ice on the mantle and the size and shape of dust grains

Worsening the problem is that Roberts et al. (2007) are forced to use observations of CO freeze-out to constrain these parameters. As a result, if there is a substantial uncertainty between the amount of depletion as found in our work, then these parameters are even more poorly constrained. They conclude that freeze-out rates are still high (98% for a quiescent core of density 10⁵) even with these effects; however, they admit a large amount of uncertainty. Aikawa et al. (2003) have also done work on surface grain interactions, which could also help in understanding desorption processes. More work is needed in this area.

Velocity structures and dynamical evolution may be tied to chemical evolution in more indirect ways, as pointed out in Lee et al. (2005). They compare chemical evolution in static and collapsing density profiles, and find that in prestellar cores, both give similar results. However, in a dynamic case, the material being pulled into the center of the core may be coming from an outer radius that has less depletion. A better understanding of these cores' dynamics and environments would be very useful. Work by Aikawa et al. (2005), Lee et al. (2004), and Lee (2005) has made much progress on this topic. However, as pointed out by M. Tafalla (2010, private communication), the line widths observed in many of these cores are inconsistent with arguments that cores such as L1521E are simply collapsing too quickly to have time for chemical freeze-out.

5. CONCLUSIONS

We find that all of our cores but one, L1521E, are best fit by models showing significant depletion, consistent with previous work in the field suggesting that unevolved cores are rare. This would imply that the depletion process happens very quickly. The lack of unevolved cores is likely an observational bias, as chemically unevolved cores may also be lower density than those in our sample set. Observations with Herschel and SCUBA2 may find many more lower density cores that may also be chemically unevolved.

We find that assuming an isothermal temperature profile can reproduce the observations, although one needs further information to know the appropriate temperature to use. This is likely due to the small temperature variations seen in dust continuum radiative transfer models. T = 10 K, however, rarely provides a good best match. We find assuming T_d = T_gas gives similar results to energetics calculations with T_d ≠ T_gas. Our work also finds no correlation between the amount of depletion and central density or central temperature, nor do we find correlations between density or temperature at the radius of depletion. It is possible that these cores are evolving at different rates that may be affected by variations in the magnetic field strengths, turbulence, and bounding pressures across Taurus.

Constraining the amount of depletion in these cores can be complicated by degeneracies between temperature and canonical abundance and difficulty in placing strong constraints on f_d. Constraining f_d more strongly is necessary to compare chemical evolution timescales to free-fall timescales, an important step in distinguishing between the relative importance of ambipolar diffusion and gravo-turbulent fragmentation in the evolution of starless cores. Ultimately, observations with better angular resolution, as well as additional understanding of the desorption processes, dynamics, and surrounding environments of these cores, will help complete our understanding of depletion and chemical evolution in starless cores.

We thank Amelia Stutz, Mario Tafalla, Jeong-Eun Lee, George Rieke, Joan Najita, Kaisa Young, Neal Evans, and Craig Kulesa for useful discussions and comments. We are additionally grateful to Floris van der Tak & Michael Hogerheijde for gracious assistance in implementing their code. We also thank Andy Marble and Megan Reiter for help with figures.

This material is based upon work supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Cooperative Agreement No. CAN-02-OSS-02 issued through the Office of Space Science.

C¹⁸O DEPLETION IN STARLESS CORES IN TAURUS

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION