A&A 489, 617-625 (2008)
DOI: 10.1051/0004-6361:200810035

Time-dependent CO depletion during the formation of protoplanetary disks

C. Brinch - R. J. van Weeren - M. R. Hogerheijde

Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

Received 23 April 2008 / Accepted 7 July 2008

Abstract
Context. Understanding the gas abundance distribution is essential when tracing star formation using molecular line observations. Variations in density and temperature conditions can cause gas to freeze-out onto dust grains, and this needs to be taken into account when modeling a collapsing molecular cloud.
Aims. This study aims to provide a realistic estimate of the CO abundance distribution throughout the collapse of a molecular cloud. We derive abundance profiles and synthetic spectral lines that can be compared with observations.
Methods. We use a two-dimensional hydrodynamical simulation of a collapsing cloud and subsequent formation of a protoplanetary disk as input to the chemical calculations. From the resulting abundances, synthetic spectra are calculated using a molecular excitation and radiation transfer code.
Results. We compare three different methods to calculate the abundance of CO. Our models also consider cosmic ray desorption and the effects of an increased CO binding energy. The resulting abundance profiles are compared with observations from the literature and are found to agree well.
Conclusions. The resulting abundance profiles agree well with analytic approximations, and the corresponding line fluxes match the observational data. Our method to calculate abundances in hydrodynamical simulations should help comparisons with observations, and can easily be generalized to include gas-phase reaction networks.

Key words: astrochemistry - hydrodynamics - stars: formation - ISM: molecules - ISM: clouds

1 Introduction

The study of low-mass young stellar objects (YSOs), from early protostellar cores to T Tauri stars surrounded by disks, involves observations of either the thermal emission from dust grains or molecular line emission. While dust emission yields information about the temperature profile (e.g. Draine & Lee 1984), and, through measurements of the shape of the Spectral Energy Distribution, the evolutionary stage (Adams et al. 1987; Lada & Wilking 1984), spectral lines are the only way to constrain the kinematical properties of YSOs. The interpretation of molecular lines is crucial to the understanding of protoplanetary disk formation and the distribution of angular momentum during these stages of star formation (Evans 1999). However, deriving the gas motion from spectral line profiles is complicated by degeneracies between the velocity field topology, inclination of the angular momentum axis, optical depth, and geometrical effects (Brinch et al. 2007).

In addition to these effects, molecular abundances can vary significant. The gas abundance is determined by the ongoing chemistry driven by the co-evolving density and temperature distributions, and to some extent also by the velocity field. The variation in molecular abundances throughout a YSO may have a huge impact on line profiles and therefore it is important to understand the chemistry in interpreting these line profiles. The effect of the chemistry is enhanced for observations of higher resolution. If observations are completed, however, at relatively low resolution ($\gtrsim$10''), the emission is broadened over a significant area of the object and a single average value of abundance is adequate. However, at higher resolution, using sub-millimeter interferometry ($\sim$1''), the emission is averaged over a smaller part of the object and variations in the abundance must be taken into account.

The problem in determining the abundance distribution was addressed before in the literature. Lee et al. (2004) adopted a semi-analytical approach using a self-similar collapse of an isothermal sphere in which they followed ``fluid-elements''. The chemistry was followed in a series of Bonnor-Ebert spheres during the pre-stellar phase and in a self-similar inside-out collapse during the protostellar phase. This model provided a good analytical description of the early stages of star formation, but, lacking any rotational velocity, no disk was formed in this scenario. Aikawa et al. (2005,2003,2001) followed the chemical evolution of contracting pre-stellar clouds using an approach similar to Lee et al. (2004) but using an isothermal cloud collapse only. Their calculations consideried only the early stages of the star formation process. Jørgensen et al. (2005,2002) introduced the ``drop model'', where chemical depletion occurs as a step function when certain temperature and density conditions of any underlying model are met. While this approach is simple yet quite successful, there may be cases where the abundance changes only gradually and a step function is no longer a good approximation. More recently, Tsamis et al. (2008) calculated the HCO⁺, CS, and N₂H⁺ abundances and spectra from a simulation of an inside-out collapsing core.

	Figure 1: The density ( left panel) and temperature ( right panel) contours of four time snapshots of the hydrodynamical simulation. Time progresses from left to right, top to bottom at $t = \{0.0(6), 0.5, 1.5, 2.5\} t_{\rm ff}$. In the right panel, contour lines refer to the dust temperature, while the color scale denotes the gas temperature.
Open with DEXTER

