© ESO, 2011

1. Introduction

Astrochemistry is a field of increasing scientific attention. Over the last three decades, a lot of effort went into understanding how interstellar chemical species are formed, destroyed, and how interstellar chemistry works in detail. For a review of recent astrochemical progress see for instance Lis et al. (2005); Herbst (2006); Lemaire & Combes (2007); van der Tak (2008). To model interstellar chemistry, reliable laboratory data on the respective reaction rate coefficients are of utmost importance. A few comprehensive databases of astrochemical reactions are available today, such as UDFA¹ (Woodall et al. 2007), the Ohio database OSU² (e.g., Garrod et al. 2008), and NIST Chemistry Webbook³. There are also efforts to pool all available reaction data into a unified database (KIDA: Kinetic Database for Astrochemistry)⁴ (Wakelam 2009). Results from astrochemical laboratory experiments are collected in these databases in an indispensable effort to make these results known to the wider community. They represent the interface between chemistry and astrophysical modeling. In the following we will mainly concentrate on UDFA06, which contains parametrized reaction rate coefficients for more than 4500 reactions, involving 13 different elements and 420 different chemical species – far more than have actually been observed in interstellar space.

Fig. 1 Two example reactions to demonstrate how extremely reactions might diverge from physically reasonable numbers if γ ≪ 0. The final fit is shown in red, the UDFA rates are shown in light blue. Left: the UDFA extrapolation to k(10   K) = 4.5  ×  104 compared to 2.5  ×  10-12 in our new fit. Right: the UDFA extrapolation to k(10   K) = 2.3 compared to 1.3  ×  10-10 in our new fit.

UDFA (and most other databases as well) tabulate two-body reactions together with a set of fit parameters α,β,γ in the form k(T) = α(T/300   K)^βexp(− γ/T) and a temperature range at which the rate coefficients hold. For some reactions, more than one entry with non-overlapping temperature ranges and respective fit parameters is present because of multiple sets of experimental data. In the given fit form, γ is equivalent to the reaction’s activation energy, the energy barrier that has to be overcome for the reaction to start. When fitting the measured reaction rates over a particular range of temperatures, a negative γ may sometimes give the best fit. Special care has to be taken when using these rate coefficients outside their given range. Especially at very low temperatures, which are relevant for many astrophysical situations, negative γ may lead to unphysically high reaction rates.

To prevent unreasonable results, Woodall et al. (2007) give a list of reactions whose rate coefficients should be set to zero at 10 K. This may be feasible in chemical models, where the gas temperature is a constant, given parameter. However, in the framework of, e.g., numerical models of photon-dominated regions (PDRs), where the temperature is determined from detailed energy balance and may cover a wide range (~5−5000 K), it is preferable to use reaction rate coefficients, which do not show discontinuities across the applied temperature range (see, e.g., Hollenbach & Tielens 1999; Röllig et al. 2007).

In Fig. 1 we show two example reactions that have a diverging reaction rate coefficient at low temperatures because of γ ≪ 0. Using these extrapolated values in a chemical network calculation may lead to artificial and very unrealistic results!

Fig. 2 Relative sulfur abundance as function of AV for a spherical cloud. The model parameters are mass: M = 1   M⊙, surface density n = 106 cm-3, and FUV intensity χ = 106 in units of Draine. Three different chemical networks were applied: a) the original UDFA rates; b) the original UDFA rates with all reactions possessing γ < 0 set to zero; and c) all reactions with γ < 0 replaced by the new fits from Table 1.

