UCLCHEM: A Gas-Grain Chemical Code

TL;DR

UCLCHEM is an open-source, time-dependent gas-grain chemical code that models molecular abundances under various physical conditions, aiding in understanding protostellar collapse and shock phenomena.

Contribution

The paper introduces UCLCHEM, a new publicly available chemical modeling tool capable of simulating complex gas-grain interactions in diverse astrophysical environments.

Findings

Reproduces CO, CS, NH₃ behavior in collapsing cores

Identifies molecular tracers for C-shocks

Classifies molecules based on shock response

Abstract

We present a publicly available, open source version of the time dependent gas-grain chemical code UCLCHEM. UCLCHEM propagates the abundances of chemical species through a large network of chemical reactions in a variety of physical conditions. The model is described in detail along with its applications. As an example of possible uses, UCLCHEM is used to explore the effect of protostellar collapse on commonly observed molecules and to study the behaviour of molecules in C-type shocks. We find the collapse of a simple Bonnor-Ebert sphere successfully reproduces most of the behaviour of CO,CS and NH$_3$ from cores observed by \citet{Tafalla2004} but cannot predict the behaviour of N$_2$H$^+$. In the C-shock application, we find that molecules can be categorized so that they can be useful observational tracers of shocks and their physical properties. Whilst many molecules are enhanced in shocked gas, we identify two groups of molecules in particular. A small number of molecules are enhanced by the sputtering of the ices as the shock propagates and then remain high in abundance throughout the shock. A second, larger set are also enhanced by sputtering but then are destroyed as the gas temperature rises. Through these applications the general applicability of UCLCHEM is demonstrated.

Tables

Table 1. 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 $\times {10}^{-17}{s}^{-1}$
Visual Extinction at Cloud Edge	1.5 mag
Cloud Radius	0.05 pc
Grain Surface Area	2.4 $\times {10}^{-22}$ cm2
$H_2$ Desorption ($ϵ$)	0.01
CR Desorption ($\varphi$)	105
UV Desorption (Y)	0.1

Table 2. Table 2: List of Shock Variable Values

Model	log(nH,0)	Vs	tsat	Tmax	Z${}_{}{}^{\ast }{}_{diss}{}^{}$
		/ km s-1	/ yr	/ K	/ cm
1	103	10	975.7	323	3.70$\times {10}^{17}$
2	104	10	97.6	323	3.70$\times {10}^{16}$
3	105	10	9.8	301	3.70$\times {10}^{15}$
4	106	10	1.0	279	3.70$\times {10}^{14}$
5	103	15	746.5	584	5.55$\times {10}^{17}$
6	104	15	74.7	584	5.55$\times {10}^{16}$
7	105	15	7.5	555	5.55$\times {10}^{15}$
8	106	15	0.8	525	5.55$\times {10}^{14}$
9	103	20	586.13	869	7.41$\times {10}^{17}$
10	104	20	58.61	869	7.41$\times {10}^{16}$
11	105	20	5.86	893	7.41$\times {10}^{15}$
12	106	20	0.59	916	7.41$\times {10}^{14}$
13	103	25	483.09	1178	9.26$\times {10}^{17}$
14	104	25	48.31	1178	9.26$\times {10}^{16}$
15	105	25	4.83	1316	9.26$\times {10}^{15}$
16	106	25	0.48	1454	9.26$\times {10}^{14}$
17	103	30	425.83	1510	1.11$\times {10}^{18}$
18	104	30	42.58	1510	1.11$\times {10}^{17}$
19	105	30	4.26	1824	1.11$\times {10}^{16}$
20	106	30	0.43	2137	1.11$\times {10}^{15}$
21	103	35	402.8	1866	1.30$\times {10}^{18}$
22	104	35	40.28	1866	1.30$\times {10}^{17}$
23	105	35	4.03	2416	1.30$\times {10}^{16}$
24	106	35	0.4	2966	1.30$\times {10}^{15}$
25	103	40	402.47	2245	1.48$\times {10}^{18}$
26	104	40	40.25	2245	1.48$\times {10}^{17}$
27	105	40	4.02	3093	1.48$\times {10}^{16}$
28	106	40	0.4	3941	1.48$\times {10}^{15}$
29	103	45	413.29	2648	1.67$\times {10}^{18}$
30	104	45	41.33	2648	1.67$\times {10}^{17}$
31	105	45	4.13	3855	1.67$\times {10}^{16}$
32	103	60	397.28	3999	2.22$\times {10}^{18}$
33	104	60	39.73	3999	2.22$\times {10}^{17}$
34	103	65	337.32	4497	2.41$\times {10}^{18}$
35	104	65	33.73	4497	2.41$\times {10}^{17}$

