UCLCHEM: A Gas-grain Chemical Code for Clouds, Cores, and C-Shocks^*

At its core, UCLCHEM is a code that sets up and solves the coupled system of ODEs that gives the fractional abundances of all the species in a parcel of gas. The ODEs are created automatically from the network (see Section 2.3). There is one ODE per species, each of which is a sum over the rates of every reaction that involves the species. At each time step, the rates of each reaction are recalculated and the third-party ODE solver DVODE (Brown et al. 1989) integrates to the end of the time step. Each type of reaction requires a different rate calculation, so UCLCHEM is limited by the kinds of reactions that have been coded. The subsections below detail each type of reaction and physical process currently included.

2.1.1. Gas-Phase Reactions

Gas-phase reactions make up the largest part of the chemical network. Two-body reactions, cosmic-ray interactions, and interaction with UV photons are all included in the code. The rate of each reaction is in cm³ s⁻¹. For two-body reactions, the rate is calculated through the Arrhenius equation,

$$\begin{eqnarray}&&{R}_{\mathrm{AB}}=\alpha {\left(\displaystyle \frac{T}{300{\rm{K}}}\right)}^{\beta }\exp \left(\displaystyle \frac{-\gamma }{T}\right),\end{eqnarray} \tag{ 1 }$$

where α, β, and γ are experimentally determined rate constants. For cosmic-ray protons, cosmic-ray induced photons, and UV photons the rates are given by

$$\begin{eqnarray}&&{R}_{\mathrm{CRProton}}={\alpha }_{\mathrm{cr}}\zeta \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{CRPhoton}}={\alpha }_{\mathrm{cr}}{\left(\displaystyle \frac{T}{300{\rm{K}}}\right)}^{\beta }\displaystyle \frac{E}{1-\omega }\zeta \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{UV}}={\alpha }_{\mathrm{rad}}{F}_{\mathrm{UV}}\exp (-{{kA}}_{v}),\end{eqnarray} \tag{ 4 }$$

where ω is the dust grain albedo, ζ is the cosmic ray ionization rate, F_UV is the UV flux, ${\alpha }_{\mathrm{rad}}$ and α_cr scale the overall radiation field and cosmic ray ionization rates for a specific reaction, β takes the same meaning as above, E is the efficiency of cosmic ray ionization events, and k is a factor increasing extinction for UV light. This formulation of the reaction rates is taken from UMIST database notes; more information can be found in McElroy et al. (2013). These rates are then used to set up and solve a series of ordinary differential equations of the form,

$$\begin{eqnarray}&&{\dot{Y}}_{\mathrm{product}}={R}_{\mathrm{AB}}{Y}_{{\rm{A}}}{Y}_{{\rm{B}}}{n}_{{\rm{H}}},\end{eqnarray} \tag{ 5 }$$

where Y is the abundance of reactants; A, B, and $\dot{Y}$ are the rate of change of the product, and n_H is the hydrogen number density. The second abundance dependency is dropped for reactions between single species and photons or cosmic rays. An equivalent negative value is added to the change in the reactants.

In addition to these reactions, the rate of dissociation of CO and H₂ due to UV is reduced by a self-shielding treatment in the code. H₂ formation is treated with the hard-coded reaction rate,

$$\begin{eqnarray}&&{R}_{{\rm{H}}2}={10}^{-17}\sqrt{T}{n}_{{\rm{H}}}{Y}_{{\rm{H}}},\end{eqnarray} \tag{ 6 }$$

where Y_H is the abundance of atomic hydrogen. This rate is taken from de Jong (1977) and performed well in chemical benchmarking tests (Röllig et al. 2007). It should be noted that, despite the self-shielding treatment and inclusion of UV photon reactions from UMIST, UCLCHEM is not a PDR code and should not be used to model regions where the UV field is sufficient to be considered a PDR. To model the chemistry in PDRs, the 3DPDR code described in Bisbas et al. (2012) is available on the UCLCHEM website, and a updated version of the 1D UCL_PDR (Bell et al. 2006) will be available soon.

2.1.2. Freeze Out

All species freeze-out at a rate given by

$$\begin{eqnarray}&&{R}_{{fr}}={\alpha }_{f}\sqrt{\displaystyle \frac{T}{m}}{a}_{g}{fr},\end{eqnarray} \tag{ 7 }$$

