SURVIVAL OF INTERSTELLAR MOLECULES TO PRESTELLAR DENSE CORE COLLAPSE AND EARLY PHASES OF DISK FORMATION

The physical structure of the collapsing core is derived using the RMHD solver of RAMSES, which integrates the equations of ideal magnetohydrodynamics (Fromang et al. 2006) and the equations of radiation-hydrodynamics under the gray flux-limited diffusion (FLD) approximation (Commerçon et al. 2011c). The initial conditions are those used in Commerçon et al. (2012a) and consist of a 1 M_☉ sphere of uniform density ρ₀ = 3.97 × 10⁻¹⁸ g cm⁻³ (corresponding to a molecular hydrogen number density $n_{\rm H_2}\approx 10^6$ cm⁻³) and temperature T₀ = 11 K. The sphere (radius r₀ ≈ 3300 AU) is in solid body rotation and threaded by a uniform magnetic field. The strength of the magnetic field is expressed in terms of the mass-to-flux to critical mass-to-flux ratio μ = (M₀/Φ)/(M₀/Φ)_c. M₀ is equal to M_☉. Φ is the magnetic flux. If (M₀/Φ) is lower than (M₀/Φ)_c (i.e., μ is lower than 1), then the core is supported against its own gravity by the magnetic field and does not collapse. In this study, we use two magnetization levels: a moderate one (μ = 10) and a low one (μ = 200). For μ = 10, we use two different initial configurations: an angle Θ = 0° between magnetic field lines and the rotation axis of the sphere, and another angle Θ = 45°. These values form three different models that we present in this paper: MU10Θ0, MU10Θ45, and MU200Θ0 (see Table 1). We computed the evolution of the collapsing dense core throughout the first collapse and first Larson core lifetime (Larson 1969), i.e., ∼3.8 × 10⁴ yr. The RMHD method used in RAMSES has been well tested, as well as the numerical convergence by Commerçon et al. (2010). Accordingly, the grid is dynamically refined as a function of the local Jeans length so that it is always resolved by at least 10 cells. In addition, it has been shown in previous studies that the gray FLD approximation is well adapted to study the early phases of star formation by comparing to more advanced moment and multifrequency models (Commerçon et al. 2011a; Vaytet et al. 2012).

Nautilus computes the abundance of species in the gas phase and at the surface of the grains as a function of time, taking into account pure gas-phase chemistry (bimolecular reactions, dissociative and radiative recombinations and associations, electron attachments, dissociations and ionizations induced by cosmic-rays and UV photons, etc.), interaction between the species in the gas phase and the grain surfaces (adsorption, thermal desorption, and non-thermal desorption by cosmic-rays), and grain surface reactions (chemical reactions, and dissociation by UV photons and cosmic-ray-induced photons). A more detailed description of this model and the chemical network can be found in Semenov et al. (2010) and Hincelin et al. (2011). The network has been updated according to the latest recommendations from the KIDA⁵ experts until 2011 October. An electronic version of our network is available at http://kida.obs.u-bordeaux1.fr/models.

To compute the initial chemical composition of the sphere, we run Nautilus for dense cloud conditions (10 K, H density⁶ of 2 × 10⁵ cm⁻³, cosmic-ray ionization rate of 1.3 × 10⁻¹⁷ s⁻¹, and a visual extinction of 30) for a time t ∼ 6 × 10⁵ yr, assuming that species are initially in the atomic form except for hydrogen, which is already molecular. We used the elemental abundances from Hincelin et al. (2011), with an oxygen elemental abundance relative to total hydrogen equal to 1.5 × 10⁻⁴ (C/O ratio equal to 1.13). For each element, we define the elemental abundance as the ratio of the number of nuclei both in the gas and on the dust grains to the total number of H nuclei. This excludes the nuclei locked in the refractory part of the grains, which is why we did not use solar elemental abundances, which take into account this refractory part. The chosen time t corresponds to the maximum of agreement between observations of molecular clouds and our simulations (see Figure 3 of Hincelin et al. 2011). Using this chemical composition as the initial condition, we then run the chemistry during the collapse using the method described in the next subsection.

