Astrochemical Bistability: Autocatalysis in Oxygen Chemistry

The occurrence of bistable solutions in models of dense cloud chemistry has been known for many years (Pineau des Forêts et al. 1992; Le Bourlot et al. 1993). These are characterized by the coexistence of two stable states connected by an unstable branch. The stable states occur in a molecular hydrogen gas and comprise a high ionization phase (the so-called HIP) in which the saturated species are under-abundant and protonation reactions are less efficient; and a low ionization phase in which saturated molecular ions are more abundant (the LIP). Several studies have shown that the occurrence of the bistable region is controlled by various combinations of the cosmic-ray ionization rates of H₂ and He, the gas density, the relative depletion of the heavy elements, interaction with dust grains, and the electron dissociative recombination rate of ${{\rm{H}}}_{3}^{+}$ (Le Bourlot et al. 1993, 1995a, 1995b; Shalabiea & Greenberg 1995; Lee et al. 1998; Forêts & Roueff 2000; Viti et al. 2001; Charnley & Markwick 2003; Boger & Sternberg 2006; Wakelam et al. 2006). However, a definitive understanding of the origin of these bistable solutions has not emerged.

Bistability is of interest because it occurs within the region of parameter space that includes the depletions and densities relevant to dense interstellar clouds. It has been proposed as the origin of the high CI/CO ratios and low O₂ abundances observed in molecular clouds (Le Bourlot et al. 1993; Viti et al. 2001), in cloud chemistry near supernova remnants (Ceccarelli et al. 2011), in Galactic spiral arm clouds (Tieftrunk et al. 1994), and in high-latitude clouds (Gerin et al. 1997). Understanding the underlying bistable solutions may provide important insight into other astronomical environments.

An early explanation of the underlying chemical mechanism concerned an ionization instability connected to the difference in the respective rates of radiative and dissociative electron recombination of H⁺ and ${{\rm{H}}}_{3}^{+}$ (Pineau des Forêts et al. 1992). Boger & Sternberg (2006) also considered an ionization instability origin in which bistability could develop in an elementary ${{\rm{H}}}_{3}^{+}$–O₂–S⁺ cycle. Bistability cannot be solely due to the presence of S⁺ since calculations show that bistable solutions occur even when S chemistry is neglected (Lee et al. 1998). It is well known in chemical kinetics that nonlinearity and feedback (i.e., autocatalysis) in the governing differential equations can lead to regions of multistability and complex chemical evolution (e.g., Epstein & Showalter 1996). Le Bourlot et al. (1993) suggested that autocatalysis could be responsible and here we demonstrate the fundamentally autocatalytic nature of astrochemical bistability.

For prescribed initial conditions, the molecular evolution of a static dark cloud is obtained by solving a nonlinear system of ordinary differential equations for the abundances of N chemicals, subject to conservation of charge and elemental nuclei (e.g., Nejad 2005). Bifurcations in the solutions are found by solving for the steady-state abundance

$$\begin{eqnarray}&&{\boldsymbol{F}}({\boldsymbol{x}};\zeta ,{n}_{{\rm{H}}},T,\beta ,{\delta }_{j},{\delta }_{{\rm{M}}})=0,\end{eqnarray} \tag{ 1 }$$

