Chempl: a playable package for modeling interstellar chemistry

In interstellar or protoplanetary conditions, where the majority of the gas under consideration has a density typically in the range of ten to the power of a few per cm³, and densities higher than ∼10¹³ cm^–3 rarely need to be considered, three-body chemical reactions are not important and are excluded in astrochemical networks.

When only first- and second-order reactions are included, Equation 1 can be written in component form as

$$\begin{eqnarray}\begin{array}{lll}\displaystyle \frac{d{y}_{i}}{dt} & = & \displaystyle \sum _{j\in { {\mathcal F} }_{i,1}}{k}_{j}{y}_{j,1}+\displaystyle \sum _{j\in { {\mathcal F} }_{i,2}}{k}_{j}{y}_{j,1}{y}_{j,2}\\ & & -\displaystyle \sum _{j\in {{\mathscr{D}}}_{i,1}}{k}_{j}{y}_{j,1}-\displaystyle \sum _{j\in {{\mathscr{D}}}_{i,2}}{k}_{j}{y}_{j,1}{y}_{j,2}{\rm{\hspace{0.17em}}},\end{array}\end{eqnarray} \tag{ 3 }$$

where ${ {\mathcal F} }_{i,1}$ and ${ {\mathcal F} }_{i,2}$ are the set of first- and second-order reactions that form species i, ${{\mathscr{D}}}_{i,1}$ and ${{\mathscr{D}}}_{i,2}$ are the set of first- and second-order reactions that destroy species i, k_j is the rate coefficient of reaction j, and y_j,1 and y_j,2 are the abundances of the first and second (when applicable) reactants of reaction j. Note that for $j\in {{\mathscr{D}}}_{i,1}$, y_j,1 is always simply y_i , and for $j\in {{\mathscr{D}}}_{i,2}$, one of y_j,1 and y_j,2 (or both) must be the same as y_i .

Examples of first-order reactions are photodissociation, cosmic-ray ionization, adsorption, and desorption reactions. These reactions are not strictly "uni-molecular", since when one of these reactions happen there must be two particles interacting, with one being the chemical species we care about and the other being a photon, a cosmic-ray particle, or a dust grain. They are treated as first order because the fluxes of photons and cosmic-rays are external inputs and are independent of local conditions, and in the case of adsorption the dust grains can be "reused" (hence effectively being constant). Examples of second-order reactions are neutral-neutral, neutral-ion, dissociative recombination, and radiative association reactions. One may consult Wakelam et al. (2012) or McElroy et al. (2013) for a list of possible types of gasphase reactions relevant for astrochemistry.

In this section we describe formulae for calculating the rates for different types of reactions. While these equations can be found in many modeling papers in astrochemistry, they are presented here for reference. More detailed pedagogical discussions can be found in Du (2012). For some reaction types, more realistic expressions for their rates may need to be used when modeling some specific problems. In CHEMPL it is possible to add new reaction types with associated handlers to accomplish this.

2.2.1. Cosmic-ray reactions

For ionization directly caused by cosmic-rays, the rate coefficient is simply a constant:

$$\begin{eqnarray}&&k=\alpha {\rm{\hspace{0.17em}}}{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 4 }$$

where α takes different value for different species. For cosmic-ray-induced photoreactions, the rate coefficient is

$$\begin{eqnarray}&&k=\alpha {(\displaystyle \frac{T}{300{\rm{\hspace{0.17em}}}{\rm{K}}})}^{\beta }\displaystyle \frac{\gamma }{1-\omega }{\rm{\hspace{0.17em}}}{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 5 }$$

where α, β, and γ are constants, and ω is the albedo of dust grains to dissociative photons, to which we assign a value of 0.6 by default. Both direct and induced reaction rates due to cosmic-rays can be rescaled by the actual cosmic-ray ionization rate (relative to the fiducial rate of ζ₀ = 1.36 × 10^–17 s^–1, see McElroy et al. 2013).

2.2.2. Photoreactions

For photoreactions (photodissociation and photoionization), the rate coefficient is

$$\begin{eqnarray}&&k=\alpha \exp (-\gamma {A}_{{\rm{V}}}){\rm{\hspace{0.17em}}}{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 6 }$$

where α is the unattenuated rate, γ is a constant to scale the visual extinction A_V to the UV band. These rates can also be rescaled by a G₀ parameter, namely the ratio between the actual UV field and the Draine interstellar field (Draine 1978).

