UNDERSTANDING THE FORMATION AND EVOLUTION OF INTERSTELLAR ICES: A BAYESIAN APPROACH

Our aim is to obtain information about the physical parameters of a molecular cloud $\boldsymbol{\theta }=(\theta _1,\theta _2,\ldots,\theta _k)$ while we measure molecular abundances $\boldsymbol{\mathcal {Y}}=(\mathcal {Y}_1,\mathcal {Y}_2, \ldots,\mathcal {Y}_n)$. These quantities are related to a (forward) function f(·) that represents the physical and chemical processes in the cloud. The main challenge is that there is no closed form function, f, mapping the parameters to the observations, which could be inverted. However, given a set of parameters, estimated abundance values for the species of interest can be computed with astrochemical models denoted here as $\mathcal {C}(\cdot)$. The addressed problem in our case is how to estimate $\boldsymbol{\theta }$ from

$$\begin{equation} \boldsymbol{\mathcal {Y}}=\mathcal {C}(\boldsymbol{\theta }) + \boldsymbol{\varepsilon } \end{equation} \tag{ 1 }$$

and, according to Idier (2008), this constitutes an inverse problem. The error term, $\boldsymbol{\varepsilon }$, represents both the observational noise and the modeling error between $\mathcal {C}(\cdot)$ and f(·).

We treat $\boldsymbol{\mathcal {Y}}$, $\boldsymbol{\theta }$, and $\boldsymbol{\varepsilon }$ as random variables and define the solution of the inverse problem to be the posterior probability distribution (PPD) of the parameters given the observations. This allows us to model the noise via its statistical properties, even though we do not know the exact instance of the noise entering our data. We can also optionally specify a priori the form of solutions that we believe to be more likely, through a prior distribution. Thereby, we can attach weights to multiple solutions that explain the data. This is the Bayesian approach to inverse problems.

We assume that we have K parameters θ_k and N solid-phase observable quantities $\mathcal {Y}_n$. The error ε_n on each observation $\mathcal {Y}_n$ is assumed to be normally distributed with variance, $\sigma _n^2$. In addition, it is assumed that the observational errors are independent. The $\sigma _n^2$ is considered to correspond to the uncertainty on $\mathcal {Y}_n$, which is solely dictated by the observation. The probability density function of the errors is given by

$$\begin{equation*} p_\varepsilon (\boldsymbol{\varepsilon }) = \prod _{n=1}^N \frac{1}{(2\pi)^\frac{1}{2}\sigma _n^2} \exp \left (\frac{\varepsilon _n^2}{2\sigma _n^2}\right). \end{equation*}$$

Using Equation (1), we can define the likelihood function, $\mathbb {L}$, of observations given a model parameterized by a set of parameters as

$$\begin{eqnarray*} \mathbb {L}(\boldsymbol{\theta };\boldsymbol{\mathcal {Y}}) &=& p_\varepsilon (\boldsymbol{\mathcal {Y}}-\mathcal {C}(\boldsymbol{\theta }))= \prod _{n=1}^N \frac{1}{(2\pi)^\frac{1}{2}\sigma _n^2} \nonumber \\ && \times \exp \left (-\frac{1}{2}\sum _{n=1}^N \left[\frac{\mathcal {C}_n(\boldsymbol{\theta })-\mathcal {Y}_n}{\sigma _n}\right]\right). \end{eqnarray*}$$

In case any prior information about the unknown parameters is available, the Bayesian approach allows for this information to be taken into account. This information can be integrated through a prior probability distribution on the parameters, say $\pi (\boldsymbol{\theta })$. Then, a parameter estimation can be performed through the PPD, using Bayes' rule:

$$\begin{equation} \pi (\boldsymbol{\theta }|\boldsymbol{\mathcal {Y}}) = \frac{\mathbb {L}(\boldsymbol{\theta };\boldsymbol{\mathcal {Y}})\pi (\boldsymbol{\theta)}}{m(\boldsymbol{\mathcal {Y}})}. \end{equation} \tag{ 2 }$$

The PPD expresses our uncertainty about the parameters after considering the observations and any prior information. The denominator is simply a normalization factor.

In reality, we are not able to access the whole PPD. Therefore, computation of parameter estimates or uncertainties is a hard task. MCMC methods are efficient methods that allow us to sample from complex probability distributions and approximate complex probability densities.