We investigate the detailed evolution of the gas-phase and solid state abundance of CO in a collapsing, rotating core. We base our study on a hydrodynamical simulation, where we follow the CO freeze-out and evaporation for a number of test particles that flow with the gas. The layout of the paper is as follows: in Sect. 2, we present our numerical simulations and equations used to solve the CO freeze-out chemistry. In Sect. 3, we present three different chemical models and compare our results with observations and the results of previous papers. Sections 4 and 5 provide discussion and a summary, respectively.

  
2 Tracing chemistry during disk formation

2.1 The physical model

Instead of using an analytic description of the gravitational collapse and subsequent disk formation, we use a grid-based two-dimensional axis-symmetric hydrodynamical scheme. The code used was described by Yorke & Bodenheimer (1999) and we adopt in particular the J-type model described in that paper.

The initial conditions of this model are a $1~M_{\odot}$ isothermal sphere with a temperature of 10 K and a power-law density slope of $\rho \propto r^{-2}$. The cloud has an outer radius of $1\times 10^{15}$ m (6667 AU), and an initial solid-body rotational perturbation of $1\times 10^{-13}$ s^-1. The model evolves under the action of gravity, while the temperature is solved self-consistently using an approximate radiation transfer method. Angular momentum is transferred by artificial viscosity using an $\alpha$-prescription (Shakura & Sunyaev 1973).

The age of the simulated system is described in terms of the initial free-fall timescale,

$$t_{\rm ff}=\sqrt{\frac{\pi^2 R^3}{8GM},}$$ (1)

where M is the mass of the cloud and R is the radius. For the initial conditions in our particular simulation one free-fall time scale is $1\times 10^5$ yr. We set the duration of the simulation to be about 2.6  $t_{\rm ff}$ at which point the simulation terminates due to extreme velocities close to the center. Figure 1 shows four time snapshots at characteristic ages. We use the same four snapshots throughout this paper, at times of 0.0, 0.5, 1.5, and 2.5 $t_{\rm ff}$ (with the exception of Fig. 1, where we use a few snapshots later than t = 0.0 to enable the temperature to evolve slightly from the isothermal initial condition).

The luminosity is given by the sum of the intrinsic stellar luminosity and the accretion luminosity

$$L_{\rm acc} = \frac{3}{4}\frac{GM_*\dot{M}}{R_*},$$ (2)

where M_* is the mass of the star, R_* is the radius of the star, and $\dot{M}$ is the mass flux onto the star. In Fig. 2, the total luminosity is plotted as a function of free-fall time. After an initial increase, the luminosity declines before increasing sharply again, and then continues to declines slowly for the remainder of the simulation. As described by Yorke & Bodenheimer (1999), the smooth increase and drop around 0.5 $t_{\rm ff}$ is controlled by the accretion luminosity, as low angular momentum material fall directly onto the star. The sharp increase in luminosity corresponds to the formation of an equilibrium disk and after then, the intrinsic stellar luminosity is dominated by the total luminosity.

The hydrodynamical simulation does not include any chemistry and we are therefore unable to relate the state of the molecules, whether in a gas phase or locked in an ice matrix, to the hydrodynamics. We therefore determine the abundance by post-processing the result of the hydrodynamical calculations. We follow the chemistry by populating the computational domain of the first snapshot with trace particles. These particles are massless and do not interact with each other in any way. The trace particles are positioned evenly in the radial direction and $9\times 10^5$ particles are used. Every trace particle represents the local molecular environment and carries information about the state of the molecules (ice or gas). The particles are then allowed to follow the flow, predetermined by the hydrodynamical calculations, and the state of the particles is updated as temperature and density conditions change throughout the hydrodynamical simulation. Finally, the state of the trace particles are mapped back onto the hydro-grid at each output time step, so that a complete history from t=0.0   $t_{\rm ff}$ to t=2.5   $t_{\rm ff}$ of the CO abundance as function of R and z is linked to the density and temperature. Taking the CO distribution into account, we can use each of these time snapshots as input models for our two-dimensional line excitation and radiative transfer code RATRAN (Hogerheijde & van der Tak 2000) to predict line profiles of CO and its optically thin isotopologue C¹⁸O.