For species such as H₂, CO, and N₂, self- and mutual- shielding can be important, and Equation (6) cannot be used. Approximate analytical formulae or tabulated rates (e.g., Draine & Bertoldi 1996; Morris & Jura 1983; Visser et al. 2009; Li et al. 2013) with dependencies on column densities and temperature are available in this case.

2.2.3. Generic two-body reaction rates

The modified Arrhenius equation is often used to calculate the two-body reaction rate coefficients:

$$\begin{eqnarray}&&k=\alpha {(\displaystyle \frac{T}{300{\rm{\hspace{0.17em}}}{\rm{K}}})}^{\beta }\exp (\displaystyle \frac{-\gamma }{T}){\rm{\hspace{0.17em}}}{{\rm{cm}}}^{3}{\rm{\hspace{0.17em}}}{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 7 }$$

where α, β, and γ are constants provided by the reaction network databases. The reaction rate, or stated more clearly, the number of occurrence of a reaction per unit time per unit volume, can be expressed as

$$\begin{eqnarray}&&k{n}_{i}{n}_{j},\end{eqnarray} \tag{ 8 }$$

where n_i and n_j are the volume density (in cm^-3) of the two reactants of a two-body reaction. The time derivative of the density of species i caused by this reaction is then (assuming i ≠ j)

$$\begin{eqnarray}&&\displaystyle \frac{d{n}_{i}}{dt}=-k{n}_{i}{n}_{j}.\end{eqnarray} \tag{ 9 }$$

In astrochemistry the abundance of a species relative to hydrogen nuclei is usually used, i.e., n_i = y_i n_H, and the time derivative of y_i is then

$$\begin{eqnarray}&&\displaystyle \frac{d{y}_{i}}{dt}=-k{n}_{{\rm{H}}}\cdot {y}_{i}{y}_{j}.\end{eqnarray} \tag{ 10 }$$

Note that the modified Arrhenius equation is an empirical formula and does not apply for all reactions for all temperature ranges (see Sect. 3.2.1).

2.2.4. Adsorption of gas phase species onto the surfaces of dust grains

The number of particles of species i adsorbed onto dust grains per unit volume per unit time is (Du & Bergin 2014)

$$\begin{eqnarray}&&s\sigma {v}_{{\rm{T}}}(i){n}_{{\rm{dust}}}{n}_{i},\end{eqnarray} \tag{ 11 }$$

where s is the sticking coefficient, σ = π a² is the dust cross section with a being the dust grain radius, n_dust and n_i are the volume density of dust grains and species i, and v_T(i) is the thermal speed of species i:

$$\begin{eqnarray}&&{v}_{{\rm{T}}}(i)={(\displaystyle \frac{8{k}_{{\rm{B}}}T}{\pi {m}_{i}})}^{1/2},\end{eqnarray} \tag{ 12 }$$

where k_B is the Boltzmann constant, and m_i is the mass of a particle of species i. Written in terms of abundance relative to hydrogen nuclei, the time derivative is

$$\begin{eqnarray}&&\displaystyle \frac{d{y}_{i}}{dt}=-s\sigma {v}_{{\rm{T}}}(i){n}_{{\rm{dust}}}{y}_{i}=-s\sigma {v}_{{\rm{T}}}(i){n}_{{\rm{H}}}\eta {y}_{i},\end{eqnarray} \tag{ 13 }$$

where η is the dust-to-hydrogen number ratio, which can be further expressed as:

$$\begin{eqnarray}\begin{array}{lll}\eta & = & {\eta }_{{\rm{M}}}\cdot \displaystyle \frac{\mu {m}_{{\rm{H}}}}{\displaystyle \frac{4\pi }{3}{a}^{3}{\rho }_{{\rm{d}}}}\\ & = & 2.8\times {10}^{-12}(\displaystyle \frac{{\eta }_{{\rm{M}}}}{0.01})(\displaystyle \frac{\mu }{1.4}){(\displaystyle \frac{a}{0.1{\rm{\hspace{0.17em}}}\mu {\rm{m}}})}^{-3}\\ & & \times {(\displaystyle \frac{{\rho }_{{\rm{d}}}}{2{{\rm{g}}{\rm{cm}}}^{-3}})}^{-1},\end{array}\end{eqnarray} \tag{ 14 }$$

