A&A 459, 813-820 (2006)
DOI: 10.1051/0004-6361:20065472
V. Wakelam^{1,2} - E. Herbst^{1,3} - F. Selsis^{4} - G. Massacrier^{4}
1 - Department of Physics, The Ohio State University, Columbus, OH 43210, USA
2 -
Observatoire Aquitain des Sciences de l'Univers, Laboratoire d'Astrodynamique, d'Astrophysique et d'Aéronomie de Bordeaux, CNRS/INSU UMR 5804, BP 89, 33270 Floirac, France
3 -
Departments of Astronomy and Chemistry, The Ohio State University, Columbus, OH 43210, USA
4 -
École Normale Supérieure de Lyon, Centre de Recherche Astronomique de Lyon,
46 allée d'Italie, 69364 Lyon Cedex 07; CNRS, UMR 5574,
Université de Lyon 1, Lyon, France
Received 21 April 2006 / Accepted 2 August 2006
Abstract
Aims. To determine whether or not gas-phase chemical models with homogeneous and time-independent physical conditions explain the many observed molecular abundances in astrophysical sources, it is crucial to estimate the uncertainties in the calculated abundances and compare them with the observed abundances and their uncertainties. Non linear amplification of the error and bifurcation may limit the applicability of chemical models. Here we study such effects on dense cloud chemistry.
Methods. Using a previously studied approach to uncertainties based on the representation of rate coefficient errors as log normal distributions, we attempted to apply our approach using as input a variety of different elemental abundances from those studied previously. In this approach, all rate coefficients are varied randomly within their log normal (Gaussian) distribution, and the time-dependent chemistry calculated anew many times so as to obtain good statistics for the uncertainties in the calculated abundances.
Results. Starting with so-called "high-metal'' elemental abundances, we found bimodal rather than Gaussian like distributions for the abundances of many species and traced these strange distributions to an extreme sensitivity of the system to changes in the ratio of the cosmic ray ionization rate
for He and that for molecular hydrogen
.
The sensitivity can be so extreme as to cause a region of bistability, which was subsequently found to be more extensive for another choice of elemental abundances. To the best of our knowledge, the bistable solutions found in this way are the same as found previously by other authors, but it is best to think of the ratio
as a control parameter perpendicular to the ``standard'' control parameter
.
Key words: astrochemistry - ISM: abundances - ISM: clouds - ISM: molecules
The equations describing the chemical evolution of quiescent cores (also known as molecular clouds) are highly non-linear and can result in extreme sensitivity to the initial conditions or to some of the parameters. For example, the elemental ratio of carbon to oxygen is well known to be an important parameter for dense cloud chemistry (Terzieva & Herbst 1998). Another well-known consequence of this non-linearity is the presence of bistable regions for some ranges of model parameters (temperature, density, comic-ray ionization rate, rate coefficients and elemental abundances). The phenomenon of bistability is characterized by two stable solutions for chemical abundances at steady-state for the same set of parameters. Bistability in dark cloud conditions, i.e. at low temperature, was first discovered by Le Bourlot et al. (1993) and subsequently studied by a large number of authors (Charnley & Markwick 2003; Pineau des Forêts & Roueff 2000; Shalabiea & Greenberg 1995; Lee et al. 1998; Boger & Sternberg 2006; Le Bourlot et al. 1995). The two solutions of this bistability are characterized by a large difference in the C/CO and O/O_{2} ratios and in the ionization fraction. For this last reason, they were named the High and Low Ionization Phases (HIP and LIP). Le Bourlot et al. (1995) showed that the region of bistability can be mapped out by variations in density, temperature, cosmic-ray ionization rate, and elemental depletions within specific ranges. Pineau des Forêts & Roueff (2000) later showed that even variations in rate coefficients can lead to the two solutions. Recently, Boger & Sternberg (2006) identified the chemical mechanisms of this phenomenon as a cycle involving the H^{+}_{3}, O_{2}, and S^{+} species.
In two previous papers, we presented a Monte Carlo method to compute the theoretical error of abundances in chemical models due to uncertainties in rate coefficients and physical conditions (Wakelam et al. 2005,2006). The errors in the rate coefficients are defined by a log-normal distribution. As a consequence, the resulting abundances at steady-state usually follow a Gaussian distribution (Wakelam et al. 2005). When applying this method, we found that the chemical abundances are very sensitive to the ratio between the ionization rates of helium and molecular hydrogen caused by cosmic-ray particles. Characterized by a variation of several orders of magnitude in some abundances for small variations of , this sensitivity can lead to bimodal rather than Gaussian distributions for the abundances. For certain ranges of the parameters, it can even lead to bistability, with the ratio appearing to be a control parameter.
This paper, which describes the sensitivity and some of its consequences, is organized as follows. In Sect. 2, we briefly present the chemical model and the uncertainty method used for this analysis. In a somewhat unorthodox way, we have decided to present, in Sect. 3, the manner in which the sensitivity was found while studying the use of so-called "high-metal'' elemental abundances. Section 4 shows the influence of the elemental abundances and the chemical network used on the sensitivity. In Sect. 5, we present and discuss the phenomenon of bistability as obtained with the variation of . In Sect. 6, the chemical reactions involved in this sensitivity/bistability are elucidated. Finally, the last section presents our conclusions.
Table 1: High- and low-metal elemental abundances with respect to H_{2}.
We used a gas-phase time-dependent chemical model with the osu.2003^{} network reported by Smith et al. (2004) (4233 reactions, 421 species and 12 elements). The model computes the evolution of the species for a fixed temperature and density. The initial conditions are all atomic except for molecular hydrogen, while the chemical model and database are the same as used in Wakelam et al. (2006). For this work, we will present the results using three different sets of elemental abundances: the high-, low- and intermediate-metal cases. The high- and low-metal elemental abundances have been defined by Graedel et al. (1982) and are listed in Table 1. The intermediate case has the same abundances as the low-metal one except that the amount of the element sulfur is raised to compared with H_{2} (also given in Table 1). The other parameters are the typical ones for quiescent cores: a kinetic temperature of 10 K, an H_{2} density of 10^{4} cm^{-3}, a visual extinction of 10 so that the photochemistry driven by the external UV photons does not occur, and a fixed cosmic-ray ionization rate of s^{-1}.
The method used in this analysis to derive uncertainties in abundance is described in detail in Wakelam et al. (2005). The uncertainties in rate coefficients derive from the UMIST database when available, but for most of the reactions, including the ionization rates of H_{2} and He, we have assumed a factor of 2 uncertainty. Note that the uncertainty in ionization rates is not an uncertainty in the parameter , but an uncertainty in the multiplicative factor that distinguishes the individual rates. Although ionization rates for H_{2} and He have been used in a large number of models, the current rates used in terms of derive back to the Ph. D. thesis of Black (1975), where the factors are 0.93 for the production of H_{2}^{+} ( ) and 0.5 for the production of He^{+} ( ), leading to a ratio of 0.54, if one neglects the minor channels for H_{2}. An older ratio of unity was used by Herbst & Klemperer (1973). New estimates, involving both direct ionization and secondary ionization by energetic electrons, are sorely needed (Black, private communication) based on the treatment of Dalgarno et al. (1999) for a mixture of atomic and molecular hydrogen and helium.
The uncertainty method consists of generating N new sets of rate coefficients by replacing each coefficient k_{i} by a random value consistent with its uncertainty factor F_{i}. We assume a normal distribution of with a standard deviation . We run the model for each set j, which produces, for each species, N values of the fractional abundances X_{j} (t) at a time t. For this work, we ignore the uncertainty in physical conditions and consider the uncertainty in the rate coefficients only. A total N of 2000 different runs was made for each set of parameters studied; this number is large enough to achieve statistical significance (see Wakelam et al. 2006).
Figure 1: Density of probability of the O_{2} abundance as a function of time ( left panel). The right panel shows the histogram of the abundance at 10^{7} yr. The elemental abundances are for the high-metal case. | |
Open with DEXTER |
Figure 2: Distributions (nomalized to the total number of runs) of the abundances of C, H_{2}O, e^{-} and H_{3}^{+} at 10^{7} yr (high-metal case). The solid lines refer to the variation of all the rate coefficients whereas the dotted lines represent the variation of and only. The vertical dashed line represents the abundance computed with the standard values of the rate coefficients. | |
Open with DEXTER |
The computation of the theoretical error due to the uncertainties in rate coefficients was first applied to the case of high-metal elemental abundances. Figure 1 shows the results of the random variation of all the rate coefficients for the O_{2} abundance as a function of time. As can be seen from this figure, the distribution of the abundance is spread into two peaks after 10^{6} yr, as steady-state approaches, with a significant number of curves in each peak and some stable solutions between the two peaks. The bimodality of the distribution at a specific time (10^{7} yr) is also shown as a histogram in Fig. 1. Such bimodal distributions are obtained for a large number of species (C, H_{2}O, SO, CS etc.) with the largest separation between the two peaks for atomic carbon (three orders of magnitude). Figure 2 shows histograms for the abundance distributions of the species C, H_{2}O, e^{-}, and H_{3}^{+} at 10^{7} yr. It can be seen that the H_{3}^{+} distribution does not show a bimodal profile but is Gaussian. The same is true for H_{2}S and OH. Since the error in the rate coefficients follows a log-normal distribution, the bimodal shapes are due to non-linear effects.
Figure 3: Histograms of the rate coefficients of the following reactions ( upper panel): He + CRP He^{+} + e^{-} ( , left box), H_{2} + CRP H_{2}^{+} + e^{-} ( , middle box) and H_{3}^{+} + e^{-} H_{2} + H (right box). The thick and thin lines represent the histograms of all the runs (hist 1) and of the runs giving the higher abundance of O_{2} (hist 2, see text) respectively. The lower panels show the ratio between hist 2 and hist 1. | |
Open with DEXTER |
We obtained the bimodal distributions by varying the rate coefficients, we then decided to identify the reactions responsible for them if at all possible. We proceeded as follows. Among the 2000 runs, we listed the sets of reactions that give one the higher of the two most probable abundances for O_{2}. Then for each rate coefficient, we plotted the histograms for all runs and the histograms for those runs giving only the higher abundance of O_{2}. In the upper panels of Fig. 3, pairs of histograms can be seen for the rate coefficient of the cosmic ray ionization of He, that of H_{2}, and the dissociative recombination of H_{3}^{+} to form H_{2} and H. The thick lines (hist1) represent histograms for all runs while the thin lines (hist2) represent those ending up with the higher O_{2} abundance. For reactions not responsible for the bimodal distributions, there should be no difference between the pairs of histograms except that the total number of runs should be lower for the second one. For most of the reactions, we indeed found no difference, as for the reaction H_{3}^{+} + e^{-} H_{2} + H in Fig. 3. In fact, there were only two exceptions: the direct ionization rate coefficients of He and H_{2} by cosmic rays ( and ), and their marked effect can be seen in Fig. 3. The high abundance of O_{2} tends to be formed with lower values of and higher values of , as can be seen most clearly in the lower panels of Fig. 3, in which the ratios of the two histograms are plotted. It is then convenient to consider the sensitivity of the abundance of O_{2} and other species to the ratio / rather than to their absolute values. To confirm their importance, we randomly varied only and , keeping the other rate coefficients constant, and we obtained the same bimodal distribution as with the fully random method at steady state. The dispersion of the abundances in this case is then almost as large as when we varied all the rate coefficients, as can be seen in the dotted lines in Fig. 2.
Figure 4: Abundance of O_{2} as a function of / at 10^{7} yr. The upper and middle panels represent the abundance computed with the random distribution of all the rate coefficients and of and only, respectively. The solid lines in the lower plot have been computed using one value of for each line and varying , while the dotted lines have been computed using one value of for each line and varying . The vertical dashed line indicates the standard value of / . | |
Open with DEXTER |
Figure 5: Abundances of C, e^{-}, H_{2}O and H_{3}^{+} as a function of / and for five different values of (solid lines). The vertical dashed line indicates the standard value of / . | |
Open with DEXTER |
The sensitivity of the chemical abundances to the ratio / can be explored in another way. In Fig. 4 (left panel), we show the steady-state abundance of O_{2}, computed by randomly varying all the rate coefficients as a function of / . For comparison, the abundance computed with the random variation of and only is shown in the middle plot of the figure. The smaller vertical dispersion is caused by not varying the other rate coefficients. Finally and for more clarity, we depict, in the lower panel, the calculated O_{2} abundance as a function of / obtained by calculations in which the ionization rates are not varied randomly. Rather, the solid lines show values of the steady-state O_{2} abundance computed for 5 values of (between 2 and /2) with varied to account for the range of the abscissa between 0.05 and 5. In the same panel, dotted lines show the O_{2} abundances computed for 5 values of with varied to account for the range of the abscissa. The abundances computed by each method are similar especially at the inflection point, confiming that the choice of the ionization ratio as a parameter of the sensitivity is a reasonable one.
The abundance of O_{2} decreases by two orders of magnitude (from 10^{-5} to 10^{-7}) when / goes from 0.3 to 0.8, and the inflection point of the curve occurs at the standard value of / used in osu.2003 (0.5/0.93). In Fig. 5 we show some other examples of this sensitivity by using five different values of , as in the solid lines of the right panel of Fig. 4. The abundance of atomic carbon has an almost equal sensitivity albeit opposite to that of O_{2}, while the H_{2}O abundance shows a similar profile to O_{2} except for a local maximum at / around 0.2. Although the electronic abundance increases towards higher / , it varies only a factor of 2.5 between / and / . The abundance of the H_{3}^{+} ion, which does not show any bimodal distribution, is linear with / . The dispersion of the abundances, for a given ratio / , due to the variation of by a factor of 2 around the standard value , depends on the species: it is smaller for H_{2}O, O_{2} and e^{-} than for C and H_{3}^{+}, for instance, but remains anyway below 0.5 dex.
Figure 6: Steady-state abundances of C, e^{-}, H_{2}O and H_{3}^{+} as functions of / for three different elemental abundances. The triangles, diamonds and crosses represent the low-, intermediate-, and high-metal cases. The vertical dashed line indicates the standard value of / . The H_{2} density is 10^{4} cm^{-3} and the temperature 10 K. The value of is fixed at its standard value while is varied. | |
Open with DEXTER |
Figure 7: Distribution (nomalized with respect to the total number of runs) of the steady-state abundance of atomic carbon for the three sets of elemental abundances. The vertical dashed line represents the abundance computed with the standard values of the rate coefficients. | |
Open with DEXTER |
The calculated sensitivity depends on the elemental abundances used for the modeling. In Fig. 6, we show the C, H_{2}O, e^{-} and H_{3}^{+} fractional abundances as a function of / for the low-, intermediate-, and high-metal elemental abundances (see Sect. 2). The sensitivity plots, as well as most others later in the paper, are obtained by fixing the helium ionization rate at its standard value so that no dispersion is seen. The H_{2} density remains at 10^{4} cm^{-3} and the temperature at 10 K. Although the sensitivity exists for the three cases, the inflection point is shifted to higher values of / for the low-metal case. For the intermediate case, on the other hand, the sensitivity seems to be stronger, leading to an abrupt jump in the abundance of C for / . A similar abrupt change can be seen for the electronic abundance. This jump is actually a manifestation of bistability (see Sect. 5). Note that for / , the abundances for the intermediate-metal case tend to be close to the low-metal ones whereas for / they are similar to the high-metal ones.
When one runs the uncertainty method with a factor of two uncertainty in and , the values of / are spread approximatively between 0.13 and 2.2. The distribution of the abundances computed with this method thus depends on the sensitivity within this range of values. The distribution of the steady-state abundances of C computed by variation of all the rate coefficients are shown in Fig. 7 for the three sets of abundances. Because the inflection point for the low-metal case is shifted to / , the uncertainty method gives Gaussian-like distributions although slightly asymmetrical towards higher C abundances. This explains why we did not see the sensitivity effect in our previous work on rate coefficient uncertainties in dark clouds (see Wakelam et al. 2006). For high-metal elemental abundances, the point of inflection occurs close to the standard value of / and the distribution for the C abundance is bimodal. When using the intermediate-metal abundances, the uncertainty method gives two distinct peaks with no solutions between them. The peak at higher abundance is the smaller since the point of inflection occurs at values of / somewhat larger than the standard. The lack of any model runs with steady-state abundances in between the two distributions stems from the sharpness of the curve around the inflection point and is a consequence of the bistability.
Figure 8: Atomic carbon abundance at steady state as a function of / for the three elemental abundances (low-, intermediate-, and high-metal). In each plot, the symbols refer to different C/O ratios. The vertical dashed line indicates the standard value of / . | |
Open with DEXTER |
In addition to its dependence on metallicity, the sensitivity to / depends strongly on the C/O elemental abundance ratio. In Fig. 8, we show the steady-state abundance of atomic carbon as a function of / computed for the three sets of elemental abundances, each with three different values of C/O - 0.2, 0.4 (the standard one), and 0.6 - obtained by maintaining a constant carbon abundance while varying the oxygen. The lower the C/O ratio is, the more the inflection point is shifted towards higher values of / . As an example, the inflection point for the high-metal case occurs around 0.2 for C/O = 0.6 whereas it rises to 1.5 for C/O = 0.2. In both cases, the inflection point is too far from the standard value of / to obtain a bimodal distribution. For the intermediate-metal case, the C/O ratio of 0.6 results in shifting the inflection point to the standard value of / , so that the analogous panel to the one in Fig. 7 would show two equal areas for the two peaks. The dispersion of the carbon abundance caused by the variation of C/O is larger for values of / smaller than the inflection point while the abundances eventually form a single plateau after the jump.
Figure 9: The left panel shows the distribution of the steady-state abundance of atomic carbon computed with a random variation of and for the high-metal elemental abundances but with a lower helium elemental abundance of 0.18 compared with H_{2}. The right panel represents the C steady-state abundance as a function of / for the three elemental abundances with the low He abundance. The vertical dashed line indicates the C abundances computed with the standard and on the left plot and the standard value of / on the right plot. | |
Open with DEXTER |
As discussed in Wakelam & Herbst (in preparation), the helium abundance used in the classic low- and high-metal elemental abundances (see Sect. 2) is quite high (0.28 compared with H_{2}) compared with the more modern value of 0.18 measured in the interstellar medium (see, for example, Baldwin et al. 1991). Naturally, the sensitivity will change if a lower value of He is used. Figure 9 shows the distribution of the steady-state C abundance computed with a random variation of and and the high-metal case with the lower helium abundance. In this case, we do not obtain a clear bimodal distribution but merely an asymmetrical profile, because the inflection point is shifted towards higher values of / , as shown in Fig. 9, and the random method barely reaches the sensitive portion of the curve.
Figure 10: Abundance of C as a function of / computed with the high- (black lines) and low- (grey lines) metal abundances. The left and right panels show the results using the Srates (see text) and rate99 networks. The vertical dashed line indicates the standard value / . | |
Open with DEXTER |
To test the influence of the chemical network utilized, we constructed the sensitivity curves at steady-state with two other databases - Srates and rate99 - using both low-metal and high-metal abundances. The rate99 network is the UMIST list of reactions (Le Teuff et al. 2000) prior to its current updating (see http://www.udfa.net) whereas Srates is the list of reactions used by Wakelam et al. (2004) to study sulphur chemistry in hot cores. This latter network contains only five elements (H, He, C, O and S) but, as discussed in Wakelam et al. (2004), the steady-state abundances computed with it are very close to the ones computed with the larger osu.2003 network at low temperature. The ionization rates for H_{2} and He are the same in all three databases.
Figure 10 is a sensitivity plot of the C abundance for five different values of between 2 and /2. The sensitivity obtained with Srates is similar to that for osu.2003, which indicates that the absence of the other elements (Si, Fe, Na etc.) does not change the sensitivity drastically. On the contrary, the results with rate99 show a much different relation between / and the atomic carbon abundance, especially for the high-metal case, where the rate99 curves show (a) a significantly lower value for the point of inflection of the sensitivity, (b) a wide dispersion due to the value of at low values of / , and (c) a much smaller jump in the abundance. We have shown in Wakelam et al. (2006) that rate99 and osu.2003 often have differences in rate coefficients larger than the estimated random errors for unstudied reactions of a factor of 2, whereas most of the rate coefficients in Srates are taken from the osu database. Since the three networks have the same and , the shape of the sensitivity can depend only on the values of other rate coefficients.
Figure 11: Illustrative scheme of hysteresis surface demonstrating bistability with two "orthogonal'' control parameters, here and / . The position of unstable solutions is indicated by dashed lines. Orthogonal cuts crossing the folding region produce hysteresis curves like the ones shown in Figs. 12 and 13. | |
Open with DEXTER |
Figure 12: Hysteresis curves of the steady-state abundance of atomic carbon as a function of / for three different elemental abundances (low-, high-, and intermediate-metal), and for five H_{2} densities between 10^{2} and 10^{6} cm^{-3}. The vertical dashed line indicates the standard value of / . Grey zones indicate regions of bistability. | |
Open with DEXTER |
In Sect. 4, we saw that the variation of / leads to a sharp change in the abundances of some species for the intermediate-metal case (see Fig. 6). This jump is indeed a manifestation of the same bistability discovered by Le Bourlot et al. (1993), as can be determined by comparing the region of the bistability in parameter space. The parameter usually used to study this bistability is the ratio between the cosmic-ray ionization rate and the density but other model parameters such as the temperature and rate coefficients have been used (Pineau des Forêts & Roueff 2000; Le Bourlot et al. 1995). What we have found here is that / is also a control parameter which explore another dimension of the bistability. If we imagine / to be an orthogonal axis to /, the bistability of Le Bourlot et al. (1993) is the same as ours when / is at its standard value. In Fig. 11, we show an illustrative scheme of this bistability as controlled by the two parameters and / . To make the bistability appear, we constructed hysteresis curves as a function of / for different H_{2} densities. These curves, shown in Fig. 12 for atomic carbon, were produced in two steps. Fist, we ran the model to compute the steady-state abundance using one value of / , then we used this solution as the initial condition and increased / to compute the next point, shown as the crosses in Figs. 12 and 13. In addition, the same principle but with decreasing / was utilized, with results shown as diamonds. In order to vary / , we fixed the helium ionization rate and changed that of hydrogen. When two solutions exist for a range of / values, a region of bistability exists.
To explore the domains of bistability using the new control parameter, we constructed hysteresis curves using the three different elemental abundances and five densities between 10^{2} and 10^{6} cm^{-3}. All the results are shown in Fig. 12 for atomic carbon. Bistability exists for all the elemental abundances but the domain of
/
where the two solutions coexist are small. To show the bistability more clearly, a magnification of this region is done for n(H_{2}) = 10^{4} cm^{-3} and the intermediate-metal case. In both the low- and high-metal cases, bistability exists only at very low density (10^{2} cm^{-3}) whereas in the intermediate-metal case, the bistability exists between 10^{2} and 10^{4} cm^{-3}. The increase in density shifts and smoothes the inflection point, making the bistability eventually disappear. Note that this is the first time that the bistability has been found with the high-metal case (Lee et al. 1998; Boger & Sternberg 2006), because it does not exist for the standard value of
/
.
We have confirmed that, for the high-metal case, using the standard value of
/
leads to no bistability with
as low as 10 cm^{-3}. Using
/
however, bistability exists in a very small density range of 70-90 cm^{-3}.
Figure 13: Hysteresis curves of the electron steady-state abundance as a function of ( left panel) and / ( right panel). Intermediate-metal elemental abundances are used for both plots. Five H_{2} densities between 10^{2} and 10^{6} cm^{-3} are used for the right panel. The vertical dashed line indicates the standard value of / . Grey zones indicate regions of bistability. | |
Open with DEXTER |
A more focused look at the two control parameters is shown in Fig. 13, where we plot the hysteresis curves of the electron abundance as a function of each parameter for the case of intermediate-metal elemental abundances. The value of is set at the standard value. For the variation of / (right panel), we used five H_{2} densities between 10^{2} and 10^{6} cm^{-3}, as in the previous figure, while for the hysteresis curve as a function of (left panel) the standard value of / (0.54) was used and the bistability found to exist in the density range 60-150 cm^{-3}. With the other control parameter and for cm^{-3}, the bistability exists for / between 0.5 and 0.7 and the domains of bistability using both parameters overlap. For higher densities, the bistable region using / shifts towards higher values and so does not exist anymore in the hysteresis curve versus . These hysteresis curves can be combined into a plot of the type shown in Fig. 11.
The sensitivity of the abundances to the ratio / is related to the non-linearity of the system. A detailed mathematical analysis of this phenomenon is in progress (Massacrier et al., in preparation). It is possible, however, to understand the difference, between the low and high / phases qualitatively from a chemical point of view. For this purpose, we have considered the main reactions of production and destruction of O, C, CO, S and SO, which are the main reservoirs of carbon, oxygen and sulphur for the intermediate- and low-metal elemental compositions in which sulfur is relatively abundant, and/or are the species most influenced by and . Let us first consider low values of / . The ionization of H_{2} leads to the formation of OH through a well-known series of reactions (see. e.g. Boger & Sternberg 2006). The radical OH then reacts with atomic oxygen to form O_{2}. The most efficient way of destroying O_{2} is to form SO by reaction with S. Since the elemental abundance of atomic sulphur is much less than that of oxygen, large amounts of O_{2} remain in the gas phase. Atomic carbon reacts with O_{2} to form CO where it is stored since CO is not efficiently destroyed. For high / , on the other hand, the destruction of CO by He^{+} produces C^{+} and O more efficiently. Carbon ion then exchanges its charge with S so large abundances of C remain in the gas phase to deplete molecular oxygen.
Since the bistability we see is an extreme manifestation of the sensitivity to / , the chemical processes involved in bistability may be related to those discussed above, although our control parameter is different from those previously studied. Pineau des Forets et al. (1992) and Le Bourlot et al. (1993,1995) identified the abundance ratio H^{+}/H_{3}^{+} as being important for bistability, with one ion identified with the HIP and the other with the LIP solution. The abundance ratio H^{+}/H_{3}^{+} depends on the abundance of electrons and as a consequence on the efficiency of the dissociative recombination of H_{3}^{+} with electrons, the rate of which was quite uncertain at the time. In our analysis however, the H_{3}^{+} ion shows little if any sensitivity to the control parameter, although the electron abundance does show a jump between the low and high / phases, which can then be thought of as analogous in some sense to the LIP and HIP solutions. Recently, Boger & Sternberg (2006) made another detailed analysis of the chemical differences between the so-called HIP and LIP solutions, arguing that bistability is due to a loop including H_{3}^{+}, O_{2} and S^{+}. Although this explanation is related in some aspects to that given for the sensitivity reported here, one must remember that our control parameter raises some new points, given the relation of the ionization rate of He to the formation of the important ion C^{+}, which is a precursor for many species, and can charge exchange to form abundant C. Moreover, S^{+}, which is prominent in the explanation of Boger & Sternberg (2006), is not the dominant charge carrier in our calculations. Note that the osu.2003 database includes the recombination of the main ions with negative grains and this does not cancel the bistability, although such a cancellation was suggested by Boger & Sternberg (2006). It may well be that the H^{+}/H_{3}^{+} abundance ratio, the H_{3}^{+}-O_{2}-S^{+} cycle, and the / ratio comprise pieces of a larger and more complex jigsaw puzzle. Since bistability also appears to be a function of the reaction network utilized, secondary reactions may play a more important role than heretofore realized. This possibility has to be clarified through the building of a simplified reaction network that still shows the main characteristics of bistability.
The sensitivity to / depends highly on the model parameters and, for this reason, it is not obvious if such effect can be detected in the interstellar medium. In an ideal situation where the cloud characteristics (at least the elemental abundances and the age) are well known and correspond to a distinct chemical regime, one could use the observations to constrain the ratio / assuming however that no local variation of this ratio is expected. In reality, the cloud characteristics are in general far from being well constrained, and an analysis of the uncertainties in the ionization rates of H_{2} and He is strongly required to understand the relevance of this sensitivity for the interstellar medium. Note that / is not sensitive to the total intensity of the cosmic-ray flux.
If the intermediate-metal elemental abundances reflect the cloud characteristics and if the ionization rate of He and H_{2} by cosmic-rays can vary by at least a factor of 2 within a cloud, the prediction of bistability over a wide range of densities suggests that strong chemical heterogeneities in dense clouds can be caused by this effect. In the absence of bistability an extreme sensitivity is predicted at most densities for low-metal and high-metal abundances and may cause some observed heterogeneities in abundance. However, we have to consider the age of the cloud. The bifurcation occurs just after 10^{5} yr, starting from our standard initial conditions, and is a maximum at steady-state (10^{7} yr). As a consequence, the chemical heterogeneity can be caused by this sensitivity only in evolved clouds ( 10^{5} - 10^{6} yr).
The large gas-phase abundances of H_{2}O and O_{2} molecules predicted by dense cloud chemical models is a long standing problem. Using high values of / , we obtain much smaller abundances for these species, closer to the upper limits in dense clouds (Bergin & Snell 2002; Pagani et al. 2003). Large uncertainties in the / ratio could then solve the problem related to the modeling of water and molecular oxygen. One has, however, to consider the overall picture, and the better agreement with H_{2}O and O_{2} should not worsen the agreement for the other species. A comparison between observations and modeling should be done in a systematic way, as done by Wakelam et al. (2006). This was not the purpose of the present study, but we can already note, as an example, that the observed abundances in L134N (North peak) of OH and C_{4}H ( and 10^{-9} respectively, Ohishi et al. 1992) are reproduced by the model when using the standard value of / but are underestimated when using the higher ratios / that reproduce H_{2}O and O_{2} abundances.
Finally, the very low H_{2} density (100 cm^{-3}), at which the bistability as a function of / is restricted to with the low- and high-metal elemental abundances, is not relevant for clouds with the assumed visual extinction ( ) since it would require a size of 26 pc, which is much larger than typical cloud sizes (Snell 1981).
Using a Monte Carlo uncertainty method to compute the theoretical error of molecular abundances in chemical models, we found a strong sensitivity of the steady-state abundances to the ratio between the ionization rate coefficient of helium and that of molecular hydrogen by cosmic rays. This sensitivity to the ratio has the consequence of changing the abundance of some species by several orders of magnitude for small variations of / . Lower values of / are characterized by large abundances of smaller molecules such as O_{2}, H_{2}O, SO, etc. whereas high values of / lead to large abundances of atoms. If the sensitivity is strong enough, it can result in the phenomenon of bistability, which we located at low densities for two sets of elemental abundances - the low- and high-metal cases, and over a wide range of densities for the case of the so-called intermediate-metal abundances, in which sulfur plays a particularly prominent role. The bistability has been determined to be the same as studied earlier by a variety of authors with other control parameters; it is perhaps simplest to consider our control parameter - / - as orthogonal to one of the others - the ratio of the overall ionization rate parameter to the density n - so that the bistability can be imagined to exist in a region inthree-dimensional space rather than in the customary two-dimensional hysteresis picture.
It is interesting to speculate that the sensitivity detected here via our uncertainty analysis is a necessary but not sufficient condition for bistability because only the most extreme sensitivity, in which the bimodal distributions detected have no solutions in between them, leads to bistability. Unfortunately, it is not guaranteed that our type of uncertainty analysis can always locate a bistable region since finding an extreme sensitivity may require large uncertainties in some parameters, such as the gas density, to take us from our standard parametric values to the bistable region.
Several questions arise when considering the dramatic effects of small / variations on the gas-phase chemistry. First, even if one assumes a given flux and energy distribution of the cosmic rays, the determination of / suffers some significant uncertainties, with the role of secondary electrons particularly important. Besides these difficulties, one might also have to consider variations of the ratio / in space, such as between one cloud and another, or between the center and the edges of a same cloud, and in time. Changes in the ratio imply important variations in the spectral distribution of the cosmic rays (Glassgold & Langer 1973). Such variations might be produced by the attenuation of the cosmic rays inside a cloud or the proximity of supernovae in star forming clouds. Whether such variations are expected or not and how they could affect the chemistry should be considered in the future. Finally,a knowledge of the helium elemental abundance appears to be crucial since the sensitivity region depends on this parameter.
Acknowledgements
V. W. and E. H. thank the National Science Foundation for its partial support of this work. G.M. acknowledges partial financial support from the CNRS/INSU programs PCMI and PNPS.