Before proceeding with the freeze-out calculations, it is interesting to consider the dynamics of the trace particles to develop an understanding of the environments that the particles traverse. By a simple radial color coding of the particles, we can effectively visualize the flow within the hydrodynamical simulation and track how material is accreted onto the disk. Four snapshots of the particle distribution are shown in Fig. 3 where the shading of the particles refers to their initial distance from the center. The disk is seen to be layered vertically rather than radially, which is somewhat counterintuitive. The vertical layering occurs because material connects to the disk from above. As the disk spreads viscously due to conservation of angular momentum, a shock front propagates outwards, pushing aside infalling material. The only possible path for the particles is therefore upwards and over the disk lobe until they can rain down onto the disk at smaller radii. This rolling motion is also reflected in the strong vertical mixing inside the outer disk. However, this behavior may be an artifact of the two-dimensional nature of the hydrodynamical code. In a 3D simulation, particles should also dissipate in the azimuthal direction, which would probably change their motion.

	Figure 2: The time evolution in the luminosity of the central source in the hydrodynamical simulation. This figure compares with Fig. 7 in Yorke & Bodenheimer (1999).
Open with DEXTER

	Figure 3: Trace particle motion in the hydrodynamical simulation. The particles are color-coded according to their initial radial position. Contour lines describe the density. The vertical color bar shows the particles initial distance from the center. The four panels ( from the top-left to bottom-right) show the particle positions at t=0.0, 0.5, 1.5, and 2.5  $t_{\rm ff}$.
Open with DEXTER

2.2 Freeze-out and evaporation

Molecular gas abundances in general depend on reaction rates and phase transitions. However, under the physical conditions present in our model, only the latter has any significant impact on the CO abundance. The state of CO is governed by rate equations. The two main mechanisms are depletion, where gas phase molecules freeze-out onto grains, and desorption, where solid state molecules evaporate into the gas phase.

The depletion rate $\lambda$ used in this paper is given by the equation (e.g. Charnley et al. 2001),

 

$$\lambda = \pi a^2 \sqrt{\frac{8kT_{\rm g}}{\pi m_{\rm CO}}} n_{\rm grain},$$ (3)

where $T_{\rm g}$ is the gas temperature, $m_{\rm CO}$ is the mass of a CO molecule, and $n_{\rm grain}$ is the grain number density, assuming a grain abundance of $1.33\times 10^{-12}$ relative to the hydrogen nucleon density. We assume a mean grain size of 0.1 $\mu$m. We also assume a sticking probability of unity.

A similar equation can be written for the desorption of molecules. The thermal desorption rate is given by Watson & Salpeter (1972),

 

$$\xi_{\rm th} = \nu(X) {\rm e}^{\left(\frac{-E_{\rm b}(X)}{k_{\rm b} T_{\rm d}}\right )},$$ (4)

where $k_{\rm b}$ is the Boltzmann constant, $T_{\rm d}$ is the dust temperature, $\nu(X)$is the vibrational frequency of X in its binding site (Hasegawa et al. 1992), and $E_{\rm b}(X)$ is the binding energy of species X onto the dust grain. In this paper we assume a default binding energy of 960 K, but other values, corresponding to mantles that do not consist of pure CO ice, are explored in Sect. 3.4.

A second desorption mechanism that is taken into account is cosmic-ray induced desorption. This mechanism was first proposed by Watson & Salpeter (1972). Energetic nuclei might eject molecules from grain surfaces by either raising the temperature of the entire grain or by spot heating close to the impact site. The formulation of Hasegawa & Herbst (1993) is used to calculate the cosmic ray desorption rate,

 

$$\xi_{\rm cr} = 3.16\times 10^{-19} \xi_{\rm th} \Big \arrowvert_{T_{\rm d} = 70~{\rm K}}.$$ (5)

Cosmic ray desorption can be added to the thermal desorption, effectively preventing a depletion of 100% under any condition.

By considering the depletion and desorption rates and defining $f_{\rm d}$ to be the normalized fractional depletion, we can write an equation for $f_{\rm d}$,

 

$$\frac{{\rm d}f_{\rm d}(\vec{r},t)}{{\rm d}t} = \lambda(\vec{r},t) f_{\rm d}(\vec{r},t)-\xi(\vec{r},t)\big (1-f_{\rm d}(\vec{r},t)\big ),$$ (6)

in which we implicitly adopt the boundary condition that the sum of CO molecules in the solid state and gas phase is constant.

Equation (6) is a general one and the rate functions can be changed readily to describe other reactions. In principle, van Weeren (2007) shows that an entire set of these equations can be solved simultaneously using gas phase reaction rates such as the UMIST database (Woodall et al. 2007), although computational demand becomes an issue for a full network.

   
3 Results