where η_M is the dust-to-grain mass ratio, with a canonical value of 0.01, m_H is the mass of a hydrogen nucleus, μ is the mean molecular mass per hydrogen nucleus, and ρ_d is the material density of dust grains. So the rate coefficient of adsorption reaction is

$$\begin{eqnarray}\begin{array}{lll}k & = & s\pi {a}^{2}{(\displaystyle \frac{8{k}_{{\rm{B}}}T}{\pi {m}_{i}})}^{1/2}{n}_{{\rm{H}}}{\eta }_{{\rm{M}}}\cdot \displaystyle \frac{\mu {m}_{{\rm{H}}}}{\displaystyle \frac{4\pi }{3}{a}^{3}{\rho }_{{\rm{d}}}}\\ & = & s\displaystyle \frac{3}{4a}{(\displaystyle \frac{8{k}_{{\rm{B}}}T}{\pi {m}_{i}})}^{1/2}{n}_{{\rm{H}}}{\eta }_{{\rm{M}}}\displaystyle \frac{\mu {m}_{{\rm{H}}}}{{\rho }_{{\rm{d}}}}\\ & = & 4.03\times {10}^{-12}{\rm{\hspace{0.17em}}}{{\rm{s}}}^{-1}{\rm{\hspace{0.17em}}}s{(\displaystyle \frac{a}{0.1\mu {\rm{m}}})}^{-1}{(\displaystyle \frac{T}{10{\rm{K}}})}^{1/2}\\ & & \times {(\displaystyle \frac{{m}_{i}}{{m}_{{\rm{H}}}})}^{-1/2}(\displaystyle \frac{{n}_{{\rm{H}}}}{{10}^{5}{{\rm{cm}}}^{-3}})(\displaystyle \frac{{\eta }_{{\rm{M}}}}{0.01})(\displaystyle \frac{\mu }{1.4})\\ & & \times {(\displaystyle \frac{{\rho }_{{\rm{d}}}}{2{{\rm{gcm}}}^{-3}})}^{-1}.\end{array}\end{eqnarray} \tag{ 15 }$$

The sticking coefficient s can be determined by experiments and is usually in the range 0.1–1. Note that for CO molecule, the above expression gives a depletion timescale of ∼4 × 10⁴ yr for the fiducial physical parameters in the parentheses.

2.2.5. Thermal desorption from dust grain surfaces

The thermal desorption rate is

$$\begin{eqnarray}&&k={\nu }_{i}\exp (-{E}_{i}/{T}_{{\rm{d}}}),\end{eqnarray} \tag{ 16 }$$

where T_d is the dust temperature (which may be different from the gas temperature), and ν_i is the vibrational frequency of species i on the dust grain surface:

$$\begin{eqnarray}\begin{array}{lll}{\nu }_{i} & = & {(\displaystyle \frac{2{\rho }_{{\rm{S}}}{E}_{i}}{{\pi }^{2}{m}_{i}})}^{1/2}\\ & = & 4.09\times {10}^{12}{\rm{\hspace{0.17em}}}{\rm{Hz}}{\rm{\hspace{0.17em}}}{(\displaystyle \frac{{\rho }_{{\rm{S}}}}{{10}^{15}{\rm{\hspace{0.17em}}}{{\rm{cm}}}^{-2}})}^{1/2}\\ & & \times {(\displaystyle \frac{{E}_{i}}{{10}^{3}{\rm{\hspace{0.17em}}}{\rm{K}}})}^{1/2}{(\displaystyle \frac{{m}_{i}}{{m}_{{\rm{H}}}})}^{-1/2}.\end{array}\end{eqnarray} \tag{ 17 }$$

E_i is the desorption energy of species i, and ρd is the surface site density typically taken to be 10¹⁵ cm^–2.

2.2.6. Cosmic-ray desorption

When energetic cosmic-ray particles hit a dust grain, the latter will be promptly heated, leading to evaporation of surface species. The evaporation only lasts a very brief period, because the dust grain cools down quickly. This cosmic-ray-induced desorption is usually implemented using the formulation of Hasegawa & Herbst (1993a):

$$\begin{eqnarray}&&{k}_{{\rm{CR}}}={k}_{i}(70{\rm{K}})f(70{\rm{K}}),\end{eqnarray} \tag{ 18 }$$