To illustrate the problem, we compare the results of three PDR model calculations using the KOSMA-τ PDR model (Röllig et al. 2006). We calculated spherical model clouds with mass M = 1 M_⊙, surface gas density n = 10⁶ cm^-3, and FUV illumination χ_Draine = 10⁶. The chemical network consists of 198 species and 3237 reactions. Except for the surface formation of molecular hydrogen, a pure gas-phase chemistry was calculated. The only difference between the models is the applied set of reaction rate coefficients. In the first model we used the original parametrization of the reaction rate coefficients⁵. In the second model we deactivated these reactions by setting their reaction rates to zero. In a third model we used our new fit to the reaction rate coefficients from Table 1. In Fig. 2 we show the relative sulfur abundance versus A_V for all three cases. For A_V > 7 the models start to deviate significantly. The relevant reaction for this comparison is (1)The original UDFA rates are valid for T > 200 K (see Fig. 1, left). The model temperature goes down to 40 K at A_V = 3.5. At this temperature, the UDFA reaction rate coefficient is a factor 670 higher than at 200 K, leading to a very efficient destruction of S. At A_V = 20 reaction 1 is three orders of magnitudes stronger than the next important destruction reaction. Switching off this reaction allows the atomic sulfur to remain in the gas up to large A_V. Our new rate coefficients give a result between these two extremes. In our comparison, the choice between the scenarios a) − c) strongly affects the whole sulfur chemistry at large A_V.

The UDFA lists the accuracies of the given fits, usually ranging between 25−50%. However, outside their given temperature ranges, negative γ may lead to errors many orders of magnitude larger. For these reactions we calculated new fit parameters to allow for continuous rate coefficients even when T becomes very low. Furthermore, we calculated a single set of fit parameters for those reactions that have multiple entries in the database. These multiple sets may or may not have overlapping temperature ranges. Our final fit is designed to be valid over the total temperature interval with only small deviations.

It is well possible that the results obtained with our new fits have large errors, especially outside the original temperature regime. The peculiar shape of the fit in Fig. 1 below T_min is a direct result of the choice of the parametrization function in UDFA and it is unlikely that the low temperature behavior of this reaction follows that particular shape. However, the divergent low temperature rates of the original UDFA parametrization are unphysical, and the new fit is our best attempt at improving this behavior while minimizing any new uncertainties.

2. Fit procedure

Below we describe the applied fit strategy. Generally, we are willing to sacrifice some of the originally given accuracy to gain continuity down to very low temperatures. To be precise, we did not calculate a new fit to previously attained experimental data. We calculated a new fit to rate coefficients calculated from the original UDFA fits over the given range. The steps are the following:

• 1.

  tabulate the reaction rate coefficients over the giventemperature range, using the given fit-accuracy as error. Theseare our artificial data;

• 2.

  if γ < 0, fix k(T = 10   K);

• 3.

  calculate new fit parameters α,β,γ⁶;

• 4.

  calculate the error of the new fit with respect to the original UDFA parametrization.

In the Appendix we show plots of all new fits presented in this paper.

2.1. Discussion of the fit strategy

The details of the fit strategy may have considerable influence on the final fit result. Particularly point 2 prevents the reaction rates from becoming unphysically high. Of course, this procedure introduces new, possibly very large errors to the rate coefficients. But these errors are usually of the same order as those from the UDFA fits. We now consider some major influences on the fit results.

2.1.1. Placement of “data” points

The number and placement of “data” points has a direct influence on the quality of the fit because it implicitly introduces a weighting of the respective temperatures. Using the Levenberg-Marquardt algorithm is in effect applying a sequence of linear least-squares fits to a nonlinear function. Regions with more “data” points contribute more to the total sum of squares and force the fit function closer to these “data” points. For instance, placing temperature points equidistant in log-space will place more grid points at lower temperatures compared to equidistant linear spacing, weighting lower temperatures stronger than higher temperatures. Optimal fit results could be achieved by adaptively changing the number and positions of grid points until the error is minimal, but to keep a consistent error estimate, we used the same temperature gridding for all fits. To fit reactions with γ < 0, we employed an equidistant temperature grid in log-space with 20 grid points, while we used linear spacing for fitting reactions with multiple entries, using 20 grid points per temperature interval.

2.1.2. Fixing k(T = 10 K)