Equations

$$R_{AB}=\alpha {\left(\frac{T}{300K}\right)}^{\beta }\exp \left(\frac{-\gamma }{T}\right),$$

$$R_{CRProton}={\alpha }_{cr}\zeta$$

$$R_{CRPhoton}={\alpha }_{cr}{\left(\frac{T}{300K}\right)}^{\beta }\frac{E}{1-\omega }\zeta$$

$$R_{UV}={\alpha }_{rad}{F}_{UV}\exp (-k{A}_v)$$

$${\dot{Y}}_{product}=R_{AB}{Y}_A{Y}_B{n}_H,$$

$$R_{H2}=10^{-17}\sqrt{T}{n}_H{Y}_H,$$

$$R_{fr}={\alpha }_f\sqrt{\frac{T}{m}}{a}_{g}fr,$$

$$C_{ion}=1+\frac{16.47\times 10^{-4}}{a_{g}T},$$

$$R_{DesH2}=R_{H2form}{ϵ}_{H2}{M}_x,$$

$$R_{DesCR}=F_{cr}A{ϵ}_{cr}{M}_x,$$

$$R_{DesUV}=F_{uv}A{ϵ}_{UV}{M}_x$$

$$\frac{dn}{dt}={\left(\frac{n_H^4}{n_0}\right)}^{\frac{1}{3}}{\left[24\pi G{m}_H{n}_0\left({\left(\frac{n_H}{n_0}\right)}^{\frac{1}{3}}-1\right)\right]}^{\frac{1}{2}}$$

Full text

Astrochemistry,ISM: abundances. ISM: molecules

\savesymbol

tablenum \restoresymbolSIXtablenum

1 Introduction

Chemistry is ubiquitous in astrophysical environments. Molecular clouds, the cold cores in which stars form, and the warm gas surrounding protostars all exhibit chemistry of varying degrees of complexity and with different dominant chemical pathways. Understanding this chemistry is vital to the study of our own origins as well as understanding the physical structure and processes involved in star formation.

Beyond this, chemistry is a useful tool for understanding the physical conditions of the region being studied. This requires well constrained chemical networks and accurate physical models so that uncertainties in the predictions of the model are much smaller than the uncertainties in the measured abundances of molecules. With the current state of the art models, networks are generally capable of putting broad constraints on the physical conditions such as maximum temperatures or minimum densities and this can be of use in poorly understood regions.

Chemical modelling is typically performed by the use of rate equations. The rates of a network of chemical reactions are calculated and used to determine the rate of change of a set of chemical species. This coupled set of ordinary differential equations (ODEs) is integrated to give the abundance of each species at any given time. These models typically centre around gas-phase reactions provided by databases such as KIDA (Wakelam et al., 2012) and UMIST (McElroy et al., 2013). Each of these databases provides a chemical code, respectively Nahoon and RATE13, which solve these networks for simple gas conditions and include other processes such as H2 formation on the dust grains and the self-shielding of CO and H2 from UV radiation.

This, of course, does not account for all astrochemical processes and many groups use more or less complicated chemical models for different purposes. In radiation-hydrodynamic models, it is common to reduce the network to gas-phased H,C and O based species to reduce chemical integration time whilst reproducing the abundances of major coolants such as CO given by more detailed chemical models. On the other hand, the modelling of dense prestellar cores or the formation of complex organic molecules requires much larger chemical networks and the addition of processes involving dust grains.