where ζ is the cosmic-ray ionization rate, n_H is the number density of hydrogen nuclei (≈2n(H₂)), T is the gas kinetic temperature, β is the photorate associated with the cosmic-ray induced radiation field (Prasad & Tarafdar 1983; Gredel et al. 1989), δ_j are the depletion factors for the major volatile elements and δ_M is a depletion factor for refractory metals (Na in this paper). The fractional abundance of any species is x(i) = n(i)/n_H and the total fractional abundance of each element, j, is given by δ_j X_j, with X_j relative to total H nuclei. An additional parameter known to influence bistability is the ${{\rm{H}}}_{3}^{+}$ electron recombination rate, α₃ (Pineau de Forêts & Roueff 2000). We initially adopted a value ten times larger than Millar et al. (1997) (Table 1), twice the value of $1.5\times {10}^{-7}{T}_{3}^{-0.5}$ ${\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$ used by Le Bourlot et al. and at the higher end of the range covered in previous studies (T₃ = T/300 K).

Table 1.  Dense Cloud Chemistry Model Parameters

Temperature	T	10 K
Visual extinction	AV	15 mag
Cosmic-ray ionization rate	ζ	5 × 10−17 s−1
${{\rm{H}}}_{3}^{+}$ electron recombination rate†	α3	$3\times {10}^{-7}{T}_{3}^{-0.5}$ ${\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$
Elemental abundances‡	XO	8.53 × 10−4
	XC	3.62 × 10−4
	XS	1.85 × 10−5
	XNa	1.5 × 10−8
	XHe	0.1

Download table as:  ASCIITypeset image

We solve the system of nonlinear Equations (1) by Newton–Raphson iteration and explore how the solutions undergo bifurcations as the control parameters are varied.

2.1. Model Calculations

To determine the underlying cause of bistability we begin with a known bistable solution and sequentially relax the δ_j control parameter. Six models are sufficient to illustrate the processes leading to astrochemical bistability. Table 1 lists the physical conditions common to each model. The relevant parameters of each model are listed in Table 2. The chemical network is derived from Millar et al. (1997). Figures 1(a)–(f) show the fractional abundances of major chemical species as a function of n_H in each of Models 1 to 6.

Zoom In Zoom Out Reset image size

Table 2.  Bifurcation Models

Model	δO	δC	δS	δNa	$\beta \ne 0$	Note
1	0.1	0.1	0.1	1.0	Yes	Reference model
2	0.1	0.1	0.01	1.0	Yes
3	0.1	0.1	0	0	Yes
4	0.1	0.1	0	0	No
5	0.1	0	0	0	No
6	0.1	0	0	0	No	As Model 5 but δHe = 0

Download table as:  ASCIITypeset image

Model 1 is taken as a reference model (see Charnley & Markwick 2003 and Le Bourlot et al. 1993). It includes the internal radiation field with carbon, oxygen, and sulfur, depleted by a factor of 0.1. Figure 1(a) shows that the HIP/LIP bistable region occurs between ${n}_{{\rm{H}}}\approx (1.8\mbox{--}4.9)\times {10}^{4}\,{\mathrm{cm}}^{-3}$. The bifurcations consist of two stable solutions connected by an unstable one (e.g., Drazin 1992). In Model 2 the sulphur depletion is increased by a factor of 10, δ_S = 0.01. Model 3 shows the effect of removing both S and Na from the network, δ_S = δ_Na = 0. In Model 4 the radiation field is turned off. Model 5 has carbon removed, δ_C = 0, so that only chemical reactions involving oxygen, hydrogen, and helium remain. Finally, Model 6 is the same as Model 5 except that it contains no helium.

Figures 1(e) and (f) (Models 5 and 6) demonstrate that the cause of astrochemical bistability lies in the oxygen chemistry, and is unconnected to the presence of metals, of sulphur, or even of carbon.

2.2. Oxygen Chemistry: Autocatalysis and Feedback

We identify the fundamental mechanism that controls the occurrence of bistable solutions as autocatalysis (e.g., Gray & Scott 1994) in the reduced Models 5 and 6. The first step is formation of the oxygen dimer

$$\begin{eqnarray}&&{\rm{O}}\,+\,\mathrm{OH}\,\longrightarrow \,{{\rm{O}}}_{2}\,+\,{\rm{H}}\end{eqnarray} \tag{ 2 }$$

followed by release of two oxygen nuclei in either of two autocatalytic steps

$$\begin{eqnarray}&&{{\rm{H}}}^{+}\,+\,{{\rm{O}}}_{2}\,\longrightarrow \,{{\rm{O}}}_{2}^{+}\,+\,{\rm{H}}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{O}}}_{2}^{+}\,+\,{{\rm{e}}}^{-}\,\longrightarrow \,{\rm{O}}\,+\,{\rm{O}}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{\mathrm{He}}^{+}\,+\,{{\rm{O}}}_{2}\,\longrightarrow \,{{\rm{O}}}^{+}\,+\,{\rm{O}}\,+\,\mathrm{He}.\end{eqnarray} \tag{ 5 }$$

The oxygen ions and atoms then react to produce OH⁺ to initiate two feedback pathways

$$\begin{eqnarray}&&{{\rm{O}}}^{+}\,+\,{{\rm{H}}}_{2}\,\longrightarrow \,{\mathrm{OH}}^{+}\,+\,{\rm{H}}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\rm{O}}\,+\,{{\rm{H}}}_{3}^{+}\,\longrightarrow \,{\mathrm{OH}}^{+}\,+\,{{\rm{H}}}_{2}\end{eqnarray} \tag{ 7 }$$

and subsequently become reincorporated into OH in the rapid sequence

