© ESO, 2011

1. Introduction

The formation of a low-mass protostar out of a cold molecular cloud is accompanied by large-scale changes in the chemical composition of the constituent gas and dust. Pre-stellar cloud cores are cold (~10 K), moderately dense (~10⁴–10⁶ cm^-3), and irradiated only by the ambient interstellar radiation field (Shu et al. 1987; Di Francesco et al. 2007; Bergin & Tafalla 2007). The material heats up as the core collapses, and the inner few hundred AU flatten out to form a circumstellar disk with densities that are orders of magnitude higher than in the parent cloud (Dullemond et al. 2007). Meanwhile, the protostar infuses the disk and remnant envelope with large fluxes of ultraviolet and X-ray photons. The chemical changes arising from these evolving physical conditions have been analysed with one-dimensional spherical models (e.g., Aikawa et al. 2008). However, two-dimensional models are required to properly describe the formation of the circumstellar disk and to analyse the full chemical evolution from core to disk.

This paper is the second in a series of publications aiming to model the chemical evolution from pre-stellar cores to circumstellar disks in two dimensions. The first paper (Visser et al. 2009b, hereafter Paper I) contained a detailed description of our semi-analytical model and an analysis of the gas and ice abundances of carbon monoxide (CO) and water (H₂O). Most CO was found to evaporate during the infall phase and to freeze out again in those parts of the disk that are colder than  ~20 K. The higher binding energy of H₂O keeps it in solid form at all times, except within  ~ 10 AU of the star. Based on the time that the infalling material spends at dust temperatures between 20 and 40 K, first-generation complex organic species were predicted to form abundantly on the grain surfaces according to the scenario of Garrod & Herbst (2006) and Garrod et al. (2008).

The current paper extends the chemical analysis to a full gas-phase network, including freeze-out onto and evaporation from dust grains, as well as basic grain-surface hydrogenation reactions. Combining semi-analytical density and velocity structures with detailed temperature profiles from full radiative transfer calculations, our aim is to bridge the gap between 1D chemical models of collapsing cores and 1+1D or 2D chemical models of T Tauri and Herbig Ae/Be disks (Bergin et al. 2007). One of the key questions is whether the chemical composition of such disks results mainly from chemical processing during the collapse or whether it is determined by in situ processing after the disk has formed.

As reviewed by Di Francesco et al. (2007) and Bergin & Tafalla (2007), the chemistry of pre-stellar cores is well understood. Because of the low temperatures and the moderately high densities, many molecules are depleted from the gas by freezing out onto the cold dust grains. The main ice constituent is H₂O, showing abundances of  ~ 10^-4relative to gas-phase H₂ (Tielens et al. 1991; Pontoppidan et al. 2004). Other abundant ices are CO₂ and CO (Gibb et al. 2004; Öberg et al. 2011a). Correspondingly, the observed gas-phase abundances of H₂O, CO and CO₂ in pre-stellar cores are low (Snell et al. 2000; Bacmann et al. 2002; Caselli et al. 2010). Nitrogen-bearing species like N₂ and NH₃ are generally less depleted than carbon- and oxygen-bearing ones (Tafalla et al. 2004), partly because they require a longer time to be formed in the gas and therefore have not yet had a chance to freeze out (Aikawa et al. 2001; Di Francesco et al. 2007), and partly because CO is no longer present as a destructive agent (Rawlings et al. 1992). The observed depletion factors are well reproduced with 1D chemical models (Bergin & Langer 1997; Lee et al. 2004).

The collapse phase is initially characterised by a gradual warm-up of the material, resulting in the evaporation of the ices according to their respective binding energies (van Dishoeck & Blake 1998; Jørgensen et al. 2004). The higher temperatures also drive a rich chemistry, especially if it gets warm enough to evaporate H₂O and organic species like CH₃OH and HCOOCH₃ (Charnley et al. 1992). Spherical models of the chemical evolution during the collapse phase (e.g., Rodgers & Charnley 2003; Lee et al. 2004; Garrod & Herbst 2006; Aikawa et al. 2008; Garrod et al. 2008) are successful at explaining the observed abundances at scales of several thousand AU, where the envelope is still close to spherically symmetric, but they cannot make the transition from the 1D spherically symmetric envelope to the 2D axisymmetric circumstellar disk. Recently, van Weeren et al. (2009) followed the chemical evolution within the framework of a 2D hydrodynamical simulation and obtained a reasonable match with observations of protostars. Even though their primary focus was still on the envelope and not on the disk, they did show how important it is to treat the chemical evolution during low-mass star formation in more than one dimension.

Once the phase of active accretion from the envelope comes to an end, the circumstellar disk settles into the comparatively static T Tauri or Herbig Ae/Be phase. So far only some simple molecules have been detected in disks (Dutrey et al. 1997; Kastner et al. 1997; Thi et al. 2004; Öberg et al. 2010, 2011b), but the inventory is expected to grow now that the Atacama Large Millimeter Array (ALMA) is becoming operational. Chemically, disks can be divided into three layers: a cold zone near the midplane, a warm molecular layer at intermediate altitudes, and a photon-dominated region at the surface (Bergin et al. 2007). However, it is unknown to what extent the chemical composition of the disk is affected by the chemical evolution in the collapsing core at earlier times, or how much chemical processing has been experienced by material ending up in planetary and cometary building blocks.

This paper aims to provide an important step towards filling in this missing link by following the chemical evolution all the way from a pre-stellar cloud core to a circumstellar disk in two spatial dimensions. The physical and chemical models are described in Sects. 2 and 3. The chemistry during the pre-collapse phase is briefly reviewed in Sect. 4, followed by an extensive discussion of the collapse-phase chemistry in Sect. 5. The abundances resulting from the collapse phase are compared to in situ processing in a static disk in Sect. 6. Sections 7 and 8 discuss the implications of our results for the origin of comets and address some caveats in the model. Conclusions are drawn in Sect. 9.

2. Collapse model

2.1. Step-wise summary

Our semi-analytical collapse model is described in detail in Paper I and summarised in Fig. 1. The model starts with a singular isothermal sphere characterised by a total mass M₀, an effective sound speed c_s, and a uniform rotation rate Ω₀. As soon as the collapse starts, at t = 0, the rotation causes the infalling material to be deflected towards the equatorial midplane. This breaks the spherical symmetry, so the entire model is run as a two-dimensional axisymmetric system. The 2D density and velocity profiles follow the solutions of Shu (1977), Cassen & Moosman (1981) and Terebey et al. (1984) for an inside-out collapse with rotation. After the disk is first formed at the midplane, it evolves by ongoing accretion from the collapsing core and by viscous spreading to conserve angular momentum (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974).

Fig. 1 Step-wise summary of our 2D axisymmetric collapse model. Steps 2 and 4 are semi-analytical, while steps 3 and 5 consist of detailed numerical simulations.

Taking the 2D density profiles from step 2, and adopting the appropriate size and luminosity for the protostar (Adams & Shu 1986; Young & Evans 2005), the next step consists of computing the dust temperature at a number of time steps. This is done with the radiative transfer code RADMC (Dullemond & Dominik 2004), which takes a 2D axisymmetric density profile but follows photons in all three dimensions. RADMC also computes the full radiation spectrum at each point in the axisymmetric disk and remnant envelope, as required for the photon-driven reactions in our chemical network (Sects. 2.3 and 3.1). The gas temperature is set equal to the dust temperature throughout the disk and the envelope. This is a poor assumption in the surface of the disk and the inner parts of the envelope (Kamp & Dullemond 2004; Jonkheid et al. 2004), the consequences of which are addressed in Sect. 8. As argued in Paper I, the disk-envelope accretion shock is not important for the chemical evolution of the material in the disk at the end of the collapse phase. It is therefore not incorporated in our model.

Given the dynamical nature of the collapse, it is easiest to solve the chemistry in a Lagrangian frame. In Paper I, the envelope was populated with several thousand parcels at t = 0 and these were followed in towards the disk and star. The method is reversed here: a regular grid of parcels is defined at the end of the collapse and the parcels are followed backwards in time to their position at t = 0. The reason for doing so is that the 2D abundance profiles at the end of the collapse, when the disk is fully formed, are of greater interest than those at the onset of collapse. In either case, step 4 of the model produces a set of infall trajectories with densities, temperatures and UV intensities as a function of time and position. These data are required for the next step: solving the time-dependent chemistry for each individual parcel. Although the parcels are followed backwards in time to get their trajectories, the chemistry is computed in the normal forward direction. The last step of our model consists of transforming the abundances from the individual parcels back into 2D axisymmetric profiles at any time step of interest.

In Paper I, the model was run for a grid of initial conditions. In the current paper, the analysis is limited to our standard set of parameters: M₀ = 1.0 M_⊙, c_s = 0.26 km s^-1 and Ω₀ = 10^-14 s^-1. The viscosity parameter α is fixed at 0.01 (Hartmann et al. 1998). Section 7 contains a brief discussion on how the results may change for other parameter values.

2.2. Differences with Paper I

The current version of the model contains several improvements over the version used in Paper I, described in detail by Visser & Dullemond (2010). Most importantly, the model now correctly treats the problem of sub-Keplerian accretion onto a 2D disk. Material falling onto the disk along an elliptic orbit has sub-Keplerian angular momentum, so it exerts a torque on the disk that results in an inward push. Several solutions are available (e.g., Cassen & Moosman 1981; Hueso & Guillot 2005), but these are not suitable for our 2D model. The ad-hoc solution from Paper I provided the appropriate qualitative physical correction – increasing the inward radial velocity of the disk material – but it did not properly conserve angular momentum. The updated model uses a new, fully consistent solution, derived directly from the equations for the conservation of mass and angular momentum (Visser & Dullemond 2010). It results in disks that are typically a factor of a few smaller than those obtained previously, but the disk masses are mostly unchanged. The new disks are a few degrees colder in the inner part and warmer in the outer part, which may further affect the chemistry.

Other changes to the model include the definition of the disk-envelope boundary and the shape of the outflow cavity. In Paper I, the disk-envelope boundary was defined as the surface where the density of the infalling envelope material equals that of the disk. The model now uses the surface where the ram pressure of the infalling material equals the thermal pressure of the disk (Visser & Dullemond 2010), providing a more physically correct description of where material becomes part of the disk. The collapse solution of Terebey et al. (1984) does not produce an outflow cavity, so one is added manually. The original model used a cavity with straight walls, but observations and theoretical predictions show that curved walls are more appropriate (Velusamy & Langer 1998; Cantó et al. 2008). The adopted description of the outflow wall is now $\mathit{z}\mathrm{=}\mathrm{(}\mathrm{0.191}\mathrm{AU}\mathrm{)}{\left(\frac{\mathit{t}}{{\mathit{t}}_{\mathrm{acc}}}\right)}^{-3}{\left(\frac{\mathit{R}}{{\mathrm{AU}}_{}}\right)}^{\mathrm{1.5}}\mathit{,}$(1)with R and z in cylindrical coordinates and t_acc = M₀ / Ṁ the time required for the entire envelope to accrete onto the star and disk. The t^-3 dependence is chosen so that the outflow starts narrow and becomes wider as the collapse proceeds. The full opening angle at t_acc is 33.6° at z = 1000 AU and 15.9° at z = 10   000 AU. There are some streamlines in the infalling envelope that hit the cavity wall, but the mass flowing along these lines is insignificant (~1%) compared to the mass accreting onto the star and the disk.