To simplify our simulations, which are already very time consuming, there is no feedback of the chemical composition calculated by Nautilus on the dynamics computed using the RAMSES code, i.e., Nautilus computations are a post-processing of RAMSES simulations (see Section 4.4). To compute the chemistry, Nautilus needs two input data: the temperature and the density as a function of time. For both the physical and chemical calculations, we assume that the object is embedded in a cloud so that direct photoprocesses can be neglected (visual extinction is fixed to 30) and that the cosmic-ray ionization rate is constant (fixed to 1.3 × 10⁻¹⁷ s⁻¹). The first step of the calculation is to obtain the time-dependent physical conditions. To this aim, we introduce 10⁶ "probe" particles, initially uniformly distributed in the sphere, which handle the temperature, density, and velocity history of the fluid through a Lagrangian approach. For any particle, at each time step, the position is registered; the temperature, density, and velocity are extracted according to the local fluid properties; and the position at the next time step is iterated. At the end of this process, we obtain the positions, the velocities, and the couples temperature/density (T(t), n(t)) as a function of time t for the million particles of the model. Figure 1 shows an example of computed trajectory and its corresponding (T(t), n(t)). Starting from the initial chemical composition of the sphere, the chemical evolution is then computed in zero-dimensional for each parcel of material ("probe" particle) using the physical conditions derived from the RMHD simulations (i.e., the 10⁶ couples (T(t), n(t))). Note that the chemistry does not have the time to reach steady state in any of our calculations. This process gives the full three-dimensional chemical structure of the collapsing dense core as a function of time. Central processing unit (CPU) time to compute the chemical evolution of 10⁶ particles (i.e., for one model) is roughly 10⁴ hr (using JADE cluster from CINES).

Zoom In Zoom Out Reset image size

Within this structure, we identify several components representing (1) the core itself (the first Larson core), (2) the outflow, (3) the rotating disk, (4) the magnetically supported pseudodisk, and (5) the envelope, based on criteria from Joos et al. (2012). More details on the extraction method will be given in a future paper (U. Hincelin et al., in preparation). The rotating disk component that we consider here is made by particles that satisfy the four following criteria:

• 1.  
  an azimuthal velocity v_ϕ more than two times larger than their radial velocity v_r,
• 2.  
  an azimuthal velocity more than two times larger than their vertical velocity v_z,
• 3.  
  a rotational support (density of rotational kinetic energy $E_k=0.5\rho v_\phi ^2$, where ρ is the mass density) at least two times larger than the thermal support (density of internal energy E_i = 1/(γ − 1) × (ρk_BT)/(μ_mm_H), where γ, k_B, μ_m, and m_H are, respectively, the adiabatic index, the Boltzmann constant, the mean molecular weight, and the mass of hydrogen nucleus),
• 4.  
  and a density above 10⁹ cm⁻³.

Figures 2 and 3 present selected particles that belong to the rotating disk for the three models presented in Table 1. The other components (central core, outflow, pseudodisk, and envelope) will be discussed in a future paper (U. Hincelin, et al. in preparation). The plotted information (in color) is the temperature (Figure 2) and the density (Figure 3) of the particles. The disks are relatively cold, around 10 K, except at low radius, where the temperature can reach about 100–200 K. Density spreads from 10⁹ to 10^{11 − 12} cm⁻³ from the outer edge of disks to their center. The temperature profile has roughly a spherical symmetry (except for model MU200Θ0 with the presence of four "hot spots" around the central core), which is due to the central core (see also Figure 2 in Commerçon et al. 2012a), while the density profile has a more complex shape, for example, with the presence of spiral arms for μ = 10 (MU10Θ0 and MU10Θ45 models). In these simulations, the disk of model MU10Θ0 is about 200 AU in diameter, whereas the one of model MU10Θ45 is roughly two times larger (∼400 AU). Due to the low magnetization level, the disk of model MU200Θ0 shows a more complex shape with fragmentation. Four "hot spots" are identified around a central disk. Each spot is characterized by an increase of both temperature and density. The central disk is about 80 AU in diameter. Table 2 shows the range of sizes, temperatures, and densities found in the disks at the end of simulations for the three models. One notes that the disk of model MU200Θ0 is fragmented due to the low magnetization level, and the disk of model MU10Θ45 is twisted because of the inclination between rotational axes and initial magnetic field lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.1. Chemical Evolution of the Matter from the Cloud to the Disk