MCMC methods are a powerful class of algorithms that produce random samples distributed according to the distribution of interest. The importance and efficiency of MCMC methods lies in the fact that these samples can be used to approximate the probability density of the distribution by only calculating it for a feasible number of parameter values. Among the several implementations of possible algorithms, we employ an MH sampling algorithm (Gilks et al. 1995). The MH algorithm will enable us to explore the parameter space and efficiently approximate the PPD. A theoretical introduction of MCMC and MH is far beyond the scope of this paper. However, in Appendix A, we briefly describe the MH algorithm and how MCMC is employed for parameter estimation in our case. Note that the tuning of the MH algorithm, as described in Appendix A, is very crucial when aiming to approximate possibly multi-modal and non-Gaussian distribution, which is the case for this study.

The chemical modeling code used in this paper, and denoted as $\mathcal {C}(\cdot)$ in Equation (1), is the UCL_CHEM time-dependent gas–grain chemical code (Viti et al. 2004) and is briefly described in Appendix B and references therein. Note that for each set of parameters, $\mathcal {C}(\cdot)$ provide us with a time series of chemical abundances. We choose to extract the chemical abundances of interest for the time points when the final density is reached and the cloud collapse has finished. Even though we ignore the previous time points, the time dependency is still taken into account and investigated through exploration of different final density values.

The parameters for our chemical modeling code create a nine-dimensional (9D) parameter space for molecular clouds, as used in our MH and described in Table 1:

$$\begin{equation*} \boldsymbol{\theta } = (n_{\rm H},\zeta,G_{\circ },C_f,fr,\epsilon, \phi, y, r). \end{equation*}$$

In a first attempt to employ a Bayesian approach for deriving branching ratios for poorly understood chemical reaction pathways, we also investigated the parameter r, which controls how much of the gas-phase oxygen turns into ice H₂O or ice OH. Parameter r reflects the percentage of O that turns into H₂O, so that 1 − r reflects the percentage of oxygen that turns into OH. Desorption efficiencies resulting from H₂ formation on grains, direct cosmic ray heating, and cosmic ray induced photodesorption are determined by parameters , ϕ, and y, as introduced and studied by Roberts et al. (2007). The freeze-out parameter in our code is effectively the sticking coefficient, a number in the range of 0%–100% that adjusts the rate per unit volume at which species deplete on the grain. For the free-collapse of a particular n_H, we used the modified formula of Rawlings et al. (1992), where parameter C_f is considered to be a retardation factor with a value less than one, to roughly mimic the magnetic and/or rotational support, or an acceleration factor with a value greater than one to simulate a collapse faster than a free-fall (e.g., due to external pressure). Table 1 lists the set of physical parameters studied in this paper along with their definition domain, $\mathbb {D}_{{\theta _k}}$. The joint definition domain, $\mathbb {D}_{{\theta }}$, represents the parameter space to explore. The selected domain limits refer to the theoretical range of possible values for molecular clouds where atoms and molecules deplete onto the dust, ensuring that extreme values are included.

Table 1. Parameter Definition Domain

Parameters $\boldsymbol{\theta }$	Unit	Definition Domain $\mathbb {D}_{{\theta }}$
ζ	10−17 s−1	1–10
G○	Habing	1–10
nH	cm−3	104–108
fr	...	0%–100%
Cf	...	0.5–3
	Yield per H2 formed	0.01–1
ϕ	Yield per cosmic ray impact	102–106
y	Yield per photon	10−3–102
r	...	0%–100%

Download table as:  ASCIITypeset image

The observational constraints of our analysis are based on data from the existing literature. Even though in this application we are primarily interested in ices, we include both gas-phase and solid-phase observations. To avoid confusion, we will use $\boldsymbol{\mathcal {Y}}$ to denote a vector containing any observed quantity and, if required, we will specify whether we refer to solid-phase or gas-phase observations.

The solid-phase observations include column densities and visual extinction data for molecular clouds in front of field stars. Such sources often provide suitable opportunities to observe and study ices in quiescent regions of the clouds (e.g., Boogert et al. 2011). We used 31 observations of H₂O, CH₃OH, CO, and CO₂ from 31 different regions of 16 different clouds found in the literature and summarized by Whittet et al. (2011). The data suggest some abundance variation, which was attributed to different evolutionary stages for different clouds. The scope of this paper lies beyond studying the behavior of a specific cloud, but rather on how to obtain statistical insight into the dynamics of common cloud classes. Therefore, the observational data is transformed into fractional abundances with respect to total H nuclei, and then the average value is computed and used for our analysis.

