ICE CHEMISTRY IN STARLESS MOLECULAR CORES

We investigate ice chemistry of four starless cores with their central density ${n}_{{\rm H}}=n({\rm H})+2n({{\rm H}}_{2})$, in cm⁻³, and interstellar extinction A_V, in mag, derived from observations. Only the core center part was considered (i.e., 0D model), where ${n}_{{\rm H}}$ and A_V are at their highest. Each simulation consists of a collapse Phase 1, when the cloud density increases, and a quiescent Phase 2 with constant physical conditions.

A single-point 0D model was chosen because the ices in the center of each core most likely do exist, despite a number of uncertain parameters. These include the efficiency of non-thermal desorption mechanisms and turbulent motions in starless cores, as observed by Redman et al. (2006), Levshakov et al. (2014), and Steinacker et al. (2014). Turbulence was not considered in the present study. Also, ices in the core center likely have experienced the full range of physical conditions of the core since its formation. The simple 0D model puts certain constraints for the interpretation of results.

Cores with a Plummer-like density profile (Plummer 1911; Whitworth & Ward-Thompson 2001), embedded in a surrounding cloud with a fixed column density ${N}_{\mathrm{out}}$, were modeled. The total column density to the core center can be calculated with

$$\begin{eqnarray}&&{N}_{{\rm H}}={\int }_{{r}_{1}}^{0}{\displaystyle \frac{{n}_{{\rm H}}}{{[1+{(r/{r}_{0})}^{2}]}^{\eta /2}}}{dr}+{N}_{\mathrm{out}},{\mathrm{cm}}^{-2},\end{eqnarray} \tag{ 1 }$$

where ${n}_{{\rm H},0}$ is the (maximum) density at the center of the core, r₀ is the radius of the central density plateau, η is the power-law slope at large radii, and r₁ is the radius of the core. The chemical modeling was done only for the center of the core, where r = 0 and the density profile of the core is necessary only for the calculation of N_H at its center. Table 1 summarizes the parameters ${n}_{{\rm H},0}$, r₀, r₁, and ${N}_{\mathrm{out}}$. These are unique for each core and were derived from Lippok et al. (2013) for the CB 17, CB 26, CB 27, and B68 cores and from Schmalzl et al. (2014) for the CB 17 core. The two "versions" of CB 17 are hereafter denoted CB 17 S and CB 17 L, for data from Schmalzl et al. (2014) and Lippok et al. (2013), respectively. This means that there are a total of five models (CB 17 S, CB 17 L, CB 26, CB 27, and B68) of four cores considered.

Table 1.  The Constant Physical Parameters of Modeled Starless Cores during Their Stable Phase 2

	Observational Data						Calculated Data
Core	${n}_{{\rm H},0}$, cm−3	r0, AU	r1, AU	η	${N}_{\mathrm{out}}$, cm−2	Reference	AV, mag	T, K	t1, kyr
CB 17 S—SMM1	2.3E+5	9.5E+3	3.0E+4	4.9	0	Sa	13.8	8.2	1465
CB 17 L—SMM	1.3E+5	1.1E+4	3.0E+4	5.0	3.6E+20	Lb	9.1	9.2	1451
CB 26—SMM	8.7E+4	8.0E+3	4.0E+4	3.0	2.2E+20	L	6.5	11.1	1437
CB 27—SMM	1.3E+5	1.3E+4	6.6E+4	6.0	3.1E+20	L	9.5	8.9	1451
B68—SMM	4.0E+5	7.0E+3	2.7E+4	5.0	3.0E+20	L	17.6	7.9	1482

Note. Assumed physical parameters: ${n}_{{\rm H},0}$—density at the center of the core; r₀—radius of the central density plateau; η—power-law slope for density; r₁—radius of the region (core) in consideration; ${N}_{\mathrm{out}}$—column density of the surrounding cloud. Calculated parameter t₁ is the length of the cloud core-collapse Phase 1.

^aSchmalzl et al. (2014). ^bLippok et al. (2013).

Download table as:  ASCIITypeset image

Table 1 data are relevant to cores in the stable Phase 2, when the contraction in Phase 1 has ended and the physical parameters of the cores are unchanging. This phase has an assumed integration time of ${t}_{2}\;\gt$ 2 Myr for each modeled core. Prior to the stable phase, a core formation phase was considered. The formation of starless cores has been described before as a free-fall (e.g., Lee et al. 2004) or "delayed free-fall" (e.g., Taquet et al. 2014). A "static," step-like contraction was employed by these authors.

In the present study, the formation of the core was described as a delayed free-fall process. The approach presented by Brown et al. (1988), Nejad et al. (1990), and Rawlings et al. (1992) was used. A Plummer-like core with an initial peak density of ${n}_{{\rm H},\mathrm{ini}}=3000$ cm⁻³ contracts to the final density ${n}_{{\rm H},0}$ in a core formation time t₁ (i.e., length of Phase 1; see Table 1). The parameter B in Equation (1) of Rawlings et al. (1992), which modifies the rate of the contraction, was taken to be 2/3. This means that the collapse to the final densities, indicated in Table 1, occurs in 1.4–1.5 Myr. Figure 1 shows an example of the evolution of physical conditions for a core (B68) during the contraction phase, up to 1.5 Myr.

Zoom In Zoom Out Reset image size

A value B = 0.7 has been used by Nejad et al. (1990), Rawlings et al. (1992), and Ruffle et al. (1999) for prestellar core models. We take $B=2/3$ as a somewhat lower value for starless cores. The exact value of B has only a limited effect on ice composition, the main subject of the present study. This is because ice accumulation was calibrated by modifying the poorly known efficiency of indirect reactive desorption to match the observed proportions of H₂O, CO, and CO₂ (Section 2.1.4).