Choosing a fixed k(T = 10   K) prevents the reaction rates from becoming unphysically high. However, the exact choice of k(T = 10   K) is completely arbitrary since no real measurement is available. Moreover, it poses a major constraint to the fit and an unfortunate choice will lead to very large deviations from the true reaction rate coefficient. Woodall et al. (2007) suggest to set k(10   K) = 0. It is possible that lim_T → 10k(T) = 0. It is also possible that k(10   K) ≫ 0, we simply do not know. To minimize the deviations from the only available “data”, we chose to set k(10   K) = k(T_min) if γ < 0 with T_min the minimal valid temperature from the UDFA parametrization. This also reduces any over/under-swinging behavior of the new fit. To illustrate how the choice of k(10   K) influences the form of the final fit function we compare in Fig. 3 the fit for k(10   K) = k(T_min) (red, solid line) and for k(10   K) = 10^-50 (red, dashed line). In this particular case, the latter leads to significantly stronger deviations from the UDFA parametrization than for k(10   K) = k(T_min).

2.1.3. Calculating new fit parameters

Fig. 3 Comparison of two fit results for different choices of k(T = 10   K). The UDFA rate coefficient, valid between 295 and 8000 K, is shown in light blue. The best fit for k(10   K) = 10-50 is shown as red, dashed line (α = 1.01(− 11),β = 0.23,γ = 913). The best fit for k(10   K) = k(Tmin) is shown as red, solid line (α = 9.95(− 12), β = 1.05, γ = −42.5).

We fitted the fit function to the list of “data” points created in steps 1 and 2. In most cases it is better to perform the fit in log-space, i.e. use log (k(T)) instead of k(T) as fit function. Otherwise, the Levenberg-Marquardt fit would have the tendency to fit k(T) to absolutely higher temperatures, while we are interested in the lower temperature regime. In some cases however, a linear fit function is more appropriate. This case is shown in Fig. 4. We try to find a simultaneous fit across the multiple temperature regimes of the reaction C + N  →  CN. The linear fit produces a continuously declining result for T → 0, while a logarithmic fit function produces the opposite behavior, which results in large deviations of the 10 K rate predicted by the UMIST parametrization and our fit result. Again, the true reaction rate behavior at very low temperatures is unknown, hence, our aim is to minimize any deviation from the UDFA fit. The logarithmic fit differs from the 10 K UDFA rate by 8 orders of magnitude, while the linear fit differs only by a factor of 40. For this reaction we chose a linear fit function, for all other fits in this paper a logarithmic fit function was used.

2.1.4. Error calculation

Strictly speaking, the term error is misleading in our context, since it suggest model deviations from real measurements⁷. Instead, whenever we are speaking of fit errors, we actually compare how two different models deviate from each other. Nevertheless, because our aim is to change the original parametrization as little as possible, we will use this deviation to quantify the quality of our fits.

Fig. 4 Comparison of fit results for a linear fit function k(T), shown as solid, red line, and a logarithmic fit function log (k(T)), shown as dashed, red line.

In Tables 1 and 2 we give the relative error brackets of our new fit, i.e.  [min(e_i),max(e_i)] , with for all T_i ∈  [T_min,T_max] . Here, k^U is the UDFA parametrization of the reaction, and k^f is our new fit.

2.2. Fitting reactions with negative γ

Fig. 5 Example reaction that can be fitted very well (maximum additional error is 11%), but has a less extreme low temperature behavior. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue. The difference in k(10 K) is 10 orders of magnitudes between the UDFA rates and our new fitted rates.

To obtain a new fit for reactions with γ ≪ 0, we used a log-spaced temperature grid with 20 grid points and a logarithmic fit function. Collider reactions, i.e. reactions of type “CL” in UDFA, were excluded from the fit. Apparently it is possible to find a very good fit for some reactions without any divergence at low temperatures and without introducing any significant additional error. The prime example for this is shown in Fig. 5. Using the UDFA reaction rates gives an unreasonably high reaction rate coefficient at 10 K of 5! The new fitted rate coefficient at 10 K is 2 × 10^-10. Even though we do not know the real reaction strength at 10 K, we consider our new result to be much more reasonable compared to the extrapolated numbers from UDFA.