In an attempt to minimize degeneracies, we introduce additional gas-phase abundances as an optional observational constraint. Due to the ill-posed nature of our problem, it is possible for our chemical model to end up with a solution space that fits the solid-phase observations perfectly, but with gas-phase abundances that are far from realistic. Hence, the addition of gas-phase observations can be considered a mathematical regularization through the introduction of additional prior information. Prior information can be naturally integrated into our Bayesian approach. The gas species observations were collected from more than one study, attempting to match the clouds, regions, or evolutionary stage of the observational sources used for the solid-phase species. If we were to fit observations of a particular source, then, ideally, every observational gas-phase constraint should be able to contribute to the regularization of our methodology. However, because here we are only attempting to explore a methodology, we found that three gas-phase species were adequate to provide insight into the efficiency of gas-phase species as a regularization factor. Abundances for NH₃ and N₂H⁺ were collected from Johnstone et al. (2010), while HCO⁺ from Schöier et al. (2002). The gas-phase observations are in the form of fractional abundances with respect to total hydrogen nuclei. Table 2 lists the average molecular abundances for all the species along with their uncertainties. We emphasize again that the error on each of the observations $\mathcal {Y}_n$ is assumed to be normally distributed with a variance of, $\sigma _n^2$, that is determined solely by the uncertainty reported in Table 2.

We run two identical sets of eight MCMC chains that differ in prior distribution information. For the first set, the prior information is noninformative and in the form of an acceptable range of possible values. Therefore, $\pi (\boldsymbol{\theta })$ is uniformly distributed on $\mathbb {D}_{{\theta }}$, as listed in Table 1. Note that the observational data, $\boldsymbol{\mathcal {Y}}$, refers only to the solid-phase molecular abundances and, in this case, the gas-phase species are ignored. In the second case, the prior information includes the observational constraints from the gas-phase species as well. Let $\boldsymbol{\mathcal {Y}}$ now include all of the observational constraints, $\boldsymbol{\mathcal {Y}}_s$ just the solid-phase, and $\boldsymbol{\mathcal {Y}}_g$ the gas-phase observational constraints. In that case, the PPD is defined as

$$\begin{equation} \pi (\boldsymbol{\theta }|\boldsymbol{\mathcal {Y}})= \pi (\boldsymbol{\theta }|\boldsymbol{\mathcal {Y}}_s,\boldsymbol{\mathcal {Y}}_g) = \frac{\pi (\boldsymbol{\mathcal {Y}}_s|\boldsymbol{\theta },\boldsymbol{\mathcal {Y}}_g)\pi (\boldsymbol{\theta }|\boldsymbol{\mathcal {Y}}_g)}{m(\boldsymbol{\mathcal {Y}})}. \end{equation} \tag{ 3 }$$

The prior information is simply the likelihood function, $\mathbb {L}(\cdot)$, of $\boldsymbol{\mathcal {Y}}_g$, given a model parameterized by $\boldsymbol{\theta }$, since

$$\begin{eqnarray*} \pi (\boldsymbol{\mathcal {Y}}_s|\boldsymbol{\theta },\boldsymbol{\mathcal {Y}}_g) &=& \mathbb {L}(\boldsymbol{\theta };\boldsymbol{\mathcal {Y}}_s), \\ \pi (\boldsymbol{\theta }|\boldsymbol{\mathcal {Y}}_g) &\propto& \mathbb {L}(\boldsymbol{\theta };\boldsymbol{\mathcal {Y}}_g)\pi (\boldsymbol{\theta)}. \end{eqnarray*}$$

Including prior information in this way is equivalent to attaching weight to the solutions that explain the gas-phase as well as the solid-phase chemistry.

In order to quantitatively investigate the effectiveness of our method to astrochemical problems, we performed a benchmark test. This benchmark test is basically our Bayesian analysis applied, this time, on synthetic observations produced by UCL_CHEM using a predefined set of parameters, $\boldsymbol{\theta _T}$. Once we have our synthetic observations, we apply our methodology and analyze the results and whether the true parameters are recovered. Knowing the solution to this test a priori allows us not only to validate the method, but also to critically perceive the nonlinear and ill-posed nature of our problem. This discussion can be found in Section 3.1. This particular selection of parameters was chosen at random and includes values that are not far from the expected parameters, or well-accepted values in the literature. The parameter values used in the test can be found in Table 3.