2.3. Radiation field

Photodissociation and photoionisation by ultraviolet (UV) radiation are important processes in the hot inner core and in the surface layers of the disk. The temperature and luminosity of the protostar change in time (Paper I), so neither the strength nor the spectral shape of the radiation field is constant. In addition, the spectral shape changes as the radiation passes through the disk and remnant envelope.

The most accurate way of obtaining the time- and location-dependent photorates is to multiply the cross section for each reaction by the UV field at each grid point. The latter can be computed from 2D radiative transfer at high spectral resolution. As this is too computationally demanding, several approximations have to be made. The first one is to assume that the wavelength-dependent attenuation of the radiation field by the dust in the disk and envelope can be represented by a single factor γ for each reaction (van Dishoeck et al. 2006). The rate coefficient for a given photoreaction at spatial coordinates r and θ and at time t can then be expressed as ${\mathit{k}}_{\mathrm{ph}}^{}\mathrm{\left(}\mathit{r,\theta ,t}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{ph}}^{\mathrm{\ast }}{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{\ast }}}\right)}^{-2}\exp \mathrm{(}\mathrm{-}\mathit{\gamma }{\mathit{A}}_{\mathrm{V}}\mathrm{)}\mathit{,}$(2)with A_V the visual extinction (see below). The unshielded rate coefficient is calculated at the stellar surface (${\mathit{k}}_{\mathrm{ph}}^{\mathrm{\ast }}$) by multiplying the cross section of the reaction by the blackbody flux at the effective temperature of the protostar, which varies between 4500 and 5500 K for the duration of the collapse (Paper I). Any excess UV from the disk-star boundary is not included. The term (r / R _∗ )^-2 accounts for the geometrical dilution of the radiation from the star across a distance r. The factor γ is discussed in Sect. 3.1.

In order to apply Eq. (2), the extinction towards each point (r,θ) is calculated with RADMC at a low spectral resolution of one frequency point per eV. The local UV flux from RADMC at (r,θ), F_RADMC (integrated from 6 to 13.6 eV), can be considered the product of the unattenuated flux and an extinction term, ${\mathit{F}}_{\mathrm{RADMC}}\mathrm{\left(}\mathit{r,\theta ,t}\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{0}}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{UV}\mathit{,}\mathrm{eff}}\mathrm{)}\mathrm{=}{\mathit{F}}_{\mathrm{\ast }}{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{\ast }}}\right)}^{-2}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{UV}\mathit{,}\mathrm{eff}}\mathrm{)}\mathit{.}$(3)F _∗  is the time-dependent UV flux at the stellar surface. The effective UV optical depth, τ_UV,eff, is an average over the many possible paths a photon can follow from the source to the point (r,θ). The large number of photons propagated through the grid by RADMC (typically 10⁵) ensures a statistical sampling of all possible trajectories. The UV optical depth is converted to the visual extinction A_V through the standard relationship A_V = τ_UV,eff / 3.02 (Bohlin et al. 1978).

3. Chemical network

The basis of our chemical network is the UMIST06 database (Woodall et al. 2007) as modified by Bruderer et al. (2009), except that X-ray chemistry is not included. The modified network includes charge-exchange reactions between positive ions and grains (Maloney et al. 1996), which can be important in the disk’s inner few AU. The cosmic-ray ionisation rate of H₂ is set to 5    × 10^-17 s^-1 (Dalgarno 2006). Attenuation of cosmic rays at column densities of more than  ~100 g cm^-2 (Umebayashi & Nakano 1981) is not taken into account, even though this can become relevant for the midplane chemistry in the inner few AU of massive disks. The network contains 162 neutral species, 251 cations and six anions, built up out of 18 elements. The initial composition is fully atomic, except that hydrogen starts as H₂. Elemental abundances are adopted from Aikawa et al. (2008). Additional values are taken from Bruderer et al. (2009), but reduced by a factor of 100 to have a consistent set of low-metal abundances. Table 1 lists the elemental abundances relative to the total hydrogen nucleus density: n_H = n(H) + 2n(H₂).

In order to set the chemical composition at the onset of collapse (t = 0), the initially atomic gas is evolved for a period of 1 Myr at n_H = 8    × 10⁴ cm^-3 and T_g = T_d = 10 K. The visual extinction during this pre-collapse phase is set to 100 mag to disable all photoprocesses, except for a minor contribution from cosmic-ray-induced photons. The resulting solid and gas-phase abundances are consistent with those observed in pre-stellar cores (e.g., Di Francesco et al. 2007). They form the initial conditions for the collapse phase for all infalling parcels. In the remainder of this paper, t = 0 always refers to the onset of collapse, following the 1 Myr pre-collapse phase here described.

3.1. Photodissociation and photoionisation

Photodissociation and photoionisation by UV radiation are important processes in the inner disk and inner envelope. Their rates are given by Eq. (2), using the effective extinction from Eq. (3). The unshielded rates at the stellar surface (${\mathit{k}}_{\mathrm{ph}}^{\mathrm{\ast }}$) are calculated for a blackbody spectrum at the time-dependent stellar temperatures from Paper I. The required cross sections are taken from our freely available database¹. Values for the extinction factor γ from Eq. (2) were tabulated by van Dishoeck et al. (2006) for three different types of radiation fields, of which the 4000 K blackbody is most appropriate for our model (Sect. 2.3). The final output of this procedure is a 3D array (two spatial coordinates and time) with rate coefficients (k_ph) for each photoreaction. When computing the infall trajectories for individual parcels (step 4 in Fig. 1), the rate coefficients at all points along each trajectory are obtained from a linear interpolation.

The photodissociation of H₂ and CO requires some special treatment. Both processes occur exclusively through discrete absorption lines, so self-shielding plays an important role. The amount of shielding for H₂ is a function of the H₂ column density; for CO, it is a function of both the CO and the H₂ column density, because some CO lines are shielded by H₂ lines. The effective UV extinction from Eq. (3) can be converted into a total hydrogen column density through A_V = τ_UV,eff / 3.02 and N_H = 1.59    × 10²¹A_V cm^-2 (Bohlin et al. 1978; Diplas & Savage 1994). Most hydrogen along each photon path is in molecular form, so the effective N(H₂) towards each spatial grid point is  ~0.5N_H. Equation (37) from Draine & Bertoldi (1996) then gives the amount of self-shielding for H₂. The unshielded dissociation rate is computed according to the one-line approximation from van Dishoeck (1987), scaled so that the rate is 4.5 × 10^-11 s^-1 in the standard Draine (1978) field. The photodissociation rate of CO is computed from our new shielding functions and cross sections (Visser et al. 2009a). Effective CO column densities for the self-shielding calculation are derived from the H₂ column densities by assuming an average N(CO) / N(H₂) ratio of 10^-5, accounting for the partial dissociation of CO along each photon path. The results are not sensitive to the exact choice of average CO abundance (van Zadelhoff et al. 2003).

3.2. Gas-grain interactions

All neutral species other than H, H₂, He, Ne and Ar are allowed to freeze out onto the dust according to Charnley et al. (2001). In cold, dense environments – such as our model cores before the onset of collapse – observations show H₂O, CO and CO₂ to be the most abundant ices (Boogert et al. 2011). As the temperature rises during the collapse, the ices evaporate according to their respective binding energies. However, the presence of non-volatile species like H₂O prevents the more volatile species like CO and CO₂ from evaporating entirely: some CO and CO₂ gets trapped in the H₂O ice (Hasegawa & Herbst 1993; Fayolle et al. 2011). The results from Paper I suggested that this is required to explain the presence of CO in solar-system comets.

Incorporating ice trapping in a chemical network is a non-trivial task (Viti et al. 2004; Fayolle et al. 2011), so the process is ignored for now. Instead, desorption of all species is treated according to the zeroth-order rate equation ${\mathit{R}}_{\mathrm{thdes}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{gr}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{gr}}^{}\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathit{\nu }\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathit{N}}_{\mathrm{ss}}\exp \left[\mathrm{-}\frac{{\mathit{E}}_{\mathrm{b}}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{\mathit{k}{\mathit{T}}_{\mathrm{d}}}\right]\mathit{,}$(4)where T_d is the dust temperature, a_gr = 0.1 μm the typical grain radius, and n_gr = 1 × 10^-12n_H the grain number density. The canonical pre-exponential factor, ν, for first-order desorption is 2    × 10¹² s^-1 (Sandford & Allamandola 1993), which is used for all ices except the four listed in Table 2. In order to get the zeroth-order pre-exponential factor, ν is multiplied by N_ss, the number of binding sites per unit grain surface. Assuming N_ss = 8 × 10¹⁴ cm^-2, the number of binding sites per grain is N_b = 1 × 10⁶ for our 0.1 μm grains. The binding energies of species other than the four in Table 2 are set to the values tabulated by Sandford & Allamandola (1993) and Aikawa et al. (1997). Species for which the binding energy is unknown are assigned the binding energy and the pre-exponential factor of H₂O. The dimensionless factor f in Eq. (4) ensures that each species desorbs according to its solid-phase abundance, and changes the overall desorption behaviour from zeroth to first order when there is less than one monolayer of ice: $\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{\max \mathrm{\left(}{\mathit{n}}_{\mathrm{ice}}\mathit{,}{\mathit{N}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{gr}}\mathrm{\right)}}\mathit{,}$(5)with n_s(X) the number density of ice species X, N_b = 10⁶ the typical number of binding sites per grain, and n_ice the total number density (per unit volume of cloud or disk) of all ice species combined. Section 8 briefly discusses how our results would change if trapping were included.

In addition to thermal desorption, our model includes desorption induced by UV photons. Laboratory experiments on the photodesorption of H₂O, CO and CO₂ all produce a yield of Y ≈ 10^-3 molecules per grain per incident UV photon (Öberg et al. 2007, 2009a), while the yield for N₂ is an order of magnitude lower (Öberg et al. 2009b). For all other ice species in our network, whose photodesorption yields have not yet been determined experimentally or theoretically, the yield is set to 10^-3. The photodesorption rate is ${\mathit{R}}_{\mathrm{phdes}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}\mathit{\pi }{\mathit{a}}_{\mathrm{gr}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{gr}}^{}\mathit{f}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathit{Y}\mathrm{\left(}\mathit{X}\mathrm{\right)}{\mathit{F}}_{\mathrm{0}}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{UV}\mathit{,}\mathrm{eff}}\mathrm{)}\mathit{,}$(6)with f the same factor as for thermal desorption. The unattenuated UV flux (F₀) and the effective UV extinction (τ_UV,eff) follow from Eq. (3). Photodesorption occurs even in strongly shielded regions because of cosmic-ray-induced photons, which is accounted for by setting a lower limit of 10⁴ cm^-2 s^-1 to the product F₀exp( − τ_UV,eff) (Shen et al. 2004).