Figure 4 shows the sum of the abundances computed in the gas phase and at the surface of the grains (ices) compared to the total hydrogen (proton) density for a selection of species, at two different steps of the simulations: (1) the initial condition, i.e., the composition of the parent molecular cloud, before the beginning of the collapse (namely, the sphere in Sections 2.2 and 2.3), and (2) the final condition, i.e., the composition in the disk. A change in this total abundance compared to the chemical content of the cloud implies that some processes of formation or destruction have occurred during the collapse in addition to just adsorption and desorption processes. As the modeled object is embedded in a dense envelope, the effects of ionization and dissociation by the external UV photons are low compared to bimolecular chemical reactions. The main result of our simulations is that the chemical content is unchanged during the collapse: there is no significant evolution of the global (ice plus gas) abundances of the species. For example, the abundances of water (H₂O), methane (CH₄), carbon monoxide (CO), and ammonia (NH₃) are the same in the parent cloud and the disk.

Zoom In Zoom Out Reset image size

There are a few exceptions: HNC, CO₂, and HNO. HNC is very sensitive to the evolution of temperature and density during the collapse: in the case of model MU200Θ0, its mean total abundance has decreased from about 5 × 10⁻⁷ (cloud composition) to 4 × 10⁻⁹ (disk composition), i.e., two orders of magnitude. The matter close to the fragments is hot so that HNC desorbs from grains and is efficiently destroyed in the gas phase by $\rm H_2COH^+ + HNC \longrightarrow HCNH^+ + H_2CO$⁷. Then, the HCN/HNC ratio (taking into account gas and ices) can exceed 100. One notes in Figure 4 that even during the cloud stage, the ratio HCN/HNC is not equal to 1, but roughly equal to 10. While the gas-phase ratio is equal to 1, HCN and HNC abundances are higher on grain surfaces than in the gas phase, so that the grain surface ratio leads to a total ratio (gas and ices) of 10. On the grain surface, HCN and HNC are photodissociated into CN and H. Then, CN and H react on the surface, but form preferentially HCN and not HNC. HNC can be formed on the surface using different pathways (C + NH, C + NH₂, or CH + NH), but these reactions are less efficient than $\rm CN + H \longrightarrow HCN$ so that the ratio HCN/HNC on grain surfaces reaches a value of 10. One notes that as a consequence, desorption of ices could enhance the gas-phase ratio HCN/HNC to a value of 10. On the contrary but to a lesser extent, the mean global abundance of CO₂ has increased from the cloud to the disk by a factor of two in model MU10Θ45. High temperature promotes CO₂ formation on the grain surface, particularly through $\rm OH + CO \longrightarrow CO_2 + H$ (Ruffle & Herbst 2001). Desorption of HNO ices is followed by gas-phase reactions mainly involving CH₃, O, CH₂, and CN that produce NO and, respectively, CH₄, OH, CH₃, and HCN.

3.2.2. Ice Composition in the Disk of Model MU10Θ0

To illustrate the chemical composition of the ices in the disks, we present in this section the distribution of the ice composition for model MU10Θ0 and for a selection of molecules observed in comets. Figure 5 shows the maximum (in blue) and minimum (in black) relative abundances compared to H₂O on the grain surfaces for these species together with the abundances observed in comets. Abundances below 10⁻³% are not shown. In some cases, the abundances are constant, like for NH₃, so that the minimum overlaps the maximum and only black is seen. On the contrary, if the species abundance shows a gradient in the disk and the smallest abundance is below 10⁻³, only the maximum value, in blue, appears in Figure 5. As an example, the abundance of methanol ice (CH₃OH) compared to water ice is about 14% all over the radius of the disk (in the figure, the black band (minimum) overlays the blue one (maximum)). The abundance of carbon monoxide ice (CO) is high from the outer edge of the disk to around 30 AU from the central first Larson core (the blue band reaches 84%). Within this radius, the CO desorbs from the grain and the ice abundance falls down to 10⁻¹¹% (this value, being under 0.001%, is not visible in Figure 5). The abundance of formic acid ice (HCOOH) is low (2 × 10⁻⁴%) all over the radius of the disk, which is below our 0.001% threshold, hence not visible in Figure 5.

Zoom In Zoom Out Reset image size