where k_i (70 K) is the thermal desorption rate of species i calculated with Equation (16) for a dust grain temperature of 70 K, and f(70 K) is the fraction of time a dust grain remains heated around 70 K by a cosmic-ray particle, estimated to be ∼3.16 × 10^–19. The number "70 K" is the peak temperature a dust grain can achieve from cosmic-ray heating balanced by evaporative cooling.

2.2.7. Photodesorption

UV photons impinging on dust grains can eject surface molecules into the gas phase. The photodesorption rates have been experimentally measured or theoretically computed for ice species such as CO, H₂O, N₂, etc. (Öberg et al. 2009a,b; Potapov et al. 2019; González Díaz et al. 2019), and empirical formulae have been used for calculating their rates.

We use a formula of the form for the rate coefficient

$$\begin{eqnarray}&&k=F{Y}_{i}\sigma \eta (1-{e}^{-{n}_{{\rm{L}}}(i)}),\end{eqnarray} \tag{ 19 }$$

where F is the UV flux (in counts cm^–2 s^–1), Y_i is the photo yield for species i, σ is the dust cross section, and η is the dust-to-gas number ratio. The number of monolayers of species i is:

$$\begin{eqnarray}&&{n}_{{\rm{L}}}(i)=\displaystyle \frac{{y}_{i}}{\eta {N}_{{\rm{S}}}},\end{eqnarray} \tag{ 20 }$$

where y_i is the abundance relative to hydrogen nucleus, and N_S is the number of surface reaction sites per grain (∼10⁶ for a dust grain radius of 0.1 μm). The yield Y may depend on the dust temperature. For ice species without experimental data, Y is assumed to be 10^–4 by default. However, Equation (19) does not accurately reflect all the complexities found in experiments and should be treated with care when photodesorption is expected to be important.

2.2.8. Chemical desorption

2.2.9. Two-body reaction rates on dust grain surfaces

For a surface reaction of the form i + j → ⋅, its rate coefficient is

$$\begin{eqnarray}&&k=-({k}_{{\rm{mig}},i}+{k}_{{\rm{mig}},j})\displaystyle \frac{{p}_{ij}}{{N}_{\eta }},\end{eqnarray} \tag{ 21 }$$

where k_mig,i is the surface migration rate of species i, p_ij is the probability for i and j to cross over the reaction barrier when they arrive at the same surface site, η is the dust-to-grain number ratio as defined in Equation (14), and N_S is the number of reactive sites per grain. The surface migration can be due to either thermal hopping:

$$\begin{eqnarray}&&{k}_{{\rm{mig}},i}={\nu }_{i}\exp (-{E}_{{\rm{diff}},i}/{T}_{{\rm{d}}})\end{eqnarray} \tag{ 22 }$$

or quantum tunneling (at very low temperatures):

$$\begin{eqnarray}&&{k}_{{\rm{mig}},i}={\nu }_{i}\exp (-\displaystyle \frac{2l}{\hslash }\sqrt{2{m}_{i}{E}_{{\rm{diff}},i}}),\end{eqnarray} \tag{ 23 }$$

where ν_i is the characteristic vibration frequency as defined in Equation (17), E_diff,i is the energy barrier against surface migration for species i, which is usually taken to be a fraction (0.3–0.5) of the desorption energy, and l is the width of the barrier. The reaction probability p_ij can be calculated similar to k_mig,i , except that the factor ν_i before the exponential function is dropped, and the diffusion barrier is replaced by the reaction barrier.

As a side note, the formation of H₂ molecules on dust grain surfaces is treated as a normal two-body surface reaction. This is different from some of the previous works, in which the H₂ formation rate is taken to be half of the adsorption rate of H atoms (see, e.g., Holdship et al. 2017).

As already described in previous sections, molecules can be adsorbed on the surface of dust grains. After being adsorbed, these molecules are able to migrate over the dust grain surface and may react with each other when two molecules meet at the same surface site. Since a dust grain has a limited surface area, when many molecules have adsorbed on its surface, newly adsorbed ones will likely sit on top of another. In this way, multiple ice layers will form on the dust grains. In such a description of dust grain structure, usually only molecules on the topmost layer are assumed to be mobile and reactive, while molecules underneath (in the mantle) are "frozen" and nonreactive (however, see Chang & Herbst 2014, 2016). When molecules in the topmost layer evaporate, the underneath layer will be revealed and becomes the new reactive layer. This is the so-called three-phase model of interstellar chemistry, in which the three phases are: gas phase, dust grain surface, and dust grain mantle.