Some of our new fits even have negative γ, which shows that γ < 0 not always has to lead to divergent behavior. If γ ≳ −100, the rate inclination for low temperatures is moderate and should not lead to unreasonable results⁸. In Table 1 we present the fit results for reactions that are listed in UDFA with γ < 0. We give the original and the new fit parameters as well as the relative error brackets of the new fit.

Some of the reactions given in Table 1 also appear in Table 2 because they have counterparts with different temperature regimes. Some of the numerical problems described above can be circumvented by choosing an alternative reaction rate coefficient with a more appropriate temperature regime. Nevertheless, we list these reactions for cases where that is not possible.

2.3. Fitting reactions with multiple entries

The fits to reactions with multiple database entries are, in most cases, much less sensitive to the details of the fit strategy, but, overall, the errors of the simultaneous fits are larger compared to the fits to reactions with γ < 0. Here, the partly very inconsistent entries in UDFA for the different temperature ranges contribute the most to the errors in the final, simultaneous fit. For cases that span a very large temperature range, the choice of T_max when setting up the artificial “data” points poses a major constraint on the fit result. This is demonstrated in Fig. 6.

In Table 2 we present the fit results for reactions with multiple entries in UDFA. The valid temperature ranges of these reaction rates may or may not overlap. Even worse, some of the multiple entries appear to be mutually exclusive. Especially the high temperature rate coefficients for the reaction He⁺+CO  →  O + C⁺ + He appear to be at odds with both theory and experiment and it may be advisable to only use the low temperature rate coefficients (see Fig. 7). However, this assessment is far beyond the scope of this paper. We apply our fit routine to all reactions in UDFA with multiple entries and leave it to the informed reader to decide whether to use the simultaneous fit or whether to select a reaction with an appropriate temperature regime.

Fig. 6 UDFA lists three C+ recombination reactions for temperatures from 10 to 41 000 K. The two reaction rates for the higher temperatures show a very different behavior compared to the reaction that is valid from 10 to 8000 K. Applying different values of Tmax during the fit leads to significantly different results.

Fig. 7 Example reaction with apparently mutually exclusive reaction rates listed in UDFA. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue. The two temperature ranges overlap and the respective fits to the reaction rates, here shown in light blue, display a very different behavior.

For the fitting, we limited the upper temperature range to 20 000 K where T_max is higher, thus neglecting higher temperatures. The errors are explicitly calculated for this range and may be significantly larger for higher temperatures. Despite the apparently large errors we consider the fits relatively good in most cases. Remember that we fit multiple, sometimes very disparate sets of “data” points, and therefore, large errors are unavoidable. This is demonstrated in Fig. 8 for two cases from Table 2.

The reader should also recall that the reaction rates from UDFA usually have errors on the order of 25−200%. One could be tempted not to use the UDFA rates to prevent the errors to rise even more because of our simultaneous fitting. However, by doing so one might implicitly accept much larger errors when moving from one temperature regime to the other. This is demonstrated in Fig. 9.

The simultaneously fitted reaction rates in Fig. 9 have a maximum relative error of 32%. This may appear large compared to the individual errors for each temperature regime. But when only using either one of the two reactions we would implicitly accept a very large error when the temperature is around 400−500 K, where both individual rates disagree by one to two orders of magnitude. Our simultaneous fit is a reasonable fit in an attempt to keep the experimental findings for both temperature regimes.

The special case of carbon recombination is shown in Fig. 6. Three different sets of fit parameters are given in UDFA for three distinct temperature ranges up to 41 000 K. Especially the high temperature rates show a very different behavior compared to the low temperature range. A simultaneous fit gives very large errors because of the large deviations from the high temperature rates. However, for most applications, temperatures up to 8000 K are sufficient and we can neglect the errors for high temperatures. In that case, the new fit has much smaller errors. For many practical purposes, one could even completely discard the two higher temperature ranges and only use the low temperature rates. For the sake of completeness, we still list the new fit temperatures up to 20 000 K.