For example, Astrochem (Maret & Bergin, 2015) includes freeze out of species onto dust grains and the non-thermal desorption of species from those grains due to UV radiation and cosmic rays making it suitable for modelling a wider range of regions than a simple gas-phase model. By further including thermal desorption, star-forming regions with high gas/dust temperature conditions can be modelled. As the temperature of the core rises, the material on the grains sublimates and proper treatment of this sublimation is required. Astrochem and UCLCHEM implement these processes, with desorption from the grains occurring according to the binding energy of each species onto the grain.

However, Collings et al. (2004) demonstrated through temperature programmed desorption experiments that multiple desorption events can occur for each species frozen on the dust grains at different temperatures dependent both on the species and on the rate of change of the temperature. UCLCHEM additionally includes these desorption events for each species on the grain surface (See Section 2.1.4). The modelling of many protostellar environments benefit from the addition of grain chemistry and proper temperature dependent sublimation is required, for instance, for massive hot cores.

Another environment where proper grain surface as well as gas-phase chemical treatment has proven useful is in shocked gas. In these systems, sublimation from the grain surface is dominated by sputtering, the collisional removal of surface material. The physical modelling of these shocks is an active area of research (eg. Van Loo et al., 2009; Anderl et al., 2013) and codes for this modelling are publicly available (Flower & Pineau des Forêts, 2015). These models have been successful in reproducing observations of shocked regions. For example, L1157-mm is a protostar driving a bipolar outflow (Gueth et al., 1997) and the bowshocks associated with that outflow are well studied. For one of these bowshocks, L1157-B1, Gusdorf et al. (2008) reproduced the SiO emission using an MHD hydrodynamical model of a continuous or C-type shock coupled with a reduced chemical network of 100 species.

As noted, reduced chemical networks often produce similar abundances to much more complex networks when simple species such as H2O and CO are considered and this is sufficient for many applications (eg. Schmalzl et al., 2014). Similarly, a simplified physical model may reproduce the shock features that are useful for chemistry and so a full chemical model can be run without fully solving the MHD equations. Jimenez-serra et al. (2008) produced a parameterized C-shock that produced a shock structure in good agreement with more complicated MHD models. This was combined with the chemical code UCLCHEM to study the chemistry in shocked environments (Viti et al., 2011). Using this parameterization with a large chemical network has been met with some success in the bowshock L1157-B1 (Viti et al., 2011; Holdship et al., 2016; Lefloch et al., 2016) despite being computationally inexpensive.

The simplicity of these models means they could be used to link observed molecular lines with the underlying chemistry to get a sense of the physical properties of the shock. Conversely, emission from a source with well constrained physical properties could be used to improve uncertain parts of the chemical network. The reaction rates of species on dust grains being a prime example. For example, Holdship et al. (2016) used observations of the L1157-B1 bowshock to constrain sulfur chemistry on ice grains by identifying likely sulfur carriers. Ideally, chemical models could be ran over large parameter spaces to find the most likely values for uncertain parts of the network when a region has large amounts of observational data available.

To this purpose UCLCHEM has been developed over many years with numerous papers discussing major updates and applications (see Viti et al., 2004; Roberts et al., 2007).There have also been benchmarking efforts in the past (Viti et al., 2001) comparing the code to other similar models. We present here the most up to date 222https://uclchem.github.io version, the source code of which is now publicly available for the first time under the MIT license. This paper describes the code in its most recent state and provides example applications of the model. In Section 2 the code is discussed, in Sections 3 and 4 two examples of the applications of UCLCHEM are given. The paper is summarized in Section 5.

2 UClCHEM

UCLCHEM is a time dependent gas-grain chemical model written in modern Fortran. It is primarily a chemical code, focusing on grain chemistry as well as gas-phase reactions. The chemistry includes freeze out, thermal and non-thermal desorption along with a gas-phase and user-provided grain surface reaction network. In addition, hard-coded treatments for H2 formation on the grains and the self-shielding of CO and H2 are contained in the chemistry module. The chemical solver calculates the rates of all the above processes to follow the fractional abundances of molecules for parcels of gas.