Species listed in Figure 5 can be divided into three groups (Table 3). The first group contains molecules with nearly uniform abundances as a function of the radius (such as water ice because of its relatively high evaporation temperature). In the second group, molecules are more abundant on ices in the outer (>25 AU), cold (∼25 K) disk (such as CO ice), while the reverse is true for molecules in the third group (such as C₂H₂). The reason for the abundance decrease in the inner parts for the second group is evaporation due to higher temperatures. In that case, the gas-phase abundance essentially reflects the initial content of the ices, with mainly one exception, HNC, which is quickly destroyed by gas-phase reactions (see Section 3.2.1) when the temperature is larger than ∼50 K (within a radius of ∼12 AU from the center).

Table 3. Groups of Molecules on Grain Surfaces According to Their Behavior in the Disk for Model MU10Θ0

Homogeneous	H2O, C2H6, CH3OH, HCOOH,
	NH2CHO, NH3, CH3CN
Outer disk	CH4, H2CO, HCN, HNCO, H2CS, HNC,
	H2S, CO, CH3CHO, S2, CO2
Inner disk	C2H2, HCOOCH3, HC3N, OCS, SO2

Download table as:  ASCIITypeset image

The projected H₂O_(ice), CO_(ice), and C₂H_2(ice) abundances (relative to total hydrogen) computed in the disk for model MU10Θ0 are shown in Figure 6. These three molecules illustrate the different behaviors presented in Table 3. The temperature is rather low, below 20 K in the majority of the disk (from the outer edge to ∼40 AU), whereas the total H density presents variations between 10⁹ and 8 × 10¹⁰ cm⁻³ from the outer edge to ∼15 AU. The inner regions (within a radius r ≈ 20 AU) present the warmest part of the disk with temperatures up to about 120 K. The abundance of water (relative to total hydrogen) trapped on interstellar grains is constant across the disk at a value of 6.5 × 10⁻⁵. The decrease of the CO ice abundance is a direct consequence of thermal evaporation in this region. When the temperature increases (as a function of the radius), electronic recombination reactions in the gas phase enhance the abundance of C₂H₂. The gas-phase C₂H₂ is then adsorbed on grains. Desorption is not efficient for this species because its adsorption energy is high enough.

Zoom In Zoom Out Reset image size

Results for model MU10Θ45 and MU200Θ0 are qualitatively similar to those of MU10Θ0, with some different distributions of molecules according to the local temperature and density conditions. Due to fragmentation and the presence of four hot spots in model MU200Θ0, ices are sublimated in these hotter regions, and molecules are eventually destroyed in the gas phase. The higher temperature of the matter next to these regions leads to a more active surface chemistry, which enhances the formation of complex organic molecules such as the methyl formate HCOOCH₃. The disk of model MU10Θ45 is on average a little bit warmer, which leads to a smaller abundance of ices but a higher abundance in the gas phase so that the total abundance is globally not changed by a significant amount (see Figure 4).

Our model stops just before the second collapse and the formation of the second Larson core, which is early in the star formation process (see Section 4.4). However, recent hydrodynamic simulations (Machida et al. 2010; Machida & Matsumoto 2011; Bate 2010, 2011) have shown that the central part (radius ≲ 1 AU) of the first Larson core will form the second Larson core, that is to say, the protostar, whereas the outer part of the first Larson core will evolve in a disk. During the protostellar phase, a small fraction of the first core material should have been transported to the outer region of the disk, and not been accreted to the protostar (see discussion of FA12).

Considering the partial interstellar origin of the molecules observed in comets (Bockelée-Morvan et al. 2000; Mumma & Charnley 2011), it is then interesting to compare the comet chemical composition and the results of our models. Due to the fact that our disk is very young, the computed ice can be considered as "pre-cometary," and so this comparison must stay qualitative.