Chempl has implemented the three-phase description of gas-grain chemistry, following the approach of Hasegawa & Herbst (1993b) and Garrod & Pauly (2011) (for a detailed description, see also Du 2012). For each surface species (except for very weakly adsorbed species such as H₂) a new mantle species is included in the network. For species i with surface population per grain i, its evolution equation is

$$\begin{eqnarray}\begin{array}{lll}\displaystyle \frac{d{n}_{{\rm{s}}}i}{dt} & = & \displaystyle \sum _{j,k}{k}_{j,k}{n}_{{\rm{S}}}(j){n}_{{\rm{S}}}k-\displaystyle \sum _{j}{k}_{ij}{n}_{{\rm{S}}}(i){n}_{{\rm{S}}}(j)\\ & & +{k}_{{\rm{ads}}}(i){n}_{{\rm{G}}}(i)-{k}_{{\rm{ads}}}{n}_{{\rm{S}}}(i)/{N}_{{\rm{S}}}\\ & & +{k}_{{\rm{des}}}{n}_{{\rm{M}}}(i)/\max ({n}_{{\rm{M}}},{n}_{{\rm{S}}})-{k}_{{\rm{des}}}(i){n}_{{\rm{S}}}(i),\end{array}\end{eqnarray} \tag{ 24 }$$

and the mantle population of i evolves according to

$$\begin{eqnarray}&&\displaystyle \frac{d{n}_{{\rm{M}}}(i)}{dt}={k}_{{\rm{ads}}}{n}_{{\rm{S}}}i/{N}_{{\rm{S}}}-{k}_{{\rm{des}}}\,{n}_{{\rm{M}}}i/\max ({n}_{{\rm{M}}},{n}_{{\rm{S}}}).\end{eqnarray} \tag{ 25 }$$

In Equation (24) the two summations are over two-body reactions producing and consuming species i. k_ads/des(i) is the adsorption/desorption rate coefficient of i, and n_G(i) is the gas phase population of i. k_ads/des is the total adsorption/desorption rate of all the species, n_S is the total population of all surface species, n_M is the total population of all mantle species, and N_S is the number of sites per grain. The quantity "population-per-grain" ${n}_{{\rm{S}}}(i)$ is related with abundance with respect to hydrogen nuclei y_i by (see also Eq. (20))

$$\begin{eqnarray}&&{y}_{i}={n}_{{\rm{S}}}(i)\eta ,\end{eqnarray} \tag{ 26 }$$

where η is the dust-to-hydrogen number ratio defined before.

Note that the three-phase formulation is statistical in nature. There are more accurate treatments such as the microscopic kinetic Monte Carlo method (Cuppen & Herbst 2007; Garrod 2013), which tracks the location of each particle on the dust grain surface.

With the inclusion of three-phase chemistry, the code will run slower, not just because of the increase in the number of variables due to addition of mantle species, but also because of the increased non-linearity of Equation (24) relative to Equation (3).

While Chempl can be used as a stand-alone program, it is mainly intended to be used as a Python module. Here we describe the basic workflow of Chempl when used in this mode. A minimal working example is shown in Table 1.

• Create the ChemModel object. This is the central Python object to work with in Chempl. During the instantiation the user may provide the file names of the reaction network, the initial abundances, and the enthalpies (needed for reaction heat calculation). These files could also be loaded later with the corresponding helper functions. The user could also add reactions one by one using the Python code, and this is the way to load reaction network stored in different formats.
• Set the physical parameters. Similarly, this could be done with an input file, or using a helper function to set the value of each parameter by its name. They could also be provided in the previous step. Optionally, the user could add time-dependency for certain parameters by providing a lookup table. Internally the code linearly interpolates the table to evaluate time-dependent parameters at arbitrary time.
• Do some "boilerplate" preparation work. These include, e.g., sorting out reactions of different types, getting the elemental composition of each species and calculating their molecular masses and quantum mobilities, and calculating the heat of each reaction. Such boilerplate work usually only need to be done once for each model, unless during some unusual investigations (e.g., changing the reaction network while the code is running).
• Setup the ODE solver. This usually amounts to setting the solver ID, which is itself optional (could be useful when running in multi-process mode). Details of the solver setup and preparation are concealed in the C++ code, which includes constructing the sparse structure of the ODE system, allocation of the needed working memories, and assignment of flags defining the working modes of the solver.
• Feed the current time and abundances, and the required output time, to the solver, then wait for the output.
• Retrieve and store the output time (which could be different from the value provided in the previous step) and abundances.
• Terminate if the output time is equal to or greater than the targeted final evolution time; otherwise go to step 5.