UCLCHEM makes use of modern Fortran’s modules to separate the chemistry and physics. Interchangeable physics modules control the number and physical conditions of gas parcels used in the chemistry to allow different scenarios to be modelled. In particular, the code is packaged with modules for molecular clouds, hot cores, C-shocks and the post-processing of hydrodynamic simulations. In this section, the chemical model is explained along with an overview of the physical models available and the python code produced to create networks automatically.

2.1 Chemical Model

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 Sec. 2.3). There is one ODE per species and each 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 and so UCLCHEM is limited by the kinds of reactions that have been coded for. 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^{3}s^{-1}$. For two body reactions, the rate is calculated through the Arrhenius equation,

[TABLE]

where $\alpha$, $\beta$ and $\gamma$ are experimentally determined rate constants. For cosmic ray protons, cosmic ray induced photons and UV photons the rates are given by,

[TABLE]

where $\omega$ is the dust grain albedo, $\zeta$ is the cosmic ray ionisation rate and FUV is the UV flux. $\alpha_{rad}$ and $\alpha_{cr}$ scale the overall radiation field and cosmic ray ionisation rates for a specific reaction, $\beta$ takes the same meaning as above, $E$ is the efficiency of cosmic ray ionisation events and $k$ is a factor increasing extinction for UV light. This formulation of the reaction rates is taken from the UMIST database notes and 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,

[TABLE]

where $Y$ is the abundance of reactants A and B and $\dot{Y}$ is the rate of change of the product. $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 H2 due to UV is reduced by a self-shielding treatment in the code. H2 formation is treated with the hard-coded reaction rate,

[TABLE]

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,

[TABLE]

where $a_{g}$ is the grain radius and $fr$ is the proportion of each species that will freeze. $\alpha_{f}$ is a branching ratio allowing the same species to freeze through different channels into several different species. Therefore, $\alpha_{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 $\alpha_{f}$ in each freeze out reaction should sum to 1 for any given species. $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,

[TABLE]

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

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

[TABLE]

Where $R_{H2form}$ is the rate of H2 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 $\epsilon$ 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).

2.1.4 Thermal Desorption

As the temperature of the dust increases, species start to desorb from the ice mantles. Laboratory experiments (eg. Ayotte et al., 2001; Burke & Brown, 2010) show that a species does not desorb in a single peak but rather in multiple desorption events. These events correspond to the removal of the pure species from the mantle surface, a strong peak at the temperature corresponding to the binding energy, and the “volcanic” desorption and co-desorption of molecules that have diffused into the water ice. Collings et al. (2004) showed that many species of astrochemical relevance could be classified as CO-like, H2O-like or intermediate, depending on what proportion of the species is removed from the mantle in each desorption event.

To the authors knowledge, UCLCHEM is the only publicly available chemical model that implements a treatment for this. Molecules are classified according to their similarity to H2O or CO, controlling their desorption behaviour. Further, the user can now set the proportion of each molecule that enters the gas phase in each event, allowing the thermal desorption treatment of each molecule to improve with the laboratory data available. The user only needs to compile a list of species, along with the binding energies and desorption fractions for each grain species, the rest of this process is set up automatically (see Sect. 2.3). Whilst any physics module could activate the desorption process at specific temperatures, this process is strongly linked to the cloud model discussed in Sect. 2.2.1.

2.1.5 Grain Surface Reactions

Proper treatment of reactions between molecules on the grain surfaces is subject of debate. As a result, UCLCHEM has implemented grain chemistry in three ways. The method used in the basic network supplied with the code, and in the network used for the applications below, is to include hydrogenation at freeze out. For example, 10% of CO freezing onto the grains is assumed to go on to form CH3OH in the network used here and so it is immediately added to the grain abundance of CH3OH.

More complicated treatments are possible without editing the code. By including reactions with modified version of the rate constants in the network, grain surface reactions can be treated in the same way as gas-phase reactions. Further, UCLCHEM calculates the rates of reactions by diffusion across the grain surfaces for reactions with a third reactant labelled ”DIFF”. In this case, the binding energy of the two actual reactants are used to calculate their diffusion rates and included in a modified rate equation according to the treatment of Occhiogrosso et al. (2014).