The ice abundances of C₂H₆, CH₃OH, NH₂CHO, and CH₃CN found in our disk for model MU10Θ0 are within a factor of two of those found in comets (see Figure 5). It is also the case for HNCO, H₂S, and H₂CS in the outer part of the disk where the temperature is low enough to maintain these species on grain surface. These species represent about one-third of the observed species in comets. On the other hand, CO₂, C₂H₂, HCOOH, HCOOCH₃, CH₃CHO, HC₃N, OCS, SO₂, and S₂ are strongly underproduced in our model, in the entire disk, compared to comets. Carbon dioxide (CO₂) and carbon sulfide (OCS) are known to be abundant on interstellar ices (Palumbo et al. 1997; Gibb et al. 2000) but not much produced by astrochemical models at low temperature because carbon, oxygen, and sulfur only diffuse slowly on surfaces (Ruffle & Herbst 2001). For example, the surface reaction $\rm HCO + O \longrightarrow CO_2 + H$ is the main pathway in our model for CO₂ production, but is in competition with $\rm HCO + H \longrightarrow H_2CO$, which is faster. In our results, CH₄ and NH₃ are reservoirs of carbon and nitrogen, whereas they are minor species in comets. The methane is probably overproduced in our model because we underestimate the formation of CO₂: the available carbon on grain surfaces is quickly hydrogenated so the carbon is rather locked into methane. For nitrogen, the reservoir is not known in comets (Bockelée-Morvan et al. 2004) and we probably lack some mechanisms to form other N-bearing species (Daranlot et al. 2012). Finally, HCOOH, CH₃CHO, and HCOOCH₃ are known to trace warmer regions where radicals can move at the surface of the grains, and are typically observed in warmer star-forming regions (Herbst & van Dishoeck 2009).

Comets seem to be made of molecules formed in different physical conditions. Table 4 illustrates this idea, proposing a formation place of molecules where the required physical conditions are present. The HNC molecule, for instance, is a strong signature of cold interstellar-like chemistry. Other molecules, on the contrary, such as HCOOCH₃, require temperatures of about 40–60 K to be formed on surfaces (Garrod & Herbst 2006), and as a consequence are not produced by our model in the disk in sufficient amount. This mixed composition may be a strong indication of radial mixing in the protoplanetary disk. This process is not included in our simulations. In fact, we find that the composition of ices is almost constant in the outer regions but that there is a gradient of abundances within 40 AU, a size consistent with the assumed region of formation of comets in our solar system (see Mumma & Charnley 2011 and references therein).

In order to assess the sensitivity of our results to the composition of the initial dense core, we computed the chemical composition for three different particles of model MU10Θ0 using several initial conditions. The three particles reach at the end of the simulation, respectively, the cold region (∼10 K), the warm region (∼30 K), and the hot region (∼100 K) of the disk.

4.2.1. A Denser Molecular Cloud

4.2.2. An Older Molecular Cloud

4.2.3. Observed Abundances in Dense Clouds as Initial Conditions

FA12 present three-dimensional modeling results of a collapsing dense core until the formation of the first hydrostatic core, that is to say, the first Larson core. They use a similar interface as ours: they use fluid particles that trace the physical conditions of the matter during the collapse, and compute the chemical evolution on these particles as a post-process. The main differences from our model are that FA12, on a physical point of view, do not compute the effect of the magnetic field on the collapse process, and on a chemical point of view, use a high-temperature gas-phase network (for a range of 100–800 K). Since we focus on the disk composition in the present paper, where the temperature does not exceed 300 K, we do not expect our results to be affected by this limitation (see Section 4.4).

Figure 7 shows the density and temperature of particles in the equatorial plane (area where |z| < 10 AU) as a function of the distance from the central core, for model MU10Θ0. At a given radius, the range of values comes from the asymmetry of the collapsing core: in particular, because of the presence of spiral arms. The figure can be compared to panel (b) of Figure 9 of FA12 (although the integration time is not exactly the same): the temperature and density profiles are qualitatively similar, except within 1–2 AU, where the density and temperature are higher in the case of FA12. In our model, an additional support against gravitational collapse exists: magnetic support, which limits contraction. As a consequence, the inner computed temperature and density are lower in our case. One of the main conclusions of FA12 is that the total gas and ice abundances of many species remain unaltered during the collapsing process until the temperature reaches about 500 K, that is to say, at a spherical radius r ≲ 3 AU. The gas-phase abundances are determined via desorption processes, depending on the adsorption energy of species on the surface of the grains. Our results show a similar conclusion, at least for the disk component: except at low radius, the chemical composition is globally unaltered during the collapse with gas-phase abundances mainly determined by desorption. In the warmer parts of the disk, that is to say, at lower radius, we also observe a formation of carbon chains and large organic molecules such as C₂H₂, HC₃N, and HCOOCH₃ (see Table 3).

Zoom In Zoom Out Reset image size