To obtain abundance distributions, we evaluate Eq. (6) using the gas temperature $T_{\rm g}$, the dust temperature $T_{\rm d}$, and the mass density, assuming that H₂ is the main mass reservoir and that the gas-to-dust ratio is 100:1, from the hydrodynamical simulation. We assume that the cloud has been in hydrostatic equilibrium prior to collapse for $3\times 10^5~{\rm yr} ~({=}1\times10^{13}$ s). The density profile before collapse is given by a spherical power-law of uniform temperature described in Sect. 2. Starting with a pure gas phase abundance of $2\times 10^{-4}$ molecules relative to H₂, we allow the abundance to evolve for $3\times 10^5$ yr with a constant temperature of 10 K and a static density profile. The new CO abundance then serves as the initial condition for the collapse, i.e. the start of the hydrodynamical simulation, which defines t = 0. In the following, we present three different ways to obtain the abundance profiles.

3.1 Model 1: drop abundance

This approach is denoted the drop abundance model, because it follows directly from the method presented by Jørgensen et al. (2005). In this approach, we do not solve Eq. (6), but consider the rate Eqs. (3) and (4). The rate equations are inverted to obtain a characteristic timescale, and, at any given time in the simulation, this timescale is evaluated based on the current density and temperature conditions. If the freeze-out timescale is shorter than the evaporation time scale, the gas phase abundance drops to a preset level (which is sufficiently low, i.e. several orders of magnitude) but if the actual age of the cloud is shorter than the freeze-out timex scale, the gas phase abundance is kept at its initial value. Using this approach, the resulting profile is a step function where all CO is either in the gas phase or in the solid state. The movement of the trace particles does not affect the result since the abundance is determined only from the current local conditions. By ensuring that the abundance decreases by only two or three orders of magnitude and not to zero, the effect of cosmic ray desorption is modeled implicitly, even though conditions for freeze-out are met.

	Figure 4: This figure shows the abundance distribution of Model 1. The top panel shows the radial abundance profile through the mid-plane of the disk. The bottom panel shows the spatial abundance distribution. The shaded area is where molecules are frozen out. The unit on the y-axis of the top panel is fractional abundance.
Open with DEXTER

The result is shown in Fig. 4. The top panel shows the radial abundance profile along the disk mid-plane (z=0), while the lower panel shows the two dimensional distribution. The result of this drop abundance model supports the idea of an evolution in the CO abundances from a pre-stellar core to a Class I object, as proposed by Jørgensen et al. (2005). The maximum size of our depleted zone is $\sim$6000 AU, just before collapse sets in, and shrinks to $\sim$500 AU (when averaged in radial and polar directions) as the cloud collapses and the disk is formed. The inner evaporated zone has a radius of $\sim$1000 AU, in agreement with values derived in Jørgensen et al. (2005).

3.2 Model 2: variable freeze-out

The second method is similar to Model 1 except that instead of using timescales, we use the actual rate equations that allow the trace particles to be in a mixed state. The values for $\lambda$ and $\xi$ are calculated for all $\vec{r}$ using the current $n(\vec{r})$ and $T(\vec{r})$, and Eq. (6) is solved to obtain the fraction of gas phase CO to solid state CO for each trace particle position. In this model, we still do not consider the history of the trace particles and the motion of the particles is therefore irrelevant and the time dependence of $\xi$ and $\lambda$ in Eq. (6) becomes negligible. Plots of the abundance distributions are shown in Fig. 5. The profiles of Model 2 are similar to those from Model 1, except for the gradual change from pure gas phase to pure solid state. The inner edge in particular is described well by the step function of Model 1 because of the exponential factor in the evaporation, described by Eq. (4). The main difference between the two models is the outer edge, and in particular at earlier times, where the gradient is more shallow.

	Figure 5: Abundance distribution of Model 2. Otherwise identical to 4.
Open with DEXTER

3.3 Model 3: dynamical evolution

Our final model takes the motion of the trace particles and their chemical history into account. Freeze-out and evaporation rates are calculated dynamically for each time step and the composition of the particles is used as the initial conditions for Eq. (6). For example, the depletion in a certain area represented by a trace particle depends on where this particle has originated, which we know because the trace particles follow the flow calculated in the hydrodynamical simulation. Two trace particles may pass by the same region at the same time and experience the same desorption and depletion conditions. However, they may have originated in different areas; for example, one particle may have originated in a warmer region and the other may have travelled from a colder region, and they will then leave the location with different compositions.

Figure 6 shows the abundance distributions for this model. Surprisingly, the profiles do not differ significantly from those of Model 2. A difference only occurs when the dynamical timescale is comparable to the chemical timescale. In this case, a gas phase trace particle may move out of a depletion zone more rapidly than depletion can make any significant change to its composition, or vice versa. The evaporation front around the disk is also sharper with respect to Model 2, which causes Model 3 to appear almost like a step function in the final part of the simulation.

	Figure 6: Abundance distribution of Model 3. Otherwise similar to 4.
