Bayesian Inference of Reaction Rates in Icy Mantles

Grain surface chemistry and its treatment in gas-grain chemical models is an area of large uncertainty. Whilst laboratory experiments are making progress, there is still much that is unknown about grain surface chemistry. Further, the results and parameters produced by experiment are often not easily translated to the rate equation approach most commonly used in astrochemical modelling. It is possible that statistical methods can reduce the uncertainty in grain surface chemical networks. In this work, a simple model of grain surface chemistry in a molecular cloud is developed and a Bayesian inference of the reactions rates is performed through MCMC sampling. Using observational data of the solid state abundances of major chemical species in molecular clouds, the posterior distributions for the rates of seven reactions producing CO, CO$_2$, CH$_3$OH and H$_2$O are calculated, in a form that is suitable for rate equation models. This represents a vital first step in the development of a method to infer reaction rates from observations of chemical abundances in astrophysical environments.


INTRODUCTION
Dust grain chemistry plays an important role in the physical processes happening deep inside dark molecular clouds during star formation (Draine 2003;Williams & Cecchi-Pestellini 2015). These dust grains are vital to every part of the star formation process and ultimately contribute to the basic matter from which icy planetesimals are formed ( van Dishoeck 2004). It is in fact evident that molecules such as water and methanol in dust grain ice mantles are primarily formed through solid state chemistry rather than accreted directly from the gas phase (Parise et al. 2005;Ceccarelli et al. 2007). In recent years, even more complex molecules have been observed in both prestellar cores and star forming regions (see reviews by Herbst & Van Dishoeck 2009;Caselli & Ceccarelli 2012), some of which can not currently be explained by pure gas phase chemistry. Therefore, chemical reactions leading to simple as well as complex molecules must occur on the surface of icy dust grains.
Experimentally, it has been known since the work of Hagen et al. (1979) and Pirronello et al. (1982) that grains can be chemical nanofactories on which surface reactions, UV photons and cosmic rays radiation can synthesize complex molecules and even prebiotic species, starting from simple atoms or molecules such as H, C, O, N, CO. Therefore, understanding the surface chemistry that takes place on dust grains is key to understanding not only the origins of stars, but also how rocky and gaseous planets are formed.
Initially, surface reaction networks in chemical models were based on chemical intuition and gas phase chemistry analogues. However, over the past two decades, laboratory astrochemists have been using experimental techniques to test and evaluate surface reactions. As a result, the efficiencies of reaction routes are being properly explored and important information on how molecules form on grain surfaces is being revealed (see Williams & Cecchi-Pestellini 2015, for a review). The first experimental work on the dust surfaces studied the formation of molecular hydrogen (Pirronello et al. 1997). Several more experiments followed studying either the formation of more complex molecules (e.g. Watanabe et al. 2005;Ioppolo et al. 2009;Minissale et al. 2015) or the ice morphology and ice mantle mechanisms (e.g. Fraser et al. 2004;Collings & McCoustra 2006). Surface reactions can be experimentally investigated within a constrained range of laboratory conditions. Typically, these conditions include different atomic fluxes, ice temperatures, ice morphologies, and mixture ratios, as well different energetic processes. The aims of the experiments are to investigate surface molecule formation, desorption and diffusion. However, experimental data for interstellar ices is limited, the main reason being that the experimentation process is neither simple nor fast. In order to make the best use of experimental resources, the chemical data that models require needs to be prioritized according to what will have most impact.
Bayesian methods are widely used in astronomy as a means of deriving posterior probability distributions for model parameters from observations (eg. Palau et al. 2014;Schmalzl et al. 2014;Bevan 2018). It is the de facto standard in the field of cosmology but is becoming more and more widely used in other areas of astronomical research (eg. Lomax et al. 2016;Testi et al. 2016). In astrochemistry, Bayesian inference has been used to derive parameters such as the gas density and cosmic ray ionization rates within a dark molecular cloud from observations of species in the gas and ices using chemical models (Makrymallis & Viti 2014).
In this work, we investigate the chemistry itself, studying the rates of reactions on the dust grain surfaces in an attempt to infer their rates and provide a list of reactions for which an accurate rate is particularly important. This is the first such work in an astrochemical context but Bayesian methods have been used to determine rate coefficients for combustion chemistry on Earth (Prager et al. 2013). This work represents a necessary first step in which we determine whether we can use a reduced chemical model and very simple observational constraints to learn more about the parameters in a grain surface chemical model. The grain chemistry model used in this work is described in Section 2. The inference process including the choice of MCMC sampler is presented in Section 3. The results of our analysis are presented in Section 4 along with additional discussion in Section 5. Finally, our conclusions are discussed in Section 6.