For the core formation epoch, r₀ was kept proportional to ${n}_{{\rm H},0}^{-1/2}$ (Keto & Caselli 2010; Taquet et al. 2014). The Plummer-like density profile of the contracting core was retained, while the density increases from ${n}_{{\rm H},\mathrm{ini}}$ to ${n}_{{\rm H},0}$. The value of r₁, which determines the physical size of the region in consideration, was assumed constant. This means that the mass of the core (i.e., the number of H atoms within a sphere of radius r₁) increases during the formation phase. ${N}_{\mathrm{out}}$ was retained constant, too.

The total integration time ${t}_{1}+{t}_{2}$ was taken to be 3.5 Myr, and the results were analyzed for an interval of 1.0–2.3 Myr (see Section 3.3). In context with ice formation, the full time consists of three periods: (1) cloud with ${n}_{{\rm H},0}\approx 3\times {10}^{3}$ cm⁻³ and no ice on the grains; (2) ice accumulation up to a maximum thickness in a rapidly contracting core and up to 200 kyr into the quiescent phase; and (3) slow processing of ice by cosmic-ray-induced photons after the freeze-out has ceased. The start of period 2 is marked with the appearance of first subsurface ice species at densities 5 × 10³ to 1 × 10⁴ cm⁻³. During the last period (3), the gas and the surface are in an approximate adsorption–desorption equilibrium. The nominal ice thickness decreases during period 3.

Ice photoprocessing leads to a declining amount of CO in the icy mantles. Because of this, quiescent Phase 2 integration times longer than 2 Myr were not considered. For longer times, the H₂O:CO:CO₂ ice abundance ratio is no longer supported by observations of interstellar ices (Gibb et al. 2004; Whittet et al. 2007; Boogert et al. 2011; Öberg et al. 2011). Mantle growth and the abundances of major ice species are considered in Section 3.

Gas and dust were assumed to be in thermal equilibrium with temperature derived from the interstellar extinction A_V, in line with Garrod & Pauly (2011). In other words, it was assumed that the cores are heated externally by the interstellar radiation field. A_V was calculated with

$$\begin{eqnarray}&&{A}_{V}={\displaystyle \frac{{N}_{{\rm H}}}{1.60\times {10}^{21}}},\end{eqnarray} \tag{ 2 }$$

Table 1 shows the calculated values for A_V and T for each of the modeled cores in the stable phase.

Table 2 shows the chemical abundances used at the start for all simulations. They are based on the abundances derived by Garrod & Herbst (2006) and Wakelam & Herbst (2008). The rate of hydrogen ionization by cosmic rays was taken to be $1.3\times {10}^{-17}$ s⁻¹. The flux of cosmic-ray-induced photons and interstellar photons was taken to be 4875 and 10⁸ cm⁻² s⁻¹, respectively, with a grain albedo of 0.5.

Table 2.  Initial Abundances of Chemical Species

Species	${n}_{i}/{n}_{{\rm H}}$
H2	0.4995
H	1.00E-03
He	9.00E-02
C	1.40E-04
N	7.50E-05
O	3.20E-04
Na	2.00E-08
Mg	2.55E-06
Si	1.95E-06
P	7.63E-08
S	1.50E-06
Cl	1.40E-08
Fe	7.40E-07

Download table as:  ASCIITypeset image

The self- and mutual shielding of H₂, CO, and N₂ was included with the use of the tabulated shielding functions (Lee et al. 1996; Li et al. 2013). Shielding was also attributed to molecules in ice. It is important to note that the surrounding molecular cloud was included in these calculations by assuming an external H₂ column density of ${N}_{\mathrm{out}}/2$. This was not included in Papers I and II, which used a physical model built on similar principles. Such a more accurate approach here is of particular importance, because more plausible values for indirect reactive desorption efficiency can now be derived (Section 2.1.4). The synthesis of hydrogen on the grain surface was considered by both Langmiur–Hinshelwood and Eley–Rideal mechanisms, although the latter is of little importance.

The grains were assumed to consist of two components. The first is inert, solid nuclei with a radius of $a=0.1\;\mu {\rm m}$ and an abundance $1.32\times {10}^{-12}$ relative to H atoms. In dark and dense conditions (${A}_{V}\gt 1.9$ mag), the second component, interstellar ice, begins to accumulate onto grain surfaces in significant amounts. The increase of grain size is calculated self-consistently, assuming a thickness of 350 pm for each ML, which is multiplied by the number of ice MLs to obtain ice thickness b. The number of adsorption sites per ML was retained constant at $1.5\times {10}^{6}$. Changes in this number likely do not influence the agreement of model results with observations (Acharyya et al. 2011; Taquet et al. 2012) and were not considered in the model. Dust grain albedo was assumed to be 0.5.

The ice layer on the grains was described with the sublayer approach, developed in Paper II. This means that a fully formed ice (accumulation from gas has ceased) consists of a surface layer (1–2 MLs) and three subsurface layers—"sublayers"—of approximately equal thickness. A sublayer nominally may contain up to several tens of MLs. Two or more sublayers in a model allow us to consider the depth-dependent composition of the ice mantles. Three sublayers were considered in the present study, similarly to Paper II. For reference, they were numbered with the first being the shallowest and the third being the sublayer adjacent to the grain nucleus. The transition from the surface to the first sublayer and, subsequently, to the second and third sublayers has been harmonized with the accretion rate, in line with Papers I and II.