Open with DEXTER

  
3.4 Cosmic ray desorption and increased binding energy

We now investigate the effects of excluding the cosmic ray desorption term (Eq. (5)) to see what effect this mechanism has on the abundances. We also increase the CO binding energy $E_{\rm b}$ in Eq. (4) to emulate different ice mantle compositions. In the following, we use Model 3 as a template.

Cosmic ray desorption ensures that depletion can never reach 100% as long as cosmic rays are able to penetrate the cloud and remove molecules from the ice matrix. If no other desorption mechanisms are present, dark cold molecular clouds will have a high degree of depletion since thermal desorption is largely ineffective at low temperatures. However, Caselli et al. (1999) found no evidence for high depletion factors and therefore some other desorption mechanism (such as for instance cosmic ray desorption) must play a certain role in these environments. However, to observe exactly how much evaporation cosmic rays account for, we calculated Model 3 without the cosmic ray desorption term. Figure 7 shows these abundance distributions. The profiles are similar to Model 3, but the amount of depletion is, as expected, about three orders of magnitude higher.

	Figure 7: Abundance distribution of Model 3 without cosmic ray desorption and the default binding energy of 960 K. Otherwise similar to 4.
Open with DEXTER

Including cosmic ray desorption obviously has a significant impact on the amount of gas phase molecules in the cloud, raising the gas-phase abundance substantially. However, the effect of cosmic ray desorption can be balanced by a higher binding energy between the molecule and grain surface. For example, raising the binding energy from 960 K to 1740 K compensates almost entirely for the effect of cosmic ray desorption, when the depletion is averaged over the entire model.

	Figure 8: Abundance distribution of Model 3, including cosmic rays, but with a CO binding energy of 1740 K. Otherwise similar to 4.
Open with DEXTER

An increased binding energy applies when the CO ice mantle is mixed with other species, typically water ice or NH₃ ice (Fraser et al. 2004). A spread in binding energies also occurs when the ice surfaces are not even, e.g. for differently layered ices, varying mantle thicknesses and rough surfaces. If the grains are more fractal, with cavities that can trap CO molecules, it becomes more difficult to desorb the molecules, which effectively raises the binding energy. Laboratory experiments of Collings et al. (2004) detected multiple peaks in the CO desorption spectrum, implying a range of binding energies. In Fig. 8, the innermost evaporation zone is smaller by a significant margin throughout the simulation, when the binding energy is increased, in comparison to when cosmic rays are not present; this implies that the two effects can be discriminated by high-resolution observations.

3.5 Comparison to observations

We compare our results with observations of CO abundances and depletion factors. The fractional depletion $f_{\rm d}$ of CO for the various models are shown in Fig. 9. The highest depletion factors occur at the end of the pre-stellar core phase just before collapse. Observations show that depletion factors decrease with decreasing envelope mass (Jørgensen et al. 2002). We observe the depletion decrease after the core begins to collapse due to CO evaporation close to the protostar. At the same time, the envelope density decreases, which also lowers the global fractional depletion defined to be the ratio of the total number of gas phase CO molecules to H₂ molecules in the entire model.

Freeze-out in the disk is most pronounced for Model 3, includes cosmic ray desorption and a high CO binding energy. The same model but with a binding energy of 960 K provides a low depletion factor at the onset of the collapse. Although low values of $f_{\rm d}$ have been measured (e.g. Bacmann et al. 2002), most observations point toward a higher value, $f_{\rm d} \geq 10$, in pre-stellar cores (e.g. Bergin et al. 2002; Jessop & Ward-Thompson 2001; Crapsi et al. 2004). The low value of $f_{\rm d}$ that we find implies that either the binding energy is higher than 960 K or that we overestimate the cosmic ray desorption rate. Shen et al. (2004) demonstrated that the importance of cosmic ray heating was overestimated by previous authors. A third possibility is that CO can be converted into other species, resulting in a lower abundance, but since we only model CO we cannot quantify this effect.

As mentioned above, a correlation was found between the amount of depletion and envelope mass. The CO abundance declines when going from a pre-stellar core, through the Class 0 phase, to a Class I object. We compare our results with the data presented by Jørgensen et al. (2002), who derived global envelope abundances for 18 sources. Their result is shown in Fig. 10, where open symbols are Class I objects, closed diamonds are Class 0 objects, and the filled squares are pre-stellar cores. Superposed on the data are curves based on our models. In the linear approximation, our models follow the same trend. The observations show no evidence of freeze-out in the Class I objects but the number of sources of this type is small and the derived envelope masses and abundances are uncertain. We also expect some intrinsic scatter due to the unique environments and characteristics of each source. The drop abundance model and Model 3, including cosmic rays and a low binding energy, have too high CO abundances (or too little depletion). Unfortunately, our simulation begins with only one solar mass and we are therefore unable to study the more massive Class 0 and pre-stellar cores.

	Figure 9: The global fractional depletion in the various models. The small bump in the curves around 0.5 $t_{\rm ff}$ is due to a change in the luminosity output of the central source as the disk is formed.