$$\begin{eqnarray}&&{\mathrm{OH}}^{+}\mathop{\longrightarrow }\limits^{{{\rm{H}}}_{2}}{{\rm{H}}}_{2}{{\rm{O}}}^{+}\mathop{\longrightarrow }\limits^{{{\rm{H}}}_{2}}{{\rm{H}}}_{3}{{\rm{O}}}^{+}\mathop{\longrightarrow }\limits^{{{\rm{e}}}^{-}}\mathrm{OH}\end{eqnarray} \tag{ 8 }$$

to complete the cycle. This is an example of indirect autocatalysis (e.g., Plasson et al. 2011).

Thus, we can identify two feedback routes to bistability: (1) an O⁺ pathway initiated by He⁺, and (2) a pure O atom pathway initiated by H⁺. In Model 6 (no He) the O atom pathway is clearly solely responsible. As one moves across the bistable region from the HIP to the LIP in Model 5, Figure 1(e) shows that, although H⁺ is the dominant ion in the HIP, close to the the bifurcation point in the LIP, He⁺ is the dominant ion. Thus, He⁺ plays an important role in producing bifurcations. Both pathways operate together in Model 5 but we find that one can come to dominate the other, depending on the value of α₃.

Figure 2(a) shows that when the rate of reaction (5) is artificially set to zero in Model 5 the removal of autocatalysis and the O⁺ pathway causes bistability to disappear. By contrast, in Model 6 (with no He), n_e ≈ n (H⁺), autocatalysis is only through reactions (3) and (4) and the higher H⁺ abundance drives reaction (3) more efficiently. As $n({{\rm{H}}}_{3}^{+})$ ∝${({\alpha }_{3}{n}_{{\rm{e}}})}^{-1}$, bistability via the O + H₃⁺ pathway can be recovered in the model of Figure 2(a) by reducing α₃(10 K) from 1.64 × 10⁻⁶ ${\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$ (see Table 1). Figure 2(b) shows the effect of reducing α₃(10 K) by a factor of 10; the resulting higher ${{\rm{H}}}_{3}^{+}$ abundance makes the O atom pathway more efficient and the bistable solution appears. Hence, in Model 5 whichever pathway produces bistability depends on the value of α₃ relative to some critical value, α_crit; by numerical experiment we find that α_crit ∼ 10⁻⁶ ${\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}$. Artificially setting the rate coefficient of reaction (4) to zero in either Model 5 or 6, removes the O atom pathway, causing bistability to disappear.

Zoom In Zoom Out Reset image size

2.3. Dense Cloud Chemistry Reconstructed

In more realistic chemical models of molecular clouds, the inclusion of carbon and sulfur means that C⁺ and S⁺ can replace or complement He⁺ and H⁺ as the major ion destroying O₂. Figures 1(b)–(d) show that C⁺ is the most abundant ion at the HIP and LIP. When carbon is added to the oxygen chemistry of Model 5, C⁺ is formed by

$$\begin{eqnarray}&&{\mathrm{He}}^{+}\,+\,{\rm{C}}\,\longrightarrow \,{{\rm{C}}}^{+}\,+\,\mathrm{He}\end{eqnarray} \tag{ 9 }$$

and by reaction with CO. In Model 4 the autocatalytic step of reaction (5) is replaced by the net effect of the two product branchings in

$$\begin{eqnarray}&&{{\rm{C}}}^{+}\,+\,{{\rm{O}}}_{2}\,\longrightarrow \,{\mathrm{CO}}^{+}\,+\,{\rm{O}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\longrightarrow \,\mathrm{CO}\,+\,{{\rm{O}}}^{+}.\end{eqnarray} \tag{ 11 }$$

Although reaction (10) followed by

$$\begin{eqnarray}&&{\mathrm{CO}}^{+}\,+\,{{\rm{e}}}^{-}\,\longrightarrow \,{\rm{C}}\,+\,{\rm{O}}\end{eqnarray} \tag{ 12 }$$

could in principle be an autocatalytic process, we find that the rapid destruction of CO⁺ by H₂ renders its contribution negligible.

The autocatalytic nature of bistable solutions can also be confirmed through artificial models, as in Section 2.2. This is illustrated in Figure 2(c) which shows that the bifurcations in Model 4 (see Figure 1(a)) disappear once the rate coefficients of reactions, (10) and (11) have been set to zero. As in the case of the artificial models of Figures 2(a) and (b), bistability can be recovered by reducing the value of α₃, as shown in Figure 2(d). In this case we find that reaction (5) is responsible for the bistable solution. We find that acting together these three reactions, (5), (10), and (11), are solely responsible for bistability in Model 4 independent of α₃. In similarly artificial versions of Model 3 (not shown), we find the same autocatalytic effects except that, at low α₃, ${{\rm{O}}}_{2}^{+}$ ions from O₂ photoionization act in concert with reaction (4). In this model, at low α₃, instead of reaction (5), reactions (3) and (4) combined with the C⁺ + O₂ autocatalytic step, are responsible for the bistable solution.

When sulphur is included at sufficiently high abundance, the reactions

$$\begin{eqnarray}&&{{\rm{C}}}^{+}\,+\,{\rm{S}}\,\longrightarrow \,{{\rm{S}}}^{+}\,+\,{\rm{C}}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\rm{H}}}^{+}\,+\,{\rm{S}}\,\longrightarrow \,{{\rm{S}}}^{+}\,+\,{\rm{H}}\end{eqnarray} \tag{ 14 }$$