THE CHEMICAL MODEL
A simple chemical model was developed that considers only the solid state chemistry in the ice mantles of dust grains in a dark molecular cloud. The simplified model is a time-dependent single-point model that generates a time series of solid phase molecular abundances as a function of the physical conditions of the molecular cloud and the chemical parameters of the defined chemical network. The chemical network consists of 23 species and 24 surface reactions that are listed in Table 1 and  Table 2 respectively.
To model the surface chemistry of a dark cloud the abundance of each solid species is derived by solving rate equations for grain-surface chemistry. The formation and destruction mechanisms for a species i are given by the following kinetic equation: (1) where k i lm is the reaction rate of all the reactions between species l and m that produce i, n i is the concentration of species i (with the subscript gas indicating the concentration of the species in the gas-phase), k r represents the reaction rates of all the reactions where species i participates as reactant, while k des i and k ads i are the desorption and adsorption rates.
The reactions in Table 2 are mainly hydrogenation reactions of common gas phase species and reactions between species that are likely to be abundant on the grains. Where possible, reactions that have been found to be efficient, or even dominant, routes to forming a species have been chosen. For example, the hydrogenation of CO to form CH 3 OH is well studied (Fuchs et al. 2009;Chuang et al. 2016) and so reactions 21-24 are the only considered route to form CH 3 OH. Similarly, the formation of CO 2 via reaction 3 is known to be efficient (Ioppolo et al. 2011) and other routes suffer from large energetic barriers. Note that in cases where H or H 2 is a product of a reaction such as in reaction 11, it is ignored and the total H abundance is not conserved. This is for simplicity and the lost H represents too small a fraction of the H abundance to affect the model.
There is no gas phase chemistry in the model and so the freeze out of species from the gas phase must be parameterized. The adsorption rate is assumed to be zero for all but the following six species: CO, CS, O, H, OH and S. To derive the adsorption rate of these species, the gas-grain chemical code uclchem (Holdship et al. 2017) was utilized. uclchem was run with a network of 220 species with gas-phase reactions from UMIST12 (McElroy et al. 2013), freeze out of gas phase species and the non-thermal desorption of grain surface species. A single point model of this full gas-grain chemistry was run in which the gas increased in density under freefall from 10 2 cm −3 to 2 × 10 4 cm −3 , which is appropriate for a dark molecular cloud. The chemistry progresses over 10 Myr at 10 K and the freeze out rates for the six species above were extracted from this model.
The freeze out rates from uclchem were inserted as source terms in the ODEs for those species in the simple grain surface model. The grain surface models starts with an abundance of zero for all species, representing bare grains. The model then progresses for 10 Myr considering only the 24 reactions in Table 2, the freeze out rates and the non-thermal desorption of each species. In this way, the grain surfaces in a dark molecular cloud are effectively modelled whilst the computation time is low  as the gas phase treatment is reduced to the incoming (freeze out) and outgoing (desorption) flux of molecules. Note that the cloud age is arbitrary, the model reaches the molecular cloud density at 6 Myr and the chemistry is then allowed to progress until a total age of 10 Myr. The exact choice of final time has only a small affect on the abundances in the model.
Whether one parameterizes the rate of surface reactions in a similar way to the Kooji-Arrhenius equation used for gas phase chemistry (Occhiogrosso et al. 2012) or considers the diffusion and reaction of species across the ice surface (Hasegawa et al. 1992;Chang et al. 2007), the rate of a reaction is constant for a given temperature and dust composition. Therefore, k in this model is treated as a constant rate of reaction in units of cm 3 s −1 as the temperature in the model is constant at 10 K. This reduces the number of parameters in the model and reflects the available data, i.e. ice phase abundances in quiescent, approximately isothermal clouds.
The result of these approximations and modifications is a model of the grain surface chemistry under the conditions of a dark molecular cloud at a constant temperature of 10 K. The freeze out rates and gas density are specific to a dark cloud and so the model is not of applicable to arbitrary ices. This model has a run time that is approximately 1000 times shorter than an equivalent run of uclchem. This reduction in run time is vital due to the number of model runs required for an MCMC inference procedure.