The chemical reactions in our model are not limited to the gas phase. As usual, the network includes the grain-surface formation of H₂ (Black & van Dishoeck 1987). Inspired by Bergin & Langer (1997) and Hollenbach et al. (2009), it also includes the grain-surface hydrogenation of C to CH₄, N to NH₃, O to H₂O, and S to H₂S. The hydrogenation is done one H atom at a time and is always in competition with thermal and photon-induced desorption. The formation of CH₄, NH₃, H₂O and H₂S does not have to start with the respective atom freezing out. For instance, CH freezing out from the gas is also subject to hydrogenation on the grain surface. The rate of each hydrogenation step is taken to be the adsorption rate of H from the gas multiplied by the probability that the H atom encounters the atom or molecule X to hydrogenate: ${\mathit{R}}_{\mathrm{hydro}}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}\mathit{\pi }{\mathit{a}}_{\mathrm{gr}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{gr}}^{}\mathit{n}\mathrm{\left(}\mathrm{H}\mathrm{\right)}{\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\sqrt{\frac{\mathrm{8}\mathit{k}{\mathit{T}}_{\mathrm{g}}}{\mathit{\pi }{\mathit{m}}_{\mathrm{p}}}}$(7)with T_g the gas temperature. The factor f′ serves a similar purpose as the factor f in Eqs. (4) and (6). Since the hydrogenation is assumed to be near-instantaneous as soon as the H atom meets X before X desorbs, X is assumed to reside always near the top layer of the ice. Hence, the abundance of solid X should not be compared to the total amount of ice (as in f), but to the combined abundance of other ice species that can be hydrogenated: ${\mathit{f}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{X}\mathrm{\right)}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{s}}\mathrm{\left(}\mathit{X}\mathrm{\right)}}{\max \mathrm{\left(}{\mathit{n}}_{\mathrm{hydro}}\mathit{,}{\mathit{N}}_{\mathrm{b}}{\mathit{n}}_{\mathrm{gr}}\mathrm{\right)}}\mathit{,}$(8)with n_hydro the sum of the solid abundances of the eleven species subject to hydrogenation: C, CH, CH₂, CH₃, N, NH, NH₂, O, OH, S and SH. The main effect of this hydrogenation scheme is to build up an ice mixture of simple saturated molecules during the pre-collapse phase, as is found observationally (e.g., Tielens et al. 1991; Gibb et al. 2004).

In reality, grain-surface chemistry is not limited to the simple hydrogenation steps included here. For example, CO can react with H to form H₂CO and CH₃OH or with OH to form CO₂ (Fuchs et al. 2009; Ioppolo et al. 2011), and O₂ can be hydrogenated to H₂O (Ioppolo et al. 2008). Grains also play an important role in the formation of more complex species (Garrod & Herbst 2006; Garrod et al. 2008). However, none of these reactions can be implemented as easily as the hydrogenation of C, N, O and S. In addition, the main focus of this paper is on simple molecules whose abundances can be well explained with conventional gas-phase chemistry. The effects of more complex grain-surface chemistry will be explored in detail in a follow-up paper.

4. Results from the pre-collapse phase

This section, together with the next two, contains the results from the gas-phase chemistry in our collapse model. The chemistry during the pre-collapse phase is briefly discussed first. The chemistry during the collapse is analysed in detail for one particular parcel in Sect. 5.1 and then generalised to others in Sect. 5.2. Finally, the collapse chemistry is compared to a static disk model in Sect. 6. The results in the current section are all consistent with available observational constraints on pre-stellar cores (e.g., Di Francesco et al. 2007).

During the 1.0 Myr pre-collapse phase, most of the initial atomic oxygen freezes out and is hydrogenated to H₂O ice. Meanwhile, H₂ is ionised by cosmic rays. The resulting H$\begin{array}{}\mathrm{+}\\ \mathrm{2}\end{array}$ reacts with H₂ to give H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$, which reacts with the remaining atomic O to ultimately form O₂: $\begin{array}{ccc}& & \mathrm{O}\mathrm{+}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\to }{\mathrm{OH}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}\\ & & {\mathrm{OH}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }\mathrm{OH}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{H}\mathit{,}\\ & & \mathrm{OH}\mathrm{+}\mathrm{O}\mathrm{\to }{\mathrm{O}}_{\mathrm{2}}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$The O₂ thus produced freezes out for the most part. At the onset of collapse, the four major oxygen reservoirs are H₂O ice (44%), CO ice (34%), O₂ ice (16%) and NO ice (3%).