can lead to S⁺ being the dominant atomic ion in the bistable regime as shown in Model 1 (see Figure 1(a)), and the autocatalytic production of oxygen atoms can occur through

$$\begin{eqnarray}&&{{\rm{S}}}^{+}\,+\,{{\rm{O}}}_{2}\,\longrightarrow \,{\mathrm{SO}}^{+}\,+\,{\rm{O}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\mathrm{SO}}^{+}\,+\,{{\rm{e}}}^{-}\,\longrightarrow \,{\rm{S}}\,+\,{\rm{O}}.\end{eqnarray} \tag{ 16 }$$

In Model 2 the δ_S is sufficiently large that S⁺ at the LIP bifurcation is not efficiently produced at the expense of C⁺ (see Figure 1(b)) and the origin of the bistable solutions are the same as in Model 3.

Returning to the reference model (Model 1), Figure 2(e) shows that, as before, artificially removing only reactions (10) and (11) causes bistability to disappear in the high α₃ case; Figure 2(f) shows that it again reappears at the lower value of α₃. For Model 1 we have found that all the four identified autocatalytic processes can combine to act in concert in realistic models of molecular clouds. With the C⁺+ O₂ autocatalytic reactions removed, as in Figure 2(f), further removing reaction (5) removes the bifurcation, as shown in Figure 2(g); this bifurcation also disappears if reaction (4) is removed instead of (5). Model 1 has the lowest S depletion and so S⁺ can play an important rôle in bistability. Figure 2(h) shows the effect of only removing reaction (16).

Thus, the occurrence of bistability sensitively depends on four autocatalytic processes acting in unison: driven by O₂ reacting with He⁺ and C⁺, and with H⁺ and S⁺.

These bistable solutions are produced in a gas-phase chemistry driven by cosmic rays. In dense molecular clouds most of the grains are negatively charged (G⁻) and only about 4% are neutral (G⁰) (e.g., Umebayashi & Nakano 1990). Atomic ion (X⁺)-grain recombination

$$\begin{eqnarray}&&{{\rm{X}}}^{+}\,\,+\,\,{{\rm{G}}}^{-}\,\,\longrightarrow \,\,{\rm{X}}\,\,+\,\,{{\rm{G}}}^{0}\end{eqnarray} \tag{ 17 }$$

will dominate electron radiative recombination, and so could affect the onset of bistable solutions. For molecular ions, ion–grain recombination is never competitive with electron dissociative recombination in our models. Published studies that have considered the gas–grain interaction disagree as to whether or not bistable solutions can persist (Le Bourlot et al. 1995b; Wakelam et al. 2006), or are suppressed (Shalabiea & Greenberg 1995; Boger & Sternberg 2006). Shalabiea & Greenberg (1995) found that while the gas–grain interaction did suppress bistability for a model adopting diffuse ISM elemental depletions, bistability still occurred for the dense cloud values of Le Bourlot et al. (1993). Le Bourlot et al. (1995b) presented a more detailed treatment of the effects of the gas–grain interaction and confirmed that it had no effect in suppressing bistability. Boger & Sternberg (2006) also adopted diffuse ISM depletions and demonstrated that atomic ion–grain neutralization could indeed suppress their bistable solutions, whereas Wakelam et al. (2006) reported no effect.

We can assess the effect of ion–grain neutralization processes on our results as follows. The steady-state number density of an atomic ion , X⁺, is