The molecules are bound to the surface or into the mantle with their adsorption (E_D) or absorption energies (E_B), respectively. E_D is the "basic value," given in the reaction network, from which all the other energies were calculated. The binding energies, used for calculating the diffusion rates, for surface (${E}_{b,S}$) and mantle (${E}_{b,M}$) species were assumed to be 0.40E_D and 0.40E_B, respectively (Garrod & Pauly 2011). In Papers I and II it was argued that ${E}_{b,M}$ has to be higher, or similar to E_D; otherwise, the movement of molecules in ice bulk would be faster than their evaporation from the surface, which fits the description of a liquid, not a solid. In line with the earlier works, it was assumed that ${E}_{B}=3{E}_{D}$, i.e., that typically ${E}_{b,M}=1.2{E}_{D}$. For CO and CO₂ molecules, the ratios ${E}_{b,S}/{E}_{D}={E}_{b,M}/{E}_{B}$ were taken to be 0.31 and 0.39, respectively (Karssemeijer & Cuppen 2014).

In dark cores, the grains may be sufficiently cold to have a significant proportion of their surface covered by molecular hydrogen. H₂ provides a much weaker binding than water-ice. The adsorption energy for H and H₂ was calculated according to Garrod & Pauly (2011). For other species, E_D was not affected by the surface coverage of H₂ because of the large difference between the hopping rate of H₂ and vibration frequency of heavy molecules (Paper II).

Accretion of neutral gaseous species onto the grains with a radius $a+b$ was calculated with sticking probabilities of 0.33 for light species (H, H₂) and 1.0 for all other species (Kalvāns & Shmeld 2010). A total of six desorption mechanisms were considered. The rate coefficients for evaporation and cosmic-ray-induced whole-grain heating (Hasegawa & Herbst 1993a) were calculated for each molecule as a thermal desorption at the equilibrium temperature and 70 K, respectively. Desorption by interstellar and cosmic-ray-induced photons was considered with uniform yields for all species, 0.003 and 0.002, respectively. These values were chosen to be in agreement with recent photodesorption experiments with interstellar ice analogs, as discussed in Paper II. Photodesorption from subsurface mantle species was permitted, if the number of overlying MLs does not exceed two. The yield for subsurface photodesorption was reduced by a factor of $1/3$. This approach was developed to be in line with the results of Andersson & van Dishoeck (2008).

A reaction-specific reactive desorption was used, with the ${\alpha }_{\mathrm{RRK}}$ parameter set to 0.03 (Paper II). This is the ratio of the surface-molecule bond frequency to the frequency at which energy is lost to the grain surface (Garrod et al. 2007).

Desorption induced by the energy released from the H+H exothermic surface reaction probably determines the relative proportions of major ice components H₂O, CO, and CO₂. In order to reproduce observations of interstellar ices, a sequence $\mathrm{CO}\gt {\mathrm{CO}}_{2}\gg {{\rm H}}_{2}{\rm O}$ must be observed for indirect reactive desorption efficiency ${f}_{{{\rm H}}_{2}\mathrm{fd}}$ (Paper II). Without indirect reactive desorption, the calculated abundance of solid carbon species in collapsing cores can be overestimated, when compared to the observational results provided by Gibb et al. (2004), Whittet et al. (2007), and Öberg et al. (2011). For example, realistic CO and CO₂ relative abundances can be achieved, if the ice contains a substantial amount of CH₄, H₂CO, and (or) CH₃OH (Garrod & Pauly 2011; Taquet et al. 2014). Of the latter three species, only methanol has been occasionally observed in interstellar ices (Boogert et al. 2011).

The parameter ${f}_{{{\rm H}}_{2}\mathrm{fd}}$ can be described as the number of desorbed molecules per H₂ formed on a grain. Exact values for the efficiency are unknown, and given the averaging and uncertainties in astrochemical numerical simulations, it is reasonable to adjust ${f}_{{{\rm H}}_{2}\mathrm{fd}}$ so that it produces an approximate fit with observations of starless cores. Such an approach can be used until a more precise value of ${f}_{{{\rm H}}_{2}\mathrm{fd}}$ is found in experiments or, for example, quantum chemistry calculations.

In order to derive ${f}_{{{\rm H}}_{2}\mathrm{fd}}$ for the present study, we employ a similar method to Paper II. A prestellar core in free-fall gravitational collapse was considered in that study. For the present research, the B68 core model (Table 1) was used. It has the highest final density and final A_V among all the starless cores in consideration. This allows the comparison of observed/modeled CO:H₂O and CO₂:H₂O ice abundance ratios for five A_V values from the data set of Whittet et al. (2007). An additional observational constraint for the model is the abundance of H₂O, CO, and CO₂ ices at their respective threshold extinctions ${A}_{\mathrm{th}}$.

The above means that, while the B68 cloud core is currently stable (Bergin et al. 2006; Redman et al. 2006), we assume that it formed in a delayed free-fall process, which ceased when the core reached its present density. This assumption has been attributed to all simulated cores in this study. While it might not be entirely correct, it is probably just as useful as other realistic assumptions, because of the poorly constrained value of ${f}_{{{\rm H}}_{2}\mathrm{fd}}$, which has to be calibrated, regardless of the exact core formation path.

Table 3 summarizes the observational and best-fit B68 model results. The desorption efficiency ${f}_{{{\rm H}}_{2}\mathrm{fd}}$ was calculated according to the empirical approach of Paper II:

Table 3.  Comparison of Observed and Calculated (B68 Core) CO and CO₂ Ice Abundances, Relative to Water

		Observations		B68 Model
Source ID	AV	$$\begin{eqnarray}&&{{\displaystyle \frac{\mathrm{CO}}{{{\rm H}}_{2}{\rm O}}}}^{{\rm a}},\%\end{eqnarray}$$	$$\begin{eqnarray}&&{\displaystyle \frac{{\mathrm{CO}}_{2}}{{{\rm H}}_{2}{\rm O}}},\%\end{eqnarray}$$	$$\begin{eqnarray}&&{\displaystyle \frac{\mathrm{CO}}{{{\rm H}}_{2}{\rm O}}},\%\end{eqnarray}$$	$$\begin{eqnarray}&&{\displaystyle \frac{{\mathrm{CO}}_{2}}{{{\rm H}}_{2}{\rm O}}},\%\end{eqnarray}$$	tb, kyr
043728.2 + 261024	6.3 ± 1.5	. . .	25.0	1.8	13.8	1417
042324.6 + 250009	10.0 ± 0.5	20.2	19.1	7.7	13.0	1451
043325.9 + 261534	11.7 ± 0.5	12.0	16.7	12.5	13.9	1459
043926.9 + 255259	15.3 ± 0.5	27.3	16.7	30.2	16.5	1471
042630.7 + 243637	17.8 ± 1.5	45.1	18.3	43.2	18.2	1476
Threshold Ice Abundances			H2O, ML	CO2, ML	CO, ML
${A}_{\mathrm{th}}$ H2O	3.2 ± 0.1	⋯	3.5	⋯	⋯	1332
${A}_{\mathrm{th}}$ CO2	4.3 ± 1.0	⋯	⋯	1.5	⋯	1375
${A}_{\mathrm{th}}$ COc	6.8 ± 1.6	⋯	⋯	⋯	0.6	1424
${A}_{\mathrm{th}}$ COd	8.1	⋯	⋯	⋯	1.3	1437

Notes. Observational data from Whittet et al. (2007), unless noted otherwise.

^aAbundance ratio for icy species. ^bIntegration time, corresponding to the A_V value. ^cWhittet et al. (2001). ^dPaper II.

Download table as:  ASCIITypeset image