2.4. Discussion

For three reactions in Tables 1 and 2 one finds updated measurements and theoretical values of the reaction rate coefficients in KIDA k^K(T). A detailed comparison of the respective reaction rate coefficient behavior is beyond the scope of this paper. We will only take a look at the updated recommended 10 K values for these reactions and compare them to our new values k^f(10   K) in Table 3.

Fig. 8 Two example reactions to demonstrate the wide range of error values. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue. Left: this reaction is listed with two very differently behaving reaction rates. Consequently, the simultaneous fit has relatively large errors especially at the interface between bow temperature regimes. Still, we consider the fit to be very good because of its overall capability to reproduce the general behavior in both regimes. Right: this reaction is listed with three distinct temperature regimes. The simultaneous fit is nevertheless a very good approximation to all three UDFA rates in the relevant temperature ranges.

For all three reactions, the 10 K rate coefficients achieved from our new fits are within a factor of two compared to the recommendations from the KIDA Kinetic Database for Astrochemistry. For the reaction O + OH, our value is only 40% higher than the recommended value. Given the overall uncertainties, the agreement between our new fits and the KIDA values is remarkable and an indication that our general fit approach is viable to achieve reasonable results.

3. Summary and conclusions

We presented an approach to overcome some numerical and practical difficulties in present-day databases of chemical reaction rates. One of these difficulties is the existence of fits to reaction rate coefficient that unphysically diverge at very low temperatures. This is especially problematic for numerical codes, which calculate chemical networks down to very low temperatures and thus have to use many of the reaction rate coefficients that are listed in chemical databases beyond their given temperature regime. We demonstrate that it is possible to find a new parametrization of these reactions that ensures a less erroneous behavior and allows calculations to temperatures of 10 K and lower.

A second complication that occurs in chemical databases is that some of the listed reactions have multiple entries because the experimental data were obtained for different temperature ranges. Any user is instructed to use the listed reaction rate coefficients only in their valid temperature regime. Unfortunately, multiple entries in these databases very often display fairly inconsistent properties. This is particularly problematic for numerical codes that continuously change the temperatures when calculating chemical networks because it introduces artificial numerical jumps in the coefficient matrix of the equation system to be solved. One solution could be to restrict oneself to only use the one reaction rate that has a temperature regime closest to the computational regime. However, in many cases this is not possible because the parameter space relevant for application spans multiple temperature regimes. Additionally, this means that one implicitely accepts quite large errors when leaving the valid temperature regime. We suggest to replace these multiple entries with a single, new rate which is fitted simultaneously to the individual ones. This minimizes the error when moving from one temperature regime to the other and additionally overcomes the problem of discontinuous reaction rate jumps.

We consider the new parametrization presented in this paper only a temporary solution until better experimental data for low temperature reaction rates are available.

Fig. 9 Simultaneously fitted reaction rates in this figure have a maximum relative error of 33%. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue.

As a side effect we demonstrate that in some cases one might find an equally good parametrization of the reaction rate coefficients without any artificial numerical side effects just by taking care not to allow too low negative values of γ. The chemical databases try to express the possibly very complicated physics of the various chemical reactions with only three fitting parameters. Of course, this is an extremely strong constraint on any parameter fit to the complex experimental data. Although, under certain circumstances, chemical reactions listed in chemical databases like UDFA, OSU, and others might reasonably be applied even outside their given scope, it is important to ensure that numerical side effects from parametrization are minimized when possible. For many numerical applications this is much more important than a smaller total error of the reaction rate coefficients in a limited range in temperature.

Any user should always take care not to blindly use reaction rate coefficients from chemical databases outside their

temperature scope. This is also true for the new parametrizations listed in this paper. In each case, one has to balance advantages against disadvantages. If you are sure that you will keep within the given temperature regime, then there is no need to use our rate coefficients because they add an additional error to the computation. Under certain circumstances, this additional error might be negligible compared to the numerical problems that come with the original rate coefficients.

Acknowledgments