We first run a model with physical conditions relevant for dark clouds. The physical parameters are listed in Table 2. The initial abundances are listed in Table 3. The network for this model is UMIST RATE12, extended by only adding the H₂ formation reaction on dust grain surface and the adsorption and evaporation of H and H₂, and no other adsorption/desorption processes and surface reactions are included.

The code takes 4–5 seconds of CPU time ² to reach 10⁸ yr. The resulting evolution of a few common species are plotted as solid curves in Figure 1. As can be seen, the system reaches steady state at ∼10⁷ yr, after which carbon mainly resides in CO, while oxygen mainly resides in CO and O₂. However, this model is not realistic at late times (≳10⁶ yr), because at the low temperature for this model, CO should freeze out onto dust grains. So this model should be only viewed as a code benchmark and it does not accurately describe the actual dark cloud chemical processes.

Zoom In Zoom Out Reset image size

So we run another model with the same initial abundances and physical parameters, but with a full three-phase chemical network. Nonthermal desorption (Garrod et al. 2007) by the heat released in surface reactions is also included. With this added complexity, the code takes ∼15 seconds of CPU time to reach 10⁸ yr. The results are plotted as dashed curves in Figure 1. For species such as CO and C, their early-time (≲10⁴ yr) evolution from the two models coincide. This is because their abundances are mainly determined by gas-phase reactions at early times. For species such as H₂O, OH, O₂, and NH₃, their evolution tracks from the two models depart away from each other from the very beginning. This is because they are efficiently formed on the dust grain surface, and the inclusion of surface reactions (together with nonthermal desorption) enhances their abundances from the very beginning.

Chemical models in which physical parameters (most frequently, density and temperature) change with time have been used in the study of warm-up phase of hot cores (Garrod & Herbst 2006). Technically, this time-dependency could be implemented with a series of models with constant physical parameters, and the parameters only gradually change from one model to another. In the series of models the final abundances of a previous model is used as initial condition of the next. However, such an implementation would be computationally inefficient, because the ODE solver would have to reinitialize itself for each individual model, which is a slow process. So it is necessary to add explicit time-dependency directly in a model, specifically, in the right-hand side of Equation (1) or Equation (3).

3.2.1. Continuity of reaction rates

In the RATE12 network each reaction has associated a temperature range in which the rate parameters are most reliable. Some reactions have two temperature ranges, which means a different set of rate parameters must be chosen at different temperatures. This could cause problems when temperature varies as a function of time: at the boundary of two temperature ranges, the reaction rates calculated from two versions of rate parameters are not necessarily continuous, which may confuse the ODE solver and the solution will become unreliable (or the solver will stop working and no solution can be obtained after the boundary of a temperature range is reached).

Chempl ensures continuity (smoothness actually) by joining the rates calculated with two sets of parameters using a weighting function of the form

$$\begin{eqnarray}&&g(x)=\displaystyle \frac{1}{1+{e}^{\pm (x-{x}_{0})/w}},\end{eqnarray} \tag{ 27 }$$

where x₀ is the interval edge, and w > 0 is the separation from one interval to another. The ± sign takes "+" when x₀ is the right edge, and takes "−" when x₀ is the left edge. The function g(x) has the property that for x < x₀, it quickly approaches unity, while for x > x₀, it quickly approaches zero. The choice of weighting function is arbitrary, in the sense that a weighting scheme is allowed as far as the transition across the boundary is smooth and without any artificial bumps.

Figure 2 shows the rate coefficients of two reactions as a function of temperature. As can be seen from the left panel, without a treatment such as the one described here, the rate coefficient would jump from one value to another when the temperature only slightly varies at the boundary (from 100 to 101 K). This may halt the ODE solver since it will have to use exceedingly small time steps (since temperature is a function of time in time-dependent models). With a continuous function, the ODE solver will be able to integrate much easier.

Zoom In Zoom Out Reset image size