The simulations from the RMHD physical model RAMSES stop just before the start of the second collapse. The formation of this second core increases the temperature in the disk at larger radii as suggested by the one-dimensional modeling of Vaytet et al. (2013). Compared to our simulations, this may mainly change the radii at which molecules evaporate. The importance of this second phase will depend on the temperature gradient obtained in the disk. This issue will be investigated in further studies.

The physical model RAMSES integrates the equations of ideal magnetohydrodynamics, assuming the medium to be a perfect conductor. The matter is then frozen on the magnetic lines. If we consider a medium that is not a perfect conductor, that is to say, in the case of nonideal (or resistive) magnetohydrodynamics, then the neutral matter will be able to cross magnetic lines. This induces a friction between ions and neutrals, called ambipolar diffusion (Mestel & Spitzer 1956; Mouschovias 1991). Ambipolar diffusion generates a net velocity difference between ions and neutral that can be sufficient to overcome activation barriers, thereby generating efficient reaction pathways. Future work could be done to couple nonideal magnetohydrodynamics with chemistry in order to take this phenomenon into account. Recent works have been done to integrate ambipolar diffusion in RAMSES (Masson et al. 2012).

As mentioned in Section 2.3, the chemical model Nautilus is used as a post-processing of the RAMSES computation. There is no chemical feedback on the physical structure of the collapsing dense core. However, a radiative cooling of the matter due to chemical species such as CO, CH₄, H₂O, C, and O can occur. The cooling rate of a given chemical species depends in particular on the temperature and the density of the medium and the abundance of the species (see, for example, Neufeld & Kaufman 1993; Dedes et al. 2011). However, at high density such as in the case of our modeling, the cooling is mostly dominated by the radiative emission of the dust.

In our modeling, each particle is independent of the others from a chemical point of view. This is valid only if mixing effects can be neglected, both molecular diffusion and turbulent mixing.

For matter collapsing onto the core and the disk, our simulations do not show any turbulent motion, and in general turbulent mixing is not expected to be efficient during this phase. The influence of molecular diffusion can be qualitatively estimated by comparing the characteristic diffusion timescale τ_diff for two nearby tracer particles, with the time τ_disk for these particles to enter the disk. τ_diff is given by

$$\begin{equation} \tau _\mathrm{diff} \simeq \frac{L^2}{c_s\lambda } \simeq L^2 n_{\rm H_2} \sigma \sqrt{\frac{\mu _m m_{\rm H}}{k_B T}}, \end{equation} \tag{ 1 }$$

where L is the distance between the two particles, c_s is the sound speed (typical velocity of molecules inside the Lagrangian particle), $n_{\rm H_2}$ is the molecular hydrogen number density, λ is the mean free path of gas-phase molecules, and σ is the collision cross section. τ_disk is roughly equal to the free fall time τ_ff:

$$\begin{equation} \tau _\mathrm{disk} \simeq \tau _\mathrm{ff}=\sqrt{\frac{3\pi }{32G\rho _0}}, \end{equation} \tag{ 2 }$$

where G is the gravitational constant. τ_ff is close to the final time given in Table 2 (τ_ff ≃ 3 × 10⁴ yr). In the collapsing region, L > 1 AU, T = 10 K, and $n_{\rm H_2}=10^6$ cm⁻³, so that τ_diff > 2 × 10⁵ yr, which is higher than τ_ff. Therefore, particles collapse faster than they can mix with each other.

When the particles reach the disk, the situation changes. Although our numerical resolution is insufficient to allow for turbulence to develop in the simulation, turbulence is expected to exist in these objects. In such a case, particles can no longer be considered as isolated; molecules are expected to steadily move from one particle to another. The detailed coupling of chemistry with turbulence being beyond the scope of the present paper, we considered only average disk abundances, as displayed in Figure 4. Note, however, that turbulent mixing is expected to homogenize the disk composition, so that the amplitude of our range of disk abundances may be overestimated.

Finally, our chemical network contains reactions whose rate coefficients are evaluated from 10 to 300 K. Above 300 K, endothermal reactions or reactions that have a large activation barrier can become efficient. Since the central part of the collapsing dense core is warmer than 300 K, a high-temperature network should be used for this region. However, the component discussed in this paper, i.e., the disk, does not experience such high temperature. Including the high-temperature gas-phase network from Harada et al. (2010) will be a future improvement of our modeling.