Open with DEXTER

	Figure 10: In this plot, we show the abundance as a function of envelope mass. Also shown in this plot are data points taken from Jørgensen et al. (2005,2002). The open symbols are Class I objects, the filled diamonds are Class 0 sources and the filled squares are pre-stellar cores.
Open with DEXTER

We can also calculate gas column densities based on our abundance distribution models. Vertically integrated column densities, $N_{\rm {CO}} = \int n_{\rm {CO}}(z) {\rm d}z$, are shown at four different free-fall times in Fig. 11. The models differ mainly in the position of the evaporation front due to the different binding energies used. The column density profiles at t = 0 can be compared with those presented by Aikawa et al. (2005,2003,2001). Although our underlying density model is different, the column densities are comparable, differing by an order of magnitude or less. Aikawa et al. (2001) measured flat profiles toward the core center or slight drops in the column within $\sim$5000 AU, which is similar to the results for our models. The only exception is our drop abundance model (Model 1), which has an increasing column density toward the center. The drop occurs at $\sim$6000 AU, and within a 6000 AU radius we have a constant abundance. The $n_{\rm {H}_2}$ density decreases with the inverse square of the radius, which means that the column density of CO increases toward the core center. This disagrees with observations (e.g. Tafalla et al. 2004) where flat or decreasing integrated intensity profiles are found at core centers and can be interpreted as a drop in the column density.

	Figure 11: CO gas column densities of the models discussed in this paper. These profiles are observed through a pencil beam, and need to be convolved with an appropriate beam to be directly compared with observations.
Open with DEXTER

Column densities decrease by about three orders of magnitude from the center toward the freeze-out zone in the disk, comparable to Aikawa et al. (1996), although column densities in our results are in general more than one order of magnitude higher. This can be explained by the differences in the adopted disk masses. A typical value adopted for disk masses is 10^-2 $~M_{\odot}$, while the disk in our simulations increases to an unrealistically high mass of about 0.4 $M_{\odot}$. Our column density distributions are also more complicated because we do not use a simple analytic description of a non-evolving disk. Aikawa et al. (1996) measured drops in column density due to freeze-out in the disk around 200-300 AU, similar to our results for 100-500 AU, depending on the model used.

In a similar way, we can calculate the column densities for the CO bound to ice. Plots of the ice column densities are shown in Fig. 12. Before collapse, the column density of CO ice increase inwards due to the shorter freeze-out timescales. When the central temperature rises, the column density flattens and drops, while in the protoplanetary disk the column density increases sharply. Overall, the column density drops in the envelope because of the decreasing density (which also results in a lower freeze-out rate). It is difficult to compare these results with observations since cold (unrelated) clouds, where CO is frozen out, often are located in front of protostars. Observed column densities of CO ice are often approximately 10¹⁸ cm^-2 in reasonable agreement with our results.

	Figure 12: Same as Fig. 11 but for CO ice column densities.
Open with DEXTER

	Figure 13: Time series of CO 3-2 spectra. The top three panels show the three different abundance models. The three lower panels show model 3, including cosmic ray desorption using three different binding energies. Time is shown by increasing line thickness and the four lines correspond to the four adopted time snapshots (0.0, 0.5, 1.5, and 2.5  $t_{\rm ff}$). Each line has been offset by 0.5 K.
Open with DEXTER

3.6 Emission lines

Simulated spectral lines are a direct way to test our models because they can be compared directly to observed spectra. As mentioned in Sect. 2, we use the hydrodynamical time snapshots including the abundance distribution mapped onto the original grid as input models for our two-dimensional molecular excitation code RATRAN. This code solves the level population using an accelerated Monte Carlo method, which is then ray-traced to create frequency dependent intensity maps. After beam convolution, spectral profiles can be extracted from these and directly compared with observations. We add a mean field of 0.2 km s^-1to the snapshots during the radiation transfer calculations to emulate the presence of micro-turbulence.