Inference Procedure
The aim of this work is to obtain information about the set of reaction rates k = (k 1 , k 2 , ..., k 23 ) of the surface chemical network, where k j is the reaction rate of reaction j. For a given set of rates, the model produces simulated molecular abundances where Y i is the abundance of species i. These quantities are related through the chemical code C(·), so that Y = C(k).
For any set of simulated abundances,the probability of the corresponding parameter values can be evaluated through the use of Bayes' rule, where d is the data, representing a set of observational constraints on Y. P (k|d) is the posterior probability distribution (PPD) of k which expresses the level of certainty about the reaction rates after considering the data and any prior information. The denominator is known as the Bayesian evidence but for the purposes of this study can simply be treated as a normalization factor. The prior probability distribution (P (k)) adopted for the reaction rates is a logarithmically uniform distribution that is non-zero when the reaction rates are between 10 −5 and 10 −30 and zero elsewhere. The limits of the prior distributions represent a larger range than that of rates typical of gas phase reactions. This is a choice to reflect the exploratory nature of the work and is expected to cover all likely rates regardless of the nature of surface reactions. The likelihood function P (d|k) must give the likelihood of having obtained the data given the assumed set of rates. Here the likelihood encodes measurement noise and is given as, where Y i are the model abundances of each species for which there is data and σ i is the Gaussian uncertainty of each observed fractional abundance.

Implementation and Data
In order to constrain the reaction rates, data is required in the form of species abundances in the ices. The solid state fractional abundances of species in quiescent gas illuminated by background stars were taken from a comprehensive recent review (Boogert et al. 2015, Table 2). Of the species in the grain surface model, this review gives constraints on the abundance of H 2 O, CO, CO 2 and CH 3 OH. The value of the median fractional abundance is provided for each along with the upper and lower quartile values. To formulate the likelihood, it is assumed these values describe a Gaussian distribution for the abundance of each species. It is also assumed that the uncertainties on the abundances are independent which is likely given that the abundances and statistics presented by Boogert et al. (2015) are combinations of different data sets for each species. The median value is taken as the mean and the upper and lower quartile values can then be assumed to be 0.68σ from that mean. The resulting abundances and uncertainties are listed in the upper half of Table 3. Due to the low number of observations, these distributions are not perfect representations of the data as the quartile values are not precisely symmetric about the median.
In order to evaluate the posterior distribution function for all values of k, a sampler must be used. The emcee python package (Foreman-Mackey et al. 2013) was chosen for this purpose. This in an implementation of the affine-invariant monte carlo sampler proposed by Goodman & Weare (2010). Rates were sampled by 128 "walkers", each producing chains of ∼10 6 samples where the frequency of the appearance a particular rate value in the chain is proportional to its likelihood. These walkers start from random positions in rate space (ie all 24 rates have a random value from 10 −30 cm 3 s −1 to 10 −5 cm 3 s −1 ).The sampling took approximately 100 hours using a single node on the DiRAC CSD3 plat- <4.0 × 10 −6 Table 3. Upper section: adopted abundances and uncertainties of species observed in the ices used as data in the parameter inference. Lower section: upper limits of the fractional abundance for other species which are used to further constrain the reactions rates. All values adapted from Boogert et al. (2015) as discussed in the text.
form's Skylake-Peta4 system utilizing emcee's built in MPI tools. In appendix A, heuristics are presented which demonstrated the chains are likely to have converged.