$$\begin{eqnarray}&&n({{\rm{X}}}^{+})\,=\,\displaystyle \frac{P({{\rm{X}}}^{+})}{{L}_{\mathrm{in}}+{k}_{{\rm{i}}}n({{\rm{O}}}_{2})+{\alpha }_{{\rm{r}}}{n}_{e}+\lambda ({{\rm{X}}}^{+})}\end{eqnarray} \tag{ 18 }$$

where P(X⁺) is the total ion production rate, k_i and α_r are the rate coefficients for reactions with O₂ and electron radiative recombination, L_in is the loss rate in all other ion–neutral reactions, and the ion–grain collision rate is λ (X⁺). Assuming a Maxwellian distribution for the ion velocities, and that all the grains are negatively charged (Umebayashi & Nakano 1990; Charnley 1997), we can write

$$\begin{eqnarray}&&\lambda ({{\rm{X}}}^{+})={\left(\displaystyle \frac{8{kT}}{\pi {M}_{{\rm{X}}}{m}_{{\rm{H}}}}\right)}^{1/2}\langle {n}_{\mathrm{gr}}(a){\sigma }_{{\rm{c}}}(a)\rangle ,\end{eqnarray} \tag{ 19 }$$

where k is the Boltzmann constant, ${m}_{{\rm{H}}}$ is the proton mass, T is the gas temperature (10 K), M_X is the ion mass (in a.m.u), n_gr(a) is the number density of grains of radius a, and σ_c(a) is the collision cross-section including the Coulomb factor. We can write

$$\begin{eqnarray}&&\lambda ({{\rm{X}}}^{+})\,=\,{C}_{\mathrm{gr}}{M}_{{\rm{X}}}^{-1/2}{n}_{{\rm{H}}},\end{eqnarray} \tag{ 20 }$$

where C_gr is the rate coefficient (${\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$) for ion–grain collisions and is given by

$$\begin{eqnarray}&&{C}_{\mathrm{gr}}\,=\,1.45\times {10}^{4}{T}^{1/2}\displaystyle \frac{\langle {n}_{\mathrm{gr}}(a){\sigma }_{{\rm{c}}}(a)\rangle }{{n}_{{\rm{H}}}}.\end{eqnarray} \tag{ 21 }$$

Grain neutralization can suppress a bistable solution produced by X⁺ reacting with O₂ when

$$\begin{eqnarray}&&\lambda ({{\rm{X}}}^{+})\,\gt \,{k}_{{\rm{i}}}n({{\rm{O}}}_{2})\end{eqnarray} \tag{ 22 }$$

or when C_gr exceeds a threshold value

$$\begin{eqnarray}&&{C}_{\mathrm{gr}}\,\gt \,{k}_{{\rm{i}}}{M}_{{\rm{X}}}^{1/2}y({{\rm{O}}}_{2}),\end{eqnarray} \tag{ 23 }$$

where y(O₂) = n(O₂)/n_H is the O₂ fractional abundance at the LIP bifurcation.

The bistable solutions are produced by S⁺ in both the Shalabiea & Greenberg model (S/H = 9 × 10⁻⁷) and in that of Boger & Sternberg (S/H ≈ 8 × 10⁻⁶). When X⁺ is S⁺, as in Model 1, ${k}_{{\rm{i}}}=1.5\times {10}^{-11}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$, $y({{\rm{O}}}_{2})\approx 8\,\times \,{10}^{-6}$, and bistability is suppressed when ${C}_{\mathrm{gr}}\gt 6.8\times {10}^{-16}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$, in good agreement with the threshold calculated by Boger & Sternberg (2006).⁴ However, in dense molecular clouds gaseous S is inferred to be depleted by a factor of ∼100 with respect to the diffuse ISM (e.g., Fuente et al. 2019), and so C⁺ rather than S⁺ drives bistability. When X⁺ is C⁺, as in Model 4, ${k}_{{\rm{i}}}=8.0\times {10}^{-10}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$, y(O₂) ≈ 8 × 10⁻⁶, and bistability is suppressed when ${C}_{\mathrm{gr}}\gt 2.2\times {10}^{-14}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$. Figure 3 shows the result of adding grain neutralization of atomic ions to Model 4, for two values of C_gr that bracket the threshold value, and confirms that even a large value of C_gr $(5.6\,\times {10}^{-15}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1})$ cannot suppress bistability.

Zoom In Zoom Out Reset image size