An in-depth analysis of spectra based on the hydrodynamical simulation is given in a separate paper (Brinch et al. 2008) and only example spectra for the abundance models derived in this paper are presented here. Figure 13 shows the CO J= 3-2 transition of the four time-snapshots used in previous figures. The top three panels correspond to the three abundance models including cosmic ray desorption, while the lower three panels show spectra based on model 3 without cosmic ray desorption and with increased binding energies. Figure 14 has a similar layout, but for the less optically thick isotopologue C¹⁸O. For comparison, CO and C¹⁸O spectra, calculated with a constant abundance profile are shown in Fig. 15. All spectra are viewed at 90$^\circ$ inclination (edge-on; similar to the orientation of the model in Fig. 1), and they were convolved with a 10'' beam, a typical value for single-dish sub-millimeter telescopes. A distance of 140 pc (Taurus star forming region) is adopted. We use an isotopic ¹⁶O:¹⁸O ratio of 560 (Wilson & Rood 1994) to calculate all C¹⁸O spectra.

	Figure 14: Similar to Fig. 13 but for C18O 3-2.
Open with DEXTER

	Figure 15: Time series of CO and C18O spectra with a constant abundance profile. Time is represented by increasing line thickness.
Open with DEXTER

We observe that the spectra in the top row of Fig. 13 are similar to each other and also to the constant abundance profiles in Fig. 15. The reason is that CO becomes optically thick before the abundance declines and is therefore independent of the underlying abundance model. We are therefore unable to use CO spectra to distinguish between the three models and cannot determine the depletion factor using CO. An exclusion of cosmic rays or increase in the binding energy lowers the intensity of the spectra, and in this case CO can be used to constrain the models.

The line profile of the earliest snapshot (the thinnest line) is to have a boxy shape with a flat line center. This is entirely due to the hydrodynamical model, which has not yet evolved from its initial isothermal state in the region where the line get optically thick. Basically we probe the same temperature through a range of velocities, which produce the flat top. The line width in the earliest snapshot in Fig. 13 is also considerably narrower than for the other lines. This is because the velocities in the model have not yet had the time to increase from the initial solid body rotation speed.

The C¹⁸O spectra however differ considerably from those derived for a constant abundance model. When the three abundance models are compared, they also differ, although to a lesser extent, from each other. This is because the lines are optically thin and a drop in abundance has a direct impact on the spectral line, no matter at what radius the drop occurs.

The spectra shown here all simulate single-dish measurements, where the emission of the source has been smeared out by a relatively large beam. If an interferometer is used instead, beam sizes smaller than 1'' can be achieved, which allows us to map the line-of-sight column densities. With these observations, it should be possible to map abundance distributions of optically thin species and constrain our models.

   
4 Discussion

The abundance profiles presented here depend on the adopted hydrodynamical model. While we believe that our results are in general applicable, details, such as the precise timescale and absolute values of the fractional depletion, depend on the initial conditions of our simulations such as mass and angular momentum and the particular hydrodynamical scheme used.

A main concern of the hydrodynamical simulation is the internal radiation transfer module used to calculate the temperature structure of the simulation. This is done with an approximate method known to overestimate the dust temperature by a few Kelvins compared with calculations by more accurate continuum radiation transfer codes (Visser, priv. comm.). While a few Kelvin may have little or no effect on the hydrodynamics and therefore the evolution of the cloud, it can certainly affect the evaporation front, which has a strong temperature dependence. While this will alter the detailed behavior of the abundances, it will not change the general trends discussed in this paper.

Although we have considered only the CO gas-phase abundance, as mentioned in Sect. 2, our method can be easily extended to include other species and even surface reactions. However, attempting to calculate the entire chemical network is considerably more computationally demanding. Where the CO models can be completed in a few hours on a standard desktop computer using $9\times 10^5$ trace particles, a test simulation with the entire chemical network indicated that it would take a factor of 100 longer to complete using a factor of 1000 fewer trace particles. Nevertheless, such large-scale chemical models present a powerful way to obtain realistic abundance profiles for hydrodynamical simulations which instruments such as the Atacama Large Millimeter Array (ALMA) will be able to test at high spatial resolution.

   
5 Conclusion

We have investigated different scenarios for the depletion of CO in a collapsing cloud. We have used a hydrodynamical code to describe the dynamical evolution of the cloud and the chemistry has been solved using a trace particle approach.

Of the various models shown in Sect. 3, that which excludes cosmic ray desorption and those with higher binding energies provide the most extreme results. The drop abundance model provides a good approximation to Model 2 and 3. Therefore, if it is possible to describe observations, using a drop profile, as achieved by Jørgensen et al. (2005), it should be possible to make a hydrodynamical simulation with customized initial conditions, so that it reproduces this profile and places an observational constraint on the hydrodynamics. In cases where an analytical model has been used, the drop abundance approximation has also proven to be good.