Upper Limits
The abundances of H 2 O, CO, CO 2 and CH 3 OH are the only strong constraints on the abundances of the species in this model. The reaction rates that are acquired as a result of performing a Bayesian parameter inference procedure with these abundances are presented and discussed in Section 4.1. However, weaker constraints do exist for other species. Boogert et al. (2015) provide upper limits on the abundances of OCS in dark clouds as well as upper limits on H 2 CO, SO 2 and H 2 S in other objects, the upper limits used are given in Table 3. Therefore, a second parameter inference procedure was performed which was identical to the first except that the likelihood function was modified to take into account the upper limits. In order to be conservative, when deriving upper limits on species that have not been detected towards background stars, the upper limits towards YSOs were increased by an order of magnitude to allow for larger abundances in molecular clouds. This is the case for the upper limits on H 2 CO and SO 2 .
To account for these upper limits, modifications were made to Equation 3. When including the upper limits, the likelihood of a model was calculated as, where δ i is 1 for observed species and 0 for species with upper limits. C is the upper-limit of a species' abundance and S(C i ) is the survival function. This modification to the likelihood is standard for left-censored data, i.e. ones where a detection limit provides only an upper limit on a quantity (Klein & Moeschberger 2003). The survival function for a Gaussian distribution is used, where erf () is the error function. σ i is assumed to be a third of the value of the upper limit. Equation 4 is equivalent to Equation 3 for detected species. However, for upper limits it takes the value of 1 − S. For model abundances much less than the upper limit, this likelihood is equal to 1 and for model abundances much larger than the limit it takes a value of zero. Thus, whilst model abundances close to the upper limit are accepted, models with much larger predicted abundances have a likelihood of zero.

Testing the Method
In order to ascertain whether this method would be able to predict reactions rates from measured abundances in the case that the model was an accurate representation of reality, it was tested using abundances obtained from the model itself. First, random rates for all 24 reactions were generated and the model was run using those rates. The abundances of the four species in Table 3 for which observed abundances are available were then stored. "Noisy" abundances were then generated by drawing randomly from a Gaussian distribution with a mean value of the model abundances and a standard deviation set to a 50% error. This produced four "observations" obtained from the model with known rates. The 50% error was chosen as it is the approximate fractional error of the real observations. The MCMC procedure was then performed to see whether the known reaction rates could be obtained.
It was found that the majority of rates could not be recovered. However, in tests where the rates for reactions that produced H 2 O, CO, CO 2 and CH 3 OH were high enough to produce observable abundances, the rates of those reactions were recovered. That is to say that intrinsically, this method appears to only be able to give information on the rates of reactions that form the species for which observational constraints are available.

Gaussian Abundance Constraints
The results are presented in the form of marginalized PPDs for the reaction rate coefficients. The density  of each marginalized PPD reveals the areas where the corresponding reaction rate is more probable based on the imposed constraints. In Figure 1, the marginalized PPDs of selected reactions are plotted. These show a large probability density only for a much smaller range of values than the prior. As expected from the tests in Section 3.4, these are generally reactions that form the constrained species. The PPDs are shown as histograms and Gaussian kernel density estimates, using the full MCMC chains from all 128 walkers. It is believed that these chains have converged and the relevant tests are discussed in Appendix A. The PPD of reactions 1, 2, 3 and 21 are well constrained and involve species directly constrained by observation. Reaction 1 provides OH required to form H 2 O and CO 2 through reactions 2 and 3. Those are in turn constrained by the observed abundances of CO 2 and H 2 O and mutual competition for OH. Reactions 21 uses up CO and so it is expected there would be an upper limit due to the observed abundance of CO and a lower limit due to competition with reaction 3. Reactions 22-24 form CH 3 OH from HCO. The competition between reactions and the correlation between the rates of reactions 21 through 24 is explored further in Section 5.2.
The other PPDs are broadly similar to the prior distributions and the implications of this should be stated. Essentially, the reactions where the rates have uniform probability distributions are reactions that do not impact the likelihood of the model. It should be noted, however, that changes in the abundance of species not included in the likelihood calculation are possible. PPDs that are similar to the priors indicate that when modelling only H 2 O, CO, CO 2 and CH 3 OH, the rates of those reactions are unimportant.