Table 3. Blind Benchmark Test

Parameters $\boldsymbol{\theta }$	Unit	Test Value
ζ	10−17s−1	2.4
G○	Habing	2.6
nH	cm−3	105
fr	...	42%
Cf	...	1.3
	Yield per H2 formed	0.02
ϕ	Yield per cosmic ray impact	150
y	Yield per photon	0.1
r	...	75%

Download table as:  ASCIITypeset image

The results of the performed test as shown in Figure 1 reveal two important insights. First of all, high probability density regions for all the parameters include, and hence recover, the true parameters. As we can see in Figure 1, all of the predefined parameter values lie under, or very close to, the highest density point of the marginal PPD. This result simply validates that both the Bayesian approaches make accurate inferences based on the given observations and the MH algorithm efficiently samples the solution space. Secondly, we can observe that, in many cases, there are additional high probability density regions. These regions prove and highlight the ill-posed nature of our problem by indicating that different parameter sets can produce similar observations. Combining the two insights, we can conclude that the Bayesian method with MCMC sampling is efficiently exploring the parameter space, revealing the solution regions that answer our ill-posed inverse problem. In addition, we can conclude that, in order to constrain our solution space, we should either introduce numerical regularization factors (e.g., gas-phase species) or scientific prior knowledge.

A visual comparison of Figures 2 and 3 reveals what we can quantitatively observe in Table 4. With noninformative uniform prior, the HDRs seem to cover large sections of the distribution, which in some cases reach 50% of the definition domain. This means that most of the parameters are not constrained enough. The most statistically straightforward parameters seem to be clearly the n_H and then the ϕ and fr parameters, presenting distinct modes and relatively low HDS. G_○ and ζ seem neither constrained nor relevant enough, while fr seems to have a clear mode, followed by a very heavy tail. The rest of the parameters present high HDS, above 40% with several disjoint HDRs and do not allow us to reach credible conclusions about the parameters. Including the prior information from the gas-phase species significantly changes the picture, as can be seen in both Figure 3 and Table 4. We can observe that the HDRs get smaller and the parameters seem more constrained. The distribution of ζ is now denser around high values (>6), while G_○ has to be low (<4). The n_H remains well-constrained with even lower HDS, while the distribution of fr now clearly constraints the parameter to low domain values. The distribution of C_f is also significantly altered: not only the HDS has dropped, but a large portion of the density has also transferred from high accelerated collapse regions to free-fall collapse regions. The nondesorption mechanisms still present a multi-modal behavior, but with significantly smaller HDRs. Their distribution clearly highlights the nonlinear way these mechanisms act together or against each other. For r, the addition of informative prior information seems to reduce the HDS as well, centralizing the density, but still slightly favoring the production of H₂O against OH. Therefore, we conclude that the addition of gas-phase species as a regularization factor outperforms the use of a noninformative uniform prior distribution alone. The HDS is reduced at an average of ∼12%, which indicates an equivalent constraint on the parameter space. Section 3.3 will discuss the statistical and numerical results of our analysis, while Section 3.4 will discuss the astrophysical implications. In both of these sections, we shall only concentrate on the results of the informative prior case.

The n_H is clearly the most constrained physical parameter. The marginal density function reveals that most of the density is between 2.2 and 5 × 10⁴ cm⁻³. The ζ is constrained to values higher than 6 × 10¹⁷ s⁻¹, while the G_○ is constrained to values lower than 4 Habing. The HDR for the fr stays between 20% and 45%, while the C_f has one distinct HDR between 0.5 and 1.55 and one long heavy tail between two and three times the default free-fall rate. The presents two modes. The first HDR is between 0.4 and 0.8 and the second between 1.2 and 1.4. The marginal distribution of ϕ, also presents two modes. One is centered around 10⁵. The second one is centered around 60. The marginal distribution for y, presents two disjoint HDRs as well. The first one indicates really low efficiency of about 10⁻⁶, while the second one indicates a slightly higher 2 × 10⁻³–8 × 10⁻². Finally, the distribution for the branching ratio parameter, r, shows high density between 40% and 70% of oxygen turning into ice water.