2.2 Physical Model

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

The “cloud” module is the main physics module of UCLCHEM. This sets up a line of parcels from the centre to the edge of a cloud of gas and controls the density, temperature and visual extinction of each effectively giving a 1D model of a cloud. The conditions at the outer edge of the cloud are set by an interstellar radiation field and a value of visual extinction at the cloud edge. The distance from the edge of the cloud and the densities of the parcels closer to the edge are used to calculate the visual extinction and hence the radiation field for each interior parcel but parcels are otherwise treated separately.

The model works in two phases. In phase 1, this module is most often used to follow the collapse of the cloud from a diffuse medium to the density required for phase 2. Starting from elemental abundances only, this produces a set of gas and grain molecular abundances that is self consistent with the network. These can be used as the starting values for phase 2 rather than assuming a set of initial abundances.

In phase 2, the temperature increases as the cloud collapses and so the envelope around a forming protostar can be modelled. As discussed above, sublimation of species from the grain is dependent both on the temperature and rate of change of temperature. In the cloud model, the temperature of the gas as a function of time and radial distance from the protostar is calculated at each time step, and compared to pre-calculated temperatures to determine when major thermal desorption events should occur. The temperature profiles and desorption peaks are from Viti et al. (2004) where the time-dependent temperature profiles and desorption temperatures for a range of high stellar masses were calculated. This has since been extended with radial dependency and lower masses. Alternatively, a module that reads and interpolates the output of hydrodynamical codes can be used to post-process the chemistry of hydrodynamical models.

Phase 1 can also be used to study cold gas. Gas in a steady state can be modelled by setting all the required parameters and turning collapse off. Further, a number of different collapse modes have been coded into the model. The standard collapse used to create a cloud for phase 2 is the freefall collapse taken from Rawlings et al. (1992). Parameterizations for the collapse of a Bonnor-Ebert sphere are also included so that the collapse of cold cores can be studied. These are described in more detail in Priestley et al. (in prep.). Each collapse model calculates the density of a parcel by following the central density of the Bonnor-Ebert sphere with time and then scaling for the radial distance of a parcel from the centre of the core. The chemistry of each parcel is evolved separately so chemical abundances as a function of radius can be studied for different collapse modes and compared to observations of collapsing cores.

2.2.2 C-shock

The “C-shock” model is the parameterized model of continuous or C-type shocks from Jimenez-serra et al. (2008) adapted for use in UCLCHEM. The parameterization calculates the evolution of the velocities of the ion and neutral fluids as a function of time and deduces the physical structure of the shock from simple principles and approximations. In addition to the changes in density and temperature, the shock model includes a treatment of sputtering. The saturation time, tsat, is defined in Jimenez-Serra et al. (2008). It corresponds to the time at which the abundances of species sputtered from the icy mantles change by less than 10% in two consecutive time steps in their calculations. This gives a measurement of the time-scales at which most of the material that should be sputtered from dust grains will have been released. As the timescales in which the mantles sputter are short in the typical conditions in molecular clouds (Jimenez-serra et al., 2008) this is treated by releasing all of the mantle material into the gas phase at tsat. The equations used to calculate the saturation time do not depend on molecular mass and so tsat is the same for all species. See Jimenez-serra et al. (2008) Appendix B for more details.

Whilst sputtering is included, there are two major grain processes that are omitted. Grain-grain collisions in the shock can cause both shattering and vaporisation. Shattering is the breaking of grains in collisions and primarily has the effect of changing the dust grain size distribution. This produces hotter and thinner shocks, particularly for shocks with higher speeds or weak magnetic fields (Anderl et al., 2013). Vaporisation is the desorption of mantle species due to grain heating from collisions. The inclusion of grain-grain processes vastly increases the complexity of the code and would not allow the code to be run so quickly with low computing power. Shattering in particular is likely to greatly affect the abundance of Si-bearing species as the destroyed dust grains would add to the gas phase abundances. It should be noted however, that in models where SiO is included in the mantle, Anderl et al. (2013) find SiO abundances that are consistent with this model. The interested reader is referred to recent work on grain-grain processes in shocks for a more detailed view of this (Guillet et al., 2010; Anderl et al., 2013).