where a_g is the grain radius and fr is the proportion of each species that will freeze. Here, α_f is a branching ratio allowing the same species to freeze through different channels into several different species. Therefore, α_f allows the user to determine what proportion of a species will freeze into different products. This is used as a proxy for the relatively uncertain surface chemistry. For example, by freezing a portion of a species as a more hydrogenated species, rather than explicitly including a hydrogenation reaction in the surface network. The values for α_f in each freeze-out reaction should sum to 1 for any given species. Variable fr allows the user to set the freeze-out efficiency; it is left at 1.0 in this work, so that non-thermal desorption alone accounts for molecules remaining in the gas-phase. An additional factor C_ion is included for positive ions, to reproduce the electrostatic attraction to the grains,

$$\begin{eqnarray}&&{C}_{\mathrm{ion}}=1+\displaystyle \frac{16.47\times {10}^{-4}}{{a}_{g}T},\end{eqnarray} \tag{ 8 }$$

with a_g and T being the average grain radius and temperature, respectively, and the constant value taken from Rawlings et al.(1992).

2.1.3. Non-thermal Desorption

Three non-thermal desorption methods from Roberts et al. (2007) have been included in UCLCHEM. Molecules can leave the grain surface due to the energy released in H₂ formation, incident cosmic rays, and by UV photons from the interstellar radiation field, as well as those produced when cosmic rays ionize the gas.

$$\begin{eqnarray}&&{R}_{\mathrm{DesH}2}={R}_{{\rm{H}}2\mathrm{form}}{\epsilon }_{{\rm{H}}2}{M}_{x},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{DesCR}}={F}_{\mathrm{cr}}A{\epsilon }_{\mathrm{cr}}{M}_{x},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{DesUV}}={F}_{\mathrm{uv}}A{\epsilon }_{\mathrm{UV}}{M}_{x},\end{eqnarray} \tag{ 11 }$$

where ${R}_{{\rm{H}}2\mathrm{form}}$ is the rate of H₂ formation, F_cr is the flux of iron nuclei cosmic rays, and F_uv is the flux of ISRF and cosmic-ray-induced UV. The values are efficiencies for each process, reflecting the number of molecules released per event. A is the total grain surface area and M_x is the proportion of the mantle that is species x. The values of these parameters are given in Table 1, and the assumptions used to obtain them are discussed in Roberts et al. (2007).

Table 1.  List of Model Parameter Values

Variable	Value
Pre-shock Temperature	10 K
Initial Density	102 cm−3
Freeze-out Efficiency	1.0
Radiation Field	1.0 Habing
Cosmic-ray Ionization Rate	1.3 × 10−17 s−1
Visual Extinction at Cloud Edge	1.5 mag
Cloud Radius	0.05 pc
Grain Surface Area	2.4 × 10−22 cm2
H2 Desorption ()	0.01
CR Desorption (ϕ)	105
UV Desorption (Y)	0.1

Download table as:  ASCIITypeset image

2.1.4. Thermal Desorption

2.1.5. Grain Surface Reactions

A number of physics modules are included with the code and can be interchanged as required. They are explained briefly below.

2.2.1. Clouds and Collapse Modes

2.2.2. C-shock

A network of all the included species and their reactions must be supplied to the model. For UCLCHEM, this consists of two parts: the gas-phase reactions taken from the UMIST database and a user-defined set of grain reactions. Grain reactions include freeze-out, non-thermal desorption, and any reactions between grain species. The UCLCHEM repository includes the UMIST12 data file (McElroy et al. 2013) and a simple set of grain reactions. Given that the code must be explicitly given the ODEs and told which species to sublimate, and in which proportions for different thermal desorption events, a Python script has been created to produce the input required to run UCLCHEM from a list of species, a UMIST file, and a user-created grain file.