The oxygen chemistry is tied closely to the carbon chemistry through CO. It is initially formed in the gas phase from CH₂, which in turn is formed from atomic C: $\begin{array}{ccc}& & \mathrm{C}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }{\mathrm{CH}}_{\mathrm{2}}\mathrm{+}\mathit{h\nu ,}\\ & & {\mathrm{CH}}_{\mathrm{2}}\mathrm{+}\mathrm{O}\mathrm{\to }\mathrm{CO}\mathrm{+}\mathrm{2}\mathrm{H}\mathit{.}\end{array}$Another early pathway from C to CO is powered by H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and goes through an HCO⁺ intermediate: $\begin{array}{ccc}& & \mathrm{C}\mathrm{+}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\to }{\mathrm{CH}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathit{,}\\ & & {\mathrm{CH}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{CH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{CH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{O}\mathrm{\to }{\mathrm{HCO}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathit{,}\\ & & {\mathrm{HCO}}^{\mathrm{+}}\mathrm{+}\mathrm{C}\mathrm{\to }\mathrm{CO}\mathrm{+}{\mathrm{CH}}^{\mathrm{+}}\mathit{.}\end{array}$The formation of CO through these two pathways accounts for most of the pre-collapse processing of carbon: at t = 0, 82% of all carbon has been converted into CO, of which 97% has frozen out onto the grains. Most of the remaining carbon is present as CH₄ ice (14% of all C), formed from the rapid grain-surface hydrogenation of atomic C.

The initial nitrogen chemistry consists mostly of converting atomic N into NH₃, N₂ and NO. The first of these is formed on the grains after freeze-out of N, in the same way that H₂O and CH₄ are formed from adsorbed O and C. The model results show two pathways leading to N₂. The first one starts with the cosmic-ray dissociation of H₂ into H⁺ and H⁻, which Cravens & Dalgarno (1978) estimated to have a 0.015% chance of happening for each cosmic-ray ionisation event: $\begin{array}{ccc}& & {\mathrm{H}}_{\mathrm{2}}\mathrm{+}\mathit{\zeta }\mathrm{\to }{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}^{\mathrm{-}}\mathit{,}\\ & & \mathrm{N}\mathrm{+}{\mathrm{H}}^{\mathrm{-}}\mathrm{\to }\mathrm{NH}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathit{,}\\ & & \mathrm{NH}\mathrm{+}\mathrm{N}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$The other pathway couples the nitrogen chemistry to the carbon chemistry. It starts with Reactions (16)–(18) to form CH$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$, followed by $\begin{array}{ccc}& & {\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }\mathrm{CH}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathit{/}\mathrm{2}\mathrm{H}\mathit{,}\\ & & \mathrm{CH}\mathrm{+}\mathrm{N}\mathrm{\to }\mathrm{CN}\mathrm{+}\mathrm{H}\mathit{,}\\ & & \mathrm{CN}\mathrm{+}\mathrm{N}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}\mathrm{+}\mathrm{C}\mathit{.}\end{array}$The nitrogen chemistry is also tied to the oxygen chemistry, forming NO out of N and OH: $\mathrm{N}\mathrm{+}\mathrm{OH}\mathrm{\to }\mathrm{NO}\mathrm{+}\mathrm{H}\mathit{,}$(27)with OH formed by Reaction (12). Nearly all of the N₂ and NO formed during the pre-collapse phase freezes out. At t = 0, solid N₂, solid NH₃ and solid NO account for 41, 32 and 22% of all nitrogen.

5. Results from the collapse phase

The collapse-phase chemistry is run for the standard set of model parameters from Paper I: M₀ = 1.0 M_⊙, c_s = 0.26 km s^-1 and Ω₀ = 10^-14 s^-1. The core collapses in t_acc = 2.52    × 10⁵ yr to form a 0.76 M_⊙ star surrounded by a compact disk (47 AU outer radius) of 0.13 M_⊙. At t = t_acc, the radius at which H₂O freezes out (the snowline) lies at 5.1 AU. The Toomre (1964)Q parameter ranges from 130 at 0.1 AU to a minimum of 1.4 at 28 AU, so the disk remains gravitationally stable by a small margin.

5.1. One single parcel

The chemical evolution is first discussed in detail for one particular infalling parcel of material. It starts near the edge of the cloud core, at r = 6710 AU from the center and θ = 48.8° degrees from the z axis (R = 5050 AU, z = 4420 AU). Its trajectory terminates at t = t_acc in the inner part of the disk, at R = 6.3 AU and z = 2.4 AU, about 0.2 AU below the surface. The physical conditions encountered along the trajectory are plotted in Fig. 2. The radiation field is characterised by the unattenuated flux and the visual extinction, where the former is expressed as a scaling factor relative to the average flux in the interstellar medium (ISM): χ = F₀ / F_ISM, with F_ISM = 8    × 10⁷ photons cm^-2 s^-1 (Draine 1978). This figure also shows the abundances of the main oxygen-, carbon- and nitrogen-bearing species. The four panels in the right column are regular plots as function of R (the coordinate along the midplane), tracing the parcel’s trajectory from t = 0 (R = 5050 AU) to t = t_acc (R = 6.3 AU). The infall velocity of the parcel increases as it gets closer to the star, so the physical conditions and chemical abundances change more rapidly at later times. Hence, the panels in the left column are plotted as a function of t_acc − t: the time remaining before the end of the collapse phase. In each individual panel, the parcel essentially moves from right to left.

Figure 3 presents a schematic overview of the parcel’s chemical evolution. It shows the infall trajectory of the parcel and the abundances of several species at four points along the trajectory. The physical conditions and the key reactions controlling those abundances are also listed. Most abundance changes for individual species are related to one specific chemical event, such as the evaporation of CO or the photodissociation of H₂O. The remainder of this subsection discusses the abundance profiles from Fig. 2 and explains them in the context of Fig. 3.

Fig. 2 Physical conditions (nH, Td, χ and AV) and abundances of the main oxygen-, carbon- and nitrogen-bearing species for the single parcel from Sect. 5.1, as function of time before the end of the collapse (left) and as function of horizontal position during the collapse (right). The grey bars correspond to the points A, B and C from Fig. 3. In each panel, the parcel essentially moves from right to left.

Fig. 3 Overview of the chemistry along the infall trajectory (solid black curve) of the single parcel from Sect. 5.1. The solid and dashed grey lines denote the surface of the disk and the outflow wall, both at t = tacc = 2.52    × 105 yr. Physical conditions, abundances (black bars: gas; grey bars: ice) and key reactions are indicated at four points (A, B, C and D) along the trajectory. The key processes governing the overall chemistry at each point are listed in the bottom right. The type of each reaction is indicated by colour, as listed in the top left.

5.1.1. Oxygen chemistry

At the onset of collapse (t = 0), the main oxygen reservoir is solid H₂O at an abundance of 8    × 10^-5 relative to n_H (Sect. 4). The abundance remains constant until the parcel gets to point C in Fig. 3, where the temperature is high enough for the H₂O ice to evaporate. The parcel is now located close to the outflow wall, so the stellar UV field is only weakly attenuated (A_V = 0.9 mag). Hence, the evaporating H₂O is immediately photodissociated into H and OH, which in turn is dissociated into O and a second H atom. At R = 17 AU (23 AU inside of point C), the dust temperature is 150 K and all solid H₂O is gone. Moving in further, the parcel enters the surface layers of the disk and is quickly shielded from the stellar radiation (A_V = 10–20 mag). The temperature decreases at the same time to 114 K, allowing some H₂O ice to reform. The final H₂O ice abundance at t_acc (point D) is 4    × 10^-8.

The dissociative recombination of H₃O⁺ (formed by Reactions (9)–(11)) initially maintains H₂O in the gas at an abundance of 7    × 10^-8 at t = 0. Following the sharp increase in the overall gas density at t = 2.1    × 10⁵ yr (Fig. 2), the freeze-out rate increases and the gas-phase H₂O abundance goes down to 3    × 10^-10 at point A in Fig. 3. Moving on towards point B, the evaporation of O₂ from the grains enables a new H₂O formation route: $\begin{array}{ccc}{\mathrm{O}}_{\mathrm{2}}\mathrm{+}{\mathrm{C}}^{\mathrm{+}}& \mathrm{\to }\mathrm{\&}\mathrm{CO}\mathrm{+}{\mathrm{O}}^{\mathrm{+}}\mathit{,}& \\ {\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}& \mathrm{\to }& {\mathrm{OH}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\end{array}$followed by Reactions (10) and (11) to give H₃O⁺, which recombines with an electron to give H₂O. The H₂O abundance thus increases to 3 × 10^-9 at R = 300 AU. Further in, at point C, solid H₂O comes off the grains as described above. However, photodissociation keeps the gas-phase abundance from growing higher than  ~10^-7. Once all H₂O ice is gone at R = 17 AU, the gas-phase abundance can no longer be sustained at 10^-7 and it drops to 3    × 10^-12. Some H₂O is eventually reformed as the parcel gets into the disk and is shielded from the stellar radiation, producing a final abundance of 9    × 10^-9 relative to n_H.

Another main oxygen reservoir at t = 0 is O₂. Its solid and gas-phase abundances in the model are 1    × 10^-5 and 4    × 10^-7, consistent with available observational constraints (Goldsmith et al. 1985, 2011; Fuente et al. 1993). Gaseous O₂ gradually continues to freeze out until it reaches a minimum gas-phase abundance of 3    × 10^-10 just inside of point A. The temperature at that time is 18 K, enough for O₂ to slowly start evaporating thermally. The gas-phase abundance is up by a factor of ten by the time the dust temperature reaches 23 K, about halfway between points A and B. The evaporation is 99% complete as the parcel reaches R = 460 AU, about 140 AU inside of point B. The gas-phase abundance remains stable at 1    × 10^-5 for the next few hundred years. Then, as the parcel gets closer to the outflow wall and into a region of lower extinction, the photodissociation of O₂ sets in and its abundance decreases to 2    × 10^-8 at point C. The evaporation and photodissociation of H₂O at that point enhances the abundances of OH and O, which react with each other to replenish some O₂. As soon as all the H₂O ice is gone, this O₂ production channel quickly disappears and the O₂ abundance drops to 1    × 10^-11. Finally, when the parcel enters the top of the disk, O₂ is no longer photodissociated and its abundance goes back up to 4    × 10^-8 at point D.

The abundance of gas-phase OH starts at 3    × 10^-7. Its main formation pathway is initially the dissociative recombination of H₃O⁺ (Reaction (12)), and its main destructors are O, N and H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$. The increase in total density at 2.1    × 10⁵ yr speeds up the destruction reactions, and the OH abundance drops to 3    × 10^-10 at point B. The evaporation of solid OH then briefly increases the gas-phase abundance to 1    × 10^-8. When all of the OH has evaporated at R = 300 AU, the gas-phase abundance goes down again to 5    × 10^-10 over the next 150 AU. As the parcel continues towards and past point C, the OH abundance is boosted to a maximum of 1    × 10^-6 by the photodissociation of H₂O. The high abundance lasts only briefly, however. As the last of the H₂O evaporates and gets photodissociated, OH can no longer be formed as efficiently, and it is dissociated itself. At the end of the collapse, the OH abundance is  ~ 10^-14.

The fifth main oxygen-bearing species is atomic O itself. Its abundance is 7    × 10^-7 at t = 0 and 1    × 10^-4 at t = t_acc, accounting for respectively 0.4 and 56% of the total amount of non-refractory oxygen. Starting from t = 0, the O abundance remains constant during the first 2.0    × 10⁵ yr of the collapse phase. The increasing overall density then speeds up the reactions with OH (forming O₂) and H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ (forming OH⁺), as well as the adsorption onto the grains, and the O abundance decreases to a low of 2    × 10^-8 just before point A. The abundance goes back up due to the evaporation of CO and, at point B, of O₂ and NO, with O formed from the following reactions: $\begin{array}{ccc}& & \mathrm{CO}\mathrm{+}{\mathrm{He}}^{\mathrm{+}}\mathrm{\to }{\mathrm{C}}^{\mathrm{+}}\mathrm{+}\mathrm{O}\mathrm{+}\mathrm{He}\mathit{,}\\ & & \mathrm{NO}\mathrm{+}\mathrm{N}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}\mathrm{+}\mathrm{O}\mathit{,}\\ & & {\mathrm{O}}_{\mathrm{2}}\mathrm{+}\mathrm{CN}\mathrm{\to }\mathrm{OCN}\mathrm{+}\mathrm{O}\mathit{.}\end{array}$Heading on towards point C, the photodissociation of O₂, NO and H₂O further drives up the amount of atomic O to the aforementioned final abundance of 1    × 10^-4.

5.1.2. Carbon chemistry

With solid and gas-phase abundances of 6    × 10^-5 and 2    × 10^-6 relative to n_H, CO is the main form of non-refractory carbon at the onset of collapse. CO is a very stable molecule and its chemistry is straightforward. The freeze-out process started during the pre-collapse phase continues up to t = 2.4    × 10⁵ yr, a few thousand years prior to reaching point A in Fig. 3, where the dust temperature of 18 K results in CO evaporating again. As the parcel continues its inward journey and is heated up further, all solid CO rapidly disappears and the gas-phase abundance goes up to 6    × 10^-5 at point B. During the remaining part of the infall trajectory, the other main carbon-bearing species are all largely converted into CO. At the end of the collapse (point D), 99.8% of all carbon and 44% of all oxygen are locked up in CO.

The second most abundant carbon-bearing ice at the onset of collapse is CH₄, at 1    × 10^-5 with respect to n_H. The gas-phase abundance of CH₄ begins at 4    × 10^-9, about a factor of 2500 lower. At point A, the evaporation of CO provides the first increase in x(CH₄) through a chain of reactions starting with the formation of C⁺ from CO and He⁺. The successive hydrogenation of C⁺ produces CH$\begin{array}{}\mathrm{+}\\ \mathrm{5}\end{array}$, which reacts with another CO molecule to form CH₄: $\begin{array}{ccc}& & \mathrm{CO}\mathrm{+}{\mathrm{He}}^{\mathrm{+}}\mathrm{\to }{\mathrm{C}}^{\mathrm{+}}\mathrm{+}\mathrm{O}\mathrm{+}\mathrm{He}\mathit{,}\\ & & {\mathrm{C}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{CH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathit{h\nu ,}\\ & & {\mathrm{CH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{CH}}_{\mathrm{5}}^{\mathrm{+}}\mathit{,}\\ & & {\mathrm{CH}}_{\mathrm{5}}^{\mathrm{+}}\mathrm{+}\mathrm{CO}\mathrm{\to }{\mathrm{CH}}_{\mathrm{4}}^{}\mathrm{+}{\mathrm{HCO}}^{\mathrm{+}}\mathit{.}\end{array}$The CH₄ ice evaporates at point B, bringing the gas-phase abundance up to 1    × 10^-5. So far, the abundances of CH₄ and CO are well coupled. The link is broken when the parcel reaches point C, where CH₄ is photodissociated, but CO is not. This difference arises from the fact that CO can only be dissociated by photons shortwards of 1076 Å, while CH₄ can be dissociated out to 1450 Å (Visser et al. 2009a). The 5300 K blackbody spectrum emitted by the protostar at this time is not powerful enough at short wavelengths to cause significant photodissociation of CO. CH₄, on the other hand, is quickly destroyed. Its final abundance at point D is 6 × 10^-10.

Neutral and ionised carbon show the same trends in their abundance profiles, with the former always more abundant by a few per cent to a few orders of magnitude. Both start the collapse phase at  ~10^-8relative to n_H. The increase in total density at 2.1    × 10⁵ yr speeds up the destruction reactions (mainly by OH and O₂ for C and by OH and H₂ for C⁺), so the abundances go down to x(C) = 4    × 10^-10 and x(C ⁺ ) = 5    × 10^-11 just outside point A. This is where CO begins to evaporate, and as a result, the C and C⁺ abundances increase again. As the parcel continues to fall in towards point B, the evaporation of O₂ and NO and the increasing total density cause a second drop in C and C⁺. Once again, though, the drop is of a temporary nature. Moving on towards point C, the parcel gets exposed to the stellar UV field. The photodissociation of CH₄ leads – via intermediate CH, CH₂ or CH₃ – to neutral C, part of which is ionised to also increase the C⁺ abundance. Finally, at point D, the photoprocesses no longer play a role, so the C and C⁺ abundances go back down. Their final values relative to n_H are 7    × 10^-12 and  ~10^-14.

5.1.3. Nitrogen chemistry

The most common nitrogen-bearing species at t = 0 is solid N₂, with an abundance of 5 × 10^-6. The corresponding gas-phase abundance is 4 × 10^-8. The evolution of N₂ parallels that of CO, because they have similar binding energies and are both very stable molecules (Bisschop et al. 2006). N₂ continues to freeze out slowly until it gets near point A in Fig. 3, where the grain temperature of 18 K causes all N₂ ice to evaporate. The gas-phase N₂ remains intact along the rest of the infall trajectory and its final abundance is 1    × 10^-5, accounting for 77% of all nitrogen.

The second largest nitrogen reservoir at the onset of collapse is NH₃, with solid and gas-phase abundances of 8    × 10^-6 and 2    × 10^-8. The gas-phase abundance receives a short boost at point A due to the evaporation of N₂, followed by $\begin{array}{ccc}& & {\mathrm{N}}_{\mathrm{2}}\mathrm{+}{\mathrm{He}}^{\mathrm{+}}\mathrm{\to }{\mathrm{N}}^{\mathrm{+}}\mathrm{+}\mathrm{N}\mathrm{+}\mathrm{He}\mathit{,}\\ & & {\mathrm{N}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }{\mathrm{NH}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{NH}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{NH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{NH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{NH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{NH}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathrm{\to }{\mathrm{NH}}_{\mathrm{4}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{NH}}_{\mathrm{4}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }{\mathrm{NH}}_{\mathrm{3}}^{}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$The binding energy of NH₃ is intermediate to that of O₂ and H₂O, so it evaporates between points B and C. Like H₂O, NH₃ is photodissociated upon evaporation. As the last of the NH₃ ice leaves the grains at R = 50 AU (10 AU outside of point C), the gas-phase reservoir is no longer replenished and x(NH₃) drops to  ~10^-14. Some NH₃ is eventually reformed as the parcel gets into the disk, and the final abundance at point D is 1    × 10^-10 relative to n_H.

With an abundance of 6 × 10^-6, solid NO is the third major initial nitrogen reservoir. Gaseous NO is a factor of twenty less abundant at t = 0: 3 × 10^-7. The NO gas is gradually destroyed prior to reaching point A by continued freeze-out and reactions with H⁺ and H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$. It experiences a brief gain at point A from the evaporation of OH and its subsequent reaction with N to give NO and H. As the parcel continues to point B, the solid NO begins to evaporate and the gas-phase abundance rises to 6    × 10^-6. Photodissociation reactions then set in around R = 100 AU and the NO abundance goes back down to 6    × 10^-9. The evaporation and photodissociation of NH₃ cause a brief spike in the NO abundance through the reactions $\begin{array}{ccc}& & {\mathrm{NH}}_{\mathrm{3}}\mathrm{+}\mathit{h\nu }\mathrm{\to }{\mathrm{NH}}_{\mathrm{2}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{NH}}_{\mathrm{2}}\mathrm{+}\mathrm{O}\mathrm{\to }\mathrm{HNO}\mathrm{+}\mathrm{H}\mathit{,}\\ & & \mathrm{HNO}\mathrm{+}\mathit{h\nu }\mathrm{\to }\mathrm{NO}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$The evaporation of the last of the NH₃ ice at R = 50 AU eliminates this channel and the NO gas abundance decreases to 3 × 10^-9 at point C. NO is now mainly sustained by the reaction between OH and N. As described above, the OH abundance drops sharply at R = 17 AU, and the NO abundance follows suit. The final abundance at point D is  ~10^-14.

The last nitrogen-bearing species from Fig. 2 is atomic N itself. It starts at an abundance of 1 × 10^-7 and slowly freezes out to an abundance of 2    × 10^-8 just before reaching point A. At point A, N₂ evaporates and is partially converted to N₂H⁺. The dissociative recombination of N₂H⁺ mostly reforms N₂, but it also produces some NH and N. The N abundance jumps back to 1    × 10^-7 and remains nearly constant at that value until the parcel reaches point B, where NO evaporates and reacts with N to produce N₂ and O. This reduces x(N) to a low of 5 × 10^-10 between points B and C. Moving in further, the parcel gets exposed to the stellar UV field, and NO and NH₃ are photodissociated to bring the N abundance to a final value of 5 × 10^-6 relative to n_H. As such, it accounts for 22% of all nitrogen at the end of the collapse.

5.2. Other parcels

At the end of the collapse (t = t_acc = 2.52 × 10⁵ yr), the parcel from Sect. 5.1 (hereafter called our reference parcel) is located at R = 6.3 AU and z = 2.4 AU, about 0.2 AU below the surface of the disk. As shown in Fig. 3, its trajectory passes close to the outflow wall, through a region of low extinction. This results in the photodissociation or photoionisation of many species. At the same time, the parcel experiences dust temperatures of up to 150 K (Fig. 2), well above the evaporation temperature of H₂O and all other non-refractory species in our network. Material that ends up in other parts of the disk encounters different physical conditions during the collapse and therefore undergoes a different chemical evolution. This subsection shows how the absence or presence of some key chemical processes, related to certain physical conditions, affects the chemical history of the entire disk. Table 3 lists the abundances of selected species at four points in the disk at t_acc.

5.2.1. Oxygen chemistry

The main oxygen reservoir at the onset of collapse is H₂O ice (Sect. 4). Its abundance remains constant at 1    × 10^-4 in our reference parcel until it gets to point C in Fig. 3, where it evaporates from the dust and is immediately photodissociated. When the parcel enters the disk, some H₂O is reformed to produce final gas-phase and solid abundances of  ~10^-8relative to n_H (Sect. 5.1.1).

Figure 4 shows the disk at t_acc, divided into seven zones according to different chemical evolutionary schemes for H₂O. The fraction of the disk mass in each zone is indicated. The material in zone 1 (61% of the disk mass) is the only material in the disk in which H₂O never evaporates during the collapse, because the temperature never gets high enough. The H₂O ice abundance is constant throughout zone 1 at t_acc at  ~1 × 10^-4. For the material ending up in the other six zones, H₂O evaporates at some point during the collapse phase. Zone 2 (1.0%) contains our reference parcel, so its H₂O history has already been described. The total H₂O abundance (gas and ice combined) at t_acc is  ~10^-8at the top of zone 2 and  ~10^-6at the bottom. The gas-ice ratio goes from  ~1 at the top to  ~10^-6at the bottom.

Fig. 4 Schematic view of the history of H2O gas and ice throughout the disk. The main oxygen reservoir at tacc is indicated for each zone; the histories are described in the text. The percentages indicate the fraction of the disk mass contained in each zone. Note the disproportionality of the R and z axes. The colours have no specific meaning other than to distinguish the different zones.

The H₂O history of zone 3 (4.7%) is the same as that of zone 2, except that it finishes with a gas-ice ratio larger than unity. In both cases H₂O evaporates and is photodissociated prior to entering the disk (point C in Fig. 3), and it partially reforms inside the disk (point D). Parcels ending up in zone 4 (0.4%) also experience the evaporation and photodissociation of H₂O. However, the low extinction against the stellar radiation in zone 4 prevents H₂O from reforming like it does in zones 2 and 3.

The material in zone 5 (0.9%) has a different history from that in zones 2–4 because it enters the disk earlier: between 1.3    × 10⁵ and 2.3    × 10⁵ yr. The material in zones 2–4 all accretes after 2.4    × 10⁵ yr. The infall trajectories terminating in zone 5 do not pass close enough to the outflow wall or the inner disk surface for photoprocesses to play a role. All H₂O in zone 5 is in the gas phase at t_acc (abundance: 1    × 10^-4) because it lies inside the disk’s snow line. The evaporation of H₂O ice does not occur until the material actually crosses the snow line. Prior to that point, the temperature never gets high enough for H₂O to leave the grains.

Fig. 5 Infall trajectory for a parcel ending up in zone 6 from Fig. 4, showing where H2O is predominantly present as ice (black) or gas (grey). The two diamonds mark the time in years after the onset of collapse.

The material ending up in zones 6 and 7 (1.6 and 30%) accretes even earlier than that ending up in zone 5: at t = 4    × 10⁴ yr. The disk at that time is only 2 AU large and several 100 K hot, so the ice mantles are completely removed. Upon reaching the disk, the zone 6 and 7 material is swept up as the disk expands to conserve angular momentum, resulting in trajectories like the one in Fig. 5. Because the zone 5 material enters the disk at a later time, its trajectories do not feature a strong outward component and terminate at smaller radii than the zone 6 and 7 trajectories. Along the outward part of the latter, material cools down and H₂O returns to the solid phase (Fig. 5). At t = 1.7    × 10⁵ yr, the parcel starts moving inwards again and comes close enough to the protostar for H₂O to evaporate a second time. All parcels ending up in zone 6 have similar trajectories and the same qualitative H₂O history. The parcels ending up in zone 7 also have H₂O evaporating during the initial infall and freezing out again during the outward part of the trajectory. However, they do not terminate close enough to the protostar for H₂O to desorb a second time. Therefore, most H₂O in zone 7 is on the grains at t_acc.

Our model is not the first one in which part of the disk contains H₂O that evaporated and readsorbed. Lunine et al. (1991), Owen et al. (1992) and Owen & Bar-Nun (1993) argued that the accretion shock at the surface of the disk is strong enough for H₂O to evaporate. However, based on the model of Neufeld & Hollenbach (1994), Paper I showed that most of the disk material does not pass through a shock that heats the dust to 100 K or more. Moreover, the material that does get shock-heated to that temperature accretes close enough to the star that the stellar radiation already heats it to more than 100 K. Hence, including the accretion shock explicitly in our model would at most result in minor changes to the chemistry.

As discussed in Sect. 5.1.1, H₂O controls part of the oxygen chemistry along the infall trajectory of our reference parcel, and it does so for other parcels, too. Figure 6 presents a schematic view of the chemical evolution of six oxygen-bearing species towards each of the seven zones from Fig. 4. The abundances are indicated qualitatively as high, intermediate or low. The horizontal axes (time, increasing from left to right) are non-linear and only indicate the order in which various events take place.

Fig. 6 Qualitative evolution of some abundances towards the seven zones with different H2O histories from Fig. 4. The horizontal axes show the time (increasing from left to right) and are non-linear. The position of points A, B, C and D from Fig. 3 is indicated for zone 2, which contains our reference parcel.

For material ending up in zone 1, H₂O never evaporates, but O₂ does. OH is inititally relatively abundant but most of it disappears when the overall density increases and the reactions with O, N and H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ become faster. The abundance of atomic O experiences a drop at the same time, but it quickly goes back up due to the evaporation of CO, O₂ and NO, followed by Reactions (30)–(32). O₂ also evaporates on its way to zones 2, 3 and 4, and because it passes through an area of low extinction, it is subsequently photodissociated. This, together with the photodissociation of OH and H₂O, causes an increase in the abundance of atomic O. Zones 2 and 3 are sufficiently shielded against the stellar UV field, so O₂ is reformed at the end of the trajectory. This does not happen in the less extincted zone 4.

Material ending up in zone 5 has the same qualitative history for gas-phase and solid O₂ as has material ending up in zone 1. The evolution of atomic O initially also shows the same pattern, but it experiences a second drop as the total density becomes even higher than it does for zone 1. The higher density also causes an additional drop in the OH abundance. En route towards zones 6 and 7, O₂ evaporates and is photodissociated during the early accretion onto the small disk. It is reformed during the outward part of the trajectory and survives until t_acc. Atomic O is also relatively abundant along most of the trajectory, although always one or two orders of magnitude below O₂. In zone 6, the O abundance decreases at the end due to the reaction with evaporating CH₃, producing H₂CO and H. Zone 7 does not get warm enough for CH₃ to evaporate, so O remains intact.

5.2.2. Carbon chemistry

The two main carbon reservoirs at the onset of collapse are CO ice and CH₄ ice (Sect. 4). Their binding energies are relatively low (855 and 1080 K), so they evaporate throughout the core soon after the collapse begins (Sect. 5.1.2). For our reference parcel, the main difference between the evolution of CO and CH₄ is the photodissociation of the latter (near point C in Fig. 3) while the former remains intact. Following the analysis for H₂O, the disk at t_acc could be divided into several zones according to different chemical evolution scenarios for CO. However, the entire disk has the same qualitative CO history: CO does not undergo any processing apart from evaporating early on in the collapse phase, and it is the main carbon reservoir throughout ( > 99%; Paper I). Hence, the disk is divided according to the evolution of CH₄ instead (Fig. 7).

Fig. 7 As Fig. 4, but for CH4 gas and ice. The main carbon reservoir at tacc is CO gas throughout the disk.

For material that ends up in zone 1 (representing 62% of the disk mass), the only major chemical event for CH₄ is the evaporation during the initial warm-up of the core. It is not photodissociated at any point during the collapse, nor does it freeze out again or react significantly with other species. Material ending up closer to the star, in zones 2 and 3 (5.3 and 0.4%), is sufficiently irradiated by the stellar UV field for CH₄ to be photodissociated. The amount of extinction increases towards the end of the trajectories terminating in these zones, allowing some CH₄ to reform. The final abundance ranges from 10^-10 in zone 3 to 10^-7 in the most shielded parts of zone 2.

Zone 4 (0.9%) contains material that accretes onto the disk several 10⁴ yr earlier than does the material in zones 2 and 3. It always remains well shielded from the stellar radiation, so the only processing of CH₄ is the evaporation during the early parts of the collapse. The CH₄ history of zones 1 and 4 is thus qualitatively the same. Zone 5 (32%) consists of material that accretes around t = 4    × 10⁴ yr and is subsequently transported outwards to conserve angular momentum (Fig. 5). CH₄ in this material evaporates before entering the young disk and is photodissociated as it gets within a few AU of the protostar. The resulting atomic C is mostly converted into CO and remains in that form for the rest of the trajectory. Hence, even though the extinction decreases again when the parcel moves outwards, no CH₄ is reformed.

Fig. 8 As Fig. 6, but for the five zones with different CH4 histories from Fig. 7.

The evolution of the abundances of CH₄ gas and ice, CO gas and ice, HCO⁺, C and C⁺ towards each of the five zones is plotted schematically in Fig. 8. The HCO⁺ abundance follows the CO abundance at a ratio that is roughly inversely proportional to the overall density. Hence, the HCO⁺ evolution is qualitatively the same towards each zone: it reaches a maximum abundance of a few 10^-10 when CO evaporates and gradually disappears as the density increases along the rest of the infall trajectories. The most complex history amongst these seven carbon-bearing species is found in C and C⁺. Towards all five zones, they are initially destroyed by reactions with H₂, O₂ and OH. Some C and C⁺ is reformed when CO evaporates, but the subsequent evaporation of O₂ and OH causes the abundances to decrease again. En route to zones 2, 3 and 5, the photodissociation of CH₄ leads to a second increase in C and C⁺, followed by a third and final decrease when the parcel moves into a more shielded area.

5.2.3. Nitrogen chemistry

Most nitrogen at the onset of collapse is in the form of solid N₂ (41%), solid NH₃ (32%) and solid NO (22%). The evolution of N₂ during the collapse is the same as that of CO, except for a minor difference in the binding energy. Both species evaporate shortly after the collapse begins and remain in the gas phase throughout the rest of the simulation. Neither one is photodissociated, because the protostar is not hot enough to provide sufficient UV flux shortwards of 1100 Å.

The evolution of NH₃ shows a lot more variation, as illustrated in Fig. 9. The disk at t_acc is divided into eight zones with different NH₃ histories. No processing occurs towards zone 1 (8.3% of the disk mass): the temperature never exceeds the 73 K required for NH₃ to evaporate, so it simply remains on the grains the whole time. Material ending up in zone 2 (25%) does get heated above 73 K, and NH₃ evaporates. However, it freezes out again at the end of the trajectory because zone 2 itself is not warm enough to sustain gaseous NH₃. The final solid NH₃ abundance in zones 1 and 2 is about 8    × 10^-6 relative to n_H (Table 3).

Fig. 9 As Fig. 4, but for NH3 gas and ice. The main nitrogen reservoir at tacc is N2 gas throughout the disk.

NH₃ ending up in zone 3 (20%) evaporates from the grains just before entering the disk, but it is immediately photodissociated or destroyed by HCO⁺. Towards the end of the trajectory, some NH₃ is reformed from the dissociative recombination of NH$\begin{array}{}\mathrm{+}\\ \mathrm{4}\end{array}$. It rapidly freezes out to produce a final solid NH₃ abundance of 10^-10–10^-8. Zone 4 (7.2%) has the same history, except that there is an additional adsorption-desorption cycle before the destruction by photons and HCO⁺.

Our standard parcel ends up in zone 5 (7.7%); its NH₃ evolution is mainly characterised by photodissociation above the disk and reformation inside it (Sect. 5.1.3). NH₃ evaporates when it gets to within about 200 AU of the star, halfway between points B and C in Fig. 3. It is immediately photodissociated, but some NH₃ is reformed once the material is shielded from the stellar radiation. The reformed NH₃ remains in the gas phase. No reformation takes place in the less extincted zone 6 (0.4%), which otherwise has the same NH₃ history as zone 5.

Material ending up in zone 7 (0.9%) does not pass close enough to the outflow wall for NH₃ to be photodissociated upon evaporation. Instead, NH₃ is mainly destroyed by HCO⁺ and attains a final abundance of  ~ 10^-9. Lastly, zone 8 (30%) contains again the material that accretes onto the disk at an early time and then moves outwards to conserve angular momentum. Its NH₃ evaporates already before reaching the disk and is subequently dissociated by the stellar UV field. As the material moves away from the star and is shielded from its radiation, some NH₃ is reformed out of NH$\begin{array}{}\mathrm{+}\\ \mathrm{4}\end{array}$ to a final abundance of  ~ 10^-10.

Fig. 10 As Fig. 6, but for the eight zones with different NH3 histories from Fig. 9.

The abundances of N₂H⁺ and atomic N are largely controlled by the evolution of N₂ and NH₃, as shown schematically in Fig. 10. In all parcels, regardless of where they end up, N₂H⁺ is mainly formed out of N₂ and H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$, so its abundance goes up when N₂ evaporates shortly after the onset of collapse. It gradually disappears again as the collapse proceeds due to the increasing density. The atomic N abundance at the onset of collapse is 1    × 10^-7. It increases when N₂ evaporates and decreases again a short while later when NO evaporates (Sect. 5.1.3). For material that ends up in zones 1 and 2, the N abundance is mostly constant for the rest of the collapse at a value of 10^-13 (inner part of zone 2) to 10^-10 (outer part of zone 1). Material that ends up in the other six zones is exposed to enough UV radiation for NO and NH₃ to be photodissociated, so there is a second increase in atomic N. Zones 3, 4, 7 and 8 are sufficiently shielded at t_acc to reform some or all NO and NH₃, and the N abundance finishes low. Less reformation is possible in zones 5 and 6, so they have a relatively large amount of atomic N at the end of the collapse phase.

5.2.4. Mixing

Given the dynamic nature of circumstellar disks, the zonal distribution presented in the preceding subsections may offer too simple a picture of the chemical composition. For example, there are as yet no first principles calculations of the processes responsible for the viscous transport in disks. The radial velocity equation used in our model (Visser & Dullemond 2010) is suitable as a zeroth-order description, but it omits the effects of mixing (Lewis & Prinn 1980; Semenov et al. 2006; Aikawa 2007) and it cannot explain important observational features like episodic accretion (Kenyon & Hartmann 1995). Indeed, with the Toomre (1964)Q parameter reaching a minimum of 1.4 at t_acc, the disk is only barely gravitationally stable. Small changes in the model or in the initial conditions would suffice to make the disk unstable at various points during the collapse phase, resulting in a redistribution of material beyond the equations currently used. Hence, in real disks, both the shapes and the locations of the zones are likely to differ from what is shown in Figs. 4, 7 and 9. A larger degree of mixing would also make the borders between the zones more diffuse than they are in our simple schematic representation. Nevertheless, the general picture from this section offers a plausible description of the chemical history towards different parts of the disk. Spectroscopic observations at 5–10 AU resolution, as possible with ALMA, are required to determine to what extent this picture holds in reality. Based on the abundances in Table 3, potential tracers are H₂O, NO and NH₃.

6. Chemical history versus local chemistry

Section 5 contains many examples of abundances increasing or decreasing on short timescales of less than a hundred years (see, e.g., Fig. 2). It appears that the abundances respond rapidly to the changing physical conditions as material falls in supersonically through the inner envelope and accretes onto the disk. However, this does not necessarily mean that the abundances are always in equilibrium. The current section explores the question of whether the disk is in chemical equilibrium at the end of the collapse, or whether its chemical composition is a non-equilibrium solution to the conditions encountered during the collapse. To that end, the chemistry is evolved for an additional 1 Myr beyond t_acc. The density, temperature, UV flux and extinction are kept constant at the values they have at t_acc, and all parcels of material are kept at the same position. Clearly, this is a purely hypothetical scenario. In reality, the disk would change in many ways after t_acc: it spreads in size, the surface layers become more strongly irradiated, the temperature changes, the dust coagulates into planetesimals, gas is photoevaporated from the surface layers, and so on. All of these processes have the potential to affect the chemical composition. However, they would also interfere with the attempt to determine whether the disk is in chemical equilibrium at t_acc. This question is most easily addressed by evolving the chemistry for an additional period of time at constant conditions. Hence, all of the physical changes that would occur during the post-collapse phase are ignored.

The 21 species discussed in Sect. 5 can be divided into two categories: those whose abundance profile changes during the post-collapse phase, and those whose abundance profile remains practically the same. Members of the “changed” category are O₂, OH, CH₄, NH₃ and NO (all gaseous), NH₃ ice, and atomic O and N. The thirteen species in the “unchanged” category are H₂O gas and ice, O₂ ice, CO gas and ice, CH₄ ice, C, C⁺, N₂ gas and ice, NO ice, HCO⁺ and N₂H⁺. The individual gas and ice abundances are summed in Fig. 11 and compared at t = t_acc and t_acc + 1 Myr. The total H₂O, CO and N₂ abundances do not change significantly during the post-collapse phase, while the total O₂, CH₄, NH₃ and NO abundances change by more than two orders of magnitude in a large part of the disk.

Fig. 11 Abundances of total H2O, O2, CO, CH4, N2, NH3 and NO (gas and ice combined) throughout the disk at the end of the collapse phase (t = tacc) and after an additional 1 Myr post-collapse phase (t = tacc + 1 Myr).

There are two areas in the disk where the abundances of all 21 species remain nearly constant: near the surface out to  ~10 AU, and at the midplane between  ~5 and  ~25 AU. The chemistry in the first area is dominated by fast photoprocesses. The second area, near the midplane, is the densest part of the disk and therefore has high collision frequencies. In both cases, the chemical timescales are short and the chemistry reaches equilibrium during the final stages of accretion.

Outside these two “equilibrium areas”, one of the key processes during the post-collapse phase is the conversion of gas-phase oxygen-bearing species into H₂O, which subsequently freezes out. Because H₂O ice is already one of the main oxygen reservoirs at t_acc throughout most of the disk (Fig. 4), its abundance only increases by about 20% during the rest of the simulation. These 20% represent 10% of the total oxygen budget and correspond to 2    × 10^-5 oxygen nuclei per hydrogen nucleus.

One of the species converted into H₂O ice is NO, which reacts with NH₂ through ${\mathrm{NH}}_{\mathrm{2}}\mathrm{+}\mathrm{NO}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathit{,}$(47)followed by adsorption of H₂O onto dust. This reaction is responsible for the post-collapse destruction of gas-phase NO outside the two “equilibrium areas”, as well as for a  ~20% increase in gas-phase N₂. NO is the main destructor of atomic N via $\mathrm{NO}\mathrm{+}\mathrm{N}\mathrm{\to }{\mathrm{N}}_{\mathrm{2}}\mathrm{+}\mathrm{O}\mathit{,}$(48)so the conversion of NO into H₂O leads to a higher N abundance. The NH₂ required for Reaction (47) is formed from NH₃: $\begin{array}{ccc}& & {\mathrm{NH}}_{\mathrm{3}}^{}\mathrm{+}{\mathrm{HCO}}^{\mathrm{+}}\mathrm{\to }{\mathrm{NH}}_{\mathrm{4}}^{\mathrm{+}}\mathrm{+}\mathrm{CO}\mathit{,}\\ & & {\mathrm{NH}}_{\mathrm{4}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }{\mathrm{NH}}_{\mathrm{2}}^{}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{}\mathit{/}\mathrm{2}\mathrm{H}\mathit{.}\end{array}$The gas-phase reservoir of NH₃ is continously fed by evaporation of NH₃ ice, so Reactions (49), (50) and (47) effectively transform all NH₃ into N₂.

Two other important species that are destroyed in the post-collapse phase are gaseous O₂ and CH₄. O₂ is ionised by cosmic rays to produce O$\begin{array}{}\mathrm{+}\\ \mathrm{2}\end{array}$, which reacts with CH₄: $\begin{array}{ccc}& & {\mathrm{CH}}_{\mathrm{4}}^{}\mathrm{+}{\mathrm{O}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\to }{\mathrm{HCOOH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}\\ & & {\mathrm{HCOOH}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }\mathrm{HCOOH}\mathrm{+}\mathrm{H}\mathit{.}\end{array}$The formic acid (HCOOH) freezes out, thus acting as a sink for both carbon and oxygen. At t = t_acc + 1 Myr, the solid HCOOH abundance is 3    × 10^-6, accounting for 3% of all oxygen and 4% of all carbon. The presence of such a sink is a common feature of disk chemistry models (e.g., Aikawa & Herbst 1999). Gas-phase species are processed by He⁺ and H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ until they form a species whose evaporation temperature is higher than the dust temperature. This is HCOOH in our case, but it could also be carbon-chain molecules like C₂H₂ or C₃H₄.

Another important post-collapse sink reaction in our model is the freeze-out of HNO. At t > t_acc, OH is primarily formed by $\mathrm{HNO}\mathrm{+}\mathrm{O}\mathrm{\to }\mathrm{NO}\mathrm{+}\mathrm{OH}$(53)throughout most of the disk. This reaction is also one of the main destruction channels for atomic O. Hence, the gradual freeze-out of HNO after t_acc leads to a decrease in OH and an increase in O. Also contributing to the higher O abundance is the fact that OH itself is an important destructor of O through Reaction (13).

The post-collapse chemistry described here is merely meant as an illustration of what might take place in a real circumstellar disk after the main accretion phase comes to an end. Depending on how the disk continues to evolve physically, other reactions may become more important than the ones listed above. However, with the abundances of several key species changing by orders of magnitude in our post-collapse experiment, the disk was clearly not in chemical equilibrium at the end of the collapse. Even the two “equilibrium areas” (near the surface out to 10 AU and from 5 to 25 AU at the midplane) are not truly in chemical equilibrium when considering less abundant species. For example, at 20 AU on the midplane, post-collapse processing increases the abundance of solid HCOOH by a factor of 13 (from 7 × 10^-8 to 9 × 10^-7 relative to n_H) and that of solid CH₃OH by a factor of 27 (from 3 × 10^-13 to 8 × 10^-12).

The lack of chemical equilibrium at the end of the collapse phase means that one has to be careful in choosing the initial abundances when running T Tauri or Herbig Ae/Be disk chemistry models. The abundances of minor species, such as most organic molecules, can be particularly sensitive to what happened before the disk was formed. Ideally, any disk chemistry model should therefore include a realistic core collapse phase to set the stage for the chemical evolution in the disk itself.

7. Comets: mixed origins?

The chemical composition of cometary ices shows many similarities to that of young protostars, but the abundance of a given species may vary by more than an order of magnitude from one comet to the next (A’Hearn et al. 1995; Mumma & Charnley 2011). Both points are illustrated in the top panel of Fig. 12, where the abundances of several species in the comets 1P/Halley, C/1995 O1 (Hale-Bopp), C/1996 B2 (Hyakutake), C/1999 H1 (Lee), C/1999 S4 (LINEAR) and 153P/Ikeya-Zhang are plotted against the gas-phase abundances in the embedded protostar IRAS 16293–2422 (warm inner envelope), the bipolar outflow L 1157, and the four hot cores W3(H₂O), G34.3+0.15, Orion HC and Orion CR. All abundances are gathered from Bockelée-Morvan et al. (2000, 2004) and Schöier et al. (2002). The cometary data are converted from abundances relative to H₂O to abundances relative to n_H by assuming n(H₂O) / n_H = 5    × 10^-5.

Fig. 12 Top: abundances observed in comets versus gas-phase abundances observed in young protostars. Middle and bottom: abundances observed in comets versus total abundances (gas and ice) from our model at t = tacc. All abundances are relative to nH. For the observational data points, the error bars indicate the spread between sources. For the model data, they indicate the range of values across the entire disk (middle) or the comet-forming region (bottom). Dotted error bars indicate the full range extends to less than 10-12. The diagonal dotted lines represent hypothetical one-to-one relationships between abundances on the x and y axes.

Each point in the top panel of Fig. 12 is characterised by the mean value from the available sources (the diamonds) and the total spread in measurements (the error bars). Uncertainties from individual measurements are not shown. The dotted line marks the one-to-one relationship where the cometary abundances equal the protostellar abundances. The data generally follow this line, suggesting that the material ending up in comets underwent little chemical processing during the evolution from protostar to main-sequence star. However, the plot also shows variations of at least an order of magnitude in the cometary abundances for CO, H₂CO, CH₃OH, HNC and H₂S, as well as smaller variations for other species. A possible explanation for the different chemical compositions is that the comets formed in different parts of the solar nebula. If that is indeed the case, there must have been some degree of chemical processing between the protostellar stage and the time where each comet formed. Another point to note is that the elemental nitrogen abundance in comets is at least a factor of three lower than it is in protostars (Wyckoff et al. 1991). Although the reason for this deficiency remains unclear, it also points at a certain degree of chemical processing.

Comets are thought to have formed at the gravitational midplane of the circumsolar disk, between 5 and 30 AU from the young Sun (Bockelée-Morvan et al. 2004). In our model disk, the 5–30 AU range is almost exactly the area containing material that accreted at an early time (4    × 10⁴ yr after the onset of collapse) and stayed in the disk for the remainder of the collapse phase (Fig. 5). As discussed in Sect. 5.2, this material undergoes a larger degree of chemical processing than any other material in the disk, due to its accreting close to the star. However, such a large degree of processing is hard to reconcile with the observed similarities between cometary and protostellar abundances, and it is therefore unlikely that comets form entirely out of material that accreted close to the star. Instead, as argued below, our model results suggest that comets consist of a mix of processed and unprocessed material.

The middle and bottom panels of Fig. 12 compare the cometary abundances from Bockelée-Morvan et al. (2000, 2004) to the volume-averaged abundances from our model at the end of the collapse phase. For all species, the gas and ice abundances from the model are summed and displayed as one. The horizontal error bars show again the spread in abundances between individual comets. In the middle panel, the vertical bars show the spread across the entire model disk; in the bottom panel, they show the spread in the 5–30 AU region at the midplane. Dotted extensions to the error bars indicate that the range of values extends to less than 10^-12.

The disk-averaged abundances tend to cluster around the hypothetical one-to-one relationship indicated by the diagonal dotted line in the middle panel of Fig. 12. Two notable exceptions are CH₃OH and HCOOCH₃. Both species are predominantly formed through grain-surface chemistry, which is not included in our model. The spread in model abundances is at least six orders of magnitude, but given the wide range of physical conditions that are sampled, such a large spread is to be expected. The 5–30 AU abundances from the bottom panel of Fig. 12 are lower for many species than the disk-wide abundances. For example, the average NH₃ abundance decreases by more than three orders of magnitude. NH₃ is photodissociated along the trajectories that terminate in the comet-forming zone, whereas other trajectories are more shielded (Sect. 5.2.3).

The abundances in the bottom panel of Fig. 12 are sampled over a smaller region, so they show less variation than those in the middle panel. However, they still take on a wider range of values than do the cometary abundances. More importantly, the 5–30 AU data show no correlation with the comet data. There is no indication from Sect. 6 that post-collapse processing brings the abundances back to near what they were in the protostellar stage. Hence, the midplane material between 5 and 30 AU in our model does not appear to be analogous to the material from which the solar-system comets formed.

How plausible is the large degree of chemical processing for material that ends up in the comet-forming zone? The amount of processing depends on the range of physical conditions encountered along a given infall trajectory, which in turn depends mostly on how close to the protostar the trajectory gets. The physics in our model are kept as simple as possible, so a trajectory similar to the one drawn in Fig. 5 may not be fully realistic. On the other hand, the back-and-forth motion results from the well-known concept of conservation of angular momentum. As the inner parts of the disk accrete onto the star, the outer parts must move out to maintain the overall angular momentum. They may be pushed inwards again at a later time if a sufficient amount of mass is accreted from the envelope at larger radii. This happens in our simple model, but also in the more physically realistic hydrodynamical simulations of Brinch et al. (2008).

Accepting some degree of back-and-forth motion as a common part in disk evolution, the question remains how close to the star the material gets before moving out to colder parts of the disk. Our revised solution to the problem of sub-Keplerian accretion (Visser & Dullemond 2010) results in radial velocity profiles that are different from the ones in Paper I. With the old profiles, all material that accretes inside of the snow line always remains inside of the snow line. The comet-forming zone therefore did not contain any material in which H₂O had evaporated and readsorbed onto the grains. Such material does exist in the comet-forming zone with the new velocity profiles (Fig. 5). However, the result depends strongly on the initial physical conditions. If the new model is run for the parameter grid from Paper I, there are three cases (out of eight) where material accretes inside of the snow line and then moves out beyond it; the standard disk discussed throughout this work is one of them. The critical question is whether the point at which the radial velocity in the disk changes sign (from inwards to outwards) lies inside or outside of the snowline. In general, the radial velocity changes sign at smaller radii for less massive disks. Such systems are therefore more likely to have material accreting inside of the snowline and then moving out. Indeed, the three cases in our parameter grid where this happens are the three with the least massive disks. Quantitatively, the dividing factor appears to be a ratio of 0.2 between the disk mass at t = t_acc and the core mass at t = 0. For larger ratios, material that accretes inside of the snow line always remains there.

The presence of crystalline silicate dust in disks provides a strong argument in favour of trajectories similar to the one from Fig. 5. Observations show crystalline fractions at R ≈ 10 AU that are significantly larger than what is found in the ISM (Bouwman et al. 2008). Crystallisation of amorphous silicates requires temperatures of 800 K. The disk is less than 100 K at 10 AU, so a possible explanation for the high crystalline fractions is that material accreted in the much hotter inner disk and then moved outwards. Although the particular trajectory from Fig. 5 does not reach the required temperature of 800 K, other trajectories do. Crystalline dust has also been detected in several comets, including C/1995 01 (Hale-Bopp) and 81P/Wild 2, providing direct evidence that part of their constituent material has been heated to temperatures well above the evaporation temperature of H₂O (Wooden et al. 1999; Keller et al. 2006).

If Hale-Bopp, Wild 2 and other comets contain crystalline silicates, they must also contain ices that underwent a large degree of chemical processing. Likewise, the presence of amorphous silicates in comets is indicative of chemically unprocessed ices. Hence, the material from which comets are formed must be of mixed origins: some of it accreted close to the star and was heated to high temperatures, while another part accreted at larger radii and remained cold. Within the context of our model, this requires that material from beyond R = 30 AU is radially mixed into the comet-forming zone between 5 and 30 AU. Alternatively, vertical mixing may bring relatively pristine material from higher up in the disk down into the comet-forming zone at the midplane. As for the chemical variations between individual comets, our model shows many examples of abundances changing with position or with time. The variations can thus be explained by having the comets form at different positions in the circumsolar disk, or at different times during the disk’s lifetime.

8. Caveats

One of the model’s three important caveats is that it sets the gas temperature equal to the dust temperature, rather than calculating it explicitly. This approach is valid for the interior of the disk (optically thick to UV radiation), but it breaks down in the surface layers and in the inner parts of the envelope (Kamp & Dullemond 2004). The chemistry in these optically thin regions is mostly controlled by photoprocesses and ion-molecule reactions, neither of which depend strongly on temperature. One aspect that would be affected is the gas-phase production of H₂O through the reactions $\begin{array}{ccc}& & \mathrm{O}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }\mathrm{OH}\mathrm{+}\mathrm{H}\mathit{,}\\ & & \mathrm{OH}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathrm{+}\mathrm{H}\mathit{,}\end{array}$which have activation barriers of 3100 and 2100 K, respectively (Natarajan & Roth 1987; Atkinson et al. 2004). Hence, our model probably underestimates the amount of H₂O in the surface of the disk.

A higher gas temperature would also change some of the physics in the model, such as increasing the pressure of the disk. This would push up the disk-envelope boundary and possibly change the spatial distribution of where material accretes onto the disk. However, the bulk of the accretion currently takes place in shielded regions, so only a small fraction of the infalling material is expected to be affected.

Another caveat is the adopted shape of the stellar spectrum, which may have larger chemical consequences than the gas temperature. Presently, the stellar radiation field is a blackbody spectrum at the effective stellar temperature. Our star never gets hotter than 5800 K (Paper I), which is not enough to produce UV photons of sufficient energy to dissociate CO and H₂. Many T Tauri stars have a UV excess (Herbig & Goodrich 1986), which would allow CO and H₂ to be photodissociated. The higher UV flux would also enhance the photoionisation of atomic C, probably resulting in C⁺ being the dominant form of carbon between points C and D in Fig. 3. In turn, this would boost the production of carbon-chain species like C₂S and HC₃N. As soon as the material enters the disk and is shielded again from the stellar radiation, C⁺ is converted back into C and CO. However, some signature of the temporary high C⁺ abundance may survive.

The third important caveat has to do with the dust acting as a possible sink for the gas-phase chemistry. In our model, once an atom or molecule freezes out, there are three possibilities: (1) it may stay on the grain if the temperature and incident UV flux remain low; (2) it may desorb again if the temperature or UV flux gets high enough; or (3) it may be hydrogenated in a few select cases (Sect. 3.2). In reality, there are at least three additional possibilities: (4) volatile species like CO, CH₄ and N₂ can be trapped in the H₂O ice matrix; (5) the atom or molecule in question may undergo grain-surface reactions not included in our model; or (6) the dust grain may grow and settle to the midplane, effectively removing the molecules in its icy mantle from the active chemistry elsewhere. The latter scenario has been invoked to explain the low H₂O abundance in the circumstellar disk of DM Tau (Bergin et al. 2010). Processes 4, 5 and 6 can all act as a sink of gas-phase CO and other abundant species, potentially diminishing their dominant role in the gas phase. A quantitative treatment of each process is required to assess how much they would alter the chemical evolution presented here.

9. Conclusions

This paper describes the two-dimensional chemical evolution during the collapse of a molecular cloud core to form a low-mass protostar and a circumstellar disk. The model, described in detail by Visser et al. (2009b) and Visser & Dullemond (2010), consists of semi-analytical density and velocity profiles throughout the core and the disk. The dust temperature and the UV radiation field – both important for the chemical evolution – are calculated with a full radiative transfer method. The chemical part of the model features a full gas-phase network, including adsorption and desorption from dust grains, as well as basic hydrogenation reactions on the grain surfaces. Starting from realistic initial conditions, the chemical evolution is computed for parcels of material terminating at selected positions in the disk. Special attention is paid to parcels ending up in the comet-forming zone. Our conclusions are as follows:

• The chemistry during the collapse phase is controlled by a smallnumber of key chemical processes, each of which is activated bychanges in the physical conditions. The evaporation of CO, CH₄and H₂O at approximately 18, 22 and 100 K is one set of such keyprocesses. Another set is the photodissociation of CH₄ and H₂O(Sect. 5.1).

• At the end of the collapse phase, the disk can be divided into several zones with different chemical histories. The different histories are related to the presence or absence of the aforementioned key processes along various infall trajectories. Spectroscopic observations at high spatial resolution are required to determine whether this zonal division really exists, or if it is smoothed out by mixing (Sect. 5.2).

• Part of the material that accretes onto the disk at early times is transported outwards to conserve angular momentum, and may remain in the disk for the rest of the collapse phase. It is heated to well above 100 K as it accretes close to the star, so H₂O and all other non-refractory species evaporate from the grains. They freeze out again when the material cools down during the subsequent outward transport (Sect. 5.2).

• When the chemistry is evolved for an additional 1 Myr at fixed physical conditions after the end of the collapse phase, the abundances of most species change throughout the disk. Hence, the disk is not in chemical equilibrium at the end of the collapse, and one should be careful in choosing the initial abundances when running stand-alone disk chemistry models (Sect. 6).

• Material that ends up in the comet-forming zone undergoes a large degree of chemical processing, including the evaporation and readsorption of H₂O and species trapped in the H₂O ice. This is consistent with the presence of crystalline silicates in comets. However, it is inconsistent with the chemical similarities observed between comets and young protostars, which are indicative of little processing. Hence, it appears that the cometary material is of mixed origins: part of it was strongly processed, and part of it was not. The chemical variations observed between individual comets suggest they were formed at different positions or times in the solar nebula. Fully pristine ices only appear in the upper and outer parts of this particular disk model, although turbulent mixing may move them into the comet-forming zone (Sect. 7).

Acknowledgments

Astrochemistry in Leiden is supported by a Spinoza Grant from the Netherlands Organization for Scientific Research (NWO) and a NOVA grant. R.V. acknowledges partial support from NASA through RSA award No. 1371476 and from the NSF through grant AST-1008800. S.D.D. acknowledges support by NASA grant NNX08AH28G.

References

All Tables

All Figures

	Fig. 1 Step-wise summary of our 2D axisymmetric collapse model. Steps 2 and 4 are semi-analytical, while steps 3 and 5 consist of detailed numerical simulations.
In the text

Fig. 2 Physical conditions (nH, Td, χ and AV) and abundances of the main oxygen-, carbon- and nitrogen-bearing species for the single parcel from Sect. 5.1, as function of time before the end of the collapse (left) and as function of horizontal position during the collapse (right). The grey bars correspond to the points A, B and C from Fig. 3. In each panel, the parcel essentially moves from right to left.

In the text

Fig. 3 Overview of the chemistry along the infall trajectory (solid black curve) of the single parcel from Sect. 5.1. The solid and dashed grey lines denote the surface of the disk and the outflow wall, both at t = tacc = 2.52    × 105 yr. Physical conditions, abundances (black bars: gas; grey bars: ice) and key reactions are indicated at four points (A, B, C and D) along the trajectory. The key processes governing the overall chemistry at each point are listed in the bottom right. The type of each reaction is indicated by colour, as listed in the top left.

In the text

	Fig. 4 Schematic view of the history of H2O gas and ice throughout the disk. The main oxygen reservoir at tacc is indicated for each zone; the histories are described in the text. The percentages indicate the fraction of the disk mass contained in each zone. Note the disproportionality of the R and z axes. The colours have no specific meaning other than to distinguish the different zones.
In the text

	Fig. 5 Infall trajectory for a parcel ending up in zone 6 from Fig. 4, showing where H2O is predominantly present as ice (black) or gas (grey). The two diamonds mark the time in years after the onset of collapse.
In the text

	Fig. 6 Qualitative evolution of some abundances towards the seven zones with different H2O histories from Fig. 4. The horizontal axes show the time (increasing from left to right) and are non-linear. The position of points A, B, C and D from Fig. 3 is indicated for zone 2, which contains our reference parcel.
In the text

	Fig. 7 As Fig. 4, but for CH4 gas and ice. The main carbon reservoir at tacc is CO gas throughout the disk.
In the text

	Fig. 8 As Fig. 6, but for the five zones with different CH4 histories from Fig. 7.
In the text

	Fig. 9 As Fig. 4, but for NH3 gas and ice. The main nitrogen reservoir at tacc is N2 gas throughout the disk.
In the text

	Fig. 10 As Fig. 6, but for the eight zones with different NH3 histories from Fig. 9.
In the text

	Fig. 11 Abundances of total H2O, O2, CO, CH4, N2, NH3 and NO (gas and ice combined) throughout the disk at the end of the collapse phase (t = tacc) and after an additional 1 Myr post-collapse phase (t = tacc + 1 Myr).
In the text

Fig. 12 Top: abundances observed in comets versus gas-phase abundances observed in young protostars. Middle and bottom: abundances observed in comets versus total abundances (gas and ice) from our model at t = tacc. All abundances are relative to nH. For the observational data points, the error bars indicate the spread between sources. For the model data, they indicate the range of values across the entire disk (middle) or the comet-forming region (bottom). Dotted error bars indicate the full range extends to less than 10-12. The diagonal dotted lines represent hypothetical one-to-one relationships between abundances on the x and y axes.

In the text