Figure 3 and Figure 4 show two example models with time-dependent temperature and density. The reaction network used for them are RATE12 with the addition of adsorption and desorption reactions, while no surface reactions beyond the formation of H₂ are included. The time evolution of temperature and density in the two models are arbitrarily designed, while keeping in mind that they may mimic realistic processes occurring in certain interstellar environments. The first model simulates a compressing and heating event, which takes about 40 CPU seconds to reach 10⁷ yr, while the second one simulates a "wavy" or oscillating temperature and density profile, and takes about 16 CPU minutes to reach 10⁷ yr. The abrupt oscillation of temperature and density, especially at late times, limits the step size used internally in the ODE solver and significantly slows down the integration speed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, the evolution curves plotted in the bottom panel of Figure 4 demonstrate two types of behavior: those which follow temperature (CO, O₂, and NH₃), and those which follow density (H₂O, OH). For species of the first type, their abundances are mostly controlled by desorption, hence their abundances jumps when the dust temperature reaches their desorption temperature, and stays flat after that. For species of the second type, either they are not able to evaporate at dust temperatures ≲40 K (the highest dust temperature in the model), and/or their abundances are mostly determined by gas phase reactions, and their abundances tend to follow the evolution of density. The behavior of C and C+ are more complex. They seem to follow both density and the rise and drop of temperature, which is because their abundances are affected by gas phase reactions which involves species that can evaporate at a dust temperature less than 40 K.

Since we largely use the UMIST RATE12 network as a starting point for Chempl, we also run a model similar to the "carbon-rich circumstellar envelope model" presented in the RATE12 paper (McElroy et al. 2013). This may be viewed as a non-strict benchmark, because a strict one is possible only if the exact parameters (including physical constants) used in their model are adopted. Here we try to use the same input as them, though an exact match is not intended for. The initial abundances and density profile are taken from table 6 of McElroy et al. (2013), and the temperature profile is taken from Cordiner & Millar (2009). The self-shielding of CO from photodissociation is treated using the approach of Morris & Jura (1983), which is also the same as in McElroy et al. (2013). While this model calculates the spatial structure of a circumstellar envelope around a carbon-rich AGB star, however, since the gas in the envelope are launched from very close to the central star with a fully-specified velocity profile, there is a one-to-one mapping between time and radius. Thus the model can be viewed as one in which physical parameters vary with time, which is how we implement it.

The results of the calculation with Chempl are shown as solid curves in Figure 5, overlaid with data from McElroy et al. (2013) shown as dashed curves. It is clear that the two calculations match very well: two versions of most of the curves are almost identical within numerical errors. This close match means our implementation of the envelope structure with Chempl is correct (although there might be some fine details that are not reproduced in our model).

Zoom In Zoom Out Reset image size

As a final example, we create a model that simulates a plane-parallel photodissociation region. The density is taken to be uniformly 10³ cm^–3, the G₀ parameter of UV field at the edge is set to 100, while the dust temperature profile is calculated using equation (9.18) of Tielens (2005), and the gas temperature is assumed to be uniformly 100 K.

The calculation has to proceed from the surface layer (A_V∼ 0) to the deeper layers to take into account the self-shielding of H₂ and CO. Here the CO self-shielding is approximated using the same approach as in 3.3, while the H₂ self-shielding is implemented following Draine & Bertoldi (1996). More accurate treatment of self and mutual shielding can also be implemented (Visser et al. 2009; Li et al. 2013; Heays et al. 2017).

The resulting abundances of different molecules are shown as a function of depth (A_V) and column density are shown in Figure 6. As is commonly known, the H/H₂ transition occurs at A_V ∼ 1 – 2, while the C⁺/C/CO transition occurs at A_V ∼ 2 – 4 (the exact transition locations depend on the density profile and the intensity of the external UV field).

Zoom In Zoom Out Reset image size

As shown before, when the physical parameters are kept constant, Chempl runs relatively fast. So it is not necessary to parallelize each individual model. However, when doing a parameter study, in which many models with different configurations need to be computed, parallelization becomes helpful. In the case of Chempl, the current parallelization approach is to make use of the Python built-in "multiprocessing" module. With this module, it is possible to run multiple independent models at the same time, which can accelerate the computational speed by a factor depending on the number of CPU cores of the system. We call this approach of parallelization "model-level" parallelism ³ .

When the physical parameters explicitly change with time, especially when the time variation is oscillatory and spiky, as shown in the example model in Figure 4, the computational speed will be significantly reduced. We do not expect an easy solution to alleviate the computation burden here, because the problem is inherently difficult. In principle, one could parallelize the matrix inversion (which is frequently carried out in the solver) in the ODE solver to speed it up, however the effectiveness of doing this remains to be tested.