We conclude that the freeze-out chemistry has little effect on optically thick spectral lines, especially when observed at low resolution. It is not unreasonable to model these lines using a constant abundance model regardless of the depletion. However, we should be careful when using optically thin lines to determine the average gas phase abundance that should be used in a constant abundance model, since the average abundancesderived from these lines depend on the freeze-out chemistry. The maximum discrepancy between the flat abundance profile model (Fig. 15) and the freeze-out models shown in Fig. 14 is a factor of three in the integrated intensity. Such a change in the average abundance does not affect the optically thick lines a lot, so depending on the accuracy required, it may not be unreasonable to use a flat abundance distribution to model data in which freeze-out is present.

Finally, the effect of chemical depletion is enhanced considerably when the resolution of the observations increases, implying that it becomes more important to model the abundance profiles accurately. This will be paticularly important for future ALMA observations, which will have a resolution of $\sim$0.01'', for which the spatial detail will be significantly higher and precise abundance models will be crucial for the accurate interpretation of the data.

Acknowledgements

The authors would like to thank the anonymous referee for thorough and useful comments. C.B. is partially supported by the European Commission through the FP6 - Marie Curie Early Stage Researcher Training programme. The research of M.R.H. is supported through a VIDI grant from the Netherlands Organization for Scientific Research.

References

• Adams, F. C., Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788 [NASA ADS] [CrossRef]
• Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684 [NASA ADS] [CrossRef] (In the text)
• Aikawa, Y., Ohashi, N., Inutsuka, S.-I., Herbst, E., & Takakuwa, S. 2001, ApJ, 552, 639 [NASA ADS] [CrossRef]
• Aikawa, Y., Ohashi, N., & Herbst, E. 2003, ApJ, 593, 906 [NASA ADS] [CrossRef]
• Aikawa, Y., Herbst, E., Roberts, H., & Caselli, P. 2005, ApJ, 620, 330 [NASA ADS] [CrossRef]
• Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2002, A&A, 389, L6 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101 [NASA ADS] [CrossRef]
• Brinch, C., Crapsi, A., Hogerheijde, M. R., & Jørgensen, J. K. 2007, A&A, 461, 1037 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Brinch, C., Hogerheijde, M. R., & Richling, S. 2008, A&A, 489, 607 [CrossRef] [EDP Sciences] (In the text)
• Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJ, 523, L165 [NASA ADS] [CrossRef] (In the text)
• Charnley, S. B., Rodgers, S. D., & Ehrenfreund, P. 2001, A&A, 378, 1024 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133 [NASA ADS] [CrossRef] (In the text)
• Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2004, A&A, 420, 957 [NASA ADS] [CrossRef] [EDP Sciences]
• Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] (In the text)
• Evans, II, N. J. 1999, ARA&A, 37, 311 [NASA ADS] [CrossRef] (In the text)
• Fraser, H. J., Collings, M. P., Dever, J. W., & McCoustra, M. R. S. 2004, MNRAS, 353, 59 [NASA ADS] [CrossRef] (In the text)
• Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] (In the text)
• Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] (In the text)
• Hogerheijde, M. R., & van der Tak, F. F. S. 2000, A&A, 362, 697 [NASA ADS] (In the text)
• Jessop, N. E., & Ward-Thompson, D. 2001, MNRAS, 323, 1025 [NASA ADS] [CrossRef]
• Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908 [NASA ADS] [CrossRef] [EDP Sciences]
• Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2005, A&A, 435, 177 [NASA ADS] [CrossRef] [EDP Sciences]
• Lada, C. J., & Wilking, B. A. 1984, ApJ, 287, 610 [NASA ADS] [CrossRef]
• Lee, J.-E., Bergin, E. A., & Evans, II, N. J. 2004, ApJ, 617, 360 [NASA ADS] [CrossRef] (In the text)
• Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
• Shen, C. J., Greenberg, J. M., Schutte, W. A., & van Dishoeck, E. F. 2004, A&A, 415, 203 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, A&A, 416, 191 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Tsamis, Y. G., Rawlings, J. M. C., Yates, J. A., & Viti, S. 2008, ArXiv e-prints, 803 (In the text)
• van Weeren, R. J. 2007, Master's thesis, Leiden University (In the text)
• Watson, W. D., & Salpeter, E. E. 1972, ApJ, 175, 659 [NASA ADS] [CrossRef] (In the text)
• Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [NASA ADS] [CrossRef] (In the text)
• Woodall, J., Agúndez, M., Markwick-Kemper, A. J., & Millar, T. J. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Yorke, H. W., & Bodenheimer, P. 1999, ApJ, 525, 330 [NASA ADS] [CrossRef] (In the text)