Inclusion of Upper Limits
The PPDs of each reaction are largely unchanged when upper limits are included indicating that the upper limits may in fact be too conservative. The only major change is that the rate of reaction 10 takes a minimum value of 10 −17 cm 3 s −1 . This is required for a models to produce a lower OCS abundance than the upper limit.
One may expect the upper limit on H 2 CO to improve the level of certainty in the reaction series 21-24. However, using the conservative value in Table 3, no change is seen in the posteriors. If the value for YSOs is utilized instead, only a small change is observed. In that case, the peaks that were apparent in the marginalized posterior of each reaction in Figure 1 become more pronounced, such that the majority of the probability density lies within them. Thus the value of H 2 CO would represent an important constraint if an appropriate value for molecular clouds could be obtained.
The most likely rates for the well constrained reactions from this MCMC analysis are presented in Table 4. These most likely rates are the same whether they are taken from the observed abundances MCMC chains or the chains that included upper limits. It is clear that the available upper limits are not sufficiently constraining to improve the parameter inference.

Model Abundances
The results of the MCMC run give the marginalized posterior distributions of the rates of each reaction. In order to understand how well these rates reproduce the observed abundances, the model must be run with rates drawn from the probability distributions. This will allow the uncertainty in the model that arises from the uncertain rates to be quantified. The model was run 1000 times, with the rates of the reactions randomly sampled from the marginalized posterior distributions derived from the upper limit MCMC procedure. By plotting the average abundance of each species and the standard deviation of those abundances, the uncertainty in the model can be evaluated.
The abundances for selected species from the model runs are plotted in Figures 2 and 3. It can immediately be seen from Figure 2 that, for the species with observational constraints, the uncertainty in the rates does not lead to a large uncertainty in the model abundances. The model appears to consistently underproduce CH 3 OH and overproduce H 2 O. However, the difference is small if the errors on the observations are accounted for. This is a good result for such a simple model. It may be that if the reduced network was expanded, these results would be improved. Equally, it  may be that the constraints are broad enough that a poor CH 3 OH abundance is not affecting the overall likelihood as much as a poor CO or CO 2 abundance would. The fractional abundances of the species with upper limits were not as well constrained and so it might be expected that they are much more varied. In Figure 3 it can be seen that this is not the case, which implies that their abundances are also strongly tied to the rates of the reactions in Figure 1. However, the average abundance of OCS is an order of magnitude higher than the . The joint probability distributions of the rates of reactions 1 and 2, darker areas represent higher probability densities. The 1,2 and 3 sigma contours are plotted. These two reactions are tightly coupled, one must take the value of ∼5.0 × 10 −18 cm 3 s −1 and act as the rate limiting step. The other must then have this value or higher in order for enough H2O to be produced in the model. observed upper limit. Examining the full abundance distribution, it appears there is a fraction of the model runs that fall outside the 67% confidence interval and give OCS abundances below the upper limits. This illustrates the problem inherent in drawing from the marginalized posteriors. Drawing from each marginalized posterior individualy gives sets of reaction rates that break the upper limit used to infer the posterior distributions.