The model discussed here has been successfully used to explain the behaviour of molecules like H2O, NH3, CS and H2S in shocks like L1157-B1 (Viti et al., 2011; Gómez-Ruiz et al., 2014; Holdship et al., 2016). The parameterization has been extensively tested against MHD C-shock simulations and performs well across a wide range of input parameters. Importantly, it is computationally simple. This allows the chemistry to be the main focus without the computational costs being too prohibitive to run large grids of models to compare quickly to observations.

2.3 Network-MakeRates

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.

3 Application I: Collapsing Cores

To provide an example of the capabilities of UCLCHEM, the results of a model run for a collapsing core of gas is presented here. In the earliest stages of star formation, protostars form from collapsing cores of gravitationally bound gas. However, the evolution of the density profiles of these cores is not well understood. UCLCHEM can be used with a variety of theoretical time-dependent density profiles. In a separate paper, Priestley et al. (in prep) will discuss these profiles and their chemical effects in detail. Here, a simple and commonly considered profile is demonstrated, that of a marginally super-critical Bonnor-Ebert (BE) sphere (Bonnor, 1956). The link to the code provided in previous sections will contain the most up to date release of UCLCHEM, a permanent record of the code used for this work is available at http://dx.doi.org/10.5281/zenodo.580044.

BE spheres are a solution to the hydrostatic equilibrium equations and have been shown to fit the observed density profiles of prestellar cores. This typically requires that the cores are close to the critical mass which would make them unstable to collapse (Kandori et al., 2005). The cloud model includes a parameterization of simulations by Foster & Chevalier (1993) who consider a BE sphere at the critical mass for stability and then increase the density throughout the sphere by 10% making it marginally super-critical. Thus, the sphere collapses and the chemistry of a prestellar core collapsing in this fashion can be modelled. Such applications are useful as it may be possible to discriminate between different collapse models from the observed chemical differentiation in the core, even if the density profile cannot be confidently described due to the lack of a sufficiently high angular resolution.

3.1 Model Setup

The model runs consider 50 gas parcels equally distributed from the core to the edge of a cloud of size 0.2 pc and with an initial density of $2\text{\times}{10}^{2}\text{\,}{\mathrm{cm}}^{-3}$. This cloud is embedded in a larger medium which adds 1.5 magnitudes to the visual extinction at the cloud edge. This cloud maintains a temperature of $10\text{\,}\mathrm{K}$ as it collapses according to the freefall equation,

[TABLE]