This work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Os 177/1–1. We acknowledge the use of UDFA (a.k.a. UMIST) (http://www.udfa.net/) chemical reaction databases. Some kinetic data we used were downloaded from the online database KIDA (Kinetic Database for Astrochemistry, http://kida.obs.u-bordeaux1.fr). The author thanks R. Simon and V. Ossenkopf for their helpful comments on the text. The author thanks the anonymous referee for his or her comments, which helped the final text a lot.

References

Appendix A: Plots of all new fits

Here, we show the plots of all new fits to allow the user an assessment of the fit quality over the entire fitted temperature range. The plots are shown in the same order as in Tables 1 and 2.

Fig. A.1 Fits to reactions with negative γ.

Fig. A.1 continued.

Fig. A.2 Fits to reactions with multiple entries.

Fig. A.2 continued.

All Tables

All Figures

	Fig. 1 Two example reactions to demonstrate how extremely reactions might diverge from physically reasonable numbers if γ ≪ 0. The final fit is shown in red, the UDFA rates are shown in light blue. Left: the UDFA extrapolation to k(10   K) = 4.5  ×  104 compared to 2.5  ×  10-12 in our new fit. Right: the UDFA extrapolation to k(10   K) = 2.3 compared to 1.3  ×  10-10 in our new fit.
In the text

Fig. 2 Relative sulfur abundance as function of AV for a spherical cloud. The model parameters are mass: M = 1   M⊙, surface density n = 106 cm-3, and FUV intensity χ = 106 in units of Draine. Three different chemical networks were applied: a) the original UDFA rates; b) the original UDFA rates with all reactions possessing γ < 0 set to zero; and c) all reactions with γ < 0 replaced by the new fits from Table 1.

In the text

	Fig. 3 Comparison of two fit results for different choices of k(T = 10   K). The UDFA rate coefficient, valid between 295 and 8000 K, is shown in light blue. The best fit for k(10   K) = 10-50 is shown as red, dashed line (α = 1.01(− 11),β = 0.23,γ = 913). The best fit for k(10   K) = k(Tmin) is shown as red, solid line (α = 9.95(− 12), β = 1.05, γ = −42.5).
In the text

	Fig. 4 Comparison of fit results for a linear fit function k(T), shown as solid, red line, and a logarithmic fit function log (k(T)), shown as dashed, red line.
In the text

	Fig. 5 Example reaction that can be fitted very well (maximum additional error is 11%), but has a less extreme low temperature behavior. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue. The difference in k(10 K) is 10 orders of magnitudes between the UDFA rates and our new fitted rates.
In the text

	Fig. 6 UDFA lists three C+ recombination reactions for temperatures from 10 to 41 000 K. The two reaction rates for the higher temperatures show a very different behavior compared to the reaction that is valid from 10 to 8000 K. Applying different values of Tmax during the fit leads to significantly different results.
In the text

	Fig. 7 Example reaction with apparently mutually exclusive reaction rates listed in UDFA. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue. The two temperature ranges overlap and the respective fits to the reaction rates, here shown in light blue, display a very different behavior.
In the text

Fig. 8 Two example reactions to demonstrate the wide range of error values. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue. Left: this reaction is listed with two very differently behaving reaction rates. Consequently, the simultaneous fit has relatively large errors especially at the interface between bow temperature regimes. Still, we consider the fit to be very good because of its overall capability to reproduce the general behavior in both regimes. Right: this reaction is listed with three distinct temperature regimes. The simultaneous fit is nevertheless a very good approximation to all three UDFA rates in the relevant temperature ranges.

In the text

	Fig. 9 Simultaneously fitted reaction rates in this figure have a maximum relative error of 33%. The final fit is shown in red, the UDFA rates for each individual temperature regime are shown in light blue.
In the text

	Fig. A.1 Fits to reactions with negative γ.
In the text

	Fig. A.1 continued.
In the text

	Fig. A.2 Fits to reactions with multiple entries.
In the text

	Fig. A.2 continued.
In the text