Network Connectivity
The PPDs presented in Figure 1 are marginalized, that is to say they represent the likelihood of a given reaction rate averaged over the values of the other rates. However, not all rates are independent and it is possible that some areas of the rate space are only likely for one reaction when a second takes a particular value. To investigate this, the joint posteriors of pairs of reactions were examined. These give the likelihood of pairs of reactions rate values so that it can be seen whether the two reactions are in some way correlated.
For example, the joint probability distribution of reactions 1 and 2 is shown in Figure 4. It can be seen that either reaction 1 or reaction 2 can take a value much higher than their respective most likely value but only Whilst a large amount of the probability density is at the location of the peaks in the marginalized posterior distributions, there is a noticeable correlation when either increases above the most likely rate. This is likely due to the fact that both reactions destroy CO and are therefore competing to produce enough CO2 and HCO respectively.
when the other is at its most likely value. This shows that, in reality, there must be certain amount of O converted to H 2 O in the model and as long as one step in that process limits the rate to the correct amount then the other can freely vary. In order to break this degeneracy, limits on the OH abundance in the ice are required. Similar joint distributions are seen for reactions 22 to 24, as a certain amount of CO must be converted to CH3OH. In Figure 5, the joint probability distribution of reactions 3 and 21 are plotted. Both reactions are less well constrained than reactions 1 or 2. This can be seen from the large area taken up by the 1σ contour. However, the high probability density areas are those where at least one reaction takes the most likely value from their respective marginalized posterior distributions. There is also a line of increased probability density where both reactions have approximately equal rates that are higher than the peak value. It is therefore likely that the reactions compete for CO and the rate of reaction 3 is poorly constrained as the availability of CO is the main factor in the amount of CO 2 produced. Tighter observational ranges on the abundances of CO, CO 2 and CH 3 OH may reduce this degeneracy and allow the rate of reaction 3 to be more clearly determined.
6. CONCLUSIONS A novel way to tackle uncertainty about surface reactions and rate coefficients using Bayesian inference was presented. To prove the efficiency of Bayesian techniques in providing insight on the chemical parameters of surface reactions, the algorithm was tested with a proof of concept example. A simple chemical code was created by parameterizing the freeze out of important species and neglecting all other adsorption to the grain surface. This left a model where only grain surface chemistry needed to be accounted for, greatly reducing the complexity and opening up the possibility of exploring a large parameter space.
The rates of the reactions in the chemical model were found through Bayesian inference. Using an MCMC sampling algorithm, the model was run with varying reaction rates and the likelihood of the model was evaluated each time. This likelihood was calculated by comparing the model to observations of ices towards background stars, which are reasonable values for a molecular cloud. It was possible to strongly constrain the rates of reactions that are involved in the production or destruction of species for which measurements exist. These rates are presented in Table 4.
Future improvements should include a more complex chemical code, including the grain surface reactions directly in a gas-grain chemical code. This would allow an improved treatment of the freeze out and non-thermal desorption amongst other effects. However, the added complexity would make this a vastly more computationally intensive procedure, initial tests with uclchem taking approximately 1000 times longer per run. Improved rates could also be achieved by including more observational data, particularly constraining species for which there are currently only upper limits. The parameter space could also be reduced by including the PPDs from the results of this work as priors in future work. It is important to understand whether an MCMC procedure has converged to a stationary distribution. The sampling is not complete if increasing the length of that chains would considerably alter the posterior distribution of the reaction rates. It is not possible to be certain that convergence has been reached but several heuristics are available and are considered in this appendix.
Most simply, the chains themselves can be inspected. Figure 6 shows every 500th step in an example chain for reaction 1. The walker repeatedly leaves the area of maximum likelihood and then returns. This cycle repeats a sufficient number of times for it to be unlikely that there is a undiscovered mode.
More rigorously, the Geweke diagnostic can be used (Geweke 1992). In this test, it is considered that if the chain has converged any two samples of the chain will have the same mean, within the variance of the samples. This is typically tested on the first 10% and the final 50% of a chain that is thought to have converged. In this work, a sample of chains were tested by breaking each chain into 10 subsamples and comparing each to the mean of the final 50% of the chain. In every case, the mean of the subsample was consistent with the mean of the larger sample. Since the value of the Geweke diagnostic should be zero, within the variance of the chains, multiple values for converged chains should follow a normal distribution. In Figure 7 the values of the diagnostic for many subsamples of the chains from this work are plotted as histogram with a normal distribution plotted for comparison. The auto-correlation time is another diagnostic that can be calculated and the tutorial provided in the emcee documentation 1 is used for this. This provides two heuristics. Firstly, once the chains reach a sufficient length that the autocorrelation time can be reliably calculated, it is likely that the chain has converged. Secondly, one use of the quantity is to calculate the sampling error in an MCMC chain. An autocorrelation time of 10 4 steps was calculated from the chains, effectively giving approximately 100 independent samples per chain. If the mean value of the chain is considered then the variance on this mean is given through the equation, where N is the number of sample and V ar[f (θ)] is the variance of the chain. In the case of an average chain in this work, the sampling error on the value of the mean is ∼1%. Finally, the posterior distribution was also evaluated using the code pyMultinest (Buchner et al. 2014) and the marginalised posteriors were consistent with those found using emcee. The consistency between these two different methods of sampling the posterior is good evidence for convergence. Ultimately, given the above heuristics and the fact that the initial MCMC runs produced approximately the same results for chains of 100,000 steps as they do in the ∼10 6 step chains used in the final work, it is assumed that the chains have, in fact, converged.