where $n_{H}$ is the hydrogen nuclei number density, $n_{0}$ is the initial central density and mH 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\text{\times}{10}^{4}\text{\,}{\mathrm{cm}}^{-3}$, the initial density used for the BE collapse. This takes approximately $3.5\text{\,}\mathrm{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 which start at $2\text{\times}{10}^{4}\text{\,}{\mathrm{cm}}^{-3}$ and collapse to $1\text{\times}{10}^{8}\text{\,}{\mathrm{cm}}^{-3}$ in either freefall collapse or according to the BE collapse parameterization.

3.2 Results

In order to test the model, it is useful to compare to observations. Tafalla et al. (2004) (T04) present observations of two starless cores, L1517B and L1498, which are well approximated by BE spheres with central densities of approximately $2.2\text{\times}{10}^{5}\text{\,}{\mathrm{cm}}^{-3}$ and $0.94\text{\times}{10}^{5}\text{\,}{\mathrm{cm}}^{-3}$ respectively. Though the reported characteristics for each core are very similar, L1498 is used as a comparison to the model here. T04 derive temperature and density profiles from the dust emission. They then derive CS, CO, NH3 and NH${}_{2}^{+}$ abundance profiles from radiative transfer fitting to the flux profile, favouring the simplest profile that recreated the data.

They find CO and CS to be depleted in the denser centres of the cores compared to the outer radii and that this is adequately represented as a step function. In contrast, NH3 is enhanced by a factor of a few in the centre; a trend they represent as a simple radial power law. N2H+ has a constant abundance for L1517B and a slight drop at large radii for L1498. Chemical models including the one presented here accurately reflect the behaviour of CS and CO. However, many struggle to reproduce the NH3 and N2H+ behaviour (Aikawa et al., 2003). Figure 1 shows results of UCLCHEM using the BE model used to fit these observations with a central density of $0.95\text{\times}{10}^{5}\text{\,}{\mathrm{cm}}^{-3}$.

From Figure 1 it can be seen that the BE sphere modelling with UCLCHEM meets with some success. For CO and CS, T04 find each can be treated as having constant abundance outside of some cut off radius, with the abundance in the core being greatly depleted. For simplicity, T04 take a central abundance of zero and whilst this is clearly not the case for UCLCHEM, there are some caveats. The T04 observations are of C17O and C18O and they report C18O abundances. For the comparison here, those values have been multiplied by 550, the C16O:C18O ratio (Wilson, 1999). Reducing the CO values by this amount gives a central abundance of 2$\times 10^{-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$\times 10^{-11}$ and therefore in agreement with zero.

T04 provide an analytic expression for NH3 which is plotted in Figure 1, the main finding is that NH3 increases by a factor of a few towards the centre of the core. The UCLCHEM model predicts the opposite, increasing with radius. UCLCHEM similarly struggles to reproduce the observations of N2H+, decreasing both towards the centre and edge of the core rather than the approximately constant abundance found in T04. It appears from the model outputs that at middle radii, N2H+ is forming through a reaction between H${}_{3}^{+}$ and N2. 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 centre reduces the amount of 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 and not the radius so there is no differentiation with distance from the core centre. Further, the freefall model reaches $1\text{\times}{10}^{5}\text{\,}{\mathrm{cm}}^{-3}$ much more quickly than the BE model and so the depletion of each species is much lower than expected as 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 short comings of the BE model could be due to a number of factors. Spherical symmetry is assumed in both the UCLCHEM model and in the radiative transfer modelling 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 behaviour of both N2H+ and NH3 are not reproduced by the BE model could indicate a problem with the nitrogen chemistry in the dense gas at core centre. However, the model is clearly superior to a simple freefall collapse and a more in depth investigation into collapse modes and initial conditions could link the observed abundance profiles and the chemical modelling. The strength of UCLCHEM is that it has been specifically designed to have easily modified, and easy to implement, physics modules. To illustrate the code’s range of physical applications, we present a second example with the chemical modelling of C-type shocks in the next section.

4 Application II: Shock Chemistry

4.1 How Can We Trace 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 modelling 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 H2O and NH3 lines were used to determine the properties of a shock. Some molecules, such as H2O trace the full shock as they are not destroyed in the hotter parts of the post shock gas, which then dominates the high-J emission. Others, such as NH3 may not and are instead destroyed in the hotter parts of the post shock gas. This behaviour 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.

4.2 Model Setup

As explained in previous sections, the model is two phased. Phase 1 starts from purely atomic gas at low density ( nH=102 cm*-3* 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 and 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 ran for this work. The phase 1 models ran here used a temperature of 10 K, a radiation field of 1 Habing and the standard cosmic ray ionization rate of 1.3 $\times$ 10*-17* s*-1*. A radius of 0.05 pc is chosen for the cloud and 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 magnitudes 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 varied 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. Note, not all of the listed variables are independent, tsat and Tmax are functions of the pre-shock density (nH,0) and shock velocity (Vs). The Tmax, Vs relationship is extracted from figures 8 and 9 of Draine et al. (1983) which give values of maximum temperature for pre-shock densities 104 and 106 cm*-3*. Since no values are available when the pre-shock density is 105 cm*-3* the maximum temperature is calculated as the average of the 104 and 106 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 Fig. 2, using model 27 as an example.

4.3 Shock Tracers

Due to the fact freezing is efficient in cold, dense clouds, many species are highly abundant in the solid phase in 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. Hydrogenated molecules in particular tend to fall into this category.

H2O and NH3 are highly abundant in the ices but not the gas phase in cold molecular clouds. Of the two, NH3 is a more useful tracer of shocks due to the difficulties in observing H2O from the ground. In the models, NH3 has a fractional abundance of 10*-14* in the pre-shock gas but can increase to 10*-5* in a shock (see Figure 3). H2S and CH3OH similarly go from extremely low abundances in the pre-shock gas to X$\sim 10^{-5}$ in a shock. Fractional abundances as a function of shock speed for different pre-shock densities are shown in Fig. 3.

An important caveat to this is that the model assumes all of the grain surface material is released into the gas phase unchanged. This may not be the case, as 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 CH3OH being released from ice mantles. They find that more than 90% of CH3OH 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 Fig 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*-1*) sputter the grains so 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, whilst bulk abundances are not sensitive to shock speed, line profiles can be used to provide further constraints.

4.4 Determining Shock Properties

The focus so far has been on average abundances across the whole shock, the observed abundance if one just 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 NH3 and H2O profiles from the L1157-B1 region and showed that when scaled to be equal at their peak emission velocity, NH3 dropped off much more quickly with increasing velocity. Viti et al. (2011) demonstrated that this is consistent with NH3 being at lower abundance in the hotter gas, with chemical models showing that a shock with v${}_{s}\textgreater$ 40 km s*-1* and pre-shock density of nH=105 cm*-3* was sufficient to destroy NH3 which, unlike H2O 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 the difference in the terminal velocity of NH3 and H2O was correlated with the shock velocity inferred from the H2O terminal velocity for a sample of protostellar outflows. Chemical modelling in the work corroborated this by showing more NH3 being lost in faster shocks. Holdship et al. (2016) showed H2S exhibits the same behaviour as NH3, dropping in abundance in high velocity gas. Given the promise of this method, molecules that exhibit this behaviour 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 modelling in Viti et al. (2011).

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

A species does not necessarily fit one category or the other for all shock conditions. The behaviour is the result of temperature and density profiles of the shock and so depends on the shock speed and initial density. In Figure 5 all modelled pairs of shock speed and initial density are shown with lists of shock-tracing and shock-destroyed species in blue and red respectively. Species which exhibit more complicated behaviour and cannot be easily categorized are simply omitted. Some species, such as H2O are almost always found to fall in one category, this allows them to be consistently used for a given purpose. ie. H2O can be used to trace the full shock in most cases. However, some species will demonstrate a certain behaviour only in very specific conditions and can be used to limit possible shock properties when this is observed. For example, H2O does not trace the full shock if the initial density is 103 cm*-3* and the shock is low speed. Such species are displayed in bold in Figure 5.

The shock-destroyed and shock-tracing molecules are both enhanced by the sputtering caused by the shock. This is important as 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 behaviour 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 and so can be reliably used to trace the full extent of the shock.

5 Summary

A publicly available, time dependent gas-grain chemical model, UCLCHEM, has been described in detail. The code solves the coupled system of ODEs that represents the chemical network when using the rate equation method of modelling chemistry. Various non-thermal and thermal desorption mechanisms are implemented along with gas phase reactions and a simple grain network. The physical models that can be used to control the gas conditions that in turn affect the chemistry are also described.

We presented the code via two applications of the new UCLCHEM: a simple prestellar collapse model which represents a typical use of the code, and an investigation into molecular tracers for C-type shocks. This code is available at https://uclchem.github.io and is actively maintained and developed by research teams at UCL and Queen Mary’s College London. Diffusion chemistry on the grain surfaces and improvements to the shock model implementation are in progress and users will benefit from continued development.

6 Acknowledgements

J.H. is funded by an STFC studentship. I.J.-S. acknowledges the financial support received from the STFC through an Ernest Rutherford Fellowship (proposal number ST/L004801/2). The authors would like to thank D. Quenard and A. Coutens for helpful discussions that contributed to the code.