How realistic are the C_gr values required to cause bistable solutions to disappear? When $\langle {n}_{\mathrm{gr}}(a)$ ${\sigma }_{{\text{}}c}(a)\rangle$ corresponds to a single-size grain distribution with a dust/gas number density ratio appropriate for dense clouds (e.g., ${n}_{\mathrm{gr}}(0.1\,\mu {\rm{m}})=1.85\,\times {10}^{-12}{{\rm{n}}}_{{\rm{H}}}$ ) C_gr$\sim 3\times {10}^{-16}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1};$ it is about seven times larger when an MRN grain size distribution (Mathis et al. 1977) is considered (Le Bourlot et al. 1995b). The latter value exceeds the S⁺ threshold but neither value will have an effect on bistability arising from C⁺. For the ion–grain collision rates⁵ considered by Boger & Sternberg (2006), ${C}_{\mathrm{gr}}\sim 2.5\,\times {10}^{-16}\mbox{--}2.5\times {10}^{-14}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$, with bistability disappearing above $\approx 7.6\times {10}^{-16}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$, the largest collision rates ($\sim {10}^{-14}\,{\mathrm{cm}}^{3}\ {{\rm{s}}}^{-1}$) are unlikely to be relevant for dense clouds. These values pertain to MRN dust distributions that are appropriate to the diffuse interstellar medium, where the size distribution extends to very small carbonaceous particles, such as large PAH molecules (e.g., Weingartner & Draine 1999). Mid-IR observations of dense clouds however show clear evidence for very large dust grains (∼1 μm) and no evidence for PAHs, suggesting that small grains have been removed from the distribution by coagulation (Steinacker et al. 2010). Thus, for realistic dense cloud elemental depletions and grain-size distributions, we conclude that the threshold for ion–grain neutralization to suppress bistability through C⁺ reactions cannot be reached.

We have demonstrated that interstellar chemistry is bistable due to the interaction of several autocatalytic processes involving molecular oxygen. By deconstructing a known bistable solution into ever simpler reduced models through omission of chemical elements, and then artificially removing selected reactions, we have identified four distinct modes of autocatalysis that can occur in dense molecular clouds.

Our calculations provide explanations of results from earlier bistability studies: the dependence on C/O ratio (Le Bourlot et al. 1995a), the He⁺ production rate (Wakelam et al. 2006), and the reasons for bistability appearing with or without sulfur (Lee et al. 1998). Apart from the previously known dependence of bistability on the ${{\rm{H}}}_{3}^{+}$ recombination rate (Forêts & Roueff 2000), we find that α₃ can also mediate between the feedback pathways (via O⁺ or O) in the autocatalytic cycle. Internally generated UV photons play a minor role in these models (although their effects do become more pronounced in artificial models, see Section 2.3). Refractory metals have no autocatalytic function and bistable solutions can occur in their absence.

The bistable solutions discovered in dark clouds are controlled by ζ/n_H, the relative elemental depletions, and the value of the ${{\rm{H}}}_{3}^{+}$ electron recombination rate. Recent measurements of the elemental ratios C/H and C/O in the dark cloud TMC-1 are consistent with those where bistability can occur (Fuente et al. 2019). The associated S/H ratio in TMC-1 is very low, as is commonly found in molecular clouds, indicating that autocatalysis involving S will not be as important as the competing processes (see Model 2). Photoprocesses have a relatively minor effect in these models but could become important at lower A_V.

We have shown that previous conclusions regarding the efficiency of atomic ion–grain neutralization in suppressing dense cloud bistability depend on the applicability of elemental depletions and dust distributions that are more appropriate for the diffuse ISM. For more realistic dense cloud model parameters (high S depletion and reduced number of small grains), grain recombination has no effect on bistable solutions.

The bistable solutions present in dense cloud chemistry are due to autocatalysis and not an ionization instability.⁶ In fact, Le Bourlot et al. (1993, 1995a), in the first explicit demonstration of bistability, presciently suggested that oxygen autocatalysis through reactions (2)–(4) could play a role in interstellar bistability, as confirmed in the calculation presented here (Model 6, Figure 1(f)). Remarkably, this insight was neglected in all subsequent studies.

Gas-phase bistability is possible when an autocatalyst can form a dimer that is subsequently destroyed in an autocatalytic step followed by reformation of the dimer. As many such reactions are fundamental to chemical models of astronomical sources, bistability can be expected to be a common phenomenon (G. Dufour & S. B. Charnley 2019, in preparation).

This research was supported by NASA's Planetary Science Division Internal Scientist Funding Program through the Fundamental Laboratory Research work package (FLaRe).