In Figure 4, we show the marginalized 2D PPD for our parameters. Note that the n_H and the fr are negatively dependent in a nearly linear way. On the other hand, G_○ and n_H seem to have a nonlinear positive correlation, hitting a plateau after a certain gas density. Similarly, the fr and the C_f may have a clear peak, but also show some evidence of a positive correlation. The relation between the cosmic ray desorption efficiency parameters, ϕ and y, reveals many distinct peaks throughout the domain space. Note that the marginalized PPD for cosmic ray ionization rate and parameter ϕ shows a clear bimodal structure. However, focusing only on the denser areas of the distribution, we can observe a potential nonlinear correlation between the cosmic rays and the efficiency of the cosmic ray related parameter ϕ. However, in general, ζ is evidently a parameter that is not sufficiently constrained. This is already obvious by the 1D marginal distribution of ζ, but the contrast of constraint between ζ and one of the most constrained parameters, such as n_H, is depicted in Figure 2(f).

Due to the nonuniqueness of our solution space, examining the joint probability distribution of the PPD provides a useful insight. The dimensionality of the distribution makes a visualization impossible; thus, we chose to extract the statistical mean and the standard deviation for each one of the parameters from the most probable mode of the joint distribution. The joint distribution was approximated using a multivariate histogram and the most probable mode was chosen in a heuristic way and corresponds to ∼35% HDR of the whole PPD. The values for the mean and standard deviation are given in Table 5. As expected, the most probable mode of the joint PPD agrees with the HDR of the marginal parameter distributions. For the unimodal 1D marginalized distributions, the most probable mode completely coincides, while, for the multi-modal cases, the most probable mode coincides with one of the modes. Hence, purely based on the statistical interpretation, we conclude that a molecular cloud that matches the observed abundances should have a low n_H, a low fr, and a low G_○. The ζ, on the other hand, is more likely to have high values, but the high standard deviation leaves room for significant variation. The collapse of the cloud may be insignificantly accelerated, while the branching ratio, r, slightly favors the branching into water, but with a high standard deviation. In terms of the nonthermal desorption efficiency parameters, we notice increased efficiency for all three of them. As a general result, we conclude that the 9D space of the joint distribution has multiple peaks. Both the marginalized distributions and the denser peak of the joint distribution indicate that some of the parameters (n_H, G_○, fr, C_f) are well-constrained, while other parameters (ζ, r, , ϕ, y) present possible variation that implies further astrophysical or statistical implications.

Here, we discuss our results for each of the parameters with regards to their astrophysical implication.

n_H. The derived credible intervals for the gas density are in very good agreement with the properties of typical collapsing dark clouds, clumps, and cores (Myers & Benson 1983; Benson & Myers 1989; Bacmann et al. 2002; Bergin & Tafalla 2007). Higher cloud densities (>10⁶ cm⁻³), that are usually expected in hot cores after the cloud has collapsed (van der Tak 2004) were explored, but showed nearly zero probability density in our analysis.

fr. Our study implies a depletion rate that is not high enough to dominate and is probably lower than 50%. Bacmann et al. (2002) suggest that freeze-out dominates when n_H exceeds ∼3 × 10⁴ cm⁻³, which is marginally the case in our study. When the freeze-out dominates and densities exceed ∼10⁵ cm³, the abundance of CO ice is found to be significantly increased to typical gaseous values (∼10⁻⁴; Pontoppidan 2006; Bergin & Tafalla 2007). Furthermore, the ice water abundance is typically 5 × 10⁻⁵–9 × 10⁻⁵ and even higher at the highest densities (Pontoppidan et al. 2005). The ice CO and H₂O abundances in our case, however, are about 0.5–1 mag lower. Therefore, along with the n_H results, the lower freeze-out rate estimated by our analysis can be explained by a different evolutionary stage of the observed clouds. According to Fontani et al. (2012), low depletion values can also imply a cloud that is going to form less massive objects.

G_○. Our analysis showed that in order to match the observed ice abundances, the G_○ is comparable to the standard interstellar radiation field of 1 Draine or ∼1.7 Habing (Draine 1978).