$$\begin{eqnarray}&&{f}_{{{\rm H}}_{2}\mathrm{fd}}=Q\mathrm{exp}(-{E}_{D}/({Q}_{1}),\end{eqnarray} \tag{ 3 }$$

where $Q=9.5\times {10}^{-5}$ and ${Q}_{1}={10}^{3}$. The adopted best-fit values of Q and O₁ are different from those in Paper II. It is possible that the values obtained here are more reliable than those of Paper II because of the more advanced and realistic physical model (Section 2.1.1).

The observational values for CO:H₂O and CO₂:H₂O ice abundance ratios were taken from Whittet et al. (2007). While newer data sets are available (Boogert et al. 2011; Öberg et al. 2011), Whittet et al. (2007) provide the A_V for each background star, which enables us to tie the observations to specific time points of the simulations with the corresponding A_V (see Paper II for more on this approach).

Several changes have been introduced into the gas-grain network, which was acquired from Laas et al. (2011). The four new gas-phase reactions of Vasyunin & Herbst (2013b) were added. All these reactions involve the methoxy radical CH₃O. Because these authors did not discern between the methoxy radical and its isomer, the hydroxymethyl radical CH₂OH, four complementary reactions were added for the latter species (1–4 in Table 4).

Table 4.  Changes in Reactions Respective to the Original Network by Laas et al. (2011)

No.	New Gas Phase Reactionsa	α	β	γ	Reference
1	${\mathrm{CH}}_{2}\mathrm{OH}+{\rm H}\;\longrightarrow \;{\mathrm{CH}}_{3}+\mathrm{OH}$	1.6E-10	0	0	NIST
2	${\mathrm{CH}}_{2}\mathrm{OH}+{\rm O}\;\longrightarrow \;{{\rm H}}_{2}\mathrm{CO}+\mathrm{OH}$	1.5E-10	0	0	NIST
3	${\mathrm{CH}}_{2}\mathrm{OH}+{\mathrm{CH}}_{3}\;\longrightarrow \;{{\rm H}}_{2}\mathrm{CO}+{\mathrm{CH}}_{4}$	1.4E-10	0	0	NIST
4	${\mathrm{CH}}_{2}\mathrm{OH}+{\mathrm{CH}}_{3}\;\longrightarrow \;{{\rm C}}_{2}{{\rm H}}_{5}\mathrm{OH}$	1.0E-15	−3.00	0	Vasyunin & Herbst (2013b)b
No.	New Surface (ice) Reactions	EA, K	⋯	⋯	Reference
5	${\rm H}+\mathrm{HCO}\;\longrightarrow \;\mathrm{CO}+{{\rm H}}_{2}$	0	⋯	⋯	This work
6	$\mathrm{OH}+\mathrm{HCO}\;\longrightarrow \;\mathrm{CO}+{{\rm H}}_{2}{\rm O}$	0	⋯	⋯	This work
7	${\mathrm{NH}}_{2}+\mathrm{HCO}\;\longrightarrow \;\mathrm{CO}+{\mathrm{NH}}_{3}$	0	⋯	⋯	This work
8	${{\rm H}}_{2}+\mathrm{HCO}\;\longrightarrow \;{\mathrm{CH}}_{3}{\rm O}$	2100c	⋯	⋯	This work
9	$\mathrm{OH}+\mathrm{SO}\;\longrightarrow \;{\mathrm{SO}}_{2}+{\rm H}$	0	⋯	⋯	Gas phase
10	${\rm O}+\mathrm{HS}\;\longrightarrow \;{\rm S}+\mathrm{OH}$	956	⋯	⋯	Gas phase
11	${\rm O}+\mathrm{HS}\;\longrightarrow \;\mathrm{SO}+{\rm H}$	0	⋯	⋯	Gas phase
12	${\rm H}+\mathrm{CL}\;\longrightarrow \;\mathrm{HCL}$	0	⋯	⋯	Gas phase
No.	Changed Surface (ice) Reactions	EA, K	⋯	⋯	Reference
13	$\mathrm{CO}+{\rm O}\;\longrightarrow \;{\mathrm{CO}}_{2}$	290	⋯	⋯	Roser et al. (2001)
14	${{\rm O}}_{2}+{\rm H}\;\longrightarrow \;{{\rm O}}_{2}{\rm H}$	600	⋯	⋯	Du et al. (2012)
15	${\rm H}+{{\rm H}}_{2}\mathrm{CO}\;\longrightarrow \;{\mathrm{CH}}_{3}{\rm O}$	2100	⋯	⋯	Woon (2002)
16	${{\rm O}}_{3}+{\rm H}\;\longrightarrow \;{{\rm O}}_{2}+\mathrm{OH}$	0	⋯	⋯	Mokrane et al. (2009)

Notes.

^aIn addition to the four reactions from Vasyunin & Herbst (2013b). Reaction rate coefficient $k=\alpha {(T/300)}^{\beta }\mathrm{exp}(-\gamma /T)$. ^bData adopted from the ${\mathrm{CH}}_{3}{\rm O}+{\mathrm{CH}}_{3}$ reaction. ^cEstimate.

Download table as:  ASCIITypeset image

The results of Paper II show that formic acid HCOOH, formaldehyde H₂CO, and formamide NH₂CHO can be overproduced in subsurface mantle conditions with the surface reaction set of Laas et al. (2011), which has been adopted from Garrod et al. (2008). In an attempt to reduce their abundance, grain surface reactions 5–8 in Table 4 have been added to the network. They tend to reduce the production efficiency of HCOOH, H₂CO, and NH₂CHO because an alternative outcome is now possible for their reactants.

Reactions 5–7 were assumed to be barrierless. Reaction 8 has an assumed activation barrier E_A = 2100 K, estimated from reactions of H₂ with compounds that have a similar structure to HCO. The barrier of the latter reaction is so high that it has little effect on the results of the current cold core model.

Reactions 9–11 were added to more adequately reproduce the chemistry of sulfur oxides, with activation barriers taken from similar gas-phase reactions. Finally, reaction 12 was added to avoid most of the chlorine atoms ending up in solid ClO, which is probably unrealistic, given the availability of atomic hydrogen, irradiation, and an aqueous medium on grain surfaces. The activation barriers for four additional surface reactions—13–16 in Table 4—were also changed, in line with available data.

The rate for reactions on grain surface was calculated according to the reaction–diffusion competition approach of Garrod & Pauly (2011). The modified rate equations method has been applied (see Paper II). As explained in Section 2.1.3, part of the surface reaction products desorb into the gas. Reaction barriers are overcome either thermally or by quantum tunneling, whichever is faster (Garrod et al. 2008).

The rate of binary reactions in bulk ice was calculated for reactants that are immobile for most of the time (Papers I and II). It was assumed that each molecule vibrates in a lattice cell with 10 neighboring molecules. Reactants have to overcome a certain energy barrier ${E}_{\mathrm{prox}}=0.1{E}_{D}$ in order to achieve sufficient proximity for a reaction to occur. Molecule diffusion in ice and reaction–diffusion competition have been taken into account in a manner similar to that for surface reactions. I refer the reader to Papers I and II for a more detailed justification and formalism for the calculation of reaction rate coefficients.

Surface and bulk-ice molecules are subjected to photodissociation by interstellar and cosmic-ray-induced photons. For surface species, the dissociation rate has been taken to be equal to that of gas-phase species. For molecules within the mantle, the photon flux is attenuated by overlying MLs. Each ML has an assumed absorption probability of 0.007 (Andersson & van Dishoeck 2008). In order to calculate the photon flux, respective to each sublayer, the value relevant to the middle ML of that sublayer was used (see Paper II).

Three major ice species—H₂O, CO, and CO₂—contain copious amounts of the eighth element, and the chemistry of oxygen in ice is closely connected to these molecules. When water and carbon oxides are excluded, the bulk-ice chemistry of oxygen for CB 17 S and B68 is initially dominated by molecular oxygen O₂, which is rather quickly converted into hydrogen peroxide H₂O₂, as shown in Figure 5. For other models, where ice photoprocessing is more intense, H₂O₂ is the more important species from the beginning. H₂O₂ is more stable in irradiated interstellar ices than O₂, which mainly arrives from the gas phase. The hydroperoxyl radical O₂H is the most abundant radical for oxygen chemistry in ice, thanks to its high ${E}_{b}$. It is generated by two mechanisms—H addition to O₂ (early times) and the photodissociation of H₂O₂ (late times). Hydroxyl radical OH begins to accumulate to a limited extent in the inner sublayer 3 only after the CO molecule has been oxidized (Figure 6).

Zoom In Zoom Out Reset image size

Interestingly, CB 26, which has the lowest A_V and highest flux of interstellar photons, has also the lowest radical content. This is because the higher temperature in CB 26 (see Table 1) allows a greater radical mobility (especially for atomic H) in the mantle and, thence, reactivity. Low radical content in CB 26 ices is characteristic also for other compound classes (Sections 4.2, 4.3, 5.1).

Because ice in CB 17 S is initially rich in O₂ (Figure 5), it has also a higher abundance for H₂O₂, which forms via reaction 14 (Table 4). Substantial amounts of H₂O₂ (relative abundance ${X}_{{{\rm H}}_{2}{{\rm O}}_{2}}$ in excess of 10⁻⁶) in sublayer 3 for CB 17 S and B68 are able to form from water photodissociation products that combine together after all CO has been oxidized in that sublayer.

Zoom In Zoom Out Reset image size

Figure 7 shows that the calculated abundance of gaseous hydrogen peroxide H₂O₂ can be reasonably high for a short period, up to 9 × 10⁻¹⁰ relative to the total hydrogen atomic number density for CB 17 S. The abundance of its associated O₂H radical, however, is lower by three orders of magnitude. Both these species have been observed in interstellar medium with roughly similar abundances of 5 × 10⁻¹¹ (Bergman et al. 2011; Parise et al. 2012). It has been suggested that reactive desorption is responsible for transferring these species from the surface to the gas (Du et al. 2012). However, reactive desorption may be too inefficient to maintain a significant amount of such heavy species in the gas (Paper II). The comparison of abundances here is indicative, only, because of the limits of the 0D model.

Zoom In Zoom Out Reset image size

Because practically all of H₂O₂ in ice resides in the sublayers, photodesorption from subsurface ice may play a role in the existence of gas-phase H₂O₂ and O₂H. More often than not, dissociative photodesorption occurs, when molecule fragments are ejected, instead of intact species (Andersson et al. 2006; Andersson & van Dishoeck 2008). This aspect has not been considered in the present model, or in the model by Du et al. (2012). It is likely that photodesorption of H₂O₂ also induces a high gas-phase abundance of O₂H.

Beginning with Goldsmith et al. (2002), molecular oxygen O₂ has been repeatedly observed in shocked gas. A few detections have been associated with the cold dense core ρ Ophiuchi A, with inferred abundances in the range $2.5\times {10}^{-7}\;\mathrm{to}\;1.25\times {10}^{-6}$ relative to total hydrogen. Such high abundances likely can be observed only during a short period in the core's evolution (Larsson et al. 2007; Liseau et al. 2010). These results agree with the calculation results shown in Figure 7, where the maximum abundances lie in the range $7\times {10}^{-7}\;\mathrm{to}\;4.5\times {10}^{-6}$. A quantitative comparison is possible in this particular case because of the centrally localized nature of the observed source (Liseau et al. 2010).

Figures 8 and 9 show that ammonia NH₃ and molecular nitrogen N₂ are the main nitrogen reservoirs in ice and gas. Nitrogen chemistry is regulated by proportions of N₂ and NH₃ in the surface and the three sublayers (Section 3.2). The relatively clear division of NH₃ (in sublayer 3) and N₂ (in sublayers 1 and 2) means that N, NH, and NH₂ radicals are not equally available in all sublayers. Interesting examples of relatively abundant minor nitrogen species in ice are formamide NH₂CHO (discussed in Section 5), hydrogen cyanide HCN, and isocyanide HNC.

Zoom In Zoom Out Reset image size

The chemistry of HCN and HNC is controlled by the proportions of CO, N₂, and NH₃ in each sublayer. HNC is more favored by ice chemistry at earlier times than cyanide HCN, while HCN is the more favored isomer during late times. The two isomers form via different pathways. HNC is mainly synthesized in the reaction NH₂+C, while HCN is the product of H addition to CN. The latter forms via C+N. This means that HNC is a product of NH₃ and CO mixture photoprocessing, while HCN is synthesized in an N₂ and CO mixture. The HNC pathway is more efficient because it requires no subsequent hydrogenation. The result is that 60% (for CB 17 S) to 90% (for CB 26) of both isomers is concentrated in the outer sublayer 1, where CO and N₂ are highly abundant, and NH₃ is available, too. The remaining 10%–40% resides in sublayer 2, while sublayer 3 is very poor on HCN and HNC (<2%), because of a lack of CO there. Because the surface closely interacts with gas-phase CO and N₂, HCN is the main isomer for the top ice layer.

Recent observed NH₃ gas-phase abundances lie in the range $9\times {10}^{-9}...7.5\times {10}^{-8}$ relative to total H (Ruoskanen et al. 2011; Levshakov et al. 2014). Crapsi et al. (2007) derives a central NH₃ abundance of $4.5\times {10}^{-9}$ for the L1544 dense core. These values qualitatively agree with calculation results for cores with ages in the range t = 1.5–1.8 Myr.

Zoom In Zoom Out Reset image size

Figure 10 shows the calculated HNC:HCN abundance ratio of the center of the respective dark cores in gas and solid phases. We note that the observed gas-phase values of this ratio are in the range of 0.54–4.5, with no distinction between YSOs and starless cores (Hirota et al. 1998). The HNC:HCN ratio in ice may exceed unity for all cores, except for CB 26, where lower abundance of CO prevents the buildup of HNC. The highest calculated HNC:HCN ratio is in the CB 27 ice phase, where it reaches 1.75 briefly at 1.54 Myr. Such a peak arises because of effective HNC production in sublayer 1, not observed in other simulations. In turn, this occurs because sublayer 1 in CB 27 is richer in NH₃ by a factor of >2 than in other models. Finally, NH₃ is more abundant because ices in CB 27 are accumulated at lower extinctions, which means a higher proportion of species that form from free atoms on the surface (H₂O, NH₃, CH₄, etc.).

Zoom In Zoom Out Reset image size

The "Standard" model here is the one described in full in Section 2. Figures 14 and 15 show the calculated abundances in the center of the starless cores, along with known observational values. The species can be categorized in several overlapping groups:

• 1.  
  those synthesized via the HCO radical, e.g., formic acid HCOOH, formaldehyde H₂CO, methyl formate HCOOCH₃, formamide NH₂CHO, and acetaldehyde CH₃CHO;
• 2.  
  those synthesized via the CH₃O radical, e.g., methanol; CH₃OH, dimethly ether CH₃OCH₃, HCOOCH₃, and as a subgroup, species whose synthesis involves CH₃, which is a daughter of methanol. This subgroup includes CH₃CHO, acetonitrile CH₃CN, and methane CH₄;
• 3.  
  carbon chain molecules, e.g., CH₃CHO, propynone CH₂CCO, ketene CH₂CO, and acetylene C₂H₂;
• 4.  
  species that are typically formed via atom addition on the surface, e.g., CH₄, H₂CO, HCO, HCOOH CH₂CCO, CH₂CO, C₂H₂, methyl imine CH₂NH, cyanamide NH₂CN, and CH₃CN. A relatively high gas-phase abundance during the freeze-out epoch is characteristic for these species.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The hydrogenation of CO and that of H₂CO are the two most important steps that regulate the general production efficiency for many oxygen-containing COMs in bulk ice. Figures 16 and 17 show the relative abundances of organic species in the four ice layers considered in the models for CB 17 L. This core can be regarded as a "median" case among the five models, without the extremes in ice composition shown by CB 17 S and B68, on one hand, and CB 26 and CB 27, on the other hand (Section 3.1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most of the species with abundances shown in Figures 16 and 17 have their peak abundances highest in the outermost sublayer 1. This indicates that even a partial ice sublimation might be sufficient to reproduce the observed gas-phase abundances of many COMs. In sublayer 1, CO is available in vast quantities at any integration time points, while the necessary H atoms are produced in H₂O or H₂ photodissociation. With the present 0D model it is not possible to determine the percentage of ice mass that should be sublimated. If the modeled ices are assumed to be representative, this proportion would be $\approx 1\%-10$%.

Figures 14 and 15 show that the model predicts relatively high gas-phase abundance for some species in the HCO group—formic acid HCOOH, formaldehyde H₂CO, and formamide NH₂CHO. Although the observed gas-phase abundances are reached for a short period of time, this occurs during the infall Phase 1 and likely cannot be treated as an agreement with observations of stable starless cores.

For the CH₃O group, the calculated abundances are significantly lower than those observed in interstellar conditions. Even assuming instant ice sublimation, observed gas-phase abundances cannot be reproduced for methyl formate HCOOCH₃, dimethyl ether CH₃OCH₃, and methanol CH₃OH. This result is similar to that of earlier prestellar core models without (Garrod & Herbst 2006; Garrod et al. 2008) and with bulk-ice photochemistry (Garrod 2013; Paper II). It is this discrepancy that prompts us to proceed with a more detailed investigation on the chemistry of COMs.

In the present surface reaction network (Garrod et al. 2008) the activation energy barriers E_A for reactions H+CO and H+H₂CO are higher than 2000 K, which means that they are inefficient for cold core conditions. Experimental evidence points to much lower values for E_A, 300–600 K (Awad et al. 2005; Fuchs et al. 2009). However, putting such low E_A values into the present model implies that almost all carbon in ice is in the form of methanol, which is unrealistic. Higher E_A is especially required for the CO+H reaction.

We present results of model calculations that use the experimental, low E_A value for formaldehyde hydrogenation. For CO hydrogenation E_A, an empirical compromise value was used, which lies between the low experimentally detected E_A of Fuchs et al. (2009) and the high (2500 K) value proposed by Garrod et al. (2008). Activation energies for the three important CO hydrogenation reactions have been summarized in Table 6. We dub this the "Compromise model" for further reference.

Table 6.  Changes in Reaction Data for the Compromise Model

		Activation Barrier EA, K
No.	Reaction	Garrod et al. (2008)a	Fuchs et al. (2009)	Compromise
17	${\rm H}+\mathrm{CO}\;\longrightarrow \;\mathrm{HCO}$	2500	390–520	1600
18	${\rm H}+{{\rm H}}_{2}\mathrm{CO}\;\longrightarrow \;{\mathrm{CH}}_{3}{\rm O}$	2100	415–490	415
19	${\rm H}+{{\rm H}}_{2}\mathrm{CO}\;\longrightarrow \;{\mathrm{CH}}_{2}\mathrm{OH}$	2500	415–490	415

Note.

^aThe "standard" values used in this paper.

Download table as:  ASCIITypeset image

Figure 18 shows species whose abundance is substantially affected by the changes in activation energies for reactions 17–19 in Table 6. Comparison with the Standard model results (Figures 14 and 15) reveals that the ice abundance of species that are formed via CH₃O is increased by up to six orders of magnitude (CH₃OCH₃ in B68). Importantly, the ice abundances of CH₃OCH₃, CH₃OH, and CH₃O are now above or similar to the detected gas-phase values.

Zoom In Zoom Out Reset image size

The synthesis of methyl formate HCOOCH₃ depends on both HCO and CH₃O radicals. We found it impossible to achieve a sufficiently high ice abundance for HCOOCH₃ with changes in E_A for reactions 17–19. Thus, the calculated relative ice abundance for this molecule remains approximately four orders of magnitude below the observed gas-phase values, which are between 10⁻¹¹ and 10⁻⁹ (Öberg et al. 2010; Bacmann et al. 2012; Cernicharo et al. 2012).

Methyl formate is also inefficiently produced in several other models considering either starless or prestellar cores (Garrod 2013; Vasyunin & Herbst 2013a, 2013b). Chang & Herbst (2014) are able to efficiently synthesize HCOOCH₃ with a Monte Carlo model, but their results also show a likely excess of CH₄ and a shortage of CO₂—discrepancies generally not observed in the above-mentioned other models. These results may indicate that the inefficient synthesis of HCOOCH₃ is not an artifact, caused by the limits of the present study. A higher temperature may be necessary for an efficient synthesis of HCOOCH₃, as suggested by Garrod & Herbst (2006).

The most significant abundance increase for the CH₃O group is observed in the B68 model. First, the lower temperatures in this core mean a low mobility for atomic H in ice. Second, the abundant CO is now more easily hydrogenated, thanks to the lower E_A, and CO soaks up much of the available atomic hydrogen in the mantle. These two factors result in an accumulation of radicals, so that they have a greater chance to react with other multi-atom species. For the present discussion on COMs, the most important radicals are CH₃O and HCO. The effect is especially visible for B68 and CB 17 S.

For species that are formed via HCO, the situation is less clear because formaldehyde H₂CO is consumed by reactions 18 and 19. In cores CB 17 L and CB 27, formaldehyde has its maximum ice abundance decreased below the detected gas-phase values. A similar situation can be observed for formic acid HCOOH, directly related to HCO.

Figure 19 shows that in the compromise model the abundances of COMs are increased in all sublayers (cf. Figures 16 and 17). The main part of their synthesis still occurs in the outer sublayer 1. Species whose synthesis involves CH₃O and CH₃ do not reach a steady abundance until the end of the simulation, mainly because they can be regarded as the end products of the CO hydrogenation sequence.

Zoom In Zoom Out Reset image size

Figures 20 and 21 show that the temperature spike positively affects the production of a number of solid species—HCOOCH₃, CH₂CCO, CH₂CO, NH₂CHO, C₂H₂, CH₄. Notably, the gas-phase abundance is substantially increased for species synthesized on the surface—CH₂CCO, H₂CO, NH₂CN, C₂H₂, CH₃CN, CH₂NH, and CH₄, thanks to reactive desorption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For most species, the changes in ice abundance are within one order of magnitude. However, for methyl formate HCOOCH₃ and formamide NH₂CHO an increase by approximately five and three orders of magnitude is observed, respectively. We therefore explore more deeply the temperature-spike-induced chemistry of these two species. Initiated at 1.6 Myr, the rise of abundances actually occurs over a period of 100 kyr after the temperature spike. Figure 22 shows the abundances of NH₂CHO and HCOOCH₃ in ice sublayers.

Zoom In Zoom Out Reset image size

Methyl formate production occurs mainly in the middle sublayer 2. This is unlike most other COMs in the Standard and Compromise models that are synthesized in the outer sublayer 1 (Figure 19). Sublayer 2 is the only one that is rich in both CO and H₂O. The former is the source species for HCO and CH₃O radicals that combine into HCOOCH₃. Meanwhile, the photodissociation of water is the principal source of atomic hydrogen, necessary for the hydrogenation sequence $\mathrm{CO}\;\longrightarrow \;\mathrm{HCO}\longrightarrow {{\rm H}}_{2}\mathrm{CO}\;\longrightarrow \;{\mathrm{CH}}_{3}{\rm O};{\mathrm{CH}}_{2}\mathrm{OH}$. Thus, CO hydrogenation by the mobile H atoms at the 20 K spike, followed by the combination of HCO and CH₃O ≈ 100 kyr after the spike, is the main mechanism for the formation of HCOOCH₃.

The case of formamide is different—NH₂CHO reaches similar abundance in all three sublayers after the 20 K spike. This is because in none of the sublayers do its main parent species—CO and NH₃—both have a high abundance. While NH₃ is mainly concentrated in sublayer 3, CO is abundant only in sublayers 1 and 2.

The ice abundance increases occur because diffusion, proximity, and activation energy barriers are more effectively overcome at 20 K. The subsequent reactions temporarily reduce ice abundances for the radical species HCO, CH₃O, and CH₂OH. The results shown in Figures 20 and 21 suggest that even such a small and short warm-up period is sufficient to qualitatively reproduce the abundance of methyl formate in ice and that the ice abundances of other organic species do not become disbalanced. The single-point 0D approach used in the model prevents a quantitative comparison between calculation results and observations.

We suggest two likely causes for a (temporal) temperature increase for a portion of ices in a starless or prestellar core. First, the gas parcel can be exposed to the interstellar radiation because of turbulent motions in the core, as suggested by Boland & de Jong (1982) and Martinell et al. (2006). Second, low-velocity transient shocks within the core (Williams & Hartquist 1984) may also result in mild heating of the gas. Both of these mechanisms also offer a hypothesis for the desorption mechanism—either photodesorption or collisional sputtering of the weakly bound, non-polar ice outer layers. These considerations are supported by the turbulent motions observed in starless cores (Redman et al. 2006; Levshakov et al. 2014; Steinacker et al. 2014).