The model runs consider 50 gas parcels equally distributed from the core to the edge of a cloud of size 0.2 pc, with an initial density of 2 × 10² cm⁻³. This cloud is embedded in a larger medium, which adds 1.5 mag to the visual extinction at the cloud edge. This cloud maintains a temperature of 10 K as it collapses, according to the freefall equation

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dt}}={\left(\displaystyle \frac{{n}_{{\rm{H}}}^{4}}{{n}_{0}}\right)}^{\tfrac{1}{3}}{\left[24\pi {{Gm}}_{{\rm{H}}}{n}_{0}\left({\left(\displaystyle \frac{{n}_{{\rm{H}}}}{{n}_{0}}\right)}^{\tfrac{1}{3}}-1\right)\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 12 }$$

where n_H is the hydrogen nuclei number density, n₀ is the initial central density, and m_H is the mass of a hydrogen nucleus. The initial molecular abundances are set to zero, and elements are set to solar abundance values (Asplund et al. 2009). The cloud then collapses until it reaches a density of 2 × 10⁴ cm⁻³, the initial density used for the BE collapse. This takes approximately 3.5 Myr. The chemical processes described in Section 2.1 occur during this collapse. This process is typically referred to as phase 1 in work using UCLCHEM, and is used to create initial abundances and physical conditions for the following phase, which is typically the process of interest. In this case, the abundances at the end of this collapse are used as the initial abundances for model runs that start at 2 × 10⁴ cm⁻³ and collapse to 1 × 10⁸ cm⁻³, either in freefall collapse or according to the BE collapse parameterization.

In order to test the model, it is useful to compare to observations. Tafalla et al. (2004) (T04, hereafter) present observations of two starless cores, L1517B and L1498, that are well-approximated by BE spheres with central densities of approximately 2.2 × 10⁵ cm⁻³ and 0.94 × 10⁵ cm⁻³, respectively. Although the reported characteristics for each core are very similar, L1498 is here used as a comparison to the model. T04 derive temperature and density profiles from the dust emission. They then derive CS, CO, NH₃, and NH₂⁺ abundance profiles from radiative transfer fitting to the flux profile, favoring the simplest profile that recreated the data.

They find CO and CS to be depleted in the denser centers of the cores compared to the outer radii, and that this is adequately represented as a step function. In contrast, NH₃ is enhanced by a factor of a few in the center; a trend they represent as a simple radial power law. N₂H⁺ has a constant abundance for L1517B, but a slight drop at large radii for L1498. Chemical models including the one presented here accurately reflect the behavior of CS and CO. However, many struggle to reproduce the NH₃ and N₂H⁺ behavior (Aikawa et al. 2003). Figure 1 shows results of UCLCHEM with the BE model used to fit these observations at a central density of 0.95 × 10⁵ cm⁻³.

Zoom In Zoom Out Reset image size

From Figure 1, it can be seen that the BE sphere modeling with UCLCHEM yields some success. For CO and CS, T04 find each can be treated as having constant abundance, outside of some cutoff radius, with the abundance in the core being greatly depleted. For simplicity, T04 take a central abundance of zero; although this is clearly not the case for UCLCHEM, there are some caveats. The T04 observations are of C¹⁷O and C¹⁸O, and they report C¹⁸O abundances. For the comparison here, those values have been multiplied by 550, the C¹⁶O:C¹⁸O ratio (Wilson 1999). Reducing the CO values by this amount gives a central abundance of 2 × 10⁻¹⁰, which is sufficiently low that the signal should be too weak to detect. Similarly, the central abundance of CS in the UCLCHEM model is 4 × 10⁻¹¹, and therefore in agreement with zero.

T04 provide an analytic expression for NH₃ which is plotted in Figure 1; the main finding is that NH₃ increases by a factor of a few toward the center of the core. The UCLCHEM model predicts the opposite, increasing with radius. UCLCHEM similarly struggles to reproduce observations of N₂H⁺, decreasing both toward the center and edge of the core, rather than the approximately constant abundance found in T04. It appears from the model outputs that, at middle radii, N₂H⁺ is forming through a reaction between ${{\rm{H}}}_{3}^{+}$ and N₂. The increased density of the middle radii makes this reaction more efficient than it is at the edges, but the even-higher density and visual extinction of the center reduces the amount of ${{\rm{H}}}_{3}^{+}$ available, and further increases the rate of freeze-out.

The freefall model is plotted on the same figure, using thinner dashed lines. The freefall model only takes into account the initial density of each parcel, not the radius, so there is no differentiation with distance from the core center. Further, the freefall model reaches 1 × 10⁵ cm⁻³ much more quickly than the BE model, so the depletion of each species is much lower than expected because freeze-out has less time to occur. In many situations, the freefall model is sufficient for single-point models that describe the densest part of a cloud of gas, but it is clearly inappropriate to model the collapse of prestellar cores.

The shortcomings of the BE model could be due to a number of factors. Spherical symmetry is assumed in both the UCLCHEM model and the radiative transfer modeling of T04, who fit to azimuthally averaged data. Alternatively, the collapse of a marginally super-critical BE sphere may not be exactly appropriate for this object, or the initial conditions for the chemical model may be incorrect. The fact the behavior of both N₂H⁺ and NH₃ are not reproduced by the BE model could indicate a problem with the nitrogen chemistry in the dense gas at core center. However, the model is clearly superior to a simple freefall collapse, and a more thorough investigation into collapse modes and initial conditions could link the observed abundance profiles and the chemical modeling. The strength of UCLCHEM is that it has been specifically designed to have physics models that are easily modified, and thus easily implemented. To illustrate the code's range of physical applications, in the next section, we present a second example with the chemical modeling of C-type shocks.

In this section, we present a second example use of UCLCHEM that focuses on the study of C-type shocks. A primary concern of chemical modeling is the identification of tracer molecules, the observation of which can provide information on a quantity that is not directly observed. We aim to find molecules that trace shocks, particularly those sensitive to the shock properties.

Here, we build on the work of Viti et al. (2011), in which H₂O and NH₃ lines were used to determine the properties of a shock. Some molecules, such as H₂O, trace the full shock because they are not destroyed in the hotter parts of the post shock gas, which then dominates the high-J emission. Others, such as NH₃, may not; they are instead destroyed in the hotter parts of the post-shock gas. This behavior depends on both the pre-shock density and the shock speed, so the determination of which molecules fall into which category in varying conditions is a rich source of information. We provide groupings of molecules under different shock conditions (pre-shock density and shock speed) so comparisons can be made between groups when observing shocked gas.

As explained in previous sections, the model is two-phased. Phase 1 starts from purely atomic gas at low density (n_H = 10² cm⁻³ in this work), which collapses according to the freefall equation, in a manner similar to the core collapse model described above. This sets the starting abundances for all molecules in phase 2, which will be self-consistent with the network, rather than being assumed values for a molecular cloud. The elemental abundances are inputs; for this work, typical elemental solar abundances from Asplund et al. (2009) were used. The cloud collapses to the initial density of phase 2 and is held there until a chosen time. The grid run for this work used a phase 1 of 6 Myr, which provided sufficient time for the required pre-shock densities to be reached.

Table 1 lists the physical parameters that were set for all models run for this work. The phase 1 models used a temperature of 10 K, a radiation field of 1 Habing, and the standard cosmic-ray ionization rate of 1.3 × 10⁻¹⁷ s⁻¹. A radius of 0.05 pc is chosen for the cloud; this sets the Av by allowing the column density between the edge of the cloud and the gas parcel to be calculated assuming a constant density. An additional 1.5 mag are added to the Av, an assumed value for the amount of extinction from the interstellar radiation field to the edge of the cloud. The final density of this phase varies, depending on the density required for phase 2 (See Table 2).

In phase 2, the C-shock occurs. Table 2 gives a full list of the C-shock models that were run. Not all of the listed variables are independent; t_sat and T_max are functions of the pre-shock density (${n}_{{\rm{H}},0}$) and shock velocity (V_s). The T_max, V_s relationship is extracted from Figures 8 and 9 of Draine et al. (1983), which give values of maximum temperature for pre-shock densities 10⁴ and 10⁶ cm⁻³. Because no values are available when the pre-shock density is 10⁵ cm⁻³, the maximum temperature is calculated as the average of the 10⁴ and 10⁶ values. More details of these variables and how they are calculated can be found in Jimenez-serra et al. (2008). The temperature and density profiles of the gas during this C-shock are shown in Figure 2, using model 27 as an example.

Zoom In Zoom Out Reset image size

Due to the fact freezing is efficient in cold, dense clouds, many species are highly abundant in the solid phase of the pre-shock cloud model used in this study. The result is that the fractional abundances of a similarly large number of species increase by orders of magnitude in shocks of any speed, and can be used to determine whether a shock has passed. These are typically molecules that do not have efficient formation mechanisms in cold gas, but are easily formed on the grain. In particular, hydrogenated molecules tend to fall into this category.

H₂O and NH₃ are highly abundant in the ices, but not the gas phase of cold molecular clouds. Of the two, NH₃ is a more useful tracer of shocks due to the difficulties in observing H₂O from the ground. In the models, NH₃ has a fractional abundance of 10⁻¹⁴ in the pre-shock gas, but can increase to 10⁻⁵ in a shock (see Figure 3). H₂S and CH₃OH similarly go from extremely low abundances in the pre-shock gas to X ∼ 10⁻⁵ in a shock. Fractional abundances as a function of shock speed for different pre-shock densities are shown in Figure 3.

Zoom In Zoom Out Reset image size

An important caveat to this is that the model assumes all of the grain surface material is released unchanged into the gas phase. This may not be the case, because the collisions that cause sputtering may have enough energy to destroy molecules as they are released. Suutarinen et al. (2014) discuss this possibility in the context of CH₃OH being released from ice mantles. They find that more than 90% of CH₃OH is destroyed either in the sputtering process itself or in gas-phase reactions after sublimation. The latter is the focus of the next section of this work, but it is not clear how much each process contributes to the destruction. However, the species presented here increase by over six orders of magnitude in gas-phase abundance, due to sputtering, if even 1% survives the process then the enhancement is observable.

An attempt was made to identify species that vary in fractional abundance by orders of magnitude depending on the shock speed. This is in contrast to those shown in Figure 3, which show a steep change from no shock to a low-velocity shock, but no variation thereafter. Such species would allow for the possibility of determining shock properties directly from the abundance of certain molecules. Despite many species increasing in shocks, it was not possible to identify any such species. In the model, even the low-velocity shocks (10 km s⁻¹) sputter the grains such that all shocks can increase gas-phase abundances in this way. Any variability due to gas temperatures reached in the shock is lost when an average over the whole shock is taken, with abundances varying by less than an order of magnitude.

It may not be possible to use the abundances directly to determine the shock velocity. However, UCLCHEM provides a useful input for the radiative transfer: a gas-phase abundance for the molecule as a function of the passage of the shock. Further, as discussed in Section 4.4, although bulk abundances are not sensitive to shock speed, line profiles can be used to provide further constraints.

The focus so far has been on average abundances across the whole shock, the observed abundance if one only resolves the shocked region. However, molecules can exhibit differences through the shock at different velocities, even when the average is unchanged. If these differences are ignored, the models can only be used to infer that a shock has passed. However, if the velocity profile of the emission is taken into account, much more can be learned from the models.

For example, Codella et al. (2010) presented NH₃ and H₂O profiles from the L1157-B1 region. They showed that, when scaled to be equal at their peak emission velocity, NH₃ dropped off much more quickly with increasing velocity. Viti et al. (2011) demonstrated that this is consistent with NH₃ being at lower abundance in the hotter gas, with chemical models showing that a shock with v_s > 40 km s⁻¹ and pre-shock density of n_H = 10⁵ cm⁻³ was sufficient to destroy NH₃ (which, unlike H₂O, does not efficiently reform in the gas phase). As such, a divergence in the velocity profiles of these molecules can indicate such a shock.

Gómez-Ruiz et al. (2016) took this further and showed that the difference in the terminal velocity of NH₃ and H₂O was correlated with the shock velocity inferred from the H₂O terminal velocity for a sample of protostellar outflows. Chemical modeling in the work corroborated this by showing more NH₃ being lost in faster shocks. Holdship et al. (2016) showed H₂S exhibits the same behavior as NH₃, dropping in abundance in high-velocity gas. Given the promise of this method, molecules that exhibit this behavior are explored so this analysis can be easily extended by observing lines of these molecules. It should be noted that the key assumption of the above studies is that the decrease in emission intensity with velocity between two molecules emitting from the same source is due to a similar change in abundance with velocity. Excitation effects will complicate this, but when comparing optically thin emission scaled to match at peak velocity, it should hold as verified with radiative transfer modeling in Viti et al. (2011).

NH₃ and H₂S are unlikely to be the only molecules enhanced in shocks and subsequently destroyed in hot gas. Species other than H₂O will remain enhanced throughout a shock and will thus be equally able to trace the whole post-shock region. H₂O can be difficult to observe with ground-based telescopes, so an alternative species showing a similar behavior could be used as a proxy. The modeling suggests two groups of molecules: shock-tracing (H₂O-like) and shock-destroyed (NH₃-like). Fractional abundances through a shock model for selected species from each categories are shown in Figure 4 to illustrate their behavior. Observing molecules from each group and comparing their velocity profiles and abundances to the predictions of the model can constrain the shock properties.

Zoom In Zoom Out Reset image size

A species does not necessarily fit one category or the other for all shock conditions. Their behavior is the result of temperature and density profiles of the shock, and so depends on the shock speed and initial density. In Table 3, all modeled pairs of shock speed and initial density are shown, with lists of shock-tracing and shock-destroyed species in blue and red, respectively. Species that exhibit more complicated behavior and cannot be easily categorized are simply omitted. Some species, such as H₂O, are almost always found to fall into one category; this allows them to be consistently used for a given purpose, i.e., H₂O can be used to trace the full shock in most cases. However, some species will demonstrate a certain behavior only in very specific conditions, and can be used to limit possible shock properties when this is observed. For example, H₂O does not trace the full shock if the initial density is 10³ cm⁻³ and the shock is low-speed. Such species are displayed in bold on Table 3.

The shock-destroyed and -tracing molecules are both enhanced by the sputtering caused by the shock. This is important because molecules that show low abundances in the pre-shock gas are unlikely to have large contributions from unshocked gas when a shocked region is observed. For example, CO shows perfect shock-tracing behavior, but is also abundant in the pre-shock gas—and its abundance is rather insensitive to shock interaction. As chemistry is the focus for this work, such species are omitted. However, note that high-J transitions of species such as CO are unlikely to be excited by ambient gas, so they can be reliably used to trace the full extent of the shock.