ζ. In dense gas ζ is measured to be in the range of 1–5 × 10⁻¹⁷ s⁻¹ (Bergin et al. 1999). However, considerable uncertainties have been reported in the literature with derived values as high as 10⁻¹⁵ s⁻¹, accounted to X-rays from a central source (Doty et al. 2004). These discrepancies may be due to whether ζ is determined via H$_3^+$ or HCO⁺ measurements. However, Dalgarno (2006) claims that given the latest evidence that the ζ range is narrow and between 10⁻¹⁶ and 10⁻¹⁵, the question should be focused not on why the estimations are different, but on why they are so similar. Our analysis confirms ζ values higher than the typical estimations and is even consistent with the 10⁻¹⁶ s⁻¹ estimations through the H$_3^+$ determination. Most importantly, our study indicates the high standard deviation of these values, highlighting that such a variation should be expected. Theoretically, this is explained considering the fact that ζ lose energy while ionizing and exciting the gas through which they travel in conjunction with the possible variation in the origin of ζ. Even though our astrochemical model does not account for the latter factors, our probabilistic approach reflects their impact.

C_f. Our study shows that the collapse of the cloud should follow the expected free-fall collapse. Higher C_f values present moderate probability density, which implies that the observational constraints could potentially also be matched with different, but also less likely, sets of parameters (e.g., higher values for both C_f and fr).

, ϕ, y. The desorption from H₂ efficiency parameter () estimates are significantly higher than the value reported by Roberts et al. (2007) ( < 0.1). The direct cosmic ray desorption efficiency (ϕ), presents two peak values. One of them agrees with Roberts et al. (2007) and is centered around 10⁵. The second one is centered around 60, which is lower than the lowest limit studied by Roberts et al. (2007). For the cosmic-ray-induced photodesorption efficiency (y), we have two probable estimates as well. The first one indicates really low efficiency. The second one presents a slightly higher efficiency that is still lower than the one estimated by Hartquist & Williams (1990) (y = 0.1), but consistent with the results of Öberg et al. (2009) for CO₂. Our analysis, in general, indicates useful credible intervals for nonthermal desorption efficiencies, highlighting that the reported nonlinearities can be tackled with further regularization factors such as molecule specific analysis and additional grain properties. Also note that our astrochemical model does not include direct UV photodesorption, which has recently be found to be efficient.

r. The branching ratio proved to be a parameter with high but anticipated variability. Its marginal probability distribution presents the most statistically normal behavior with a mean that implies a shared branching ratio of oxygen freezing into ice H₂O and ice OH, favoring slightly solid H₂O. The first laboratory experiment to reproduce the ice H₂O formation (Dulieu et al. 2010) implied that the hydrogenation of oxygen is an important route for water formation. Furthermore, Cazaux et al. (2010) state that species such as OH are only transitory and quickly turn into ice water. However, they also state that ∼30% of the O coming on the grain is released in the gas phase as OH, which can freeze back as H₂O. A pathway that is included in our model and can explain both the high water abundance on the grains and the shared branching ratio, r. At last the high branching ratio toward ice OH highlights the importance of ice OH for the production of ice CO₂.

We now look at the correlation between our parameters as presented in Figure 4. The relation between the n_H and depletion or freeze-out has been the subject of many studies (Bacmann et al. 2002; Christie et al. 2012; Fontani et al. 2012; Hocuk et al. 2014). They all conclude that the amount of depletion, the ice abundances, and the density of the cloud should all scale together (Rawlings et al. 1992). Even though our analysis suggests a clear anti-correlation between n_H and fr (Figure 4(a)), this result is completely in line with the literature, since we are not analyzing the time evolution of the cloud, but instead focus on parameter fitting at specific time points. This negative correlation suggests that the less gas density we have, the higher the depletion should be in order to match the observed ice abundances. Our results are also in line with the negative correlation between the depletion factor and n_H, derived by Fontani et al. (2012) from CO observations. The positive correlation between n_H and G_○ depicted in Figure 4(b) confirms that the denser a cloud is, the higher ζ values are needed to match the observations. The plateau after a density value (∼7 × 10⁴ cm⁻³) indicates that the explored radiation field domain space is not high enough to penetrate the cloud after a density threshold. When C_f is increased, the freeze-out timescale needs to be decreased since the final n_H is reached quicker. This reduced timescale requires higher fr values in order to simulate the observed ice abundances, and this relation is depicted in Figure 4(c). Even though it is not very straightforward, the relation between the cosmic ray desorption efficiency parameters, ϕ and y, is very interesting (Figure 4(d)). In most cases, the cosmic ray photodesorption efficiency is either low or high for both direct cosmic ray heating and cosmic ray induced cases. However, there is a significant peak when the direct cosmic ray impact is very efficient, while the cosmic ray induced impact is very inefficient.